1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18 #ifndef clitkImage2DicomDoseGenericFilter_txx
19 #define clitkImage2DicomDoseGenericFilter_txx
21 /* =================================================
22 * @file clitkImage2DicomDosemGenericFilter.txx
28 ===================================================*/
32 #include "clitkImage2DicomDoseGenericFilter.h"
33 #include "clitkCommon.h"
34 #include "itkMinimumMaximumImageCalculator.h"
35 #include "itkImageRegionIterator.h"
36 #include "itkImageRegionConstIterator.h"
37 #if GDCM_MAJOR_VERSION >= 2
38 #include "gdcmUIDGenerator.h"
39 #include <gdcmImageHelper.h>
40 #include <gdcmAttribute.h>
41 #include <gdcmReader.h>
42 #include <gdcmWriter.h>
48 //#include "gdcmBase.h"
49 //#include "gdcmDocEntry.h"
57 //-----------------------------------------------------------
59 //-----------------------------------------------------------
60 template<class args_info_type>
61 Image2DicomDoseGenericFilter<args_info_type>::Image2DicomDoseGenericFilter()
68 //-----------------------------------------------------------
70 //-----------------------------------------------------------
71 template<class args_info_type>
72 void Image2DicomDoseGenericFilter<args_info_type>::Update()
74 // Read the Dimension and PixelType
76 std::string PixelType;
77 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
81 if(Dimension==2) UpdateWithDim<2>(PixelType);
82 else if(Dimension==3) UpdateWithDim<3>(PixelType);
83 // else if (Dimension==4)UpdateWithDim<4>(PixelType);
85 std::cout<<"Error, Only for 2 or 3 Dimensions!!!"<<std::endl ;
90 //-------------------------------------------------------------------
91 // Update with the number of dimensions
92 //-------------------------------------------------------------------
93 template<class args_info_type>
94 template<unsigned int Dimension>
96 Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
98 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
100 if(PixelType == "short") {
101 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
102 UpdateWithDimAndPixelType<Dimension, signed short>();
104 else if(PixelType == "unsigned_short"){
105 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
106 UpdateWithDimAndPixelType<Dimension, unsigned short>();
109 else if (PixelType == "unsigned_char") {
110 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
111 UpdateWithDimAndPixelType<Dimension, unsigned char>();
114 // else if (PixelType == "char"){
115 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
116 // UpdateWithDimAndPixelType<Dimension, signed char>();
118 else if (PixelType == "double") {
119 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
120 UpdateWithDimAndPixelType<Dimension, double>();
123 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
124 UpdateWithDimAndPixelType<Dimension, float>();
128 //-------------------------------------------------------------------
129 // Update with the number of dimensions and the pixeltype read from
130 // the dicom files. The MHD files may be resampled to match the
131 // dicom spacing (and number of slices). Rounding errors in resampling
132 // are handled by removing files when generating the output dicom
134 //-------------------------------------------------------------------
135 template<class args_info_type>
136 template <unsigned int Dimension, class PixelType>
138 Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
141 #if GDCM_MAJOR_VERSION >= 2
143 typedef itk::Image<PixelType, Dimension> InputImageType;
144 typedef unsigned short int OutputPixelType;
145 typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
146 typedef itk::ImageFileReader< InputImageType > ReaderType;
147 typedef itk::ImageSeriesReader< InputImageType > ReaderSeriesType;
148 typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType> WriterSerieType;
149 typedef itk::ImageRegionIterator< InputImageType > IteratorType;
150 typedef itk::MinimumMaximumImageCalculator <InputImageType> ImageCalculatorFilterType;
151 typedef itk::GDCMImageIO ImageIOType;
153 //-----------------------------------------------------------------------------
154 // opening image input file
155 typename ReaderType::Pointer reader = ReaderType::New();
156 const char * filename = m_ArgsInfo.input_arg;
157 reader->SetFileName( filename );
159 typename InputImageType::Pointer image = reader->GetOutput();
160 std::ostringstream value;
162 // Read Dicom model file
163 typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
164 ImageIOType::Pointer gdcmIO = ImageIOType::New();
165 std::string filename_out = m_ArgsInfo.output_arg;
166 gdcmIO->LoadPrivateTagsOn();
167 gdcmIO->KeepOriginalUIDOn();
168 typename ReaderSeriesType::FileNamesContainer fileNames;
169 fileNames.push_back(m_ArgsInfo.DicomInputFile_arg);
170 readerSeries->SetImageIO( gdcmIO );
171 readerSeries->SetFileNames( fileNames );
173 readerSeries->Update();
174 } catch (itk::ExceptionObject &excp) {
175 std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
176 std::cerr << excp << std::endl;
179 // update output dicom keys/tags
180 typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
181 typename ReaderSeriesType::DictionaryArrayType outputArray;
182 typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
183 CopyDictionary (*inputDict, *outputDict);
186 typename InputImageType::PointType origin = image->GetOrigin();
188 value<<origin[0]<<'\\'<<origin[1]<<'\\'<<origin[2];
189 itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0032", value.str());
193 typename InputImageType::DirectionType direction = image->GetDirection();
195 value<<direction[0][0]<<'\\'<<direction[0][1]<<'\\'<<direction[0][2]<<'\\'<<direction[1][0]<<'\\'<<direction[1][1]<<'\\'<<direction[1][2];
196 itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0037", value.str());
200 typename InputImageType::SizeType imageSize = image->GetLargestPossibleRegion().GetSize();
202 int NbCols=imageSize[0]; // col
203 int NbRows=imageSize[1]; // row
204 int NbFrames=imageSize[2]; // frame
205 itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0008", NumberToString(NbFrames));
206 itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0010", NumberToString(NbRows));
207 itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0011", NumberToString(NbCols));
213 typename InputImageType::SpacingType Spacing = image->GetSpacing();
215 value<<Spacing[0]<<'\\'<<Spacing[1];
216 itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0030", value.str());
219 itk::EncapsulateMetaData<std::string>(*outputDict, "0018|0050", value.str());
226 for (int i=1; i<NbFrames ; i++){
231 itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000c", value.str());
235 typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
236 imageCalculatorFilter->SetImage(image);
237 imageCalculatorFilter->ComputeMaximum();
238 float highestValue=imageCalculatorFilter->GetMaximum();
239 double doseScaling = highestValue/(pow(2,16)-1);
242 itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000e", value.str());
246 typename OutputImageType::Pointer imageOutput = OutputImageType::New();
247 imageOutput->SetRegions(image->GetLargestPossibleRegion());
248 imageOutput->SetSpacing(image->GetSpacing());
249 imageOutput->SetOrigin(image->GetOrigin());
250 imageOutput->SetDirection(image->GetDirection());
251 imageOutput->Allocate();
253 typename itk::ImageRegionConstIterator<InputImageType> imageIterator(image,image->GetLargestPossibleRegion());
254 typename itk::ImageRegionIterator<OutputImageType> imageOutputIterator(imageOutput,imageOutput->GetLargestPossibleRegion());
255 while(!imageIterator.IsAtEnd()) {
256 // Set the current pixel to white
257 imageOutputIterator.Set((signed short int)(imageIterator.Get()/doseScaling));
260 ++imageOutputIterator;
263 // Output directory and filenames
264 typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
265 outputArray.push_back(outputDict);
266 writerSerie->SetInput( imageOutput );
267 writerSerie->SetImageIO( gdcmIO );
268 typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
269 fileNamesOutput.push_back(filename_out);
270 writerSerie->SetFileNames( fileNamesOutput );
271 writerSerie->SetMetaDataDictionaryArray(&outputArray);
275 if (m_ArgsInfo.verbose_flag)
276 std::cout << writerSerie << std::endl;
277 writerSerie->Update();
278 } catch( itk::ExceptionObject & excp ) {
279 std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
280 std::cerr << excp << std::endl;
283 //Read sequence dicom tag with gdcm
284 gdcm::Reader readerTemplateGDCM;
285 readerTemplateGDCM.SetFileName( fileNames[0].c_str() );
286 readerTemplateGDCM.Read();
287 gdcm::File &fileTemplate = readerTemplateGDCM.GetFile();
288 gdcm::DataSet &dsTemplate = fileTemplate.GetDataSet();
289 const unsigned int ptr_len = 42;
290 char *ptrTemplate = new char[ptr_len];
291 memset(ptrTemplate,0,ptr_len);
293 const gdcm::DataElement &referenceRTPlanSq = dsTemplate.GetDataElement(gdcm::Tag(0x300c, 0x02));
295 //Create the Dose Grid Scaling data element (ITK 4.13 do not take into account - it works well with ITK 4.5.1)
296 gdcm::DataElement deDoseGridScaling( gdcm::Tag(0x3004,0x0e) );
297 deDoseGridScaling.SetVR( gdcm::VR::DS );
298 deDoseGridScaling.SetByteValue(NumberToString(doseScaling).c_str(), (uint32_t)strlen(NumberToString(doseScaling).c_str()));
300 //Copy/Write sequence dicom tag with gdcm
301 gdcm::Reader readerOutputGDCM;
302 readerOutputGDCM.SetFileName( fileNamesOutput[0].c_str() );
303 readerOutputGDCM.Read();
304 gdcm::File &file = readerOutputGDCM.GetFile();
305 gdcm::DataSet &dsOutput = file.GetDataSet();
307 dsOutput.Insert(referenceRTPlanSq);
308 dsOutput.Replace(deDoseGridScaling);
311 w.SetFileName( fileNamesOutput[0].c_str() );
314 //---------------------------------------------------------------------------------------
316 // The previous way of writting DICOM-RT-DOSE works only for ITK
317 // and resulting RT-DOSE files are not readable by commercial systems.
318 // The nex step is to copy again the output file with another syntax, allowing to make RT-DOSE files readable by commercial systems.
319 // see Jean-Pierre Roux comment below.
321 /*gdcm::FileHelper *fh = new gdcm::FileHelper(args_info.OutputFile_arg);
325 dataSize = fh->GetImageDataRawSize();
326 imageData = fh->GetImageDataRaw();// somewhat important : Loads the Pixels in memory !
328 fh->SetWriteModeToRaw();
329 fh->SetWriteTypeToDcmExplVR();
331 bool res = fh->Write(args_info.OutputFile_arg);
334 std::cout <<"Fail to write [" << args_info.OutputFile_arg << "]" <<std::endl;
335 else std::cout<<"\n DICOM File re-written, using the FileHelper syntax, in order to be processed by commercial systems !"<<std::endl;
340 w.SetFileName(m_ArgsInfo.output_arg);
344 //---------------------------------------------------------------------------------------
345 //---------------------------------------------------------------------------------------
346 // Jean-Pierre Roux help for DICOM-RT-DOSE writting
347 //---------------------------------------------------------------------------------------
348 //---------------------------------------------------------------------------------------
350 // de Jean-Pierre Roux <jpr@creatis.insa-lyon.fr>
351 // répondre à jpr@creatis.insa-lyon.fr
352 // à Loic Grevillot <loic.grevillot@gmail.com>
353 // cc jpr@creatis.insa-lyon.fr,
354 // Joël Schaerer <joel.schaerer@gmail.com>,
355 // David Sarrut <David.Sarrut@creatis.insa-lyon.fr>,
356 // Joel Schaerer <joel.schaerer@creatis.insa-lyon.fr>
357 // date 12 juillet 2010 12:23
358 // objet Re: DICOM RT DOSE
360 // masquer les détails 12:23 (Il y a 21 heures)
364 // J'aurais écrit à peut prèt la même chose.
365 // (Ci après un extrait de mon code -Example/ReWrite.cxx-)
367 // L'utilisation d'un FileHelper (qui ne changera rien dans ce cas précis) est une mesure de précaution, car, l'élément 7FE0|0010 peut être compressé (ce qui n'est pas le cas pour tes images), que la manière de stocker les pixels ainsi compressés était parfois un peu ... curieuse.
368 // Je décompresse, et réécrit non compressé.
370 // ============================
371 // gdcm::File *f = new gdcm::File();
372 // f->SetMaxSizeLoadEntry(0x7fffffff);
373 // f->SetFileName( fileName );
374 // bool res = f->Load();
381 // if (!f->IsReadable())
383 // std::cerr << "Sorry, not a Readable DICOM / ACR File" <<std::endl;
388 // gdcm::FileHelper *fh = new gdcm::FileHelper(f);
392 // dataSize = fh->GetImageDataRawSize();
393 // imageData = fh->GetImageDataRaw();// somewhat important : Loads the Pixels in memory !
395 // fh->SetWriteModeToRaw();
396 // fh->SetWriteTypeToDcmExplVR();
397 // res = fh->Write(outputFileName);
400 // std::cout <<"Fail to write [" << outputFileName << "]" <<std::endl;
405 //---------------------------------------------------------------------------------------
406 //---------------------------------------------------------------------------------------
407 // pour supprimer des tags:
408 //---------------------------------------------------------------------------------------
409 //---------------------------------------------------------------------------------------
410 // SOLUCE 1 QUI MARCHE POUR UN SQItem SIMPLE
413 //a= ((gdcm::SQItem*)mDCMFile)->GetDocEntry(0xfffe,0xe00d);
414 a= ((gdcm::SQItem*)mDCMFile)->GetDocEntry(0x3004,0x000e);
415 ((gdcm::SQItem*)mDCMFile)->RemoveEntry(a);
416 mDCMFile->Write (args_info.OutputFile_arg, type);
419 // SOLUCE 2 QUI MARCHE POUR UNE SeqEntry->SQItem
421 std::cout<<"\ntest correction fichier apres ecriture\n"<<std::endl;
422 gdcm::SeqEntry *seqEntry = mDCMFile->GetSeqEntry(0x300c,0x0002);
423 gdcm::SQItem* currentItem = seqEntry->GetFirstSQItem();
425 //a= currentItem->GetDocEntry(0x0008,0x1155);
426 a= currentItem->GetDocEntry(0xfffe,0xe00d);
427 currentItem->RemoveEntry(a);
428 mDCMFile->Write (args_info.OutputFile_arg, type);
432 //a=GetDocEntry(0x7fe0,0x0000);
433 //((gdcm::SQItem*)mDCMFile)->RemoveEntry(a);
435 //-----------------------------------------------------------------------------------
436 std::cout <<"\n## DICOM Image to RT DOSE application is ended..."<<std::endl;
437 std::cout <<"#########################################################\n" << std::endl;
440 std::cout << "Use GDCM2" << std::endl;
443 //---------------------------------------------------------------------------
446 //---------------------------------------------------------------------------
447 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
449 typedef itk::MetaDataDictionary DictionaryType;
451 DictionaryType::ConstIterator itr = fromDict.Begin();
452 DictionaryType::ConstIterator end = fromDict.End();
453 typedef itk::MetaDataObject< std::string > MetaDataStringType;
457 itk::MetaDataObjectBase::Pointer entry = itr->second;
459 MetaDataStringType::Pointer entryvalue =
460 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
463 std::string tagkey = itr->first;
464 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
465 itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
470 //---------------------------------------------------------------------------
472 //---------------------------------------------------------------------------
473 template <typename T> std::string NumberToString ( T Number )
475 std::ostringstream ss;
479 //---------------------------------------------------------------------------
484 #endif //#define clitkImage2DicomDoseGenericFilter_txx