00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "ioDICOM.hpp"
00033 #include <vlCore/LoadWriterManager.hpp>
00034 #include <vlCore/FileSystem.hpp>
00035 #include <vlCore/glsl_math.hpp>
00036
00037 #include <gdcmReader.h>
00038 #include <gdcmWriter.h>
00039 #include <gdcmAttribute.h>
00040 #include <gdcmImageReader.h>
00041 #include <gdcmImageWriter.h>
00042 #include <gdcmImage.h>
00043 #include <gdcmPhotometricInterpretation.h>
00044 #include <gdcmImage.h>
00045 #include <gdcmImageWriter.h>
00046
00047 #include <memory>
00048
00049 using namespace vl;
00050
00051
00052 static void to8bits(int bits_used, void* ptr, int px_count, bool reverse)
00053 {
00054 unsigned char* px = (unsigned char*)ptr;
00055 if (bits_used >= 8)
00056 return;
00057 unsigned int max1 = (1<<bits_used)-1;
00058 unsigned int max2 = 0xFF;
00059 for(int i=0; i<px_count; ++i)
00060 if (reverse)
00061 px[i] = 0xFF - (unsigned char)( (float)px[i]/max1*max2 );
00062 else
00063 px[i] = (unsigned char)( (float)px[i]/max1*max2 );
00064 }
00065 static void to16bits(int bits_used, void* ptr, int px_count, bool reverse)
00066 {
00067 unsigned short* px = (unsigned short*)ptr;
00068
00069 if (bits_used >= 16)
00070 return;
00071 unsigned int max1 = (1<<bits_used)-1;
00072 unsigned int max2 = 0xFFFF;
00073 for(int i=0; i<px_count; ++i)
00074 {
00075 if (reverse)
00076 px[i] = 0xFFFF - (unsigned short)( (float)px[i]/max1*max2 );
00077 else
00078 px[i] = (unsigned short)( (float)px[i]/max1*max2 );
00079 }
00080 }
00081 static void to32bits(int bits_used, void* ptr, int px_count, bool reverse)
00082 {
00083 unsigned int* px = (unsigned int*)ptr;
00084 if (bits_used >= 32)
00085 return;
00086 unsigned int max1 = (1<<bits_used)-1;
00087 unsigned int max2 = 0xFFFFFFFF;
00088 for(int i=0; i<px_count; ++i)
00089 if (reverse)
00090 px[i] = 0xFFFFFFFF - (unsigned int)( (float)px[i]/max1*max2 );
00091 else
00092 px[i] = (unsigned int)( (float)px[i]/max1*max2 );
00093 }
00094
00095 ref<Image> vl::loadDICOM(const String& path)
00096 {
00097 ref<VirtualFile> file = defFileSystem()->locateFile(path);
00098 if ( !file )
00099 {
00100 Log::error( Say("File '%s' not found.\n") << path );
00101 return NULL;
00102 }
00103 else
00104 return loadDICOM(file.get());
00105 }
00106
00108 ref<Image> vl::loadDICOM(VirtualFile* vfile)
00109 {
00110 gdcm::ImageReader reader;
00111
00112 if (!vfile->open(OM_ReadOnly))
00113 {
00114 Log::error( Say("loadDICOM: cannot open file '%s'\n") << vfile->path() );
00115 return NULL;
00116 }
00117
00118 const int bufsize = 128*1024;
00119 char buffer[bufsize];
00120 long long count = 0;
00121 std::stringstream strstr;
00122 while( (count=vfile->read(buffer, bufsize)) )
00123 strstr.write(buffer,(int)count);
00124 std::istream* istr = &strstr;
00125 reader.SetStream( *istr );
00126 if( !reader.Read() )
00127 {
00128 std::cerr << "Could not read: " << vfile->path().toStdString() << std::endl;
00129 vfile->close();
00130 return NULL;
00131 }
00132
00133 gdcm::File &file = reader.GetFile();
00134 const gdcm::Image &image = reader.GetImage();
00135
00136 #if 0
00137 printf("GDCM --- --- --- --- ---\n");
00138 image.Print(std::cout);
00139 #endif
00140
00141 unsigned int ndim = image.GetNumberOfDimensions();
00142 const unsigned int *dims = image.GetDimensions();
00143 gdcm::PixelFormat pf = image.GetPixelFormat();
00144 const gdcm::PhotometricInterpretation &pi = image.GetPhotometricInterpretation();
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 ref<KeyValues> tags = new KeyValues;
00157 tags->set("Origin") = Say("%n %n %n") << image.GetOrigin()[0] << image.GetOrigin()[1] << image.GetOrigin()[2];
00158 tags->set("Spacing") = Say("%n %n %n") << image.GetSpacing()[0] << image.GetSpacing()[1] << image.GetSpacing()[2];
00159 tags->set("Intercept") = Say("%n") << image.GetIntercept();
00160 tags->set("Slope") = Say("%n") << image.GetSlope();
00161 tags->set("DirectionCosines") = Say("%n %n %n %n %n %n")
00162 << image.GetDirectionCosines()[0] << image.GetDirectionCosines()[1] << image.GetDirectionCosines()[2]
00163 << image.GetDirectionCosines()[3] << image.GetDirectionCosines()[4] << image.GetDirectionCosines()[5];
00164 tags->set("BitsStored") = Say("%n") << pf.GetBitsStored();
00165
00166 gdcm::DataSet& ds = file.GetDataSet();
00167
00168 {
00169 gdcm::Attribute<0x28,0x1050> win_center;
00170 const gdcm::DataElement& de = ds.GetDataElement( win_center.GetTag() );
00171 if(!de.IsEmpty())
00172 {
00173 win_center.SetFromDataElement( ds.GetDataElement( win_center.GetTag() ) );
00174 tags->set("WindowCenter") = Say("%n") << win_center.GetValue();
00175 }
00176 }
00177
00178 {
00179 gdcm::Attribute<0x28,0x1051> win_width;
00180 const gdcm::DataElement& de = ds.GetDataElement( win_width.GetTag() );
00181 if(!de.IsEmpty())
00182 {
00183 win_width.SetFromDataElement( ds.GetDataElement( win_width.GetTag() ) );
00184 tags->set("WindowWidth") = Say("%n") << win_width.GetValue();
00185 }
00186 }
00187
00188 {
00189 gdcm::Attribute<0x28,0x1052> rescale_intercept;
00190 const gdcm::DataElement& de = ds.GetDataElement( rescale_intercept.GetTag() );
00191 if(!de.IsEmpty())
00192 {
00193 rescale_intercept.SetFromDataElement( ds.GetDataElement( rescale_intercept.GetTag() ) );
00194 tags->set("RescaleIntercept") = Say("%n") << rescale_intercept.GetValue();
00195 }
00196 else
00197 tags->set("RescaleIntercept") = Say("%n") << 0;
00198 }
00199
00200 {
00201 gdcm::Attribute<0x28,0x1053> rescale_slope;
00202 const gdcm::DataElement& de = ds.GetDataElement( rescale_slope.GetTag() );
00203 if(!de.IsEmpty())
00204 {
00205 rescale_slope.SetFromDataElement( ds.GetDataElement( rescale_slope.GetTag() ) );
00206 tags->set("RescaleSlope") = Say("%n") << rescale_slope.GetValue();
00207 }
00208 else
00209 tags->set("RescaleIRescaleSlopentercept") = Say("%n") << 1;
00210 }
00211
00212 #if 0
00213 printf("TAGS --- --- --- --- ---\n");
00214 tags->print();
00215 #endif
00216
00217 int w=1,h=1,d=1;
00218 if (ndim>=1)
00219 w = dims[0];
00220 if (ndim>=2)
00221 h = dims[1];
00222 if (ndim>=3)
00223 d = dims[2];
00224 h = h * d;
00225
00226 ref<Image> img = new Image;
00227 img->setObjectName( vfile->path().toStdString().c_str() );
00228 if (pf.GetSamplesPerPixel() == 1 && pi == gdcm::PhotometricInterpretation::PALETTE_COLOR)
00229 {
00230 if (pf.GetBitsStored() <= 8)
00231 {
00232 img->allocate2D( w, h, 1, vl::IF_RGBA, vl::IT_UNSIGNED_BYTE);
00233 image.GetBuffer((char*)img->pixels());
00234
00235 const gdcm::LookupTable& lut = image.GetLUT();
00236 std::auto_ptr<char> rgba( new char[(1<<lut.GetBitSample())*4] );
00237 lut.GetBufferAsRGBA((unsigned char*)rgba.get());
00238 for(int i=w*h; i--; )
00239 {
00240 int idx = img->pixels()[i];
00241 img->pixels()[i*4+0] = rgba.get()[idx*4+0];
00242 img->pixels()[i*4+1] = rgba.get()[idx*4+1];
00243 img->pixels()[i*4+2] = rgba.get()[idx*4+2];
00244 img->pixels()[i*4+3] = rgba.get()[idx*4+3];
00245 }
00246 }
00247 else
00248 if (pf.GetBitsStored() <= 16)
00249 {
00250 img->allocate2D( w, h, 1, vl::IF_RGBA, vl::IT_UNSIGNED_BYTE);
00251 image.GetBuffer((char*)img->pixels());
00252
00253 const gdcm::LookupTable& lut = image.GetLUT();
00254 std::auto_ptr<char> rgba( new char[(1<<lut.GetBitSample())*4] );
00255 lut.GetBufferAsRGBA((unsigned char*)rgba.get());
00256 for(int i=w*h; i--; )
00257 {
00258 int idx = ((unsigned short*)img->pixels())[i];
00259 img->pixels()[i*4+0] = rgba.get()[idx*4+0];
00260 img->pixels()[i*4+1] = rgba.get()[idx*4+1];
00261 img->pixels()[i*4+2] = rgba.get()[idx*4+2];
00262 img->pixels()[i*4+3] = rgba.get()[idx*4+3];
00263 }
00264 }
00265 }
00266 else
00267 if (pf.GetSamplesPerPixel() == 1 && (pi == gdcm::PhotometricInterpretation::MONOCHROME2 || pi == gdcm::PhotometricInterpretation::MONOCHROME1))
00268 {
00269 if (pf.GetBitsStored() <= 8)
00270 {
00271 img->allocate2D( w, h, 1, vl::IF_LUMINANCE, vl::IT_UNSIGNED_BYTE);
00272 image.GetBuffer((char*)img->pixels());
00273 if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
00274 to8bits(pf.GetBitsStored(), img->pixels(), w*h, true);
00275 else
00276 to8bits(pf.GetBitsStored(), img->pixels(), w*h, false);
00277 }
00278 else
00279 if (pf.GetBitsStored() <= 16)
00280 {
00281 img->allocate2D( w, h, 1, vl::IF_LUMINANCE, vl::IT_UNSIGNED_SHORT);
00282 image.GetBuffer((char*)img->pixels());
00283 if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
00284 to16bits(pf.GetBitsStored(), img->pixels(), w*h, true);
00285 else
00286 to16bits(pf.GetBitsStored(), img->pixels(), w*h, false);
00287 }
00288 else
00289 if (pf.GetBitsStored() <= 32)
00290 {
00291 img->allocate2D( w, h, 1, vl::IF_LUMINANCE, vl::IT_UNSIGNED_INT);
00292 image.GetBuffer((char*)img->pixels());
00293 if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
00294 to32bits(pf.GetBitsStored(), img->pixels(), w*h, true);
00295 else
00296 to32bits(pf.GetBitsStored(), img->pixels(), w*h, false);
00297 }
00298 }
00299 else
00300 if (pf.GetSamplesPerPixel() == 3)
00301 {
00302 if (pf.GetBitsStored() <= 8)
00303 {
00304 img->allocate2D( w, h, 1, vl::IF_RGB, vl::IT_UNSIGNED_BYTE);
00305 image.GetBuffer((char*)img->pixels());
00306 if (image.GetPlanarConfiguration())
00307 {
00308 int slice_pix_count = w*h/d;
00309 std::auto_ptr<char> px ( new char[slice_pix_count*3] );
00310 for(int slice=0; slice<d; ++slice)
00311 {
00312 char* red = (char* )img->pixels() + slice*slice_pix_count*3 + 0;
00313 char* green = (char* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count;
00314 char* blue = (char* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count*2;
00315 for(int i=0; i<slice_pix_count*3; i+=3, ++red, ++green, ++blue)
00316 {
00317 px.get()[i+0] = *red;
00318 px.get()[i+1] = *green;
00319 px.get()[i+2] = *blue;
00320 }
00321 memcpy(img->pixels(), px.get(), img->requiredMemory());
00322 }
00323 }
00324 }
00325 else
00326 if (pf.GetBitsStored() <= 16)
00327 {
00328 img->allocate2D( w, h, 1, vl::IF_RGB, vl::IT_UNSIGNED_SHORT);
00329 image.GetBuffer((char*)img->pixels());
00330 if (image.GetPlanarConfiguration())
00331 {
00332 int slice_pix_count = w*h/d;
00333 std::auto_ptr<short> px ( new short[slice_pix_count*3] );
00334 for(int slice=0; slice<d; ++slice)
00335 {
00336 short* red = (short* )img->pixels() + slice*slice_pix_count*3 + 0;
00337 short* green = (short* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count;
00338 short* blue = (short* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count*2;
00339 for(int i=0; i<slice_pix_count*3; i+=3, ++red, ++green, ++blue)
00340 {
00341 px.get()[i+0] = *red;
00342 px.get()[i+1] = *green;
00343 px.get()[i+2] = *blue;
00344 }
00345 memcpy(img->pixels(), px.get(), img->requiredMemory());
00346 }
00347 }
00348 }
00349 else
00350 {
00351 vl::Log::error("DICOM format not supported: SamplesPerPixel = 3, BitsStored > 16.\n\n");
00352 image.Print(std::cout);
00353 vfile->close();
00354 return NULL;
00355 }
00356 }
00357 else
00358 {
00359 vl::Log::error("DICOM format not supported.\n");
00360 image.Print(std::cout);
00361 vfile->close();
00362 return NULL;
00363 }
00364 vfile->close();
00365 img->setHeight(h/d);
00366 img->setDepth(d>1?d:0);
00367 img->setTags(tags.get());
00368 img->flipVertically();
00369 return img;
00370 }
00371
00372 bool vl::isDICOM(VirtualFile* file)
00373 {
00374 file->open(OM_ReadOnly);
00375 file->seekSet(128);
00376 char signature[] = { 'D', 'I', 'C', 'M' };
00377 char buf[] = { 0, 0, 0, 0 };
00378 file->read(buf,4);
00379 bool ok = memcmp(buf,signature,4) == 0;
00380 file->close();
00381 return ok;
00382 }
00383
00384 bool vl::saveDICOM(const Image* src, const String& path)
00385 {
00386 return saveDICOM(src, new DiskFile(path));
00387 }
00388
00389 bool vl::saveDICOM(const Image* src, VirtualFile* fout)
00390 {
00391 if (src->dimension()<1 || src->dimension()>3)
00392 {
00393 vl::Log::error("saveDICOM(): uncompatible image dimension().\n");
00394 return false;
00395 }
00396
00397 ref<Image> img = new Image;
00398 *img = *src;
00399
00400
00401 unsigned short bps = 0;
00402 switch(img->type())
00403 {
00404 case vl::IT_UNSIGNED_BYTE: bps = 1; break;
00405 case vl::IT_UNSIGNED_SHORT: bps = 2; break;
00406 default:
00407 vl::Log::error("saveDICOM(): unsupported image type().\n");
00408 return false;
00409 }
00410
00411
00412 unsigned short spp = 0;
00413 switch(img->format())
00414 {
00415 case vl::IF_ALPHA: spp = 1; break;
00416 case vl::IF_LUMINANCE: spp = 1; break;
00417 case vl::IF_LUMINANCE_ALPHA: spp = 1; break;
00418 case vl::IF_RGB: spp = 3; break;
00419 case vl::IF_BGR: spp = 3; break;
00420 case vl::IF_RGBA: spp = 3; break;
00421 case vl::IF_BGRA: spp = 3; break;
00422 default:
00423 vl::Log::error("saveDICOM(): unsupported image format().\n");
00424 return false;
00425 }
00426
00427 if (img->format() == vl::IF_LUMINANCE_ALPHA)
00428 img = img->convertFormat(vl::IF_LUMINANCE);
00429 if (img->format() == vl::IF_BGR)
00430 img = img->convertFormat(vl::IF_RGB);
00431 if (img->format() == vl::IF_BGRA || img->format() == vl::IF_RGBA)
00432 img = img->convertFormat(vl::IF_RGB);
00433
00434 gdcm::SmartPointer<gdcm::Image> im = new gdcm::Image;
00435 im->SetNumberOfDimensions( img->dimension() );
00436 int dims[] = { img->width(), img->height()?img->height():1, img->depth()?img->depth():1 };
00437 for(int i=0; i<img->dimension(); ++i)
00438 im->SetDimension( i, dims[i] );
00439
00440 im->GetPixelFormat().SetScalarType(bps==8?gdcm::PixelFormat::INT8:gdcm::PixelFormat::INT16);
00441 im->GetPixelFormat().SetBitsAllocated(bps*8);
00442 im->GetPixelFormat().SetBitsStored(bps*8);
00443 im->GetPixelFormat().SetHighBit(bps*8-1);
00444 im->GetPixelFormat().SetSamplesPerPixel(spp);
00445 gdcm::PhotometricInterpretation pi = spp==1?gdcm::PhotometricInterpretation::MONOCHROME2:gdcm::PhotometricInterpretation::RGB;
00446 im->SetPhotometricInterpretation(pi);
00447
00448
00449 unsigned long len = im->GetBufferLength();
00450 unsigned long req = img->requiredMemory();
00451 if( len != req )
00452 {
00453
00454 vl::Log::error("saveDICOM(): buffer size computation error.\n");
00455 VL_TRAP()
00456 return false;
00457 }
00458
00459 img->flipVertically();
00460
00461 gdcm::DataElement pixeldata( gdcm::Tag(0x7fe0,0x0010) );
00462 pixeldata.SetByteValue( (const char*)img->pixels(), len );
00463 im->SetDataElement( pixeldata );
00464
00465 gdcm::ImageWriter w;
00466 w.SetImage( *im );
00467 std::ostringstream ostr;
00468 w.SetStream(ostr);
00469 if( !w.Write() )
00470 {
00471 vl::Log::error("saveDICOM(): error writing stream.\n");
00472 return false;
00473 }
00474 fout->open(OM_WriteOnly);
00475 int bcount = (int)fout->write(ostr.str().c_str(), ostr.str().length());
00476 fout->close();
00477 if( bcount != (int)ostr.str().length() )
00478 {
00479 vl::Log::error("saveDICOM(): error writing file.\n");
00480 return false;
00481 }
00482
00483 return true;
00484 }