FGFS-4.1/OpenSceneGraph
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Various improvements to the dicom loader to be able to handle a broader range of dicom files

Browse Source
This commit is contained in:
Robert Osfield 2008-10-02 15:45:08 +00:00
parent 66c857b645
commit 3b4184295e
7 changed files with 426 additions and 182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
examples/osgvolume/osgvolume.cpp
View File

@ -37,6 +37,7 @@
#include <osg/BlendFunc> #include <osg/BlendFunc>
#include <osg/BlendEquation> #include <osg/BlendEquation>
#include <osg/TransferFunction> #include <osg/TransferFunction>
#include <osg/MatrixTransform>
#include <osgDB/Registry> #include <osgDB/Registry>
#include <osgDB/ReadFile> #include <osgDB/ReadFile>
@ -2019,7 +2020,12 @@ int main( int argc, char **argv )
std::string filename = arguments[pos]; std::string filename = arguments[pos];
if (osgDB::getLowerCaseFileExtension(filename)=="dicom") 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 else
{ {
@ -2032,32 +2038,11 @@ int main( int argc, char **argv )
if (fileType == osgDB::DIRECTORY) if (fileType == osgDB::DIRECTORY)
{ {
osgDB::DirectoryContents contents = osgDB::getDirectoryContents(filename); osg::Image *image = osgDB::readImageFile(filename+".dicom");
if(image)
std::sort(contents.begin(), contents.end());
ImageList imageList;
for(osgDB::DirectoryContents::iterator itr = contents.begin();
itr != contents.end();
++itr)
{ {
std::string localFile = filename + "/" + *itr; images.push_back(image);
std::cout<<"contents = "<<localFile<<std::endl; }
if (osgDB::fileType(localFile) == osgDB::REGULAR_FILE)
{
// not an option so assume string is a filename.
osg::Image *image = osgDB::readImageFile(localFile);
if(image)
{
imageList.push_back(image);
}
}
}
// pack the textures into a single texture.
ProcessRow processRow;
images.push_back(createTexture3D(imageList, processRow, numComponentsDesired, s_maximumTextureSize, t_maximumTextureSize, r_maximumTextureSize, resizeToPowerOfTwo));
} }
else if (fileType == osgDB::REGULAR_FILE) else if (fileType == osgDB::REGULAR_FILE)
{ {
@ -2098,8 +2083,8 @@ int main( int argc, char **argv )
} }
#if 1
osg::RefMatrix* matrix = dynamic_cast<osg::RefMatrix*>(images.front()->getUserData()); osg::RefMatrix* matrix = dynamic_cast<osg::RefMatrix*>(images.front()->getUserData());
#if 0
if (matrix) if (matrix)
{ {
osg::notify(osg::NOTICE)<<"Image has Matrix = "<<*matrix<<std::endl; osg::notify(osg::NOTICE)<<"Image has Matrix = "<<*matrix<<std::endl;
@ -2109,23 +2094,32 @@ int main( int argc, char **argv )
} }
#endif #endif
osg::Vec4 minValue, maxValue; osg::Vec4 minValue(FLT_MAX, FLT_MAX, FLT_MAX, FLT_MAX);
osg::Vec4 maxValue(-FLT_MAX, -FLT_MAX, -FLT_MAX, -FLT_MAX);
bool computeMinMax = false; bool computeMinMax = false;
for(Images::iterator itr = images.begin(); for(Images::iterator itr = images.begin();
itr != images.end(); itr != images.end();
++itr) ++itr)
{ {
#if 0 osg::Vec4 localMinValue, localMaxValue;
osg::RefMatrix* matrix = dynamic_cast<osg::RefMatrix*>((*itr)->getUserData()); if (osgVolume::computeMinMax(itr->get(), localMinValue, localMaxValue))
if (matrix)
{ {
std::cout<<"matrix = "<<*matrix<<std::endl; if (localMinValue.r()<minValue.r()) minValue.r() = localMinValue.r();
if (localMinValue.g()<minValue.g()) minValue.g() = localMinValue.g();
if (localMinValue.b()<minValue.b()) minValue.b() = localMinValue.b();
if (localMinValue.a()<minValue.a()) minValue.a() = localMinValue.a();
if (localMaxValue.r()>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)<<" ("<<localMinValue<<") ("<<localMaxValue<<") "<<(*itr)->getFileName()<<std::endl;
computeMinMax = true;
} }
#endif
if (osgVolume::computeMinMax(itr->get(), minValue, maxValue)) computeMinMax = true;
} }
if (computeMinMax) if (computeMinMax)
{ {
osg::notify(osg::NOTICE)<<"Min value "<<minValue<<std::endl; osg::notify(osg::NOTICE)<<"Min value "<<minValue<<std::endl;
@ -2188,6 +2182,7 @@ int main( int argc, char **argv )
osg::notify(osg::NOTICE)<<"Creating sequence of "<<images.size()<<" volumes."<<std::endl; osg::notify(osg::NOTICE)<<"Creating sequence of "<<images.size()<<" volumes."<<std::endl;
osg::ref_ptr<osg::ImageSequence> imageSequence = new osg::ImageSequence; osg::ref_ptr<osg::ImageSequence> imageSequence = new osg::ImageSequence;
imageSequence->setLength(10.0);
image_3d = imageSequence.get(); image_3d = imageSequence.get();
for(Images::iterator itr = images.begin(); for(Images::iterator itr = images.begin();
itr != images.end(); itr != images.end();
@ -2242,6 +2237,15 @@ int main( int argc, char **argv )
numSlices, sliceEnd, alphaFunc); numSlices, sliceEnd, alphaFunc);
} }
if (matrix && rootNode)
{
osg::MatrixTransform* mt = new osg::MatrixTransform;
mt->setMatrix(*matrix);
mt->addChild(rootNode);
rootNode = mt;
}
if (!outputFile.empty()) if (!outputFile.empty())
{ {
std::string ext = osgDB::getFileExtension(outputFile); std::string ext = osgDB::getFileExtension(outputFile);

2
examples/osgvolume/volume_iso_frag.cpp
View File

@ -102,7 +102,7 @@ char volume_iso_frag[] = "uniform sampler3D baseTexture;\n"
" vec3 grad = vec3(px-nx, py-ny, pz-nz);\n" " vec3 grad = vec3(px-nx, py-ny, pz-nz);\n"
" vec3 normal = normalize(grad);\n" " vec3 normal = normalize(grad);\n"
"\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"
" \n" " \n"
"#if 0\n" "#if 0\n"

2
examples/osgvolume/volume_tf_iso_frag.cpp
View File

@ -104,7 +104,7 @@ char volume_tf_iso_frag[] = "uniform sampler3D baseTexture;\n"
" vec3 grad = vec3(px-nx, py-ny, pz-nz);\n" " vec3 grad = vec3(px-nx, py-ny, pz-nz);\n"
" vec3 normal = normalize(grad);\n" " vec3 normal = normalize(grad);\n"
"\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"
" color.x *= lightScale;\n" " color.x *= lightScale;\n"
" color.y *= lightScale;\n" " color.y *= lightScale;\n"

3
include/osgVolume/ImageUtils
View File

@ -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.*/ /** 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); 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);
} }

3
src/osgDB/Registry.cpp
View File

@ -217,6 +217,9 @@ Registry::Registry()
addFileExtensionAlias("ivz", "gz"); addFileExtensionAlias("ivz", "gz");
addFileExtensionAlias("ozg", "gz"); addFileExtensionAlias("ozg", "gz");
addFileExtensionAlias("mag", "dicom");
addFileExtensionAlias("ph", "dicom");
addFileExtensionAlias("ima", "dicom");
addFileExtensionAlias("dcm", "dicom"); addFileExtensionAlias("dcm", "dicom");
addFileExtensionAlias("dic", "dicom"); addFileExtensionAlias("dic", "dicom");

377
src/osgPlugins/dicom/ReaderWriterDICOM.cpp
View File

@ -15,6 +15,7 @@
#include <osgVolume/Volume> #include <osgVolume/Volume>
#include <osgVolume/Brick> #include <osgVolume/Brick>
#include <osgVolume/ImageUtils>
#ifdef USE_DCMTK #ifdef USE_DCMTK
#define HAVE_CONFIG_H #define HAVE_CONFIG_H
@ -44,11 +45,18 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
ReaderWriterDICOM() ReaderWriterDICOM()
{ {
supportsExtension("mag","dicom image format");
supportsExtension("ph","dicom image format");
supportsExtension("ima","dicom image format");
supportsExtension("dic","dicom image format"); supportsExtension("dic","dicom image format");
supportsExtension("dcm","dicom image format"); supportsExtension("dcm","dicom image format");
supportsExtension("dicom","dicom image format"); supportsExtension("dicom","dicom image format");
// supportsExtension("*","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<typename T> template<typename T>
T* readData(std::istream& fin, unsigned int length, unsigned int& numComponents) const 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); result.getImage()->setUserData(0);
osg::notify(osg::NOTICE)<<"Locator "<<*matrix<<std::endl; notice()<<"Locator "<<*matrix<<std::endl;
} }
volume->addChild(brick.get()); volume->addChild(brick.get());
@ -152,7 +160,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
std::string fileName = osgDB::findDataFile( file, options ); std::string fileName = osgDB::findDataFile( file, options );
if (fileName.empty()) return ReadResult::FILE_NOT_FOUND; if (fileName.empty()) return ReadResult::FILE_NOT_FOUND;
osg::notify(osg::NOTICE)<<"Reading DICOM file "<<fileName<<std::endl; notice()<<"Reading DICOM file "<<fileName<<std::endl;
typedef unsigned short PixelType; typedef unsigned short PixelType;
const unsigned int Dimension = 3; const unsigned int Dimension = 3;
@ -194,7 +202,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
osg::RefMatrix* matrix = new osg::RefMatrix; osg::RefMatrix* matrix = new osg::RefMatrix;
osg::notify(osg::NOTICE)<<"width = "<<width<<" height = "<<height<<" depth = "<<depth<<std::endl; notice()<<"width = "<<width<<" height = "<<height<<" depth = "<<depth<<std::endl;
for(unsigned int i=0; i<Dimension; ++i) for(unsigned int i=0; i<Dimension; ++i)
{ {
(*matrix)(i,i) = inputImage->GetSpacing()[i]; (*matrix)(i,i) = inputImage->GetSpacing()[i];
@ -224,6 +232,66 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
#ifdef USE_DCMTK #ifdef USE_DCMTK
void convertPixelTypes(const DiPixel* pixelData,
EP_Representation& pixelRep, unsigned int& numPlanes,
GLenum& dataType, GLenum& pixelFormat, unsigned int& pixelSize) const
{
dataType = GL_UNSIGNED_BYTE;
pixelRep = pixelData->getRepresentation();
switch(pixelRep)
{
case(EPR_Uint8):
dataType = GL_UNSIGNED_BYTE;
pixelSize = 1;
break;
case(EPR_Sint8):
dataType = GL_BYTE;
pixelSize = 1;
break;
case(EPR_Uint16):
dataType = GL_UNSIGNED_SHORT;
pixelSize = 2;
break;
case(EPR_Sint16):
dataType = GL_SHORT;
pixelSize = 2;
break;
case(EPR_Uint32):
dataType = GL_UNSIGNED_INT;
pixelSize = 4;
break;
case(EPR_Sint32):
dataType = GL_INT;
pixelSize = 4;
break;
default:
dataType = 0;
break;
}
pixelFormat = GL_INTENSITY;
numPlanes = pixelData->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;
}
}
typedef std::vector<std::string> Files; typedef std::vector<std::string> Files;
bool getDicomFilesInDirectory(const std::string& path, Files& files) const bool getDicomFilesInDirectory(const std::string& path, Files& files) const
{ {
@ -236,8 +304,8 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
++itr) ++itr)
{ {
std::string localFile = path + "/" + *itr; std::string localFile = path + "/" + *itr;
osg::notify(osg::INFO)<<"contents = "<<localFile<<std::endl; info()<<"contents = "<<localFile<<std::endl;
if (osgDB::getLowerCaseFileExtension(localFile)=="dcm" && if (acceptsExtension(osgDB::getLowerCaseFileExtension(localFile)) &&
osgDB::fileType(localFile) == osgDB::REGULAR_FILE) osgDB::fileType(localFile) == osgDB::REGULAR_FILE)
{ {
files.push_back(localFile); files.push_back(localFile);
@ -250,7 +318,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
virtual ReadResult readImage(const std::string& file, const osgDB::ReaderWriter::Options* options) const virtual ReadResult readImage(const std::string& file, const osgDB::ReaderWriter::Options* options) const
{ {
osg::notify(osg::NOTICE)<<"Reading DICOM file "<<file<<" using DCMTK"<<std::endl; notice()<<"Reading DICOM file "<<file<<" using DCMTK"<<std::endl;
std::string ext = osgDB::getLowerCaseFileExtension(file); std::string ext = osgDB::getLowerCaseFileExtension(file);
if (!acceptsExtension(ext)) return ReadResult::FILE_NOT_HANDLED; if (!acceptsExtension(ext)) return ReadResult::FILE_NOT_HANDLED;
@ -304,6 +372,8 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
typedef std::map<double, FileInfo> DistanceFileInfoMap; typedef std::map<double, FileInfo> DistanceFileInfoMap;
typedef std::map<osg::Vec3d, DistanceFileInfoMap> OrientationFileInfoMap; typedef std::map<osg::Vec3d, DistanceFileInfoMap> OrientationFileInfoMap;
OrientationFileInfoMap orientationFileInfoMap; OrientationFileInfoMap orientationFileInfoMap;
unsigned int totalNumSlices = 0;
for(Files::iterator itr = files.begin(); for(Files::iterator itr = files.begin();
itr != files.end(); itr != files.end();
@ -317,55 +387,95 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
fileInfo.filename = *itr; fileInfo.filename = *itr;
double pixelSize_y = 1.0; 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; fileInfo.matrix(1,1) = pixelSize_y;
} }
double pixelSize_x = 1.0; if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, value,1).good())
if (fileformat.getDataset()->findAndGetFloat64(DCM_PixelSpacing, pixelSize_x,1).good())
{ {
pixelSize_x = value;
fileInfo.matrix(0,0) = pixelSize_x; fileInfo.matrix(0,0) = pixelSize_x;
} }
// Get slice thickness // Get slice thickness
double sliceThickness = 1.0; if (fileformat.getDataset()->findAndGetFloat64(DCM_SliceThickness, value).good())
if (fileformat.getDataset()->findAndGetFloat64(DCM_SliceThickness, sliceThickness).good())
{ {
fileInfo.matrix(2,2) = sliceThickness; sliceThickness = value;
notice()<<"sliceThickness = "<<sliceThickness<<std::endl;
fileInfo.sliceThickness = sliceThickness;
}
notice()<<"tagExistsWithValue(DCM_NumberOfFrames)="<<fileformat.getDataset()->tagExistsWithValue(DCM_NumberOfFrames)<<std::endl;
notice()<<"tagExistsWithValue(DCM_NumberOfSlices)="<<fileformat.getDataset()->tagExistsWithValue(DCM_NumberOfSlices)<<std::endl;
Uint32 numFrames;
if (fileformat.getDataset()->findAndGetUint32(DCM_NumberOfFrames, numFrames).good())
{
fileInfo.numSlices = numFrames;
notice()<<"Read number of frames = "<<numFrames<<std::endl;
} }
double imagePositionPatient[3] = {0, 0, 0};
OFString numFramesStr;
if (fileformat.getDataset()->findAndGetOFString(DCM_NumberOfFrames, numFramesStr).good())
{
fileInfo.numSlices = atoi(numFramesStr.c_str());
notice()<<"Read number of frames = "<<numFramesStr<<std::endl;
}
if (fileformat.getDataset()->findAndGetUint16(DCM_NumberOfFrames, numOfSlices).good())
{
fileInfo.numSlices = numOfSlices;
notice()<<"Read number of frames = "<<numOfSlices<<std::endl;
}
if (fileformat.getDataset()->findAndGetUint16(DCM_NumberOfSlices, numOfSlices).good())
{
fileInfo.numSlices = numOfSlices;
notice()<<"Read number of slices = "<<numOfSlices<<std::endl;
}
// patient position // patient position
for(int i=0; i<3; ++i) for(int i=0; i<3; ++i)
{ {
if (fileformat.getDataset()->findAndGetFloat64(DCM_ImagePositionPatient, imagePositionPatient[i],i).good()) if (fileformat.getDataset()->findAndGetFloat64(DCM_ImagePositionPatient, imagePositionPatient[i],i).good())
{ {
osg::notify(osg::INFO)<<"Read DCM_ImagePositionPatient["<<i<<"], "<<imagePositionPatient[i]<<std::endl; notice()<<"Read DCM_ImagePositionPatient["<<i<<"], "<<imagePositionPatient[i]<<std::endl;
} }
else else
{ {
osg::notify(osg::INFO)<<"Have not read DCM_ImagePositionPatient["<<i<<"]"<<std::endl; notice()<<"Have not read DCM_ImagePositionPatient["<<i<<"]"<<std::endl;
} }
} }
//notice()<<"imagePositionPatient[2]="<<imagePositionPatient[2]<<std::endl;
fileInfo.matrix.setTrans(imagePositionPatient[0],imagePositionPatient[1],imagePositionPatient[2]); fileInfo.matrix.setTrans(imagePositionPatient[0],imagePositionPatient[1],imagePositionPatient[2]);
double imageOrientationPatient[6] = { 1.0, 0.0, 0.0,
0.0, -1.0, 0.0 };
for(int i=0; i<6; ++i) for(int i=0; i<6; ++i)
{ {
if (fileformat.getDataset()->findAndGetFloat64(DCM_ImageOrientationPatient, imageOrientationPatient[i],i).good()) double value = 0.0;
if (fileformat.getDataset()->findAndGetFloat64(DCM_ImageOrientationPatient, value,i).good())
{ {
osg::notify(osg::INFO)<<"Read imageOrientationPatient["<<i<<"], "<<imageOrientationPatient[i]<<std::endl; imageOrientationPatient[i] = value;
notice()<<"Read imageOrientationPatient["<<i<<"], "<<imageOrientationPatient[i]<<std::endl;
} }
else else
{ {
osg::notify(osg::INFO)<<"Have not read imageOrientationPatient["<<i<<"]"<<std::endl; notice()<<"Have not read imageOrientationPatient["<<i<<"]"<<std::endl;
} }
} }
osg::Vec3d dirX(imageOrientationPatient[0],imageOrientationPatient[1],imageOrientationPatient[2]); osg::Vec3d dirX(imageOrientationPatient[0],imageOrientationPatient[1],imageOrientationPatient[2]);
osg::Vec3d dirY(imageOrientationPatient[3],imageOrientationPatient[4],imageOrientationPatient[5]); osg::Vec3d dirY(imageOrientationPatient[3],imageOrientationPatient[4],imageOrientationPatient[5]);
osg::Vec3d dirZ = dirX ^ dirY; osg::Vec3d dirZ = dirX ^ dirY;
@ -388,16 +498,18 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
fileInfo.distance = dirZ * (osg::Vec3d(0.0,0.0,0.0)*fileInfo.matrix); fileInfo.distance = dirZ * (osg::Vec3d(0.0,0.0,0.0)*fileInfo.matrix);
osg::notify(osg::INFO)<<"dirX = "<<dirX<<std::endl; notice()<<"dirX = "<<dirX<<std::endl;
osg::notify(osg::INFO)<<"dirY = "<<dirY<<std::endl; notice()<<"dirY = "<<dirY<<std::endl;
osg::notify(osg::INFO)<<"dirZ = "<<dirZ<<std::endl; notice()<<"dirZ = "<<dirZ<<std::endl;
osg::notify(osg::INFO)<<"matrix = "<<fileInfo.matrix<<std::endl; notice()<<"matrix = "<<fileInfo.matrix<<std::endl;
osg::notify(osg::INFO)<<"pos = "<<osg::Vec3d(0.0,0.0,0.0)*fileInfo.matrix<<std::endl; notice()<<"pos = "<<osg::Vec3d(0.0,0.0,0.0)*fileInfo.matrix<<std::endl;
osg::notify(osg::INFO)<<"dist = "<<fileInfo.distance<<std::endl; notice()<<"dist = "<<fileInfo.distance<<std::endl;
osg::notify(osg::INFO)<<std::endl; notice()<<std::endl;
(orientationFileInfoMap[dirZ])[fileInfo.distance] = fileInfo; (orientationFileInfoMap[dirZ])[fileInfo.distance] = fileInfo;
totalNumSlices += fileInfo.numSlices;
} }
if (orientationFileInfoMap.empty()) return 0; if (orientationFileInfoMap.empty()) return 0;
@ -408,14 +520,14 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
itr != orientationFileInfoMap.end(); itr != orientationFileInfoMap.end();
++itr) ++itr)
{ {
osg::notify(osg::INFO)<<"Orientation = "<<itr->first<<std::endl; notice()<<"Orientation = "<<itr->first<<std::endl;
DistanceFileInfoMap& dfim = itr->second; DistanceFileInfoMap& dfim = itr->second;
for(DistanceFileInfoMap::iterator ditr = dfim.begin(); for(DistanceFileInfoMap::iterator ditr = dfim.begin();
ditr != dfim.end(); ditr != dfim.end();
++ditr) ++ditr)
{ {
FileInfo& fileInfo = ditr->second; FileInfo& fileInfo = ditr->second;
osg::notify(osg::INFO)<<" d = "<<fileInfo.distance<<" "<<fileInfo.filename<<std::endl; notice()<<" d = "<<fileInfo.distance<<" "<<fileInfo.filename<<std::endl;
} }
} }
@ -423,104 +535,68 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
DistanceFileInfoMap& dfim = orientationFileInfoMap.begin()->second; DistanceFileInfoMap& dfim = orientationFileInfoMap.begin()->second;
if (dfim.empty()) return 0; 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 "<<totalDistance<<std::endl;
double averageThickness = totalNumSlices<=1 ? 1.0 : totalDistance / double(totalNumSlices-1);
notice()<<"Average thickness "<<averageThickness<<std::endl;
for(DistanceFileInfoMap::iterator ditr = dfim.begin(); for(DistanceFileInfoMap::iterator ditr = dfim.begin();
ditr != dfim.end(); ditr != dfim.end();
++ditr) ++ditr)
{ {
FileInfo& fileInfo = ditr->second; FileInfo& fileInfo = ditr->second;
fileInfoList.push_back(fileInfo);
osg::notify(osg::INFO)<<" d = "<<fileInfo.distance<<" "<<fileInfo.filename<<std::endl;
}
double totalDistance = dfim.rbegin()->first - dfim.begin()->first;
double averageThickness = dfim.size()<=1 ? 1.0 : totalDistance / double(dfim.size()-1);
osg::notify(osg::INFO)<<"Average thickness "<<averageThickness<<std::endl;
for(FileInfoList::iterator itr = fileInfoList.begin();
itr != fileInfoList.end();
++itr)
{
FileInfo& fileInfo = *itr;
std::auto_ptr<DicomImage> dcmImage(new DicomImage(fileInfo.filename.c_str())); std::auto_ptr<DicomImage> dcmImage(new DicomImage(fileInfo.filename.c_str()));
if (dcmImage.get()) if (dcmImage.get())
{ {
if (dcmImage->getStatus()==EIS_Normal) if (dcmImage->getStatus()==EIS_Normal)
{ {
// get the pixel data // get the pixel data
const DiPixel* pixelData = dcmImage->getInterData(); const DiPixel* pixelData = dcmImage->getInterData();
if(!pixelData) return ReadResult::ERROR_IN_READING_FILE; if(!pixelData)
{
warning()<<"Error: no data in DicomImage object."<<std::endl;
return ReadResult::ERROR_IN_READING_FILE;
}
osg::ref_ptr<osg::Image> 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) if (!image)
{ {
osg::notify(osg::NOTICE)<<"dcmImage->getWidth() = "<<dcmImage->getWidth()<<std::endl; pixelRep = curr_pixelRep;
osg::notify(osg::NOTICE)<<"dcmImage->getHeight() = "<<dcmImage->getHeight()<<std::endl; numPlanes = curr_numPlanes;
osg::notify(osg::NOTICE)<<"dcmImage->getFrameCount() = "<<dcmImage->getFrameCount()<<std::endl; dataType = curr_dataType;
osg::notify(osg::NOTICE)<<"pixelData->getCount() = "<<pixelData->getCount()<<std::endl; pixelFormat = curr_pixelFormat;
osg::notify(osg::NOTICE)<<"pixelFormat = "; pixelSize = curr_pixelSize;
dataType = GL_UNSIGNED_BYTE;
pixelRep = pixelData->getRepresentation();
switch(pixelRep)
{
case(EPR_Uint8):
dataType = GL_UNSIGNED_BYTE;
pixelSize = 1;
osg::notify(osg::NOTICE)<<"unsigned char"<<std::endl;
break;
case(EPR_Sint8):
dataType = GL_BYTE;
pixelSize = 1;
osg::notify(osg::NOTICE)<<"char"<<std::endl;
break;
case(EPR_Uint16):
dataType = GL_UNSIGNED_SHORT;
pixelSize = 2;
osg::notify(osg::NOTICE)<<"unsigned short"<<std::endl;
break;
case(EPR_Sint16):
dataType = GL_SHORT;
pixelSize = 2;
osg::notify(osg::NOTICE)<<"short"<<std::endl;
break;
case(EPR_Uint32):
dataType = GL_UNSIGNED_INT;
pixelSize = 4;
osg::notify(osg::NOTICE)<<"unsigned int"<<std::endl;
break;
case(EPR_Sint32):
osg::notify(osg::NOTICE)<<"int"<<std::endl;
dataType = GL_INT;
pixelSize = 4;
break;
default:
osg::notify(osg::NOTICE)<<"unidentified"<<std::endl;
dataType = 0;
break;
}
pixelFormat = GL_INTENSITY;
numPlanes = pixelData->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;
}
(*matrix)(0,0) = fileInfo.matrix(0,0); (*matrix)(0,0) = fileInfo.matrix(0,0);
(*matrix)(1,0) = fileInfo.matrix(1,0); (*matrix)(1,0) = fileInfo.matrix(1,0);
@ -535,38 +611,45 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
image = new osg::Image; image = new osg::Image;
image->setUserData(matrix.get()); image->setUserData(matrix.get());
image->setFileName(fileName.c_str()); image->setFileName(fileName.c_str());
image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), files.size() * dcmImage->getFrameCount(), image->allocateImage(dcmImage->getWidth(), dcmImage->getHeight(), totalNumSlices,
pixelFormat, dataType); pixelFormat, dataType);
osg::notify(osg::NOTICE)<<"Image dimensions = "<<image->s()<<", "<<image->t()<<", "<<image->r()<<std::endl; notice()<<"Image dimensions = "<<image->s()<<", "<<image->t()<<", "<<image->r()<<std::endl;
}
else if (pixelData->getPlanes()>numPlanes || pixelData->getRepresentation()>pixelRep)
{
notice()<<"Need to reallocated "<<image->s()<<", "<<image->t()<<", "<<image->r()<<std::endl;
// record the previous image settings to use when we copy back the content.
osg::ref_ptr<osg::Image> 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 && osgVolume::copyImage(imageAdapter.get(), 0,0,0, imageAdapter->s(), imageAdapter->t(), imageAdapter->r(),
pixelData->getRepresentation()==pixelRep && image.get(), 0, 0, imageNum,
dcmImage->getWidth()==image->s() && false);
dcmImage->getHeight()==image->t())
{ imageNum += dcmImage->getFrameCount();
int numFramesToCopy = std::min(static_cast<unsigned int>(image->r()-imageNum),
static_cast<unsigned int>(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;
}
} }
else else
{ {
osg::notify(osg::NOTICE)<<"Error in reading dicom file "<<fileName.c_str()<<", error = "<<DicomImage::getString(dcmImage->getStatus())<<std::endl; warning()<<"Error in reading dicom file "<<fileName.c_str()<<", error = "<<DicomImage::getString(dcmImage->getStatus())<<std::endl;
} }
} }
} }
@ -576,7 +659,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
return ReadResult::ERROR_IN_READING_FILE; return ReadResult::ERROR_IN_READING_FILE;
} }
osg::notify(osg::NOTICE)<<"Spacing = "<<*matrix<<std::endl; notice()<<"Spacing = "<<*matrix<<std::endl;
return image.get(); return image.get();
} }
@ -587,7 +670,8 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
FileInfo(): FileInfo():
numX(0), numX(0),
numY(0), numY(0),
numSlices(0), numSlices(1),
sliceThickness(0.0),
distance(0.0) {} distance(0.0) {}
FileInfo(const FileInfo& rhs): FileInfo(const FileInfo& rhs):
@ -596,6 +680,7 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
numX(rhs.numX), numX(rhs.numX),
numY(rhs.numY), numY(rhs.numY),
numSlices(rhs.numSlices), numSlices(rhs.numSlices),
sliceThickness(rhs.sliceThickness),
distance(distance) {} distance(distance) {}
FileInfo& operator = (const FileInfo& rhs) FileInfo& operator = (const FileInfo& rhs)
@ -606,16 +691,18 @@ class ReaderWriterDICOM : public osgDB::ReaderWriter
matrix = rhs.matrix; matrix = rhs.matrix;
numX = rhs.numX; numX = rhs.numX;
numY = rhs.numY; numY = rhs.numY;
sliceThickness = rhs.sliceThickness;
numSlices = rhs.numSlices; numSlices = rhs.numSlices;
distance = rhs.distance; distance = rhs.distance;
} }
std::string filename; std::string filename;
osg::Matrixd matrix; osg::Matrixd matrix;
unsigned int numX; unsigned int numX;
unsigned int numY; unsigned int numY;
unsigned int numSlices; unsigned int numSlices;
double distance; double sliceThickness;
double distance;
}; };
}; };

147
src/osgVolume/ImageUtils.cpp
View File

@ -119,5 +119,152 @@ bool osgVolume::offsetAndScaleImage(osg::Image* image, const osg::Vec4& offset,
return true; return true;
} }
template<typename SRC, typename DEST>
void _copyRowAndScale(const SRC* src, DEST* dest, int num, float scale)
{
if (scale==1.0)
{
for(int i=0; i<num; ++i)
{
*dest = DEST(*src);
++dest; ++src;
}
}
else
{
for(int i=0; i<num; ++i)
{
*dest = DEST(float(*src)*scale);
++dest; ++src;
}
}
}
template<typename DEST>
void _copyRowAndScale(const unsigned char* src, GLenum srcDataType, DEST* dest, int num, float scale)
{
switch(srcDataType)
{
case(GL_BYTE): _copyRowAndScale((char*)src, dest, num, scale); break;
case(GL_UNSIGNED_BYTE): _copyRowAndScale((unsigned char*)src, dest, num, scale); break;
case(GL_SHORT): _copyRowAndScale((short*)src, dest, num, scale); break;
case(GL_UNSIGNED_SHORT): _copyRowAndScale((unsigned short*)src, dest, num, scale); break;
case(GL_INT): _copyRowAndScale((int*)src, dest, num, scale); break;
case(GL_UNSIGNED_INT): _copyRowAndScale((unsigned int*)src, dest, num, scale); break;
case(GL_FLOAT): _copyRowAndScale((float*)src, dest, num, scale); break;
}
}
void _copyRowAndScale(const unsigned char* src, GLenum srcDataType, unsigned char* dest, GLenum dstDataType, int num, float scale)
{
switch(dstDataType)
{
case(GL_BYTE): _copyRowAndScale(src, srcDataType, (char*)dest, num, scale); break;
case(GL_UNSIGNED_BYTE): _copyRowAndScale(src, srcDataType, (unsigned char*)dest, num, scale); break;
case(GL_SHORT): _copyRowAndScale(src, srcDataType, (short*)dest, num, scale); break;
case(GL_UNSIGNED_SHORT): _copyRowAndScale(src, srcDataType, (unsigned short*)dest, num, scale); break;
case(GL_INT): _copyRowAndScale(src, srcDataType, (int*)dest, num, scale); break;
case(GL_UNSIGNED_INT): _copyRowAndScale(src, srcDataType, (unsigned int*)dest, num, scale); break;
case(GL_FLOAT): _copyRowAndScale(src, srcDataType, (float*)dest, num, scale); break;
}
}
bool osgVolume::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)
{
if ((src_s+width) > (dest_s + destImage->s()))
{
osg::notify(osg::NOTICE)<<"copyImage("<<srcImage<<", "<<src_s<<", "<< src_t<<", "<<src_r<<", "<<width<<", "<<height<<", "<<depth<<std::endl;
osg::notify(osg::NOTICE)<<" "<<destImage<<", "<<dest_s<<", "<< dest_t<<", "<<dest_r<<", "<<doRescale<<")"<<std::endl;
osg::notify(osg::NOTICE)<<" input width too large."<<std::endl;
return false;
}
if ((src_t+height) > (dest_t + destImage->t()))
{
osg::notify(osg::NOTICE)<<"copyImage("<<srcImage<<", "<<src_s<<", "<< src_t<<", "<<src_r<<", "<<width<<", "<<height<<", "<<depth<<std::endl;
osg::notify(osg::NOTICE)<<" "<<destImage<<", "<<dest_s<<", "<< dest_t<<", "<<dest_r<<", "<<doRescale<<")"<<std::endl;
osg::notify(osg::NOTICE)<<" input height too large."<<std::endl;
return false;
}
if ((src_r+depth) > (dest_r + destImage->r()))
{
osg::notify(osg::NOTICE)<<"copyImage("<<srcImage<<", "<<src_s<<", "<< src_t<<", "<<src_r<<", "<<width<<", "<<height<<", "<<depth<<std::endl;
osg::notify(osg::NOTICE)<<" "<<destImage<<", "<<dest_s<<", "<< dest_t<<", "<<dest_r<<", "<<doRescale<<")"<<std::endl;
osg::notify(osg::NOTICE)<<" input depth too large."<<std::endl;
return false;
}
float scale = 1.0f;
if (doRescale && srcImage->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("<<srcImage<<", "<<src_s<<", "<< src_t<<", "<<src_r<<", "<<width<<", "<<height<<", "<<depth<<std::endl;
//osg::notify(osg::NOTICE)<<" "<<destImage<<", "<<dest_s<<", "<< dest_t<<", "<<dest_r<<", "<<doRescale<<")"<<std::endl;
if (srcImage->getDataType() == destImage->getDataType() && !doRescale)
{
//osg::notify(osg::NOTICE)<<" Compatible pixelFormat and dataType."<<std::endl;
for(int slice = 0; slice<depth; ++slice)
{
for(int row = 0; row<height; ++row)
{
const unsigned char* srcData = srcImage->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."<<std::endl;
for(int slice = 0; slice<depth; ++slice)
{
for(int row = 0; row<height; ++row)
{
const unsigned char* srcData = srcImage->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("<<srcImage<<", "<<src_s<<", "<< src_t<<", "<<src_r<<", "<<width<<", "<<height<<", "<<depth<<std::endl;
osg::notify(osg::NOTICE)<<" "<<destImage<<", "<<dest_s<<", "<< dest_t<<", "<<dest_r<<", "<<doRescale<<")"<<std::endl;
osg::notify(osg::NOTICE)<<" Incompatible pixelFormat and dataType."<<std::endl;
return false;
}
}