]> Creatis software - clitk.git/blob - common/clitkDicomRTDoseIO.cxx
Add code to write dicom sequence tag in gateSimulation2Dicom
[clitk.git] / 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