common/clitkDicomRTDoseIO.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18
  19 // clitk include
  20 #include "clitkDicomRTDoseIO.h"
  21 #include "clitkCommon.h"
  22
  23 // itk include
  24 #if GDCM_MAJOR_VERSION == 2
  25 // IMPLEMENTATION NOTE:
  26 // The following has been done without the use of vtkGDCMImageReader which directly
  27 // handle RTDOSE image. Another approach would have been to use gdcm::ImageReader
  28 // which also properly recognize RTDOSE
  29   #include <gdcmReader.h>
  30   #include <gdcmAttribute.h>
  31 #else
  32   #include <gdcmFileHelper.h>
  33 #endif
  34
  35 #define GFOV_SPACING_TOL (1e-1)
  36
  37 //--------------------------------------------------------------------
  38 // Read Image Information
  39 void clitk::DicomRTDoseIO::ReadImageInformation()
  40 {
  41 #if GDCM_MAJOR_VERSION < 2
  42   if(m_GdcmFile)
  43     delete m_GdcmFile;
  44   m_GdcmFile = new gdcm::File;
  45 #endif
  46
  47   int i;
  48   std::string tmp;
  49   float ipp[3];
  50   int dim[3];
  51   float spacing[3];
  52   double *gfov;    /* gfov = GridFrameOffsetVector */
  53   int gfov_len;
  54 #if GDCM_MAJOR_VERSION < 2
  55   const char *gfov_str;
  56   int rc;
  57 #endif
  58
  59 #if GDCM_MAJOR_VERSION == 2
  60   m_GdcmImageReader.SetFileName(m_FileName.c_str());
  61   m_GdcmImageReader.Read();
  62   gdcm::File* m_GdcmFile = &m_GdcmImageReader.GetFile();
  63 #else
  64   m_GdcmFile->SetMaxSizeLoadEntry (0xffff);
  65   m_GdcmFile->SetFileName (m_FileName.c_str());
  66   m_GdcmFile->SetLoadMode (0);
  67   m_GdcmFile->Load();
  68 #endif
  69
  70   /* Modality -- better be RTSTRUCT */
  71 #if GDCM_MAJOR_VERSION == 2
  72   gdcm::DataSet &ds = m_GdcmFile->GetDataSet();
  73
  74   gdcm::Attribute<0x8,0x60> at1;
  75   at1.SetFromDataSet(ds);
  76   tmp = at1.GetValue();
  77 #else
  78   tmp = m_GdcmFile->GetEntryValue (0x0008, 0x0060);
  79 #endif
  80   if (strncmp (tmp.c_str(), "RTDOSE", strlen("RTDOSE"))) {
  81     itkExceptionMacro(<< "Error.  Input file not RTDOSE: " << m_FileName);
  82   }
  83
  84   /* ImagePositionPatient */
  85 #if GDCM_MAJOR_VERSION == 2
  86   gdcm::Attribute<0x20,0x32> at2;
  87   at2.SetFromDataSet(ds);
  88   ipp[0] = at2.GetValue(0);
  89   ipp[1] = at2.GetValue(1);
  90   ipp[2] = at2.GetValue(2);
  91 #else
  92   tmp = m_GdcmFile->GetEntryValue (0x0020, 0x0032);
  93   rc = sscanf (tmp.c_str(), "%f\\%f\\%f", &ipp[0], &ipp[1], &ipp[2]);
  94   if (rc != 3) {
  95     itkExceptionMacro(<< "Error parsing RTDOSE ipp.");
  96   }
  97 #endif
  98
  99   /* Rows */
 100 #if GDCM_MAJOR_VERSION == 2
 101   gdcm::Attribute<0x28,0x10> at3;
 102   at3.SetFromDataSet(ds);
 103   dim[1] = at3.GetValue();
 104 #else
 105   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0010);
 106   rc = sscanf (tmp.c_str(), "%d", &dim[1]);
 107   if (rc != 1) {
 108     itkExceptionMacro(<< "Error parsing RTDOSE rows.");
 109   }
 110 #endif
 111
 112   /* Columns */
 113 #if GDCM_MAJOR_VERSION == 2
 114   gdcm::Attribute<0x28,0x11> at4;
 115   at4.SetFromDataSet(ds);
 116   dim[0] = at4.GetValue();
 117 #else
 118   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0011);
 119   rc = sscanf (tmp.c_str(), "%d", &dim[0]);
 120   if (rc != 1) {
 121     itkExceptionMacro(<< "Error parsing RTDOSE columns.");
 122   }
 123 #endif
 124
 125   /* PixelSpacing */
 126 #if GDCM_MAJOR_VERSION == 2
 127   gdcm::Attribute<0x28,0x30> at5;
 128   at5.SetFromDataSet(ds);
 129   spacing[0] = at5.GetValue(0);
 130   spacing[1] = at5.GetValue(1);
 131 #else
 132   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0030);
 133   rc = sscanf (tmp.c_str(), "%g\\%g", &spacing[0], &spacing[1]);
 134   if (rc != 2) {
 135     itkExceptionMacro(<< "Error parsing RTDOSE pixel spacing.");
 136   }
 137 #endif
 138
 139   /* GridFrameOffsetVector */
 140 #if GDCM_MAJOR_VERSION == 2
 141   gdcm::Attribute<0x3004,0x000C> at6;
 142   const gdcm::DataElement& de6 = ds.GetDataElement( at6.GetTag() );
 143   at6.SetFromDataElement(de6);
 144   gfov = (double*) at6.GetValues();
 145   gfov_len = at6.GetNumberOfValues();
 146 #else
 147   tmp = m_GdcmFile->GetEntryValue (0x3004, 0x000C);
 148   gfov = 0;
 149   gfov_len = 0;
 150   gfov_str = tmp.c_str();
 151   while (1) {
 152     int len;
 153     gfov = (double*) realloc (gfov, (gfov_len + 1) * sizeof(double));
 154     rc = sscanf (gfov_str, "%lf%n", &gfov[gfov_len], &len);
 155     if (rc != 1) {
 156       break;
 157     }
 158     gfov_len ++;
 159     gfov_str += len;
 160     if (gfov_str[0] == '\\') {
 161       gfov_str ++;
 162     }
 163   }
 164 #endif
 165
 166   dim[2] = gfov_len;
 167   if (gfov_len == 0) {
 168     itkExceptionMacro(<< "Error parsing RTDOSE gfov.");
 169   }
 170
 171   /* --- Analyze GridFrameOffsetVector --- */
 172
 173   /* (1) Make sure first element is 0. */
 174   if (gfov[0] != 0.) {
 175     itkExceptionMacro(<< "Error RTDOSE gfov[0] is not 0.");
 176   }
 177
 178   /* (2) Handle case where gfov_len == 1 (only one slice). */
 179   if (gfov_len == 1) {
 180     spacing[2] = spacing[0];
 181   }
 182
 183   /* (3) Check to make sure spacing is regular. */
 184   for (i = 1; i < gfov_len; i++) {
 185     if (i == 1) {
 186       spacing[2] = gfov[1] - gfov[0];
 187     } else {
 188       double sp = gfov[i] - gfov[i-1];
 189       if (fabs(sp - spacing[2]) > GFOV_SPACING_TOL) {
 190         itkExceptionMacro(<< "Error RTDOSE grid has irregular spacing:"
 191                           << sp << " vs " << spacing[2]);
 192       }
 193     }
 194   }
 195
 196 #if GDCM_MAJOR_VERSION < 2
 197   free (gfov);
 198 #endif
 199
 200   /* DoseGridScaling */
 201 #if GDCM_MAJOR_VERSION == 2
 202   gdcm::Attribute<0x3004,0x000E> at7 = { 1. } ;
 203   at7.SetFromDataSet(ds);
 204   m_DoseScaling = at7.GetValue();
 205 #else
 206   m_DoseScaling = 1.0;
 207   tmp = m_GdcmFile->GetEntryValue (0x3004, 0x000E);
 208   rc = sscanf (tmp.c_str(), "%f", &m_DoseScaling);
 209 #endif
 210   /* If element doesn't exist, scaling is 1.0 */
 211
 212   // set dimension values
 213   SetNumberOfDimensions(3);
 214   for(int i=0; i<3; i++) {
 215     SetDimensions(i, dim[i]);
 216     SetSpacing(i, spacing[i]);
 217     SetOrigin(i, ipp[i]);
 218   }
 219
 220   // set other information
 221   SetByteOrderToLittleEndian();
 222   SetPixelType(itk::ImageIOBase::SCALAR);
 223   SetNumberOfComponents(1);
 224   SetComponentType(itk::ImageIOBase::FLOAT);
 225 }
 226
 227 //--------------------------------------------------------------------
 228 // Read Image Information
 229 bool clitk::DicomRTDoseIO::CanReadFile(const char* FileNameToRead)
 230 {
 231 #if GDCM_MAJOR_VERSION == 2
 232   gdcm::Reader creader;
 233   creader.SetFileName(FileNameToRead);
 234   if (!creader.Read())
 235     return false;
 236
 237   gdcm::DataSet &ds = creader.GetFile().GetDataSet();
 238   gdcm::Attribute<0x8,0x60> at1;
 239   at1.SetFromDataSet(ds);
 240   std::string tmp = at1.GetValue();
 241 #else
 242   // opening dicom input file
 243   gdcm::File dcmFile;
 244   dcmFile.SetFileName(FileNameToRead);
 245   if (!dcmFile.OpenFile())
 246     return false;
 247
 248   //Taken from plastimatch
 249   dcmFile.SetMaxSizeLoadEntry (0xffff);
 250   dcmFile.SetLoadMode (0);
 251   dcmFile.Load();
 252
 253   /* Modality -- better be RTDOSE */
 254   std::string tmp = dcmFile.GetEntryValue (0x0008, 0x0060);
 255 #endif
 256   if (!strncmp (tmp.c_str(), "RTDOSE", strlen("RTDOSE"))) {
 257     return true;
 258   }
 259
 260   return false;
 261 }
 262
 263 //--------------------------------------------------------------------
 264 template <class T>
 265 void
 266 clitk::DicomRTDoseIO::dose_copy_raw (float *img_out, T *img_in, int nvox, float scale)
 267 {
 268   for (int i = 0; i < nvox; i++) {
 269     img_out[i] = img_in[i] * scale;
 270   }
 271 }
 272
 273 //--------------------------------------------------------------------
 274 // Read Image Content
 275 void clitk::DicomRTDoseIO::Read(void * buffer)
 276 {
 277   unsigned long npix = 1;
 278   for(unsigned int i=0; i<GetNumberOfDimensions(); i++)
 279     npix *= GetDimensions(i);
 280
 281   /* PixelData */
 282 #if GDCM_MAJOR_VERSION == 2
 283   gdcm::Image &i = m_GdcmImageReader.GetImage();
 284
 285   char* image_data = new char[i.GetBufferLength()];
 286   i.GetBuffer(image_data);
 287
 288   gdcm::DataSet& ds = m_GdcmImageReader.GetFile().GetDataSet();
 289   float *img = (float*) buffer;
 290
 291   gdcm::Attribute<0x28, 0x100> pixel_size;
 292   pixel_size.SetFromDataSet(ds);
 293   if (pixel_size.GetValue() == 16)  {
 294     unsigned short* image_data2 = (unsigned short*) image_data ;
 295     dose_copy_raw (img, image_data2, npix, m_DoseScaling);
 296   }
 297   else if (pixel_size.GetValue() == 32)  {
 298     unsigned long* image_data2 = (unsigned long*) image_data;
 299     dose_copy_raw (img, image_data2, npix, m_DoseScaling);
 300   } else {
 301     itkExceptionMacro(<< "Error RTDOSE not type 16U and 32U (type="
 302                       << pixel_size.GetValue() << ")");
 303   }
 304
 305   delete [] image_data;
 306 #else
 307   float *img = (float*) buffer;
 308
 309   gdcm::FileHelper m_GdcmFile_helper (m_GdcmFile);
 310
 311   //size_t image_data_size = m_GdcmFile_helper.GetImageDataSize();
 312   if (strcmp (m_GdcmFile->GetPixelType().c_str(), "16U")==0) {
 313     unsigned short* image_data = (unsigned short*) m_GdcmFile_helper.GetImageData();
 314     dose_copy_raw (img, image_data, npix, m_DoseScaling);
 315   } else if (strcmp(m_GdcmFile->GetPixelType().c_str(),"32U")==0) {
 316     unsigned long* image_data = (unsigned long*) m_GdcmFile_helper.GetImageData ();
 317     dose_copy_raw (img, image_data, npix, m_DoseScaling);
 318   } else {
 319     itkExceptionMacro(<< "Error RTDOSE not type 16U and 32U (type="
 320                       << m_GdcmFile->GetPixelType().c_str() << ")");
 321   }
 322   delete m_GdcmFile;
 323 #endif
 324
 325   /* GCS FIX: Do I need to do something about endian-ness? */
 326
 327
 328
 329   /* Copy data to volume */
 330 //   float *img = (float*) vol->img;
 331 //   for (i = 0; i < vol->npix; i++) {
 332   //img[i] = image_data[i] * m_DoseScaling;
 333 //   }
 334 }
 335
 336 //--------------------------------------------------------------------
 337 // Write Image Information
 338 void clitk::DicomRTDoseIO::WriteImageInformation(bool keepOfStream)
 339 {
 340 }
 341
 342 //--------------------------------------------------------------------
 343 // Write Image Information
 344 bool clitk::DicomRTDoseIO::CanWriteFile(const char* FileNameToWrite)
 345 {
 346   return false;
 347 }
 348
 349 //--------------------------------------------------------------------
 350 // Write Image
 351 void clitk::DicomRTDoseIO::Write(const void * buffer)
 352 {
 353 }
 354