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 typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
286 typename ReaderSeriesType::DictionaryArrayType outputArray;
287 typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
288 CopyDictionary (*inputDict, *outputDict);
290 std::string entryId, value;
291 entryId = "0054|0011";
292 value = NumberToString(energyNumber);
293 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", value );
294 entryId = "0054|0021";
295 value = NumberToString(headNumber);
296 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", value );
298 outputArray.push_back(outputDict);
299 // Output directory and filenames
300 //itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputFilename_arg ); // create if it doesn't exist
301 typedef itk::ImageFileWriter<OutputImageType> WriterType;
302 typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType> WriterSerieType;
303 typename WriterType::Pointer writer = WriterType::New();
304 typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
306 writer->SetInput( mhdCorrectOrder );
307 writer->SetImageIO( gdcmIO );
309 writer->SetFileName( filename_out );
310 //writer->SetMetaDataDictionary(outputDict);
311 writerSerie->SetInput( mhdCorrectOrder );
312 writerSerie->SetImageIO( gdcmIO );
313 typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
314 fileNamesOutput.push_back(filename_out);
315 writerSerie->SetFileNames( fileNamesOutput );
316 writerSerie->SetMetaDataDictionaryArray(&outputArray);
320 if (m_ArgsInfo.verbose_flag)
321 std::cout << writer << std::endl;
323 writerSerie->Update();
324 } catch( itk::ExceptionObject & excp ) {
325 std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
326 std::cerr << excp << std::endl;
329 std::cout << "Use GDCM2" << std::endl;
333 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
335 typedef itk::MetaDataDictionary DictionaryType;
337 DictionaryType::ConstIterator itr = fromDict.Begin();
338 DictionaryType::ConstIterator end = fromDict.End();
339 typedef itk::MetaDataObject< std::string > MetaDataStringType;
343 itk::MetaDataObjectBase::Pointer entry = itr->second;
345 MetaDataStringType::Pointer entryvalue =
346 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
349 std::string tagkey = itr->first;
350 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
351 itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
357 template <typename T>
358 std::string NumberToString ( T Number )
360 std::ostringstream ss;
366 #endif //#define clitkGateSimulation2DicomGenericFilter_txx