Visualization Library v1.0.3A lightweight C++ OpenGL middleware for 2D/3D graphics |
[Download] [Tutorials] [All Classes] [Grouped Classes] |
00001 /**************************************************************************************/ 00002 /* */ 00003 /* Visualization Library */ 00004 /* http://visualizationlibrary.org */ 00005 /* */ 00006 /* Copyright (c) 2005-2010, Michele Bosi */ 00007 /* All rights reserved. */ 00008 /* */ 00009 /* Redistribution and use in source and binary forms, with or without modification, */ 00010 /* are permitted provided that the following conditions are met: */ 00011 /* */ 00012 /* - Redistributions of source code must retain the above copyright notice, this */ 00013 /* list of conditions and the following disclaimer. */ 00014 /* */ 00015 /* - Redistributions in binary form must reproduce the above copyright notice, this */ 00016 /* list of conditions and the following disclaimer in the documentation and/or */ 00017 /* other materials provided with the distribution. */ 00018 /* */ 00019 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND */ 00020 /* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED */ 00021 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */ 00022 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR */ 00023 /* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */ 00024 /* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */ 00025 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON */ 00026 /* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */ 00027 /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS */ 00028 /* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ 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 /* for debugging purposes only 00147 const double *origin = image.GetOrigin(); 00148 unsigned int planar_conf = image.GetPlanarConfiguration(); 00149 unsigned int rows = image.GetRows(); 00150 unsigned int cols = image.GetColumns(); 00151 unsigned int buflen = image.GetBufferLength(); 00152 int swap = image.GetNeedByteSwap(); 00153 int overlays = image.GetNumberOfOverlays(); 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 // bytes per sample 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 // sample count; 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; // converted to IF_LUMINANCE 00418 case vl::IF_RGB: spp = 3; break; 00419 case vl::IF_BGR: spp = 3; break; // converted to IF_RGB 00420 case vl::IF_RGBA: spp = 3; break; // converted to IF_RGB 00421 case vl::IF_BGRA: spp = 3; break; // converted to IF_RGB 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 // im->SetTransferSyntax(gdcm::TransferSyntax::JPEGBaselineProcess1); 00448 00449 unsigned long len = im->GetBufferLength(); 00450 unsigned long req = img->requiredMemory(); 00451 if( len != req ) 00452 { 00453 // is a img->byteAlignment() issue? 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 }