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