root/OpenSceneGraph/trunk/src/osgPlugins/dicom/ReaderWriterDICOM.cpp @ 13041

Revision 13041, 42.4 kB (checked in by robert, 3 years ago)

Ran script to remove trailing spaces and tabs

  • Property svn:eol-style set to native
Line 
1// Released under the OSGPL license, as part of the OpenSceneGraph distribution.
2//
3// ReaderWriter for sgi's .rgb format.
4// specification can be found at http://local.wasp.uwa.edu.au/~pbourke/dataformats/sgirgb/sgiversion.html
5
6#include <osg/Image>
7#include <osg/Notify>
8#include <osg/Geode>
9#include <osg/GL>
10#include <osg/io_utils>
11#include <osg/ImageUtils>
12
13#include <osgDB/FileNameUtils>
14#include <osgDB/FileUtils>
15#include <osgDB/Registry>
16
17#include <osgVolume/Volume>
18#include <osgVolume/VolumeTile>
19#include <osgVolume/RayTracedTechnique>
20
21#ifdef  USE_DCMTK
22    #ifndef _MSC_VER
23        #define HAVE_CONFIG_H
24    #endif
25
26    #include <dcmtk/config/osconfig.h>
27    #include <dcmtk/dcmdata/dcfilefo.h>
28    #include <dcmtk/dcmdata/dcdeftag.h>
29    #include <dcmtk/dcmdata/dcuid.h>
30    #include <dcmtk/dcmimgle/dcmimage.h>
31    #include <dcmtk/dcmimgle/dimopx.h>
32    #include <dcmtk/dcmimage/dicopx.h>
33    #include "dcmtk/dcmimage/diregist.h"
34#endif
35
36#ifdef USE_ITK
37    #include <itkImageFileReader.h>
38    #include <itkImageFileWriter.h>
39    #include <itkImage.h>
40    #include <itkImageRegionConstIterator.h>
41    #include <itkMetaDataDictionary.h>
42    #include <itkMetaDataObject.h>
43    #include <itkGDCMImageIO.h>
44#endif
45
46#include <stdio.h>
47#include <stdlib.h>
48#include <string.h>
49
50#include <memory>
51
52class ReaderWriterDICOM : public osgDB::ReaderWriter
53{
54    public:
55
56        ReaderWriterDICOM()
57        {
58            supportsExtension("mag","dicom image format");
59            supportsExtension("ph","dicom image format");
60            supportsExtension("ima","dicom image format");
61            supportsExtension("dic","dicom image format");
62            supportsExtension("dcm","dicom image format");
63            supportsExtension("dicom","dicom image format");
64            // supportsExtension("*","dicom image format");
65        }
66
67        std::ostream& warning() const { return osg::notify(osg::WARN); }
68        std::ostream& notice() const { return osg::notify(osg::NOTICE); }
69        std::ostream& info() const { return osg::notify(osg::INFO); }
70
71        template<typename T>
72        T* readData(std::istream& fin, unsigned int length, unsigned int& numComponents) const
73        {
74            numComponents = length/sizeof(T);
75            T* data = new T[numComponents];
76            fin.read((char*)data, numComponents*sizeof(T));
77
78            // read over any padding
79            length -= numComponents*sizeof(T);
80            while(fin && length>0) { fin.get(); --length; }
81
82            return data;
83        }
84
85        template<typename T>
86        void printData(std::ostream& out, T* data, unsigned int numComponents) const
87        {
88            if (sizeof(T)==1)
89            {
90                for(unsigned int i=0; i<numComponents; ++i)
91                {
92                    if (data[i]>32) out<<data[i];
93                    else out<<".";
94                }
95            }
96            else
97            {
98                for(unsigned int i=0; i<numComponents; ++i)
99                {
100                    if (i==0) out<<data[i];
101                    else out<<", "<<data[i];
102                }
103            }
104        }
105
106        virtual const char* className() const { return "DICOM Image Reader/Writer"; }
107
108        bool isFileADicom(const std::string& filename) const
109        {
110            std::ifstream fin(filename.c_str(), std::ios::in | std::ios::binary);
111            if (!fin) return false;
112
113            char str[133];
114            str[128]=str[129]=str[130]=str[131]=0;
115            fin.getline(str, sizeof(str));
116            return (str[128]=='D' && str[129]=='I' && str[130]=='C' && str[131]=='M');
117        }
118
119        typedef std::vector<std::string> Files;
120        bool getDicomFilesInDirectory(const std::string& path, Files& files) const
121        {
122            osgDB::DirectoryContents contents = osgDB::getDirectoryContents(path);
123
124            std::sort(contents.begin(), contents.end());
125
126            for(osgDB::DirectoryContents::iterator itr = contents.begin();
127                itr != contents.end();
128                ++itr)
129            {
130                if ((*itr).empty()) continue;
131
132                if ((*itr)[0]=='.')
133                {
134                    info()<<"Ignoring tempory file "<<*itr<<std::endl;
135                    continue;
136                }
137
138                std::string localFile = path + "/" + *itr;
139
140                if (isFileADicom(localFile))
141                {
142                    files.push_back(localFile);
143                }
144            }
145
146            return !files.empty();
147        }
148
149        virtual ReadResult readObject(std::istream& fin,const osgDB::ReaderWriter::Options* options =NULL) const
150        {
151            return readImage(fin, options);
152        }
153
154        virtual ReadResult readObject(const std::string& file, const osgDB::ReaderWriter::Options* options =NULL) const
155        {
156            return readImage(file, options);
157        }
158
159        virtual ReadResult readNode(std::istream& fin,const osgDB::ReaderWriter::Options* options =NULL) const
160        {
161            return readImage(fin, options);
162        }
163
164        virtual ReadResult readNode(const std::string& file, const osgDB::ReaderWriter::Options* options =NULL) const
165        {
166            ReadResult result = readImage(file, options);
167            if (!result.validImage()) return result;
168
169            osg::ref_ptr<osgVolume::VolumeTile> tile = new osgVolume::VolumeTile;
170            tile->setVolumeTechnique(new osgVolume::RayTracedTechnique());
171
172            osg::ref_ptr<osgVolume::ImageLayer> layer= new osgVolume::ImageLayer(result.getImage());
173            layer->rescaleToZeroToOneRange();
174
175            osgVolume::SwitchProperty* sp = new osgVolume::SwitchProperty;
176            sp->setActiveProperty(0);
177
178            float alphaFunc = 0.1f;
179
180            osgVolume::AlphaFuncProperty* ap = new osgVolume::AlphaFuncProperty(alphaFunc);
181            osgVolume::SampleDensityProperty* sd = new osgVolume::SampleDensityProperty(0.005);
182            osgVolume::TransparencyProperty* tp = new osgVolume::TransparencyProperty(1.0);
183
184            {
185                // Standard
186                osgVolume::CompositeProperty* cp = new osgVolume::CompositeProperty;
187                cp->addProperty(ap);
188                cp->addProperty(sd);
189                cp->addProperty(tp);
190
191                sp->addProperty(cp);
192            }
193
194            {
195                // Light
196                osgVolume::CompositeProperty* cp = new osgVolume::CompositeProperty;
197                cp->addProperty(ap);
198                cp->addProperty(sd);
199                cp->addProperty(tp);
200                cp->addProperty(new osgVolume::LightingProperty);
201
202                sp->addProperty(cp);
203            }
204
205            {
206                // Isosurface
207                osgVolume::CompositeProperty* cp = new osgVolume::CompositeProperty;
208                cp->addProperty(sd);
209                cp->addProperty(tp);
210                cp->addProperty(new osgVolume::IsoSurfaceProperty(alphaFunc));
211
212                sp->addProperty(cp);
213            }
214
215            {
216                // MaximumIntensityProjection
217                osgVolume::CompositeProperty* cp = new osgVolume::CompositeProperty;
218                cp->addProperty(ap);
219                cp->addProperty(sd);
220                cp->addProperty(tp);
221                cp->addProperty(new osgVolume::MaximumIntensityProjectionProperty);
222
223                sp->addProperty(cp);
224            }
225
226            layer->addProperty(sp);
227
228            tile->setLayer(layer.get());
229
230            // get matrix providing size of texels (in mm)
231            osgVolume::ImageDetails* details = dynamic_cast<osgVolume::ImageDetails*>(result.getImage()->getUserData());
232            osg::RefMatrix* matrix = details ? details->getMatrix() : 0;
233
234            if (details)
235            {
236                layer->setTexelOffset(details->getTexelOffset());
237                layer->setTexelScale(details->getTexelScale());
238            }
239
240            if (matrix)
241            {
242                osgVolume::Locator* locator = new osgVolume::Locator(*matrix);
243
244                tile->setLocator(locator);
245                layer->setLocator(locator);
246
247                // result.getImage()->setUserData(0);
248
249                info()<<"Locator "<<*matrix<<std::endl;
250            }
251            else
252            {
253                info()<<"No Locator found on osg::Image"<<std::endl;
254            }
255
256            return tile.release();
257        }
258
259
260        virtual ReadResult readImage(std::istream& fin,const osgDB::ReaderWriter::Options*) const
261        {
262            return 0;
263        }
264
265
266#ifdef USE_ITK
267        virtual ReadResult readImage(const std::string& file, const osgDB::ReaderWriter::Options* options) const
268        {
269            std::string ext = osgDB::getLowerCaseFileExtension(file);
270            std::string fileName = file;
271            if (ext=="dicom")
272            {
273                fileName = osgDB::getNameLessExtension(file);
274            }
275
276            fileName = osgDB::findDataFile( fileName, options );
277            if (fileName.empty()) return ReadResult::FILE_NOT_FOUND;
278
279            Files files;
280
281            osgDB::FileType fileType = osgDB::fileType(fileName);
282            if (fileType==osgDB::DIRECTORY)
283            {
284                getDicomFilesInDirectory(fileName, files);
285            }
286            else if (isFileADicom(fileName))
287            {
288                files.push_back(fileName);
289            }
290
291            if (files.empty())
292            {
293                return ReadResult::FILE_NOT_FOUND;
294            }
295
296
297            typedef std::vector< osg::ref_ptr<osg::Image> > Images;
298            Images images;
299            for(Files::iterator itr = files.begin();
300                itr != files.end();
301                ++itr)
302            {
303                ReadResult result = readSingleITKImage(*itr, options);
304                if (result.success()) images.push_back(result.getImage());
305                else return result;
306            }
307
308            if (images.empty()) return ReadResult::ERROR_IN_READING_FILE;
309
310            if (images.size()==1) return images[0].get();
311
312
313            typedef std::map<float, osg::ref_ptr<osg::Image> > DistanceImageMap;
314            typedef std::map<osg::Vec3, DistanceImageMap> OrientationDistanceImageMap;
315            OrientationDistanceImageMap orientationDistanceImageMap;
316            for(Images::iterator itr = images.begin();
317                itr != images.end();
318                ++itr)
319            {
320                osg::Image* image = itr->get();
321                osgVolume::ImageDetails* details = dynamic_cast<osgVolume::ImageDetails*>(image->getUserData());
322                osg::RefMatrix* matrix = details ? details->getMatrix() : 0;
323                if (matrix)
324                {
325                    osg::Vec3 p0 = osg::Vec3(0.0, 0.0, 0.0) * (*matrix);
326                    osg::Vec3 p1 = osg::Vec3(0.0, 0.0, 1.0) * (*matrix);
327                    osg::Vec3 direction = p1-p0;
328                    direction.normalize();
329                    float distance = p0 * direction;
330                    info()<<" direction="<<direction<<"  distance = "<<distance<<std::endl;
331                    orientationDistanceImageMap[direction][distance] = image;
332                }
333            }
334
335            if (orientationDistanceImageMap.empty()) return ReadResult::ERROR_IN_READING_FILE;
336
337            DistanceImageMap& dim = orientationDistanceImageMap.begin()->second;
338
339            double totalThickness = dim.rbegin()->first - dim.begin()->first;
340
341            int width = 0;
342            int height = 0;
343            int depth = 0;
344            for(DistanceImageMap::iterator itr = dim.begin();
345                itr != dim.end();
346                ++itr)
347            {
348                osg::Image* image = itr->second.get();
349                if (image->s() > width) width = image->s();
350                if (image->t() > height) height = image->t();
351                depth += image->r();
352            }
353
354
355            osg::ref_ptr<osg::Image> image3D = new osg::Image;
356            image3D->allocateImage(width, height, depth, GL_LUMINANCE, GL_UNSIGNED_BYTE, 1);
357            int r = 0;
358            for(DistanceImageMap::iterator itr = dim.begin();
359                itr != dim.end();
360                ++itr)
361            {
362                osg::Image* image = itr->second.get();
363                osg::copyImage(image, 0,0,0, image->s(), image->t(), image->r(),
364                            image3D.get(), 0, 0, r,
365                            false);
366                r += image->r();
367            }
368
369            osg::Image* firstImage = dim.begin()->second.get();
370            osgVolume::ImageDetails* details = dynamic_cast<osgVolume::ImageDetails*>(firstImage->getUserData());
371            osg::RefMatrix* matrix = details ? details->getMatrix() : 0;
372            if (matrix)
373            {
374                osgVolume::ImageDetails* details3D = new osgVolume::ImageDetails(*details);
375                details3D->getMatrix()->preMult(osg::Matrix::scale(1.0,1.0,totalThickness));
376                image3D->setUserData(details3D);
377            }
378
379            return image3D.get();
380        }
381
382        virtual ReadResult readSingleITKImage(const std::string& fileName, const osgDB::ReaderWriter::Options* options) const
383        {
384
385            typedef unsigned short PixelType;
386            const unsigned int Dimension = 3;
387            typedef itk::Image< PixelType, Dimension > ImageType;
388
389
390            typedef itk::ImageFileReader< ImageType > ReaderType;
391            ReaderType::Pointer reader = ReaderType::New();
392            reader->SetFileName( fileName.c_str() );
393
394            typedef itk::GDCMImageIO           ImageIOType;
395            ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
396            reader->SetImageIO( gdcmImageIO );
397
398            try
399            {
400                reader->Update();
401            }
402            catch (itk::ExceptionObject & e)
403            {
404                std::cerr << "exception in file reader " << std::endl;
405                std::cerr << e.GetDescription() << std::endl;
406                std::cerr << e.GetLocation() << std::endl;
407                return ReadResult::ERROR_IN_READING_FILE;
408            }
409
410            ImageType::Pointer inputImage = reader->GetOutput();
411
412            ImageType::RegionType region = inputImage->GetBufferedRegion();
413            ImageType::SizeType size = region.GetSize();
414            ImageType::IndexType start = region.GetIndex();
415
416            //inputImage->GetSpacing();
417            //inputImage->GetOrigin();
418
419            unsigned int width = size[0];
420            unsigned int height = size[1];
421            unsigned int depth = size[2];
422
423            osg::RefMatrix* matrix = new osg::RefMatrix;
424
425            info()<<"width = "<<width<<" height = "<<height<<" depth = "<<depth<<std::endl;
426            for(unsigned int i=0; i<Dimension; ++i)
427            {
428                (*matrix)(i,i) = inputImage->GetSpacing()[i];
429                (*matrix)(3,i) = inputImage->GetOrigin()[i];
430            }
431
432            osg::Image* image = new osg::Image;
433            image->allocateImage(width, height, depth, GL_LUMINANCE, GL_UNSIGNED_BYTE, 1);
434
435            unsigned char* data = image->data();
436            typedef itk::ImageRegionConstIterator< ImageType > IteratorType;
437            IteratorType it(inputImage, region);
438
439            it.GoToBegin();
440            while (!it.IsAtEnd())
441            {
442                *data = it.Get();
443                ++data;
444                ++it;
445            }
446
447            osgVolume::ImageDetails* details = new osgVolume::ImageDetails;
448            details->setMatrix(matrix);
449
450            image->setUserData(details);
451
452            matrix->preMult(osg::Matrix::scale(double(image->s()), double(image->t()), double(image->r())));
453
454            return image;
455        }
456#endif
457
458#ifdef USE_DCMTK
459
460        void convertPixelTypes(const DiPixel* pixelData,
461                            EP_Representation& pixelRep, int& numPlanes,
462                            GLenum& dataType, GLenum& pixelFormat, unsigned int& pixelSize) const
463        {
464            dataType = GL_UNSIGNED_BYTE;
465            pixelRep = pixelData->getRepresentation();
466            switch(pixelRep)
467            {
468                case(EPR_Uint8):
469                    dataType = GL_UNSIGNED_BYTE;
470                    pixelSize = 1;
471                    break;
472                case(EPR_Sint8):
473                    dataType = GL_BYTE;
474                    pixelSize = 1;
475                    break;
476                case(EPR_Uint16):
477                    dataType = GL_UNSIGNED_SHORT;
478                    pixelSize = 2;
479                    break;
480                case(EPR_Sint16):
481                    dataType = GL_SHORT;
482                    pixelSize = 2;
483                    break;
484                case(EPR_Uint32):
485                    dataType = GL_UNSIGNED_INT;
486                    pixelSize = 4;
487                    break;
488                case(EPR_Sint32):
489                    dataType = GL_INT;
490                    pixelSize = 4;
491                    break;
492                default:
493                    dataType = 0;
494                    pixelSize = 1;
495                    break;
496            }
497
498            pixelFormat = GL_INTENSITY;
499            numPlanes = pixelData->getPlanes();
500            switch(numPlanes)
501            {
502                case(1):
503                    pixelFormat = GL_LUMINANCE;
504                    break;
505                case(2):
506                    pixelFormat = GL_LUMINANCE_ALPHA;
507                    pixelSize *= 2;
508                    break;
509                case(3):
510                    pixelFormat = GL_RGB;
511                    pixelSize *= 3;
512                    break;
513                case(4):
514                    pixelFormat = GL_RGBA;
515                    pixelSize *= 4;
516                    break;
517            }
518        }
519
520
521        struct SeriesIdentifier
522        {
523            std::string SeriesInstanceUID;
524            std::string SeriesDescription;
525            double Orientation[6];
526
527            SeriesIdentifier()
528            {
529                Orientation[0] = 1.0;
530                Orientation[1] = 0.0;
531                Orientation[2] = 0.0;
532                Orientation[3] = 0.0;
533                Orientation[4] = 1.0;
534                Orientation[5] = 0.0;
535            }
536
537            void set(DcmDataset* dataset)
538            {
539                OFString seriesInstanceUIDStr;
540                if (dataset->findAndGetOFString(DCM_SeriesInstanceUID, seriesInstanceUIDStr).good())
541                {
542                    SeriesInstanceUID = seriesInstanceUIDStr.c_str();
543                }
544
545                OFString seriesDescriptionStr;
546                if (dataset->findAndGetOFString(DCM_SeriesDescription, seriesDescriptionStr).good())
547                {
548                    SeriesDescription = seriesDescriptionStr.c_str();
549                }
550
551                for(int i=0; i<6; ++i)
552                {
553                    double value = 0.0;
554                    if (dataset->findAndGetFloat64(DCM_ImageOrientationPatient, value,i).good())
555                    {
556                        Orientation[i] = value;
557                    }
558                }
559            }
560
561            bool operator == (const SeriesIdentifier& rhs) const
562            {
563                if (SeriesInstanceUID != rhs.SeriesInstanceUID) return false;
564                if (SeriesDescription != rhs.SeriesDescription) return false;
565
566                for(unsigned int i=0; i<6; ++i)
567                {
568                    if (Orientation[i] != rhs.Orientation[i]) return false;
569                }
570                return true;
571            }
572
573            bool operator < (const SeriesIdentifier& rhs) const
574            {
575                if (SeriesInstanceUID < rhs.SeriesInstanceUID) return true;
576                if (rhs.SeriesInstanceUID < SeriesInstanceUID) return false;
577
578                if (SeriesDescription < rhs.SeriesDescription) return true;
579                if (rhs.SeriesDescription < SeriesDescription) return false;
580
581                for(unsigned int i=0; i<6; ++i)
582                {
583                    if (Orientation[i] >= rhs.Orientation[i]) return false;
584                }
585                return true;
586            }
587        };
588
589        virtual ReadResult readImage(const std::string& file, const osgDB::ReaderWriter::Options* options) const
590        {
591            std::string ext = osgDB::getLowerCaseFileExtension(file);
592            std::string fileName = file;
593            if (ext=="dicom")
594            {
595                fileName = osgDB::getNameLessExtension(file);
596            }
597
598            fileName = osgDB::findDataFile( fileName, options );
599            if (fileName.empty()) return ReadResult::FILE_NOT_FOUND;
600
601            Files files;
602
603            osgDB::FileType fileType = osgDB::fileType(fileName);
604            if (fileType==osgDB::DIRECTORY)
605            {
606                getDicomFilesInDirectory(fileName, files);
607            }
608            else if (isFileADicom(fileName))
609            {
610                files.push_back(fileName);
611            }
612            else
613            {
614                return ReadResult::FILE_NOT_HANDLED;
615            }
616
617            if (files.empty())
618            {
619                return ReadResult::FILE_NOT_FOUND;
620            }
621
622            info()<<"Reading DICOM file "<<file<<" using DCMTK"<<std::endl;
623
624
625            osg::ref_ptr<osgVolume::ImageDetails> details = new osgVolume::ImageDetails;
626            details->setMatrix(new osg::RefMatrix);
627
628            unsigned int imageNum = 0;
629            EP_Representation pixelRep = EPR_Uint8;
630            int numPlanes = 0;
631            GLenum pixelFormat = 0;
632            GLenum dataType = 0;
633            unsigned int pixelSize = 0;
634
635            typedef std::list<FileInfo> FileInfoList;
636            FileInfoList fileInfoList;
637
638            SeriesIdentifier seriesIdentifier;
639
640            typedef std::map<double, FileInfo> DistanceFileInfoMap;
641            typedef std::map<SeriesIdentifier, DistanceFileInfoMap> SeriesFileInfoMap;
642            SeriesFileInfoMap seriesFileInfoMap;
643
644            typedef std::map<std::string, ReadResult> ErrorMap;
645            ErrorMap errorMap;
646
647            for(Files::iterator itr = files.begin();
648                itr != files.end();
649                ++itr)
650            {
651                DcmFileFormat fileformat;
652                const std::string& dicom_filename = *itr;
653                OFCondition status = fileformat.loadFile(dicom_filename.c_str());
654                if(!status.good())
655                {
656                    errorMap[dicom_filename] = ReadResult::ERROR_IN_READING_FILE;
657                    continue;
658                }
659
660                FileInfo fileInfo;
661                fileInfo.filename = *itr;
662
663                SeriesIdentifier seriesIdentifier;
664                seriesIdentifier.set(fileformat.getDataset());
665
666                // code for reading the intercept and scale that is required to convert to Hounsfield units.
667                bool rescaling = false;
668                double rescaleIntercept = 0.0;
669                double rescaleSlope = 1.0;
670                const char *classUID = NULL;
671                if (fileformat.getDataset()->findAndGetString(DCM_SOPClassUID, classUID).good())
672                {
673                    info()<<" classUID = "<<classUID<<std::endl;
674                    if (0 == strcmp(classUID, UID_CTImageStorage))
675                    {
676                        info()<<" is a UID_CTImageStorage "<<std::endl;
677                    }
678
679                }
680
681
682
683                rescaling = fileformat.getDataset()->findAndGetFloat64(DCM_RescaleIntercept, rescaleIntercept).good();
684                rescaling &= fileformat.getDataset()->findAndGetFloat64(DCM_RescaleSlope, rescaleSlope).good();
685                if (rescaling)
686                {
687                    fileInfo.rescaleIntercept = rescaleIntercept;
688                    fileInfo.rescaleSlope = rescaleSlope;
689                    info()<<" rescaleIntercept = "<<rescaleIntercept<<std::endl;
690                    info()<<" rescaleSlope = "<<rescaleSlope<<std::endl;
691                }
692
693
694                double value = 0.0;
695                if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,0).good())
696                {
697                    fileInfo.pixelSize_x = value;
698                }
699
700                if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,1).good())
701                {
702                    fileInfo.pixelSize_y = value;
703                }
704
705                if (fileformat.getDataset()->findAndGetFloat64(DCM_SpacingBetweenSlices, value,0).good())
706                {
707                    info()<<"DCM_SpacingBetweenSlices = "<<value<<std::endl;
708                    fileInfo.sliceThickness = value;
709                }
710
711
712                // Get slice thickness
713                if (fileformat.getDataset()->findAndGetFloat64(DCM_SliceThickness, value).good())
714                {
715                    info()<<"DCM_SliceThickness = "<<value<<std::endl;
716                    fileInfo.sliceThickness = value;
717                }
718
719                info()<<"tagExistsWithValue(DCM_NumberOfFrames)="<<fileformat.getDataset()->tagExistsWithValue(DCM_NumberOfFrames)<<std::endl;
720                info()<<"tagExistsWithValue(DCM_NumberOfSlices)="<<fileformat.getDataset()->tagExistsWithValue(DCM_NumberOfSlices)<<std::endl;
721
722                Uint16 numOfSlices = 1;
723                Uint32 numFrames = 1;
724                if (fileformat.getDataset()->findAndGetUint32(DCM_NumberOfFrames, numFrames).good())
725                {
726                    fileInfo.numSlices = numFrames;
727                    info()<<"Read number of frames = "<<numFrames<<std::endl;
728                }
729
730
731                OFString numFramesStr;
732                if (fileformat.getDataset()->findAndGetOFString(DCM_NumberOfFrames, numFramesStr).good())
733                {
734                    fileInfo.numSlices = atoi(numFramesStr.c_str());
735                    info()<<"Read number of frames = "<<numFramesStr<<std::endl;
736                }
737
738                if (fileformat.getDataset()->findAndGetUint16(DCM_NumberOfFrames, numOfSlices).good())
739                {
740                    fileInfo.numSlices = numOfSlices;
741                    info()<<"Read number of frames = "<<numOfSlices<<std::endl;
742                }
743
744                if (fileformat.getDataset()->findAndGetUint16(DCM_NumberOfSlices, numOfSlices).good())
745                {
746                    //fileInfo.numSlices = numOfSlices;
747                    info()<<"Read number of slices = "<<numOfSlices<<std::endl;
748                }
749
750
751                // patient position
752                double imagePositionPatient[3] = {0.0, 0.0, 0.0};
753                for(int i=0; i<3; ++i)
754                {
755                    if (fileformat.getDataset()->findAndGetFloat64(DCM_ImagePositionPatient, imagePositionPatient[i],i).good())
756                    {
757                        info()<<"Read DCM_ImagePositionPatient["<<i<<"], "<<imagePositionPatient[i]<<std::endl;
758                    }
759                    else
760                    {
761                        info()<<"Have not read DCM_ImagePositionPatient["<<i<<"]"<<std::endl;
762                    }
763                }
764                fileInfo.position.set(imagePositionPatient[0],imagePositionPatient[1],imagePositionPatient[2]);
765
766                double imageOrientationPatient[6] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0 };
767                for(int i=0; i<6; ++i)
768                {
769                    double value = 0.0;
770                    if (fileformat.getDataset()->findAndGetFloat64(DCM_ImageOrientationPatient, value,i).good())
771                    {
772                        imageOrientationPatient[i] = value;
773                        info()<<"Read imageOrientationPatient["<<i<<"], "<<imageOrientationPatient[i]<<std::endl;
774                    }
775                    else
776                    {
777                        info()<<"Have not read imageOrientationPatient["<<i<<"]"<<std::endl;
778                    }
779                }
780
781                fileInfo.dirX.set(imageOrientationPatient[0],imageOrientationPatient[1],imageOrientationPatient[2]);
782                fileInfo.dirY.set(imageOrientationPatient[3],imageOrientationPatient[4],imageOrientationPatient[5]);
783                fileInfo.dirZ = fileInfo.dirX ^ fileInfo.dirY;
784                fileInfo.dirZ.normalize();
785                fileInfo.distance = fileInfo.dirZ * fileInfo.position;
786
787
788                info()<<"pixelSize_x="<<fileInfo.pixelSize_x<<std::endl;
789                info()<<"pixelSize_y="<<fileInfo.pixelSize_x<<std::endl;
790
791                info()<<"dirX.length() = "<<fileInfo.dirX.length()<<std::endl;
792                info()<<"dirY.length() = "<<fileInfo.dirY.length()<<std::endl;
793                info()<<"dot_product = "<<fileInfo.dirX*fileInfo.dirY<<std::endl;
794                info()<<"dirX = "<<fileInfo.dirX<<std::endl;
795                info()<<"dirY = "<<fileInfo.dirY<<std::endl;
796                info()<<"dirZ = "<<fileInfo.dirZ<<std::endl;
797                info()<<"pos = "<<fileInfo.position<<std::endl;
798                info()<<"dist = "<<fileInfo.distance<<std::endl;
799                info()<<std::endl;
800
801                (seriesFileInfoMap[seriesIdentifier])[fileInfo.distance] = fileInfo;
802
803            }
804
805            if (seriesFileInfoMap.empty()) return 0;
806
807            for(SeriesFileInfoMap::iterator itr = seriesFileInfoMap.begin();
808                itr != seriesFileInfoMap.end();
809                ++itr)
810            {
811                notice()<<"Description = "<<itr->first.SeriesDescription<<", Orientation = "<<itr->first.Orientation<<std::endl;
812
813                unsigned int totalNumSlices = 0;
814
815                DistanceFileInfoMap& dfim = itr->second;
816                for(DistanceFileInfoMap::iterator ditr = dfim.begin();
817                    ditr != dfim.end();
818                    ++ditr)
819                {
820                    FileInfo& fileInfo = ditr->second;
821                    totalNumSlices += fileInfo.numSlices;
822                    info()<<"   d = "<<fileInfo.distance<<" "<<fileInfo.filename<<" fileInfo.numSlices="<<fileInfo.numSlices<<std::endl;
823                }
824
825                if (dfim.empty()) continue;
826
827                osg::ref_ptr<osg::Image> image;
828
829                double totalDistance = 0.0;
830                if (dfim.size()>1)
831                {
832                    totalDistance = fabs(dfim.rbegin()->first - dfim.begin()->first);
833                }
834                else
835                {
836                    totalDistance = dfim.begin()->second.sliceThickness * double(dfim.begin()->second.numSlices);
837                }
838
839                info()<<"Total Number slices "<<totalNumSlices<<std::endl;
840                info()<<"Total Distance including ends "<<totalDistance<<std::endl;
841
842                double averageThickness = totalNumSlices<=1 ? 1.0 : totalDistance / double(totalNumSlices-1);
843
844                info()<<"Average thickness "<<averageThickness<<std::endl;
845
846                for(DistanceFileInfoMap::iterator ditr = dfim.begin();
847                    ditr != dfim.end();
848                    ++ditr)
849                {
850                    FileInfo& fileInfo = ditr->second;
851
852                    std::auto_ptr<DicomImage> dcmImage(new DicomImage(fileInfo.filename.c_str()));
853                    if (dcmImage.get())
854                    {
855                        if (dcmImage->getStatus()==EIS_Normal)
856                        {
857
858                            EP_Representation curr_pixelRep;
859                            int curr_numPlanes;
860                            GLenum curr_pixelFormat;
861                            GLenum curr_dataType;
862                            unsigned int curr_pixelSize;
863
864                            // get the pixel data
865                            const DiPixel* pixelData = dcmImage->getInterData();
866                            if(!pixelData)
867                            {
868                                warning()<<"Error: no data in DicomImage object."<<std::endl;
869                                return ReadResult::ERROR_IN_READING_FILE;
870                            }
871
872                            // create the new image
873                            convertPixelTypes(pixelData,
874                                            curr_pixelRep, curr_numPlanes,
875                                            curr_dataType, curr_pixelFormat, curr_pixelSize);
876
877                            // dcmImage->getFrameCount()
878
879                            osg::ref_ptr<osg::Image> imageAdapter = new osg::Image;
880
881                            if (dcmImage->isMonochrome())
882                            {
883                                imageAdapter->setImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(),
884                                                    curr_pixelFormat,
885                                                    curr_pixelFormat,
886                                                    curr_dataType,
887                                                    (unsigned char*)(pixelData->getData()),
888                                                    osg::Image::NO_DELETE);
889
890                            }
891                            else
892                            {
893                                imageAdapter->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(),
894                                                curr_pixelFormat, curr_dataType);
895
896                                void* data = imageAdapter->data(0,0,0);
897                                unsigned long size = dcmImage->createWindowsDIB( data,
898                                                                                imageAdapter->getTotalDataSize(),
899                                                                                0,
900                                                                                imageAdapter->getPixelSizeInBits(),
901                                                                                0,
902                                                                                0);
903
904                                if (size==0)
905                                {
906                                    info()<<"  dcmImage->createWindowsDIB() failed to create required imagery."<<std::endl;
907                                    continue;
908                                }
909                            }
910
911                            if (!image)
912                            {
913                                pixelRep = curr_pixelRep;
914                                numPlanes = curr_numPlanes;
915                                dataType = curr_dataType;
916                                pixelFormat = curr_pixelFormat;
917                                pixelSize = curr_pixelSize;
918
919                                osg::RefMatrix* matrix = details->getMatrix();
920
921                                (*matrix)(0,0) = fileInfo.dirX.x();
922                                (*matrix)(1,0) = fileInfo.dirX.y();
923                                (*matrix)(2,0) = fileInfo.dirX.z();
924
925                                (*matrix)(0,1) = fileInfo.dirY.x();
926                                (*matrix)(1,1) = fileInfo.dirY.y();
927                                (*matrix)(2,1) = fileInfo.dirY.z();
928
929                                (*matrix)(0,2) = fileInfo.dirZ.x();
930                                (*matrix)(1,2) = fileInfo.dirZ.y();
931                                (*matrix)(2,2) = fileInfo.dirZ.z();
932
933                                matrix->preMultScale(osg::Vec3d(
934                                    fileInfo.pixelSize_x * dcmImage->getWidth(),
935                                    fileInfo.pixelSize_y * dcmImage->getHeight(),
936                                    averageThickness * totalNumSlices));
937
938                                (*matrix)(3,0) = fileInfo.position.x();
939                                (*matrix)(3,1) = fileInfo.position.y();
940                                (*matrix)(3,2) = fileInfo.position.z();
941
942                                (*matrix)(3,3) = 1.0;
943
944                                // note from Robert Osfield, testing various dicom files I have found that the rescaleIntercept
945                                // for CT data doesn't look to be applicable as an straight value offset, so we'll ignore for now.
946                                // details->setTexelOffset(fileInfo.rescaleIntercept);
947                                double s = fileInfo.rescaleSlope;
948                                switch(dataType)
949                                {
950                                    case(GL_BYTE): s *= 128.0; break;
951                                    case(GL_UNSIGNED_BYTE): s *= 255.0; break;
952                                    case(GL_SHORT): s *= 32768.0; break;
953                                    case(GL_UNSIGNED_SHORT): s *= 65535.0; break;
954                                    case(GL_INT): s *= 2147483648.0; break;
955                                    case(GL_UNSIGNED_INT): s *= 4294967295.0; break;
956                                    default: break;
957                                }
958
959                                details->setTexelScale(osg::Vec4(s,s,s,s));
960
961                                image = new osg::Image;
962                                image->setUserData(details.get());
963                                image->setFileName(fileName.c_str());
964                                image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices,
965                                                    pixelFormat, dataType);
966
967
968                                //matrix->preMult(osg::Matrix::scale(double(image->s()), double(image->t()), double(image->r())));
969
970                                info()<<"Image dimensions = "<<image->s()<<", "<<image->t()<<", "<<image->r()<<" pixelFormat=0x"<<std::hex<<pixelFormat<<" dataType=0x"<<std::hex<<dataType<<std::dec<<std::endl;
971                            }
972                            else if (pixelData->getPlanes()>numPlanes ||
973                                    pixelData->getRepresentation()>pixelRep)
974                            {
975                                info()<<"Need to reallocated "<<image->s()<<", "<<image->t()<<", "<<image->r()<<std::endl;
976
977                                // record the previous image settings to use when we copy back the content.
978                                osg::ref_ptr<osg::Image> previous_image = image;
979
980                                // create the new image
981                                convertPixelTypes(pixelData,
982                                                pixelRep, numPlanes,
983                                                dataType, pixelFormat, pixelSize);
984
985                                image = new osg::Image;
986                                image->setUserData(previous_image->getUserData());
987                                image->setFileName(fileName.c_str());
988                                image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices,
989                                                    pixelFormat, dataType);
990                                osg::copyImage(previous_image.get(), 0,0,0, previous_image->s(), previous_image->t(), imageNum,
991                                            image.get(), 0, 0, 0,
992                                            false);
993                            }
994
995                            info()<<"copyImage(, fileInfo.distance"<<fileInfo.distance<<", imageNum="<<imageNum<<std::endl;
996
997                            osg::copyImage(imageAdapter.get(), 0,0,0, imageAdapter->s(), imageAdapter->t(), imageAdapter->r(),
998                                        image.get(), 0, 0, imageNum,
999                                        false);
1000
1001                            imageNum += dcmImage->getFrameCount();
1002                        }
1003                        else
1004                        {
1005                            warning()<<"Error in reading dicom file "<<fileInfo.filename<<", error = "<<DicomImage::getString(dcmImage->getStatus())<<std::endl;
1006                            info()<<"    dcmImage->getPhotometricInterpretation()="<<DicomImage::getString(dcmImage->getPhotometricInterpretation())<<std::endl;
1007                            info()<<"    dcmImage->width="<<dcmImage->getWidth()<<", height="<<dcmImage->getHeight()<<" FrameCount="<< dcmImage->getFrameCount()<<std::endl;
1008                        }
1009                    }
1010                }
1011
1012                info()<<"Image matrix = "<<*(details->getMatrix())<<std::endl;
1013
1014                return image.get();
1015            }
1016
1017            if (!errorMap.empty())
1018            {
1019                for(ErrorMap::iterator itr = errorMap.begin();
1020                    itr != errorMap.end();
1021                    ++itr)
1022                {
1023                    warning()<<"Error in reading file "<<itr->first<<std::endl;
1024                }
1025            }
1026
1027            return ReadResult::ERROR_IN_READING_FILE;
1028        }
1029#endif
1030
1031        struct FileInfo
1032        {
1033            FileInfo():
1034                rescaleIntercept(0.0),
1035                rescaleSlope(1.0),
1036                numX(0),
1037                numY(0),
1038                numSlices(1),
1039                pixelSize_x(0.0),
1040                pixelSize_y(0.0),
1041                sliceThickness(0.0),
1042                distance(0.0),
1043                position(0.0,0.0,0.0),
1044                dirX(1.0,0.0,0.0),
1045                dirY(0.0,1.0,0.0),
1046                dirZ(0.0,0.0,1.0) {}
1047
1048            FileInfo(const FileInfo& rhs):
1049                filename(rhs.filename),
1050                rescaleIntercept(rhs.rescaleIntercept),
1051                rescaleSlope(rhs.rescaleSlope),
1052                numX(rhs.numX),
1053                numY(rhs.numY),
1054                numSlices(rhs.numSlices),
1055                pixelSize_x(rhs.pixelSize_x),
1056                pixelSize_y(rhs.pixelSize_y),
1057                sliceThickness(rhs.sliceThickness),
1058                distance(distance),
1059                position(rhs.position),
1060                dirX(rhs.dirX),
1061                dirY(rhs.dirY),
1062                dirZ(rhs.dirZ) {}
1063
1064            FileInfo& operator = (const FileInfo& rhs)
1065            {
1066                if (&rhs == this) return *this;
1067
1068                filename = rhs.filename;
1069                rescaleIntercept = rhs.rescaleIntercept;
1070                rescaleSlope = rhs.rescaleSlope;
1071                numX = rhs.numX;
1072                numY = rhs.numY;
1073                pixelSize_x = rhs.pixelSize_x;
1074                pixelSize_y = rhs.pixelSize_y;
1075                sliceThickness = rhs.sliceThickness;
1076                numSlices = rhs.numSlices;
1077                distance = rhs.distance;
1078                position = rhs.position;
1079                dirX = rhs.dirX;
1080                dirY = rhs.dirY;
1081                dirZ = rhs.dirZ;
1082
1083                return *this;
1084            }
1085
1086            std::string     filename;
1087            double          rescaleIntercept;
1088            double          rescaleSlope;
1089            unsigned int    numX;
1090            unsigned int    numY;
1091            unsigned int    numSlices;
1092            double          pixelSize_x;
1093            double          pixelSize_y;
1094            double          sliceThickness;
1095            double          distance;
1096            osg::Vec3d      position;
1097            osg::Vec3d      dirX;
1098            osg::Vec3d      dirY;
1099            osg::Vec3d      dirZ;
1100        };
1101
1102};
1103
1104// now register with Registry to instantiate the above
1105// reader/writer.
1106REGISTER_OSGPLUGIN(dicom, ReaderWriterDICOM)
Note: See TracBrowser for help on using the browser.