diff --git a/src/osgPlugins/dicom/ReaderWriterDICOM.cpp b/src/osgPlugins/dicom/ReaderWriterDICOM.cpp index 500cb7ec8..8ecb8c1c0 100644 --- a/src/osgPlugins/dicom/ReaderWriterDICOM.cpp +++ b/src/osgPlugins/dicom/ReaderWriterDICOM.cpp @@ -559,8 +559,6 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter osg::ref_ptr details = new osgVolume::ImageDetails; details->setMatrix(new osg::RefMatrix); - osg::ref_ptr image; - unsigned int imageNum = 0; EP_Representation pixelRep = EPR_Uint8; int numPlanes = 0; @@ -575,8 +573,6 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter typedef std::map OrientationFileInfoMap; OrientationFileInfoMap orientationFileInfoMap; - unsigned int totalNumSlices = 0; - typedef std::map ErrorMap; ErrorMap errorMap; @@ -596,13 +592,6 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter FileInfo fileInfo; fileInfo.filename = *itr; - double pixelSize_y = 1.0; - double pixelSize_x = 1.0; - double sliceThickness = 1.0; - double imagePositionPatient[3] = {0, 0, 0}; - double imageOrientationPatient[6] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0 }; - Uint16 numOfSlices = 1; - // code for reading the intercept and scale that is required to convert to Hounsfield units. bool rescaling = false; double rescaleIntercept = 0.0; @@ -632,28 +621,33 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter double value = 0.0; if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,0).good()) { - pixelSize_y = value; - fileInfo.matrix(1,1) = pixelSize_y; + fileInfo.pixelSize_x = value; } if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,1).good()) { - pixelSize_x = value; - fileInfo.matrix(0,0) = pixelSize_x; + fileInfo.pixelSize_y = value; } + if (fileformat.getDataset()->findAndGetFloat64(DCM_SpacingBetweenSlices, value,0).good()) + { + info()<<"DCM_SpacingBetweenSlices = "<findAndGetFloat64(DCM_SliceThickness, value).good()) { - sliceThickness = value; - info()<<"sliceThickness = "<second; - if (dfim.empty()) return 0; + osg::ref_ptr image; - double totalDistance = 0.0; - if (dfim.size()>1) - { - totalDistance = fabs(dfim.rbegin()->first - dfim.begin()->first); - } - else - { - totalDistance = dfim.begin()->second.sliceThickness * double(dfim.begin()->second.numSlices); - } - - info()<<"Total Distance including ends "<second; - - std::auto_ptr dcmImage(new DicomImage(fileInfo.filename.c_str())); - if (dcmImage.get()) + double totalDistance = 0.0; + if (dfim.size()>1) { - if (dcmImage->getStatus()==EIS_Normal) + totalDistance = fabs(dfim.rbegin()->first - dfim.begin()->first); + } + else + { + totalDistance = dfim.begin()->second.sliceThickness * double(dfim.begin()->second.numSlices); + } + + info()<<"Total Number slices "<second; + + std::auto_ptr dcmImage(new DicomImage(fileInfo.filename.c_str())); + if (dcmImage.get()) { - - EP_Representation curr_pixelRep; - int curr_numPlanes; - GLenum curr_pixelFormat; - GLenum curr_dataType; - unsigned int curr_pixelSize; - - // get the pixel data - const DiPixel* pixelData = dcmImage->getInterData(); - if(!pixelData) + if (dcmImage->getStatus()==EIS_Normal) { - warning()<<"Error: no data in DicomImage object."<getFrameCount() - - osg::ref_ptr imageAdapter = new osg::Image; - - if (dcmImage->isMonochrome()) - { - imageAdapter->setImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(), - curr_pixelFormat, - curr_pixelFormat, - curr_dataType, - (unsigned char*)(pixelData->getData()), - osg::Image::NO_DELETE); - - } - else - { - imageAdapter->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(), - curr_pixelFormat, curr_dataType); - - void* data = imageAdapter->data(0,0,0); - unsigned long size = dcmImage->createWindowsDIB( data, - imageAdapter->getTotalDataSize(), - 0, - imageAdapter->getPixelSizeInBits(), - 0, - 0); - - if (size==0) + // get the pixel data + const DiPixel* pixelData = dcmImage->getInterData(); + if(!pixelData) { - info()<<" dcmImage->createWindowsDIB() failed to create required imagery."<getMatrix(); - - (*matrix)(0,0) = fileInfo.matrix(0,0); - (*matrix)(1,0) = fileInfo.matrix(1,0); - (*matrix)(2,0) = fileInfo.matrix(2,0); - (*matrix)(0,1) = fileInfo.matrix(0,1); - (*matrix)(1,1) = fileInfo.matrix(1,1); - (*matrix)(2,1) = fileInfo.matrix(2,1); - (*matrix)(0,2) = fileInfo.matrix(0,2) * averageThickness; - (*matrix)(1,2) = fileInfo.matrix(1,2) * averageThickness; - (*matrix)(2,2) = fileInfo.matrix(2,2) * averageThickness; - - // note from Robert Osfield, testing various dicom files I have found that the rescaleIntercept - // for CT data doesn't look to be applicable as an straight value offset, so we'll ignore for now. - // details->setTexelOffset(fileInfo.rescaleIntercept); - double s = fileInfo.rescaleSlope; - switch(dataType) - { - case(GL_BYTE): s *= 128.0; break; - case(GL_UNSIGNED_BYTE): s *= 255.0; break; - case(GL_SHORT): s *= 32768.0; break; - case(GL_UNSIGNED_SHORT): s *= 65535.0; break; - case(GL_INT): s *= 2147483648.0; break; - case(GL_UNSIGNED_INT): s *= 4294967295.0; break; - default: break; - } - - details->setTexelScale(osg::Vec4(s,s,s,s)); - - image = new osg::Image; - image->setUserData(details.get()); - image->setFileName(fileName.c_str()); - image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices, - pixelFormat, dataType); - - - matrix->preMult(osg::Matrix::scale(double(image->s()), double(image->t()), double(image->r()))); - - info()<<"Image dimensions = "<s()<<", "<t()<<", "<r()<<" pixelFormat=0x"<getPlanes()>numPlanes || - pixelData->getRepresentation()>pixelRep) - { - info()<<"Need to reallocated "<s()<<", "<t()<<", "<r()< previous_image = image; // create the new image convertPixelTypes(pixelData, - pixelRep, numPlanes, - dataType, pixelFormat, pixelSize); + curr_pixelRep, curr_numPlanes, + curr_dataType, curr_pixelFormat, curr_pixelSize); - image = new osg::Image; - image->setUserData(previous_image->getUserData()); - image->setFileName(fileName.c_str()); - image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices, - pixelFormat, dataType); + // dcmImage->getFrameCount() - osg::copyImage(previous_image.get(), 0,0,0, previous_image->s(), previous_image->t(), imageNum, - image.get(), 0, 0, 0, + osg::ref_ptr imageAdapter = new osg::Image; + + if (dcmImage->isMonochrome()) + { + imageAdapter->setImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(), + curr_pixelFormat, + curr_pixelFormat, + curr_dataType, + (unsigned char*)(pixelData->getData()), + osg::Image::NO_DELETE); + + } + else + { + imageAdapter->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(), + curr_pixelFormat, curr_dataType); + + void* data = imageAdapter->data(0,0,0); + unsigned long size = dcmImage->createWindowsDIB( data, + imageAdapter->getTotalDataSize(), + 0, + imageAdapter->getPixelSizeInBits(), + 0, + 0); + + if (size==0) + { + info()<<" dcmImage->createWindowsDIB() failed to create required imagery."<getMatrix(); + + (*matrix)(0,0) = fileInfo.dirX.x(); + (*matrix)(1,0) = fileInfo.dirX.y(); + (*matrix)(2,0) = fileInfo.dirX.z(); + + (*matrix)(0,1) = fileInfo.dirY.x(); + (*matrix)(1,1) = fileInfo.dirY.y(); + (*matrix)(2,1) = fileInfo.dirY.z(); + + (*matrix)(0,2) = fileInfo.dirZ.x(); + (*matrix)(1,2) = fileInfo.dirZ.y(); + (*matrix)(2,2) = fileInfo.dirZ.z(); + + matrix->preMultScale(osg::Vec3d( + fileInfo.pixelSize_x * dcmImage->getWidth(), + fileInfo.pixelSize_y * dcmImage->getHeight(), + averageThickness * totalNumSlices)); + + (*matrix)(3,0) = fileInfo.position.x(); + (*matrix)(3,1) = fileInfo.position.y(); + (*matrix)(3,2) = fileInfo.position.z(); + + (*matrix)(3,3) = 1.0; + + // note from Robert Osfield, testing various dicom files I have found that the rescaleIntercept + // for CT data doesn't look to be applicable as an straight value offset, so we'll ignore for now. + // details->setTexelOffset(fileInfo.rescaleIntercept); + double s = fileInfo.rescaleSlope; + switch(dataType) + { + case(GL_BYTE): s *= 128.0; break; + case(GL_UNSIGNED_BYTE): s *= 255.0; break; + case(GL_SHORT): s *= 32768.0; break; + case(GL_UNSIGNED_SHORT): s *= 65535.0; break; + case(GL_INT): s *= 2147483648.0; break; + case(GL_UNSIGNED_INT): s *= 4294967295.0; break; + default: break; + } + + details->setTexelScale(osg::Vec4(s,s,s,s)); + + image = new osg::Image; + image->setUserData(details.get()); + image->setFileName(fileName.c_str()); + image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices, + pixelFormat, dataType); + + + //matrix->preMult(osg::Matrix::scale(double(image->s()), double(image->t()), double(image->r()))); + + info()<<"Image dimensions = "<s()<<", "<t()<<", "<r()<<" pixelFormat=0x"<getPlanes()>numPlanes || + pixelData->getRepresentation()>pixelRep) + { + info()<<"Need to reallocated "<s()<<", "<t()<<", "<r()< previous_image = image; + + // create the new image + convertPixelTypes(pixelData, + pixelRep, numPlanes, + dataType, pixelFormat, pixelSize); + + image = new osg::Image; + image->setUserData(previous_image->getUserData()); + image->setFileName(fileName.c_str()); + image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices, + pixelFormat, dataType); + osg::copyImage(previous_image.get(), 0,0,0, previous_image->s(), previous_image->t(), imageNum, + image.get(), 0, 0, 0, + false); + } + + info()<<"copyImage(, fileInfo.distance"<getPhotometricInterpretation()="<getPhotometricInterpretation())<width="<getWidth()<<", height="<getHeight()<<" FrameCount="<< dcmImage->getFrameCount()<s(), imageAdapter->t(), imageAdapter->r(), - image.get(), 0, 0, imageNum, - false); - - imageNum += dcmImage->getFrameCount(); - } - else - { - warning()<<"Error in reading dicom file "<getPhotometricInterpretation()="<getPhotometricInterpretation())<width="<getWidth()<<", height="<getHeight()<<" FrameCount="<< dcmImage->getFrameCount()<