Visualization Library 2.1.0

A lightweight C++ OpenGL middleware for 2D/3D graphics

VL     Star     Watch     Fork     Issue

[Download] [Tutorials] [All Classes] [Grouped Classes]
ioDICOM.cpp
Go to the documentation of this file.
1 /**************************************************************************************/
2 /* */
3 /* Visualization Library */
4 /* http://visualizationlibrary.org */
5 /* */
6 /* Copyright (c) 2005-2020, Michele Bosi */
7 /* All rights reserved. */
8 /* */
9 /* Redistribution and use in source and binary forms, with or without modification, */
10 /* are permitted provided that the following conditions are met: */
11 /* */
12 /* - Redistributions of source code must retain the above copyright notice, this */
13 /* list of conditions and the following disclaimer. */
14 /* */
15 /* - Redistributions in binary form must reproduce the above copyright notice, this */
16 /* list of conditions and the following disclaimer in the documentation and/or */
17 /* other materials provided with the distribution. */
18 /* */
19 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND */
20 /* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED */
21 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
22 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR */
23 /* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */
24 /* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
25 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON */
26 /* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
27 /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS */
28 /* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
29 /* */
30 /**************************************************************************************/
31 
32 #include "ioDICOM.hpp"
34 #include <vlCore/FileSystem.hpp>
35 #include <vlCore/glsl_math.hpp>
36 
37 #include <gdcmReader.h>
38 #include <gdcmWriter.h>
39 #include <gdcmAttribute.h>
40 #include <gdcmImageReader.h>
41 #include <gdcmImageWriter.h>
42 #include <gdcmImage.h>
43 #include <gdcmPhotometricInterpretation.h>
44 #include <gdcmImage.h>
45 #include <gdcmImageWriter.h>
46 
47 #include <memory>
48 
49 using namespace vl;
50 
51 //---------------------------------------------------------------------------
52 static void to8bits(int bits_used, void* ptr, int px_count, bool reverse)
53 {
54  unsigned char* px = (unsigned char*)ptr;
55  if (bits_used >= 8)
56  return;
57  unsigned int max1 = (1<<bits_used)-1;
58  unsigned int max2 = 0xFF;
59  for(int i=0; i<px_count; ++i)
60  if (reverse)
61  px[i] = 0xFF - (unsigned char)( (float)px[i]/max1*max2 );
62  else
63  px[i] = (unsigned char)( (float)px[i]/max1*max2 );
64 }
65 static void to16bits(int bits_used, void* ptr, int px_count, bool reverse)
66 {
67  unsigned short* px = (unsigned short*)ptr;
68 
69  if (bits_used >= 16)
70  return;
71  unsigned int max1 = (1<<bits_used)-1;
72  unsigned int max2 = 0xFFFF;
73  for(int i=0; i<px_count; ++i)
74  {
75  if (reverse)
76  px[i] = 0xFFFF - (unsigned short)( (float)px[i]/max1*max2 );
77  else
78  px[i] = (unsigned short)( (float)px[i]/max1*max2 );
79  }
80 }
81 static void to32bits(int bits_used, void* ptr, int px_count, bool reverse)
82 {
83  unsigned int* px = (unsigned int*)ptr;
84  if (bits_used >= 32)
85  return;
86  unsigned int max1 = (1<<bits_used)-1;
87  unsigned int max2 = 0xFFFFFFFF;
88  for(int i=0; i<px_count; ++i)
89  if (reverse)
90  px[i] = 0xFFFFFFFF - (unsigned int)( (float)px[i]/max1*max2 );
91  else
92  px[i] = (unsigned int)( (float)px[i]/max1*max2 );
93 }
94 //-----------------------------------------------------------------------------
96 {
98  if ( !file )
99  {
100  Log::error( Say("File '%s' not found.\n") << path );
101  return NULL;
102  }
103  else
104  return loadDICOM(file.get());
105 }
106 //-----------------------------------------------------------------------------
109 {
110  gdcm::ImageReader reader;
111 
112  if (!vfile->open(OM_ReadOnly))
113  {
114  Log::error( Say("loadDICOM: cannot open file '%s'\n") << vfile->path() );
115  return NULL;
116  }
117 
118  const int bufsize = 128*1024;
119  char buffer[bufsize];
120  long long count = 0;
121  std::stringstream strstr;
122  while( (count=vfile->read(buffer, bufsize)) )
123  strstr.write(buffer,(int)count);
124  std::istream* istr = &strstr;
125  reader.SetStream( *istr );
126  if( !reader.Read() )
127  {
128  std::cerr << "Could not read: " << vfile->path().toStdString() << std::endl;
129  vfile->close();
130  return NULL;
131  }
132 
133  gdcm::File &file = reader.GetFile();
134  const gdcm::Image &image = reader.GetImage();
135 
136  #if 0
137  printf("GDCM --- --- --- --- ---\n");
138  image.Print(std::cout);
139  #endif
140 
141  unsigned int ndim = image.GetNumberOfDimensions();
142  const unsigned int *dims = image.GetDimensions();
143  gdcm::PixelFormat pf = image.GetPixelFormat();
144  const gdcm::PhotometricInterpretation &pi = image.GetPhotometricInterpretation();
145 
146  /* for debugging purposes only
147  const double *origin = image.GetOrigin();
148  unsigned int planar_conf = image.GetPlanarConfiguration();
149  unsigned int rows = image.GetRows();
150  unsigned int cols = image.GetColumns();
151  unsigned int buflen = image.GetBufferLength();
152  int swap = image.GetNeedByteSwap();
153  int overlays = image.GetNumberOfOverlays();
154  */
155 
156  ref<KeyValues> tags = new KeyValues;
157  tags->set("Origin") = Say("%n %n %n") << image.GetOrigin()[0] << image.GetOrigin()[1] << image.GetOrigin()[2];
158  tags->set("Spacing") = Say("%n %n %n") << image.GetSpacing()[0] << image.GetSpacing()[1] << image.GetSpacing()[2];
159  tags->set("Intercept") = Say("%n") << image.GetIntercept();
160  tags->set("Slope") = Say("%n") << image.GetSlope();
161  tags->set("DirectionCosines") = Say("%n %n %n %n %n %n")
162  << image.GetDirectionCosines()[0] << image.GetDirectionCosines()[1] << image.GetDirectionCosines()[2]
163  << image.GetDirectionCosines()[3] << image.GetDirectionCosines()[4] << image.GetDirectionCosines()[5];
164  tags->set("BitsStored") = Say("%n") << pf.GetBitsStored();
165 
166  gdcm::DataSet& ds = file.GetDataSet();
167 
168  {
169  gdcm::Attribute<0x28,0x1050> win_center;
170  const gdcm::DataElement& de = ds.GetDataElement( win_center.GetTag() );
171  if(!de.IsEmpty())
172  {
173  win_center.SetFromDataElement( ds.GetDataElement( win_center.GetTag() ) );
174  tags->set("WindowCenter") = Say("%n") << win_center.GetValue();
175  }
176  }
177 
178  {
179  gdcm::Attribute<0x28,0x1051> win_width;
180  const gdcm::DataElement& de = ds.GetDataElement( win_width.GetTag() );
181  if(!de.IsEmpty())
182  {
183  win_width.SetFromDataElement( ds.GetDataElement( win_width.GetTag() ) );
184  tags->set("WindowWidth") = Say("%n") << win_width.GetValue();
185  }
186  }
187 
188  {
189  gdcm::Attribute<0x28,0x1052> rescale_intercept;
190  const gdcm::DataElement& de = ds.GetDataElement( rescale_intercept.GetTag() );
191  if(!de.IsEmpty())
192  {
193  rescale_intercept.SetFromDataElement( ds.GetDataElement( rescale_intercept.GetTag() ) );
194  tags->set("RescaleIntercept") = Say("%n") << rescale_intercept.GetValue();
195  }
196  else
197  tags->set("RescaleIntercept") = Say("%n") << 0;
198  }
199 
200  {
201  gdcm::Attribute<0x28,0x1053> rescale_slope;
202  const gdcm::DataElement& de = ds.GetDataElement( rescale_slope.GetTag() );
203  if(!de.IsEmpty())
204  {
205  rescale_slope.SetFromDataElement( ds.GetDataElement( rescale_slope.GetTag() ) );
206  tags->set("RescaleSlope") = Say("%n") << rescale_slope.GetValue();
207  }
208  else
209  tags->set("RescaleIRescaleSlopentercept") = Say("%n") << 1;
210  }
211 
212  #if 0
213  printf("TAGS --- --- --- --- ---\n");
214  tags->print();
215  #endif
216 
217  int w=1,h=1,d=1;
218  if (ndim>=1)
219  w = dims[0];
220  if (ndim>=2)
221  h = dims[1];
222  if (ndim>=3)
223  d = dims[2];
224  h = h * d;
225 
226  ref<Image> img = new Image;
227  img->setObjectName( vfile->path().toStdString().c_str() );
228  if (pf.GetSamplesPerPixel() == 1 && pi == gdcm::PhotometricInterpretation::PALETTE_COLOR)
229  {
230  if (pf.GetBitsStored() <= 8)
231  {
233  image.GetBuffer((char*)img->pixels());
234 
235  const gdcm::LookupTable& lut = image.GetLUT();
236  std::auto_ptr<char> rgba( new char[(1<<lut.GetBitSample())*4] );
237  lut.GetBufferAsRGBA((unsigned char*)rgba.get());
238  for(int i=w*h; i--; )
239  {
240  int idx = img->pixels()[i];
241  img->pixels()[i*4+0] = rgba.get()[idx*4+0];
242  img->pixels()[i*4+1] = rgba.get()[idx*4+1];
243  img->pixels()[i*4+2] = rgba.get()[idx*4+2];
244  img->pixels()[i*4+3] = rgba.get()[idx*4+3];
245  }
246  }
247  else
248  if (pf.GetBitsStored() <= 16)
249  {
251  image.GetBuffer((char*)img->pixels());
252 
253  const gdcm::LookupTable& lut = image.GetLUT();
254  std::auto_ptr<char> rgba( new char[(1<<lut.GetBitSample())*4] );
255  lut.GetBufferAsRGBA((unsigned char*)rgba.get());
256  for(int i=w*h; i--; )
257  {
258  int idx = ((unsigned short*)img->pixels())[i];
259  img->pixels()[i*4+0] = rgba.get()[idx*4+0];
260  img->pixels()[i*4+1] = rgba.get()[idx*4+1];
261  img->pixels()[i*4+2] = rgba.get()[idx*4+2];
262  img->pixels()[i*4+3] = rgba.get()[idx*4+3];
263  }
264  }
265  }
266  else
267  if (pf.GetSamplesPerPixel() == 1 && (pi == gdcm::PhotometricInterpretation::MONOCHROME2 || pi == gdcm::PhotometricInterpretation::MONOCHROME1))
268  {
269  if (pf.GetBitsStored() <= 8)
270  {
272  image.GetBuffer((char*)img->pixels());
273  if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
274  to8bits(pf.GetBitsStored(), img->pixels(), w*h, true);
275  else
276  to8bits(pf.GetBitsStored(), img->pixels(), w*h, false);
277  }
278  else
279  if (pf.GetBitsStored() <= 16)
280  {
282  image.GetBuffer((char*)img->pixels());
283  if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
284  to16bits(pf.GetBitsStored(), img->pixels(), w*h, true);
285  else
286  to16bits(pf.GetBitsStored(), img->pixels(), w*h, false);
287  }
288  else
289  if (pf.GetBitsStored() <= 32)
290  {
292  image.GetBuffer((char*)img->pixels());
293  if (pi == gdcm::PhotometricInterpretation::MONOCHROME1)
294  to32bits(pf.GetBitsStored(), img->pixels(), w*h, true);
295  else
296  to32bits(pf.GetBitsStored(), img->pixels(), w*h, false);
297  }
298  }
299  else
300  if (pf.GetSamplesPerPixel() == 3)
301  {
302  if (pf.GetBitsStored() <= 8)
303  {
304  img->allocate2D( w, h, 1, vl::IF_RGB, vl::IT_UNSIGNED_BYTE);
305  image.GetBuffer((char*)img->pixels());
306  if (image.GetPlanarConfiguration())
307  {
308  int slice_pix_count = w*h/d;
309  std::auto_ptr<char> px ( new char[slice_pix_count*3] );
310  for(int slice=0; slice<d; ++slice)
311  {
312  char* red = (char* )img->pixels() + slice*slice_pix_count*3 + 0;
313  char* green = (char* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count;
314  char* blue = (char* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count*2;
315  for(int i=0; i<slice_pix_count*3; i+=3, ++red, ++green, ++blue)
316  {
317  px.get()[i+0] = *red;
318  px.get()[i+1] = *green;
319  px.get()[i+2] = *blue;
320  }
321  memcpy(img->pixels(), px.get(), img->requiredMemory());
322  }
323  }
324  }
325  else
326  if (pf.GetBitsStored() <= 16)
327  {
329  image.GetBuffer((char*)img->pixels());
330  if (image.GetPlanarConfiguration())
331  {
332  int slice_pix_count = w*h/d;
333  std::auto_ptr<short> px ( new short[slice_pix_count*3] );
334  for(int slice=0; slice<d; ++slice)
335  {
336  short* red = (short* )img->pixels() + slice*slice_pix_count*3 + 0;
337  short* green = (short* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count;
338  short* blue = (short* )img->pixels() + slice*slice_pix_count*3 + slice_pix_count*2;
339  for(int i=0; i<slice_pix_count*3; i+=3, ++red, ++green, ++blue)
340  {
341  px.get()[i+0] = *red;
342  px.get()[i+1] = *green;
343  px.get()[i+2] = *blue;
344  }
345  memcpy(img->pixels(), px.get(), img->requiredMemory());
346  }
347  }
348  }
349  else
350  {
351  vl::Log::error("DICOM format not supported: SamplesPerPixel = 3, BitsStored > 16.\n\n");
352  image.Print(std::cout);
353  vfile->close();
354  return NULL;
355  }
356  }
357  else
358  {
359  vl::Log::error("DICOM format not supported.\n");
360  image.Print(std::cout);
361  vfile->close();
362  return NULL;
363  }
364  vfile->close();
365  img->setHeight(h/d);
366  img->setDepth(d>1?d:0);
367  img->setTags(tags.get());
368  img->flipVertically();
369  return img;
370 }
371 //---------------------------------------------------------------------------
373 {
374  file->open(OM_ReadOnly);
375  file->seekSet(128);
376  char signature[] = { 'D', 'I', 'C', 'M' };
377  char buf[] = { 0, 0, 0, 0 };
378  file->read(buf,4);
379  bool ok = memcmp(buf,signature,4) == 0;
380  file->close();
381  return ok;
382 }
383 //---------------------------------------------------------------------------
384 bool vl::saveDICOM(const Image* src, const String& path)
385 {
386  return saveDICOM(src, new DiskFile(path));
387 }
388 //-----------------------------------------------------------------------------
389 bool vl::saveDICOM(const Image* src, VirtualFile* fout)
390 {
391  if (src->dimension()<1 || src->dimension()>3)
392  {
393  vl::Log::error("saveDICOM(): uncompatible image dimension().\n");
394  return false;
395  }
396 
397  ref<Image> img = new Image;
398  *img = *src;
399 
400  // bytes per sample
401  unsigned short bps = 0;
402  switch(img->type())
403  {
404  case vl::IT_UNSIGNED_BYTE: bps = 1; break;
405  case vl::IT_UNSIGNED_SHORT: bps = 2; break;
406  default:
407  vl::Log::error("saveDICOM(): unsupported image type().\n");
408  return false;
409  }
410 
411  // sample count;
412  unsigned short spp = 0;
413  switch(img->format())
414  {
415  case vl::IF_ALPHA: spp = 1; break;
416  case vl::IF_LUMINANCE: spp = 1; break;
417  case vl::IF_LUMINANCE_ALPHA: spp = 1; break; // converted to IF_LUMINANCE
418  case vl::IF_RGB: spp = 3; break;
419  case vl::IF_BGR: spp = 3; break; // converted to IF_RGB
420  case vl::IF_RGBA: spp = 3; break; // converted to IF_RGB
421  case vl::IF_BGRA: spp = 3; break; // converted to IF_RGB
422  default:
423  vl::Log::error("saveDICOM(): unsupported image format().\n");
424  return false;
425  }
426 
427  if (img->format() == vl::IF_LUMINANCE_ALPHA)
428  img = img->convertFormat(vl::IF_LUMINANCE);
429  if (img->format() == vl::IF_BGR)
430  img = img->convertFormat(vl::IF_RGB);
431  if (img->format() == vl::IF_BGRA || img->format() == vl::IF_RGBA)
432  img = img->convertFormat(vl::IF_RGB);
433 
434  gdcm::SmartPointer<gdcm::Image> im = new gdcm::Image;
435  im->SetNumberOfDimensions( img->dimension() );
436  int dims[] = { img->width(), img->height()?img->height():1, img->depth()?img->depth():1 };
437  for(int i=0; i<img->dimension(); ++i)
438  im->SetDimension( i, dims[i] );
439 
440  im->GetPixelFormat().SetScalarType(bps==8?gdcm::PixelFormat::INT8:gdcm::PixelFormat::INT16);
441  im->GetPixelFormat().SetBitsAllocated(bps*8);
442  im->GetPixelFormat().SetBitsStored(bps*8);
443  im->GetPixelFormat().SetHighBit(bps*8-1);
444  im->GetPixelFormat().SetSamplesPerPixel(spp);
445  gdcm::PhotometricInterpretation pi = spp==1?gdcm::PhotometricInterpretation::MONOCHROME2:gdcm::PhotometricInterpretation::RGB;
446  im->SetPhotometricInterpretation(pi);
447  // im->SetTransferSyntax(gdcm::TransferSyntax::JPEGBaselineProcess1);
448 
449  unsigned long len = im->GetBufferLength();
450  unsigned long req = img->requiredMemory();
451  if( len != req )
452  {
453  // is a img->byteAlignment() issue?
454  vl::Log::error("saveDICOM(): buffer size computation error.\n");
455  VL_TRAP()
456  return false;
457  }
458 
459  img->flipVertically();
460 
461  gdcm::DataElement pixeldata( gdcm::Tag(0x7fe0,0x0010) );
462  pixeldata.SetByteValue( (const char*)img->pixels(), len );
463  im->SetDataElement( pixeldata );
464 
465  gdcm::ImageWriter w;
466  w.SetImage( *im );
467  std::ostringstream ostr;
468  w.SetStream(ostr);
469  if( !w.Write() )
470  {
471  vl::Log::error("saveDICOM(): error writing stream.\n");
472  return false;
473  }
474  fout->open(OM_WriteOnly);
475  int bcount = (int)fout->write(ostr.str().c_str(), ostr.str().length());
476  fout->close();
477  if( bcount != (int)ostr.str().length() )
478  {
479  vl::Log::error("saveDICOM(): error writing file.\n");
480  return false;
481  }
482 
483  return true;
484 }
long long read(void *buffer, long long byte_count)
Reads byte_count bytes from a file. Returns the number of bytes actually read.
Definition: VirtualFile.cpp:82
VLCORE_EXPORT FileSystem * defFileSystem()
Returns the default FileSystem used by VisualizationLibrary.
Definition: pimpl.cpp:97
int depth() const
Definition: Image.hpp:211
const T * get() const
Definition: Object.hpp:128
An abstract class representing a file.
Definition: VirtualFile.hpp:60
A simple String formatting class.
Definition: Say.hpp:124
const unsigned char * pixels() const
Raw pointer to pixels.
Definition: Image.hpp:170
void setObjectName(const char *name)
The name of the object, by default set to the object&#39;s class name in debug builds.
Definition: Object.hpp:220
VLCORE_EXPORT ref< Image > loadDICOM(VirtualFile *file)
Loads a DICOM file.
Definition: ioDICOM.cpp:108
void setHeight(int y)
Definition: Image.hpp:126
The String class implements an advanced UTF16 (Unicode BMP) string manipulation engine.
Definition: String.hpp:62
static void error(const String &message)
Use this function to provide information about run-time errors: file not found, out of memory...
Definition: Log.cpp:165
virtual ref< VirtualFile > locateFile(const String &full_path, const String &alternate_path=String()) const
Looks for a VirtualFile on the disk and in the currently active FileSystem.
Definition: FileSystem.cpp:61
int requiredMemory() const
Returns the number of bytes requested to store the image.
Definition: Image.cpp:532
virtual void close()=0
Closes the file.
const String & path() const
Returns the path of the file.
Definition: VirtualFile.hpp:98
void setDepth(int z)
Definition: Image.hpp:128
VLCORE_EXPORT bool isDICOM(VirtualFile *file)
Checks if the given file is a DICOM file.
Definition: ioDICOM.cpp:372
Visualization Library main namespace.
VLCORE_EXPORT bool saveDICOM(const Image *src, const String &path)
Writes a DICOM file.
Definition: ioDICOM.cpp:384
int height() const
Definition: Image.hpp:209
#define VL_TRAP()
Definition: checks.hpp:70
EImageDimension dimension() const
Definition: Image.cpp:372
Implements the OpenGL Shading Language convenience functions for scalar and vector operations...
int width() const
Definition: Image.hpp:207
long long write(const void *buffer, long long byte_count)
Writes byte_count bytes to a file. Returns the number of bytes actually written.
Definition: VirtualFile.cpp:90
virtual bool open(EOpenMode mode)=0
Opens the file in the specified mode.
#define NULL
Definition: OpenGLDefs.hpp:81
bool seekSet(long long offset)
Changes the current read/write position of a file.
void setTags(KeyValues *tags)
A set of key/value couples that can be used to attach extra information to an image like DICOM inform...
Definition: Image.hpp:158
ref< Image > convertFormat(EImageFormat new_format) const
Converts the format() of an image.
Definition: Image.cpp:1719
void print()
Definition: KeyValues.cpp:63
A VirtualFile that operates on regular disk files.
Definition: DiskFile.hpp:64
EImageType type() const
Definition: Image.hpp:217
Implements a generic 1d, 2d, 3d and cubemap image that can have mipmaps.
Definition: Image.hpp:54
The ref<> class is used to reference-count an Object.
Definition: Object.hpp:55
void allocate2D(int x, int y, int bytealign, EImageFormat format, EImageType type)
Definition: Image.cpp:608
std::string toStdString() const
Returns a UTF8 encoded std::string.
Definition: String.cpp:1156
String & set(const String &key)
Definition: KeyValues.hpp:53
void flipVertically()
Definition: Image.cpp:942
A set of key/value pairs usually used to associate generic information, tags, attributes etc...
Definition: KeyValues.hpp:42
EImageFormat format() const
Definition: Image.hpp:215