From 3b4184295e14e0d6758579154508db56a755355d Mon Sep 17 00:00:00 2001 From: Robert Osfield Date: Thu, 2 Oct 2008 15:45:08 +0000 Subject: [PATCH] Various improvements to the dicom loader to be able to handle a broader range of dicom files --- examples/osgvolume/osgvolume.cpp | 74 ++-- examples/osgvolume/volume_iso_frag.cpp | 2 +- examples/osgvolume/volume_tf_iso_frag.cpp | 2 +- include/osgVolume/ImageUtils | 3 + src/osgDB/Registry.cpp | 3 + src/osgPlugins/dicom/ReaderWriterDICOM.cpp | 377 +++++++++++++-------- src/osgVolume/ImageUtils.cpp | 147 ++++++++ 7 files changed, 426 insertions(+), 182 deletions(-) diff --git a/examples/osgvolume/osgvolume.cpp b/examples/osgvolume/osgvolume.cpp index 1eef6035f..9e2ae80a5 100644 --- a/examples/osgvolume/osgvolume.cpp +++ b/examples/osgvolume/osgvolume.cpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -2019,7 +2020,12 @@ int main( int argc, char **argv ) std::string filename = arguments[pos]; if (osgDB::getLowerCaseFileExtension(filename)=="dicom") { - images.push_back(osgDB::readImageFile( filename )); + // not an option so assume string is a filename. + osg::Image *image = osgDB::readImageFile(filename); + if(image) + { + images.push_back(image); + } } else { @@ -2032,32 +2038,11 @@ int main( int argc, char **argv ) if (fileType == osgDB::DIRECTORY) { - osgDB::DirectoryContents contents = osgDB::getDirectoryContents(filename); - - std::sort(contents.begin(), contents.end()); - - ImageList imageList; - for(osgDB::DirectoryContents::iterator itr = contents.begin(); - itr != contents.end(); - ++itr) + osg::Image *image = osgDB::readImageFile(filename+".dicom"); + if(image) { - std::string localFile = filename + "/" + *itr; - std::cout<<"contents = "<(images.front()->getUserData()); +#if 0 if (matrix) { osg::notify(osg::NOTICE)<<"Image has Matrix = "<<*matrix<((*itr)->getUserData()); - if (matrix) + osg::Vec4 localMinValue, localMaxValue; + if (osgVolume::computeMinMax(itr->get(), localMinValue, localMaxValue)) { - std::cout<<"matrix = "<<*matrix<maxValue.r()) maxValue.r() = localMaxValue.r(); + if (localMaxValue.g()>maxValue.g()) maxValue.g() = localMaxValue.g(); + if (localMaxValue.b()>maxValue.b()) maxValue.b() = localMaxValue.b(); + if (localMaxValue.a()>maxValue.a()) maxValue.a() = localMaxValue.a(); + + osg::notify(osg::NOTICE)<<" ("<getFileName()<get(), minValue, maxValue)) computeMinMax = true; } - if (computeMinMax) { osg::notify(osg::NOTICE)<<"Min value "< imageSequence = new osg::ImageSequence; + imageSequence->setLength(10.0); image_3d = imageSequence.get(); for(Images::iterator itr = images.begin(); itr != images.end(); @@ -2242,6 +2237,15 @@ int main( int argc, char **argv ) numSlices, sliceEnd, alphaFunc); } + if (matrix && rootNode) + { + osg::MatrixTransform* mt = new osg::MatrixTransform; + mt->setMatrix(*matrix); + mt->addChild(rootNode); + + rootNode = mt; + } + if (!outputFile.empty()) { std::string ext = osgDB::getFileExtension(outputFile); diff --git a/examples/osgvolume/volume_iso_frag.cpp b/examples/osgvolume/volume_iso_frag.cpp index afa32401d..462910a99 100644 --- a/examples/osgvolume/volume_iso_frag.cpp +++ b/examples/osgvolume/volume_iso_frag.cpp @@ -102,7 +102,7 @@ char volume_iso_frag[] = "uniform sampler3D baseTexture;\n" " vec3 grad = vec3(px-nx, py-ny, pz-nz);\n" " vec3 normal = normalize(grad);\n" "\n" - " float lightScale = 0.1 + abs(dot(normal.xyz, eyeDirection));\n" + " float lightScale = 0.1 + abs(dot(normal.xyz, eyeDirection))*0.9;\n" " \n" " \n" "#if 0\n" diff --git a/examples/osgvolume/volume_tf_iso_frag.cpp b/examples/osgvolume/volume_tf_iso_frag.cpp index f2d652a32..2f59bf065 100644 --- a/examples/osgvolume/volume_tf_iso_frag.cpp +++ b/examples/osgvolume/volume_tf_iso_frag.cpp @@ -104,7 +104,7 @@ char volume_tf_iso_frag[] = "uniform sampler3D baseTexture;\n" " vec3 grad = vec3(px-nx, py-ny, pz-nz);\n" " vec3 normal = normalize(grad);\n" "\n" - " float lightScale = 0.1 + abs(dot(normal.xyz, eyeDirection));\n" + " float lightScale = 0.1 + abs(dot(normal.xyz, eyeDirection))*0.9;\n" " \n" " color.x *= lightScale;\n" " color.y *= lightScale;\n" diff --git a/include/osgVolume/ImageUtils b/include/osgVolume/ImageUtils index 5468fdec9..9ed63c99e 100644 --- a/include/osgVolume/ImageUtils +++ b/include/osgVolume/ImageUtils @@ -126,6 +126,9 @@ extern OSGVOLUME_EXPORT bool computeMinMax(const osg::Image* image, osg::Vec4& m /** Compute the min max colour values in the image.*/ extern OSGVOLUME_EXPORT bool offsetAndScaleImage(osg::Image* image, const osg::Vec4& offset, const osg::Vec4& scale); +/** Compute source image to destination image.*/ +extern OSGVOLUME_EXPORT bool copyImage(const osg::Image* srcImage, int src_s, int src_t, int src_r, int width, int height, int depth, + osg::Image* destImage, int dest_s, int dest_t, int dest_r, bool doRescale = false); } diff --git a/src/osgDB/Registry.cpp b/src/osgDB/Registry.cpp index 1339ee342..aeac5e463 100644 --- a/src/osgDB/Registry.cpp +++ b/src/osgDB/Registry.cpp @@ -217,6 +217,9 @@ Registry::Registry() addFileExtensionAlias("ivz", "gz"); addFileExtensionAlias("ozg", "gz"); + addFileExtensionAlias("mag", "dicom"); + addFileExtensionAlias("ph", "dicom"); + addFileExtensionAlias("ima", "dicom"); addFileExtensionAlias("dcm", "dicom"); addFileExtensionAlias("dic", "dicom"); diff --git a/src/osgPlugins/dicom/ReaderWriterDICOM.cpp b/src/osgPlugins/dicom/ReaderWriterDICOM.cpp index f9253c3e4..9e76bc441 100644 --- a/src/osgPlugins/dicom/ReaderWriterDICOM.cpp +++ b/src/osgPlugins/dicom/ReaderWriterDICOM.cpp @@ -15,6 +15,7 @@ #include #include +#include #ifdef USE_DCMTK #define HAVE_CONFIG_H @@ -44,11 +45,18 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter ReaderWriterDICOM() { + supportsExtension("mag","dicom image format"); + supportsExtension("ph","dicom image format"); + supportsExtension("ima","dicom image format"); supportsExtension("dic","dicom image format"); supportsExtension("dcm","dicom image format"); supportsExtension("dicom","dicom image format"); // supportsExtension("*","dicom image format"); } + + std::ostream& warning() const { return osg::notify(osg::WARN); } + std::ostream& notice() const { return osg::notify(osg::NOTICE); } + std::ostream& info() const { return osg::notify(osg::NOTICE); } template T* readData(std::istream& fin, unsigned int length, unsigned int& numComponents) const @@ -128,7 +136,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter result.getImage()->setUserData(0); - osg::notify(osg::NOTICE)<<"Locator "<<*matrix<addChild(brick.get()); @@ -152,7 +160,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter std::string fileName = osgDB::findDataFile( file, options ); if (fileName.empty()) return ReadResult::FILE_NOT_FOUND; - osg::notify(osg::NOTICE)<<"Reading DICOM file "< DistanceFileInfoMap; typedef std::map OrientationFileInfoMap; OrientationFileInfoMap orientationFileInfoMap; + + unsigned int totalNumSlices = 0; for(Files::iterator itr = files.begin(); itr != files.end(); @@ -317,55 +387,95 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter fileInfo.filename = *itr; double pixelSize_y = 1.0; - if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, pixelSize_y,0).good()) + 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; + + double value = 0.0; + if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,0).good()) { + pixelSize_y = value; fileInfo.matrix(1,1) = pixelSize_y; } - double pixelSize_x = 1.0; - if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, pixelSize_x,1).good()) + if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,1).good()) { + pixelSize_x = value; fileInfo.matrix(0,0) = pixelSize_x; } // Get slice thickness - double sliceThickness = 1.0; - if (fileformat.getDataset()->findAndGetFloat64(DCM_SliceThickness, sliceThickness).good()) + if (fileformat.getDataset()->findAndGetFloat64(DCM_SliceThickness, value).good()) { - fileInfo.matrix(2,2) = sliceThickness; + sliceThickness = value; + notice()<<"sliceThickness = "<second; if (dfim.empty()) return 0; + 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); + } + + notice()<<"Total Distance including ends "<second; - fileInfoList.push_back(fileInfo); - osg::notify(osg::INFO)<<" d = "<first - dfim.begin()->first; - double averageThickness = dfim.size()<=1 ? 1.0 : totalDistance / double(dfim.size()-1); - - osg::notify(osg::INFO)<<"Average thickness "< dcmImage(new DicomImage(fileInfo.filename.c_str())); if (dcmImage.get()) { if (dcmImage->getStatus()==EIS_Normal) { + // get the pixel data const DiPixel* pixelData = dcmImage->getInterData(); - if(!pixelData) return ReadResult::ERROR_IN_READING_FILE; + if(!pixelData) + { + warning()<<"Error: no data in DicomImage object."< imageAdapter = new osg::Image; + + EP_Representation curr_pixelRep; + unsigned int curr_numPlanes; + GLenum curr_pixelFormat; + GLenum curr_dataType; + unsigned int curr_pixelSize; + // create the new image + convertPixelTypes(pixelData, + curr_pixelRep, curr_numPlanes, + curr_dataType, curr_pixelFormat, curr_pixelSize); + + imageAdapter->setImage(dcmImage->getWidth(), dcmImage->getHeight(), dcmImage->getFrameCount(), + curr_pixelFormat, + curr_pixelFormat, + curr_dataType, + (unsigned char*)(pixelData->getData()), + osg::Image::NO_DELETE); if (!image) { - osg::notify(osg::NOTICE)<<"dcmImage->getWidth() = "<getWidth()<getHeight() = "<getHeight()<getFrameCount() = "<getFrameCount()<getCount() = "<getCount()<getPlanes(); - switch(numPlanes) - { - case(1): - pixelFormat = GL_LUMINANCE; - break; - case(2): - pixelFormat = GL_LUMINANCE_ALPHA; - pixelSize *= 2; - break; - case(3): - pixelFormat = GL_RGB; - pixelSize *= 3; - break; - case(4): - pixelFormat = GL_RGBA; - pixelSize *= 4; - break; - } + pixelRep = curr_pixelRep; + numPlanes = curr_numPlanes; + dataType = curr_dataType; + pixelFormat = curr_pixelFormat; + pixelSize = curr_pixelSize; (*matrix)(0,0) = fileInfo.matrix(0,0); (*matrix)(1,0) = fileInfo.matrix(1,0); @@ -535,38 +611,45 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter image = new osg::Image; image->setUserData(matrix.get()); image->setFileName(fileName.c_str()); - image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), files.size() * dcmImage->getFrameCount(), + image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices, pixelFormat, dataType); - osg::notify(osg::NOTICE)<<"Image dimensions = "<s()<<", "<t()<<", "<r()<t()<<", "<r()<getPlanes()>numPlanes || pixelData->getRepresentation()>pixelRep) + { + notice()<<"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); + + osgVolume::copyImage(previous_image.get(), 0,0,0, previous_image->s(), previous_image->t(), imageNum, + image.get(), 0, 0, 0, + false); + } - if (pixelData->getPlanes()==numPlanes && - pixelData->getRepresentation()==pixelRep && - dcmImage->getWidth()==image->s() && - dcmImage->getHeight()==image->t()) - { - int numFramesToCopy = std::min(static_cast(image->r()-imageNum), - static_cast(dcmImage->getFrameCount())); - unsigned int numPixels = dcmImage->getWidth() * dcmImage->getHeight() * numFramesToCopy; - unsigned int dataSize = numPixels * pixelSize; - - if (invertOrigiantion) - { - memcpy(image->data(0,0,image->r()-imageNum-1), pixelData->getData(), dataSize); - } - else - { - memcpy(image->data(0,0,imageNum), pixelData->getData(), dataSize); - } - - imageNum += numFramesToCopy; - } + osgVolume::copyImage(imageAdapter.get(), 0,0,0, imageAdapter->s(), imageAdapter->t(), imageAdapter->r(), + image.get(), 0, 0, imageNum, + false); + + imageNum += dcmImage->getFrameCount(); } else { - osg::notify(osg::NOTICE)<<"Error in reading dicom file "< (dest_t + destImage->t())) + { + osg::notify(osg::NOTICE)<<"copyImage("< (dest_r + destImage->r())) + { + osg::notify(osg::NOTICE)<<"copyImage("<getDataType() != destImage->getDataType()) + { + switch(srcImage->getDataType()) + { + case(GL_BYTE): scale = 1.0f/128.0f ; break; + case(GL_UNSIGNED_BYTE): scale = 1.0f/255.0f; break; + case(GL_SHORT): scale = 1.0f/32768.0f; break; + case(GL_UNSIGNED_SHORT): scale = 1.0f/65535.0f; break; + case(GL_INT): scale = 1.0f/2147483648.0f; break; + case(GL_UNSIGNED_INT): scale = 1.0f/4294967295.0f; break; + case(GL_FLOAT): scale = 1.0f; break; + } + switch(destImage->getDataType()) + { + case(GL_BYTE): scale *= 128.0f ; break; + case(GL_UNSIGNED_BYTE): scale *= 255.0f; break; + case(GL_SHORT): scale *= 32768.0f; break; + case(GL_UNSIGNED_SHORT): scale *= 65535.0f; break; + case(GL_INT): scale *= 2147483648.0f; break; + case(GL_UNSIGNED_INT): scale *= 4294967295.0f; break; + case(GL_FLOAT): scale *= 1.0f; break; + } + } + + if (srcImage->getPixelFormat() == destImage->getPixelFormat()) + { + //osg::notify(osg::NOTICE)<<"copyImage("<getDataType() == destImage->getDataType() && !doRescale) + { + //osg::notify(osg::NOTICE)<<" Compatible pixelFormat and dataType."<data(src_s, src_t+row, src_r+slice); + unsigned char* destData = destImage->data(dest_s, dest_t+row, dest_r+slice); + memcpy(destData, srcData, (width*destImage->getPixelSizeInBits())/8); + } + } + return true; + } + else + { + //osg::notify(osg::NOTICE)<<" Compatible pixelFormat and incompatible dataType."<data(src_s, src_t+row, src_r+slice); + unsigned char* destData = destImage->data(dest_s, dest_t+row, dest_r+slice); + unsigned int numComponents = osg::Image::computeNumComponents(destImage->getPixelFormat()); + + _copyRowAndScale(srcData, srcImage->getDataType(), destData, destImage->getDataType(), (width*numComponents), scale); + } + } + + return true; + } + } + else + { + osg::notify(osg::NOTICE)<<"copyImage("<