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>
38 #include <gdcmWriter.h>
44 #include "itkImageRegionIterator.h"
45 #include "itkMetaImageIO.h"
46 #include "itkMetaDataDictionary.h"
53 //-----------------------------------------------------------
55 //-----------------------------------------------------------
56 template<class args_info_type>
57 GateSimulation2DicomGenericFilter<args_info_type>::GateSimulation2DicomGenericFilter()
64 //-----------------------------------------------------------
66 //-----------------------------------------------------------
67 template<class args_info_type>
68 void GateSimulation2DicomGenericFilter<args_info_type>::Update()
70 // Read the Dimension and PixelType
72 std::string PixelType;
73 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
77 if(Dimension==2) UpdateWithDim<2>(PixelType);
78 else if(Dimension==3) UpdateWithDim<3>(PixelType);
79 // else if (Dimension==4)UpdateWithDim<4>(PixelType);
81 std::cout<<"Error, Only for 2 or 3 Dimensions!!!"<<std::endl ;
86 //-------------------------------------------------------------------
87 // Update with the number of dimensions
88 //-------------------------------------------------------------------
89 template<class args_info_type>
90 template<unsigned int Dimension>
92 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
94 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
96 if(PixelType == "short") {
97 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
98 UpdateWithDimAndPixelType<Dimension, signed short>();
100 else if(PixelType == "unsigned_short"){
101 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
102 UpdateWithDimAndPixelType<Dimension, unsigned short>();
105 else if (PixelType == "unsigned_char") {
106 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
107 UpdateWithDimAndPixelType<Dimension, unsigned char>();
110 // else if (PixelType == "char"){
111 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
112 // UpdateWithDimAndPixelType<Dimension, signed char>();
114 else if (PixelType == "double") {
115 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
116 UpdateWithDimAndPixelType<Dimension, double>();
119 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
120 UpdateWithDimAndPixelType<Dimension, float>();
124 //-------------------------------------------------------------------
125 // Update with the number of dimensions and the pixeltype read from
126 // the dicom files. The MHD files may be resampled to match the
127 // dicom spacing (and number of slices). Rounding errors in resampling
128 // are handled by removing files when generating the output dicom
130 //-------------------------------------------------------------------
131 template<class args_info_type>
132 template <unsigned int Dimension, class PixelType>
134 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
137 #if GDCM_MAJOR_VERSION >= 2
139 typedef itk::Image<PixelType, Dimension> InputImageType;
140 typedef itk::Image<PixelType, Dimension> OutputImageType;
141 typedef itk::ImageFileReader< InputImageType > ReaderType;
142 typedef itk::ImageSeriesReader< InputImageType > ReaderSeriesType;
143 typedef itk::GDCMImageIO ImageIOType;
146 // Read Dicom model file
147 typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
148 ImageIOType::Pointer gdcmIO = ImageIOType::New();
149 std::string filename_out = m_ArgsInfo.outputFilename_arg;
150 gdcmIO->LoadPrivateTagsOn();
151 gdcmIO->KeepOriginalUIDOn();
152 typename ReaderSeriesType::FileNamesContainer fileNames;
153 fileNames.push_back(m_ArgsInfo.inputModelFilename_arg);
154 readerSeries->SetImageIO( gdcmIO );
155 readerSeries->SetFileNames( fileNames );
157 readerSeries->Update();
158 } catch (itk::ExceptionObject &excp) {
159 std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
160 std::cerr << excp << std::endl;
164 // Read the input (MHD file)
165 typedef typename InputImageType::RegionType RegionType;
166 typedef typename RegionType::SizeType SizeType;
167 typedef itk::ImageFileReader<InputImageType> InputReaderType;
168 typename InputReaderType::Pointer volumeReader = InputReaderType::New();
169 volumeReader->SetFileName( m_InputFileName);
170 volumeReader->Update();
171 typename InputImageType::Pointer input = volumeReader->GetOutput();
174 //Read the metadata informations in the mhd
176 tObj.AddUserField("NumberOfEnergyWindows", MET_STRING);
177 tObj.AddUserField("NbOfHeads", MET_STRING);
178 tObj.AddUserField("NbOfProjections", MET_STRING);
179 tObj.Read(m_InputFileName.c_str());
180 char *p_energyNumber, *p_headNumber, *p_rotationNumber;
181 int energyNumber, headNumber, rotationNumber;
182 p_energyNumber = static_cast<char*>(tObj.GetUserField("NumberOfEnergyWindows"));
183 p_headNumber = static_cast<char*>(tObj.GetUserField("NbOfHeads"));
184 p_rotationNumber = static_cast<char*>(tObj.GetUserField("NbOfProjections"));
185 energyNumber = atoi(p_energyNumber);
186 headNumber = atoi(p_headNumber);
187 rotationNumber = atoi(p_rotationNumber);
188 for (unsigned int i=0; i<energyNumber; ++i) {
189 std::string tagEnergyName("EnergyWindow_");
190 tagEnergyName += NumberToString(i).c_str();
191 std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
192 tagEnergyThreshold += "_threshold";
193 tagEnergyUphold += "_uphold";
194 tObj.AddUserField(tagEnergyName.c_str(), MET_STRING);
195 tObj.AddUserField(tagEnergyThreshold.c_str(), MET_STRING);
196 tObj.AddUserField(tagEnergyUphold.c_str(), MET_STRING);
198 tObj.Read(m_InputFileName.c_str());
199 std::vector<char*> p_EnergyWindowName(energyNumber), p_EnergyWindowThreshold(energyNumber), p_EnergyWindowUphold(energyNumber);
200 std::vector<int> energyWindowThreshold(energyNumber), energyWindowUphold(energyNumber);
201 for (unsigned int i=0; i<energyNumber; ++i) {
202 std::string tagEnergyName("EnergyWindow_");
203 tagEnergyName += NumberToString(i).c_str();
204 std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
205 tagEnergyThreshold += "_threshold";
206 tagEnergyUphold += "_uphold";
207 p_EnergyWindowName[i] = static_cast<char*>(tObj.GetUserField(tagEnergyName.c_str()));
208 p_EnergyWindowThreshold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyThreshold.c_str()));
209 p_EnergyWindowUphold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyUphold.c_str()));
210 energyWindowThreshold[i] = atoi(p_EnergyWindowThreshold[i]);
211 energyWindowUphold[i] = atoi(p_EnergyWindowUphold[i]);
214 //TODO if the input size/spacing and dicom model ones are different
216 // Create a new mhd image with the correct dicom order slices
217 typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
218 mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
219 mhdCorrectOrder->Allocate();
220 unsigned int zAxis(0); //z value for the input mhd image
221 for (unsigned int energy = 0; energy < energyNumber; ++energy) {
222 for (unsigned int head = 0; head < headNumber; ++head) {
223 for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
224 std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
226 typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
227 startIteratorIndexCorrectOrder[0] = 0;
228 startIteratorIndexCorrectOrder[1] = 0;
229 startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
231 typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
232 startIteratorIndexOriginalOrder[0] = 0;
233 startIteratorIndexOriginalOrder[1] = 0;
234 startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
236 typename InputImageType::SizeType regionSizeIterator;
237 regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
238 regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
239 regionSizeIterator[2] = 1;
241 typename InputImageType::RegionType regionIteratorCorrectOrder;
242 regionIteratorCorrectOrder.SetSize(regionSizeIterator);
243 regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
245 typename InputImageType::RegionType regionIteratorOriginalOrder;
246 regionIteratorOriginalOrder.SetSize(regionSizeIterator);
247 regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
249 itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
250 itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
251 while(!CorrectOrderIterator.IsAtEnd()) {
252 CorrectOrderIterator.Set(OriginalOrderIterator.Get());
253 ++CorrectOrderIterator;
254 ++OriginalOrderIterator;
263 // update output dicom keys/tags
264 // string for distinguishing items inside sequence:
265 const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE");
266 std::string tempString = ITEM_ENCAPSULATE_STRING + "01";
267 typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
268 typename ReaderSeriesType::DictionaryArrayType outputArray;
269 typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
270 CopyDictionary (*inputDict, *outputDict);
271 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(energyNumber));
272 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
273 outputArray.push_back(outputDict);
276 // Output directory and filenames
277 typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType> WriterSerieType;
278 typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
279 writerSerie->SetInput( mhdCorrectOrder );
280 writerSerie->SetImageIO( gdcmIO );
281 typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
282 fileNamesOutput.push_back(filename_out);
283 writerSerie->SetFileNames( fileNamesOutput );
284 writerSerie->SetMetaDataDictionaryArray(&outputArray);
289 if (m_ArgsInfo.verbose_flag)
290 std::cout << writerSerie << std::endl;
291 writerSerie->Update();
292 } catch( itk::ExceptionObject & excp ) {
293 std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
294 std::cerr << excp << std::endl;
298 //Write sequence dicom tag with gdcm
300 reader.SetFileName( fileNamesOutput[0].c_str() );
302 gdcm::File &file = reader.GetFile();
303 gdcm::DataSet &ds2 = file.GetDataSet();
304 const unsigned int ptr_len = 42;
305 char *ptr = new char[ptr_len];
306 memset(ptr,0,ptr_len);
310 gdcm::SmartPointer<gdcm::SequenceOfItems> rotationSq = new gdcm::SequenceOfItems();
311 rotationSq->SetLengthToUndefined();
312 // Create a dataelement
313 gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
314 rotationDE.SetVR( gdcm::VR::US );
315 char essai = (char)rotationNumber;
316 char *p_essai = &essai;
317 rotationDE.SetByteValue(p_essai, 1);
319 gdcm::Item rotationIt;
320 rotationIt.SetVLToUndefined();
321 gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
322 rotationDS.Insert(rotationDE);
323 rotationSq->AddItem(rotationIt);
324 // Insert sequence into data set
325 gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) );
326 rotationDEParent.SetVR(gdcm::VR::SQ);
327 rotationDEParent.SetValue(*rotationSq);
328 rotationDEParent.SetVLToUndefined();
331 gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) );
332 energyDEParent.SetVR(gdcm::VR::SQ);
334 gdcm::SmartPointer<gdcm::SequenceOfItems> energySq = new gdcm::SequenceOfItems();
335 energySq->SetLengthToUndefined();
336 for (unsigned int i=0; i<energyNumber; ++i) {
337 gdcm::SmartPointer<gdcm::SequenceOfItems> energyThresholdSq = new gdcm::SequenceOfItems();
338 energyThresholdSq->SetLengthToUndefined();
339 // Create a dataelement
340 gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) );
341 gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) );
342 energyThresholdDE.SetVR( gdcm::VR::DS );
343 energyUpholdDE.SetVR( gdcm::VR::DS );
344 energyThresholdDE.SetByteValue(p_EnergyWindowThreshold[i], (uint32_t)strlen(p_EnergyWindowThreshold[i]));
345 energyUpholdDE.SetByteValue(p_EnergyWindowUphold[i], (uint32_t)strlen(p_EnergyWindowUphold[i]));
347 gdcm::Item energyThresholdIt;
348 energyThresholdIt.SetVLToUndefined();
349 gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
350 energyThresholdDS.Insert(energyThresholdDE);
351 energyThresholdDS.Insert(energyUpholdDE);
352 energyThresholdSq->AddItem(energyThresholdIt);
353 // Insert sequence into data set
354 gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) );
355 energyThresholdDEParent.SetVR(gdcm::VR::SQ);
356 energyThresholdDEParent.SetValue(*energyThresholdSq);
357 energyThresholdDEParent.SetVLToUndefined();
358 // Create a dataelement
359 gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) );
360 energyNameDE.SetVR( gdcm::VR::SH );
361 energyNameDE.SetByteValue(p_EnergyWindowName[i], (uint32_t)strlen(p_EnergyWindowName[i]));
364 energyIt.SetVLToUndefined();
365 gdcm::DataSet &energyDS = energyIt.GetNestedDataSet();
366 energyDS.Insert(energyNameDE);
367 energyDS.Insert(energyThresholdDEParent);
368 energySq->AddItem(energyIt);
370 // Insert sequence into data set
371 energyDEParent.SetValue(*energySq);
372 energyDEParent.SetVLToUndefined();
373 ds2.Insert(energyDEParent);
374 ds2.Insert(rotationDEParent);
378 w.SetFileName( fileNamesOutput[0].c_str() );
382 std::cout << "Use GDCM2" << std::endl;
385 //---------------------------------------------------------------------------
388 //---------------------------------------------------------------------------
389 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
391 typedef itk::MetaDataDictionary DictionaryType;
393 DictionaryType::ConstIterator itr = fromDict.Begin();
394 DictionaryType::ConstIterator end = fromDict.End();
395 typedef itk::MetaDataObject< std::string > MetaDataStringType;
399 itk::MetaDataObjectBase::Pointer entry = itr->second;
401 MetaDataStringType::Pointer entryvalue =
402 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
405 std::string tagkey = itr->first;
406 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
407 itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
412 //---------------------------------------------------------------------------
415 //---------------------------------------------------------------------------
416 template <typename T> std::string NumberToString ( T Number )
418 std::ostringstream ss;
422 //---------------------------------------------------------------------------
426 #endif //#define clitkGateSimulation2DicomGenericFilter_txx