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 clitkGateSimulation2DicomGenericFilter_txx
19 #define clitkGateSimulation2DicomGenericFilter_txx
21 /* =================================================
22 * @file clitkGateSimulation2DicomGenericFilter.txx
28 ===================================================*/
32 #include "clitkResampleImageWithOptionsFilter.h"
33 #if GDCM_MAJOR_VERSION >= 2
34 #include "gdcmUIDGenerator.h"
35 #include <gdcmImageHelper.h>
36 #include <gdcmAttribute.h>
37 #include <gdcmReader.h>
43 #include "itkImageRegionIterator.h"
50 //-----------------------------------------------------------
52 //-----------------------------------------------------------
53 template<class args_info_type>
54 GateSimulation2DicomGenericFilter<args_info_type>::GateSimulation2DicomGenericFilter()
61 //-----------------------------------------------------------
63 //-----------------------------------------------------------
64 template<class args_info_type>
65 void GateSimulation2DicomGenericFilter<args_info_type>::Update()
67 // Read the Dimension and PixelType
69 std::string PixelType;
70 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
74 if(Dimension==2) UpdateWithDim<2>(PixelType);
75 else if(Dimension==3) UpdateWithDim<3>(PixelType);
76 // else if (Dimension==4)UpdateWithDim<4>(PixelType);
78 std::cout<<"Error, Only for 2 or 3 Dimensions!!!"<<std::endl ;
83 //-------------------------------------------------------------------
84 // Update with the number of dimensions
85 //-------------------------------------------------------------------
86 template<class args_info_type>
87 template<unsigned int Dimension>
89 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
91 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
93 if(PixelType == "short") {
94 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
95 UpdateWithDimAndPixelType<Dimension, signed short>();
97 else if(PixelType == "unsigned_short"){
98 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
99 UpdateWithDimAndPixelType<Dimension, unsigned short>();
102 else if (PixelType == "unsigned_char") {
103 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
104 UpdateWithDimAndPixelType<Dimension, unsigned char>();
107 // else if (PixelType == "char"){
108 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
109 // UpdateWithDimAndPixelType<Dimension, signed char>();
111 else if (PixelType == "double") {
112 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
113 UpdateWithDimAndPixelType<Dimension, double>();
116 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
117 UpdateWithDimAndPixelType<Dimension, float>();
121 //-------------------------------------------------------------------
122 // Update with the number of dimensions and the pixeltype read from
123 // the dicom files. The MHD files may be resampled to match the
124 // dicom spacing (and number of slices). Rounding errors in resampling
125 // are handled by removing files when generating the output dicom
127 //-------------------------------------------------------------------
128 template<class args_info_type>
129 template <unsigned int Dimension, class PixelType>
131 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
134 #if GDCM_MAJOR_VERSION == 2
136 typedef itk::Image<PixelType, Dimension> InputImageType;
137 typedef itk::Image<PixelType, Dimension> OutputImageType;
138 typedef itk::ImageFileReader< InputImageType > ReaderType;
139 typedef itk::ImageSeriesReader< InputImageType > ReaderSeriesType;
140 typedef itk::GDCMImageIO ImageIOType;
143 // Read Dicom model file
144 typename ReaderType::Pointer reader = ReaderType::New();
145 typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
146 ImageIOType::Pointer gdcmIO = ImageIOType::New();
147 std::string filename_out = m_ArgsInfo.outputFilename_arg;
148 gdcmIO->LoadPrivateTagsOn();
149 gdcmIO->KeepOriginalUIDOn();
150 reader->SetImageIO( gdcmIO );
151 reader->SetFileName( m_ArgsInfo.inputModelFilename_arg );
152 typename ReaderSeriesType::FileNamesContainer fileNames;
153 fileNames.push_back(m_ArgsInfo.inputModelFilename_arg);
154 readerSeries->SetImageIO( gdcmIO );
155 readerSeries->SetFileNames( fileNames );
158 readerSeries->Update();
159 } catch (itk::ExceptionObject &excp) {
160 std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
161 std::cerr << excp << std::endl;
165 // Read the input (MHD file)
166 typedef typename InputImageType::RegionType RegionType;
167 typedef typename RegionType::SizeType SizeType;
168 typedef itk::ImageFileReader<InputImageType> InputReaderType;
169 typename InputReaderType::Pointer volumeReader = InputReaderType::New();
170 volumeReader->SetFileName( m_InputFileName);
171 volumeReader->Update();
172 typename InputImageType::Pointer input = volumeReader->GetOutput();
175 //TODO if the input size/spacing and dicom model ones are different
176 /* if (input->GetLargestPossibleRegion().GetSize() != reader->GetOutput()->GetLargestPossibleRegion().GetSize()) {
178 // resampling is carried out on the fly if resolution or size between
179 // the input mhd and input dicom series is different
182 typedef clitk::ResampleImageWithOptionsFilter<InputImageType, InputImageType> ResampleImageFilterType;
183 typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
184 filter->SetInput(input);
185 filter->SetVerboseOptions(m_Verbose);
186 filter->SetGaussianFilteringEnabled(false);
187 filter->SetDefaultPixelValue(0);
189 const SizeType& input_size = input->GetLargestPossibleRegion().GetSize();
190 SizeType output_size;
191 for (unsigned int i = 0; i < Dimension; i++)
192 output_size[i] = input_size[i];
193 filter->SetOutputSize(output_size);
195 std::cout << "Warning: The image size differs between the MHD file and the input dicom series. Performing resampling with default options using mhd size as reference (for advanced options, use clitkResampleImage)." << std::endl;
196 std::cout << "MHD -> " << input->GetLargestPossibleRegion().GetSize() << std::endl;
197 std::cout << "dicom -> " << reader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
202 input = filter->GetOutput();
203 } catch( itk::ExceptionObject & excp ) {
204 std::cerr << "Error: Exception thrown while resampling!!" << std::endl;
205 std::cerr << excp << std::endl;
211 //Read the dicom file to find the energy, the head & the steps (TODO: do it on the mhd filetext)
212 gdcm::Reader hreader;
213 hreader.SetFileName(m_ArgsInfo.inputModelFilename_arg);
215 gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
216 int energyNumber, headNumber, rotationNumber; //Read the number of energy, of head and rotation steps
217 gdcm::Attribute<0x54,0x11> energyNumber_att;
218 energyNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
219 energyNumber = energyNumber_att.GetValue();
221 gdcm::Attribute<0x54,0x21> headNumber_att;
222 headNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
223 headNumber = headNumber_att.GetValue();
225 gdcm::Tag squenceTag(0x54,0x52); //rotation, read a sequence first
226 const gdcm::DataElement &sequence = hreader.GetFile().GetDataSet().GetDataElement( squenceTag );
227 gdcm::DataSet &sequenceDataSet = sequence.GetValueAsSQ()->GetItem(1).GetNestedDataSet();
228 gdcm::Tag rotationNumber_Tag(0x54,0x53);
229 const gdcm::DataElement &rotationNumber_Element = sequenceDataSet.GetDataElement( rotationNumber_Tag );
230 gdcm::Attribute<0x54,0x53> rotationNumber_att;
231 rotationNumber_att.SetFromDataElement(rotationNumber_Element);
232 rotationNumber = rotationNumber_att.GetValue();
237 // Create a new mhd image with the dicom order slices
238 typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
239 mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
240 mhdCorrectOrder->Allocate();
241 unsigned int zAxis(0); //z value for the input mhd image
242 for (unsigned int energy = 0; energy < energyNumber; ++energy) {
243 for (unsigned int head = 0; head < headNumber; ++head) {
244 for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
245 std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
247 typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
248 startIteratorIndexCorrectOrder[0] = 0;
249 startIteratorIndexCorrectOrder[1] = 0;
250 startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
252 typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
253 startIteratorIndexOriginalOrder[0] = 0;
254 startIteratorIndexOriginalOrder[1] = 0;
255 startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
257 typename InputImageType::SizeType regionSizeIterator;
258 regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
259 regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
260 regionSizeIterator[2] = 1;
262 typename InputImageType::RegionType regionIteratorCorrectOrder;
263 regionIteratorCorrectOrder.SetSize(regionSizeIterator);
264 regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
266 typename InputImageType::RegionType regionIteratorOriginalOrder;
267 regionIteratorOriginalOrder.SetSize(regionSizeIterator);
268 regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
270 itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
271 itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
272 while(!CorrectOrderIterator.IsAtEnd()) {
273 CorrectOrderIterator.Set(OriginalOrderIterator.Get());
274 ++CorrectOrderIterator;
275 ++OriginalOrderIterator;
284 // update output dicom keys/tags
285 // string for distinguishing items inside sequence:
286 const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE");
287 std::string tempString = ITEM_ENCAPSULATE_STRING + "01";
288 typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
289 typename ReaderSeriesType::DictionaryArrayType outputArray;
290 typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
291 typename ReaderSeriesType::DictionaryRawPointer sequenceDict = new typename ReaderSeriesType::DictionaryType;
292 // Dictionary for encapsulating the sequences.
293 itk::GDCMImageIO::Pointer gdcmIOStructureSetRoi = itk::GDCMImageIO::New();
294 typename ReaderSeriesType::DictionaryType &structureSetRoiSeq = gdcmIOStructureSetRoi->GetMetaDataDictionary();
295 CopyDictionary (*inputDict, *outputDict);
296 //CopyDictionary (*inputDict, *sequenceDict);
298 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(energyNumber));
299 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
301 itk::EncapsulateMetaData<std::string>(*sequenceDict, "0054|0053", NumberToString(rotationNumber));
302 //itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0053", NumberToString(rotationNumber));
303 itk::EncapsulateMetaData<typename ReaderSeriesType::DictionaryType>(structureSetRoiSeq, tempString, *sequenceDict);
304 itk::EncapsulateMetaData<typename ReaderSeriesType::DictionaryType>(*outputDict, "0054|0052",structureSetRoiSeq);
306 outputArray.push_back(outputDict);
307 // Output directory and filenames
308 //itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputFilename_arg ); // create if it doesn't exist
309 typedef itk::ImageFileWriter<OutputImageType> WriterType;
310 typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType> WriterSerieType;
311 typename WriterType::Pointer writer = WriterType::New();
312 typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
314 writer->SetInput( mhdCorrectOrder );
315 writer->SetImageIO( gdcmIO );
317 writer->SetFileName( filename_out );
318 //writer->SetMetaDataDictionary(outputDict);
319 writerSerie->SetInput( mhdCorrectOrder );
320 writerSerie->SetImageIO( gdcmIO );
321 typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
322 fileNamesOutput.push_back(filename_out);
323 writerSerie->SetFileNames( fileNamesOutput );
324 writerSerie->SetMetaDataDictionaryArray(&outputArray);
328 if (m_ArgsInfo.verbose_flag)
329 std::cout << writer << std::endl;
331 writerSerie->Update();
332 } catch( itk::ExceptionObject & excp ) {
333 std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
334 std::cerr << excp << std::endl;
337 std::cout << "Use GDCM2" << std::endl;
341 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
343 typedef itk::MetaDataDictionary DictionaryType;
345 DictionaryType::ConstIterator itr = fromDict.Begin();
346 DictionaryType::ConstIterator end = fromDict.End();
347 typedef itk::MetaDataObject< std::string > MetaDataStringType;
351 itk::MetaDataObjectBase::Pointer entry = itr->second;
353 MetaDataStringType::Pointer entryvalue =
354 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
357 std::string tagkey = itr->first;
358 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
359 itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
365 template <typename T>
366 std::string NumberToString ( T Number )
368 std::ostringstream ss;
374 #endif //#define clitkGateSimulation2DicomGenericFilter_txx