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 #include <gdcmFileHelper.h>
  25
  26 #define GFOV_SPACING_TOL (1e-1)
  27
  28 //--------------------------------------------------------------------
  29 // Read Image Information
  30 void clitk::DicomRTDoseIO::ReadImageInformation()
  31 {
  32   if(m_GdcmFile)
  33     delete m_GdcmFile;
  34   m_GdcmFile = new gdcm::File;
  35
  36   int i, rc;
  37   std::string tmp;
  38   float ipp[3];
  39   int dim[3];
  40   float spacing[3];
  41   float *gfov;    /* gfov = GridFrameOffsetVector */
  42   int gfov_len;
  43   const char *gfov_str;
  44
  45   m_GdcmFile->SetMaxSizeLoadEntry (0xffff);
  46   m_GdcmFile->SetFileName (m_FileName.c_str());
  47   m_GdcmFile->SetLoadMode (0);
  48   m_GdcmFile->Load();
  49
  50   /* Modality -- better be RTSTRUCT */
  51   tmp = m_GdcmFile->GetEntryValue (0x0008, 0x0060);
  52   if (strncmp (tmp.c_str(), "RTDOSE", strlen("RTDOSE"))) {
  53     itkExceptionMacro(<< "Error.  Input file not RTDOSE: " << m_FileName);
  54   }
  55
  56   /* ImagePositionPatient */
  57   tmp = m_GdcmFile->GetEntryValue (0x0020, 0x0032);
  58   rc = sscanf (tmp.c_str(), "%f\\%f\\%f", &ipp[0], &ipp[1], &ipp[2]);
  59   if (rc != 3) {
  60     itkExceptionMacro(<< "Error parsing RTDOSE ipp.");
  61   }
  62
  63   /* Rows */
  64   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0010);
  65   rc = sscanf (tmp.c_str(), "%d", &dim[1]);
  66   if (rc != 1) {
  67     itkExceptionMacro(<< "Error parsing RTDOSE rows.");
  68   }
  69
  70   /* Columns */
  71   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0011);
  72   rc = sscanf (tmp.c_str(), "%d", &dim[0]);
  73   if (rc != 1) {
  74     itkExceptionMacro(<< "Error parsing RTDOSE columns.");
  75   }
  76
  77   /* PixelSpacing */
  78   tmp = m_GdcmFile->GetEntryValue (0x0028, 0x0030);
  79   rc = sscanf (tmp.c_str(), "%g\\%g", &spacing[0], &spacing[1]);
  80
  81   if (rc != 2) {
  82     itkExceptionMacro(<< "Error parsing RTDOSE pixel spacing.");
  83   }
  84
  85   /* GridFrameOffsetVector */
  86   tmp = m_GdcmFile->GetEntryValue (0x3004, 0x000C);
  87   gfov = 0;
  88   gfov_len = 0;
  89   gfov_str = tmp.c_str();
  90   while (1) {
  91     int len;
  92     gfov = (float*) realloc (gfov, (gfov_len + 1) * sizeof(float));
  93     rc = sscanf (gfov_str, "%g%n", &gfov[gfov_len], &len);
  94     if (rc != 1) {
  95       break;
  96     }
  97     gfov_len ++;
  98     gfov_str += len;
  99     if (gfov_str[0] == '\\') {
 100       gfov_str ++;
 101     }
 102   }
 103   dim[2] = gfov_len;
 104   if (gfov_len == 0) {
 105     itkExceptionMacro(<< "Error parsing RTDOSE gfov.");
 106   }
 107
 108   /* --- Analyze GridFrameOffsetVector --- */
 109
 110   /* (1) Make sure first element is 0. */
 111   if (gfov[0] != 0.) {
 112     itkExceptionMacro(<< "Error RTDOSE gfov[0] is not 0.");
 113   }
 114
 115   /* (2) Handle case where gfov_len == 1 (only one slice). */
 116   if (gfov_len == 1) {
 117     spacing[2] = spacing[0];
 118   }
 119
 120   /* (3) Check to make sure spacing is regular. */
 121   for (i = 1; i < gfov_len; i++) {
 122     if (i == 1) {
 123       spacing[2] = gfov[1] - gfov[0];
 124     } else {
 125       float sp = gfov[i] - gfov[i-1];
 126       if (fabs(sp - spacing[2]) > GFOV_SPACING_TOL) {
 127         itkExceptionMacro(<< "Error RTDOSE grid has irregular spacing:"
 128                           << sp << " vs " << spacing[2]);
 129       }
 130     }
 131   }
 132
 133   free (gfov);
 134
 135   /* DoseGridScaling */
 136   m_DoseScaling = 1.0;
 137   tmp = m_GdcmFile->GetEntryValue (0x3004, 0x000E);
 138   rc = sscanf (tmp.c_str(), "%f", &m_DoseScaling);
 139   /* If element doesn't exist, scaling is 1.0 */
 140
 141   // set dimension values
 142   SetNumberOfDimensions(3);
 143   for(int i=0; i<3; i++) {
 144     SetDimensions(i, dim[i]);
 145     SetSpacing(i, spacing[i]);
 146     SetOrigin(i, ipp[i]);
 147   }
 148
 149   // set other information
 150   SetByteOrderToLittleEndian();
 151   SetPixelType(itk::ImageIOBase::SCALAR);
 152   SetNumberOfComponents(1);
 153   SetComponentType(itk::ImageIOBase::FLOAT);
 154 }
 155
 156 //--------------------------------------------------------------------
 157 // Read Image Information
 158 bool clitk::DicomRTDoseIO::CanReadFile(const char* FileNameToRead)
 159 {
 160   // opening dicom input file
 161   gdcm::File dcmFile;
 162   dcmFile.SetFileName(FileNameToRead);
 163   if (!dcmFile.OpenFile())
 164     return false;
 165
 166   //Taken from plastimatch
 167   dcmFile.SetMaxSizeLoadEntry (0xffff);
 168   dcmFile.SetLoadMode (0);
 169   dcmFile.Load();
 170
 171   /* Modality -- better be RTDOSE */
 172   std::string tmp = dcmFile.GetEntryValue (0x0008, 0x0060);
 173   if (!strncmp (tmp.c_str(), "RTDOSE", strlen("RTDOSE"))) {
 174     return true;
 175   }
 176
 177   return false;
 178 }
 179
 180 //--------------------------------------------------------------------
 181 template <class T>
 182 void
 183 clitk::DicomRTDoseIO::dose_copy_raw (float *img_out, T *img_in, int nvox, float scale)
 184 {
 185   for (int i = 0; i < nvox; i++) {
 186     img_out[i] = img_in[i] * scale;
 187   }
 188 }
 189
 190 //--------------------------------------------------------------------
 191 // Read Image Content
 192 void clitk::DicomRTDoseIO::Read(void * buffer)
 193 {
 194   float *img = (float*) buffer;
 195
 196   unsigned long npix = 1;
 197   for(unsigned int i=0; i<GetNumberOfDimensions(); i++)
 198     npix *= GetDimensions(i);
 199
 200   /* PixelData */
 201   gdcm::FileHelper m_GdcmFile_helper (m_GdcmFile);
 202
 203   //size_t image_data_size = m_GdcmFile_helper.GetImageDataSize();
 204   if (strcmp (m_GdcmFile->GetPixelType().c_str(), "16U")==0) {
 205     unsigned short* image_data = (unsigned short*) m_GdcmFile_helper.GetImageData();
 206     dose_copy_raw (img, image_data, npix, m_DoseScaling);
 207   } else if (strcmp(m_GdcmFile->GetPixelType().c_str(),"32U")==0) {
 208     unsigned long* image_data = (unsigned long*) m_GdcmFile_helper.GetImageData ();
 209     dose_copy_raw (img, image_data, npix, m_DoseScaling);
 210   } else {
 211     itkExceptionMacro(<< "Error RTDOSE not type 16U and 32U (type="
 212                       << m_GdcmFile->GetPixelType().c_str() << ")");
 213   }
 214
 215   /* GCS FIX: Do I need to do something about endian-ness? */
 216
 217
 218
 219   /* Copy data to volume */
 220 //   float *img = (float*) vol->img;
 221 //   for (i = 0; i < vol->npix; i++) {
 222   //img[i] = image_data[i] * m_DoseScaling;
 223 //   }
 224
 225   delete m_GdcmFile;
 226 }
 227
 228 //--------------------------------------------------------------------
 229 // Write Image Information
 230 void clitk::DicomRTDoseIO::WriteImageInformation(bool keepOfStream)
 231 {
 232 }
 233
 234 //--------------------------------------------------------------------
 235 // Write Image Information
 236 bool clitk::DicomRTDoseIO::CanWriteFile(const char* FileNameToWrite)
 237 {
 238   return false;
 239 }
 240
 241 //--------------------------------------------------------------------
 242 // Write Image
 243 void clitk::DicomRTDoseIO::Write(const void * buffer)
 244 {
 245 }
 246