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 clitkPartitionEnergyWindowDicomGenericFilter_txx
19 #define clitkPartitionEnergyWindowDicomGenericFilter_txx
21 /* =================================================
22 * @file clitkPartitionEnergyWindowDicomGenericFilter.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>
39 #include <gdcmDataElement.h>
46 #include "itkImageRegionIterator.h"
47 #include "itkMetaImageIO.h"
48 #include "itkMetaDataDictionary.h"
55 //-----------------------------------------------------------
57 //-----------------------------------------------------------
58 template<class args_info_type>
59 PartitionEnergyWindowDicomGenericFilter<args_info_type>::PartitionEnergyWindowDicomGenericFilter()
66 //-----------------------------------------------------------
68 //-----------------------------------------------------------
69 template<class args_info_type>
70 void PartitionEnergyWindowDicomGenericFilter<args_info_type>::Update()
72 // Read the Dimension and PixelType
74 std::string PixelType;
75 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
79 if(Dimension==2) UpdateWithDim<2>(PixelType);
80 else if(Dimension==3) UpdateWithDim<3>(PixelType);
81 // else if (Dimension==4)UpdateWithDim<4>(PixelType);
83 std::cout<<"Error, Only for 2 or 3 Dimensions!!!"<<std::endl ;
88 //-------------------------------------------------------------------
89 // Update with the number of dimensions
90 //-------------------------------------------------------------------
91 template<class args_info_type>
92 template<unsigned int Dimension>
94 PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
96 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
98 if(PixelType == "short") {
99 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
100 UpdateWithDimAndPixelType<Dimension, signed short>();
102 else if(PixelType == "unsigned_short"){
103 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
104 UpdateWithDimAndPixelType<Dimension, unsigned short>();
107 else if (PixelType == "unsigned_char") {
108 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
109 UpdateWithDimAndPixelType<Dimension, unsigned char>();
112 // else if (PixelType == "char"){
113 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
114 // UpdateWithDimAndPixelType<Dimension, signed char>();
116 else if (PixelType == "double") {
117 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
118 UpdateWithDimAndPixelType<Dimension, double>();
121 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
122 UpdateWithDimAndPixelType<Dimension, float>();
126 //-------------------------------------------------------------------
127 // Update with the number of dimensions and the pixeltype read from
128 // the dicom files. The MHD files may be resampled to match the
129 // dicom spacing (and number of slices). Rounding errors in resampling
130 // are handled by removing files when generating the output dicom
132 //-------------------------------------------------------------------
133 template<class args_info_type>
134 template <unsigned int Dimension, class PixelType>
136 PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
139 #if GDCM_MAJOR_VERSION >= 2
141 typedef itk::Image<PixelType, Dimension> InputImageType;
142 typedef itk::Image<PixelType, Dimension> OutputImageType;
143 typedef itk::ImageSeriesReader< InputImageType > ReaderType;
144 typedef itk::ImageSeriesWriter< InputImageType, InputImageType > WriterType;
145 typedef itk::GDCMImageIO ImageIOType;
146 typedef typename InputImageType::RegionType RegionType;
147 typedef typename RegionType::SizeType SizeType;
150 // Read Dicom model file
151 typename ReaderType::Pointer reader = ReaderType::New();
152 ImageIOType::Pointer gdcmIO = ImageIOType::New();
153 gdcmIO->LoadPrivateTagsOn();
154 gdcmIO->KeepOriginalUIDOn();
155 reader->SetImageIO( gdcmIO );
156 typename ReaderType::FileNamesContainer fileNamesInput;
157 fileNamesInput.push_back(m_ArgsInfo.input_arg);
158 reader->SetFileNames( fileNamesInput );
161 } catch (itk::ExceptionObject &excp) {
162 std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
163 std::cerr << excp << std::endl;
165 typename InputImageType::Pointer input = reader->GetOutput();
167 typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
168 typename ReaderType::DictionaryArrayType outputArray;
169 typename ReaderType::DictionaryRawPointer outputDict = new typename ReaderType::DictionaryType;
170 CopyDictionary (*inputDict, *outputDict);
171 itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(1));
172 outputArray.push_back(outputDict);
175 //Read the number of slices and energy window, rotation
176 int nbSlice = input->GetLargestPossibleRegion().GetSize()[2];
177 gdcm::Reader hreader;
178 hreader.SetFileName(m_ArgsInfo.input_arg);
180 gdcm::DataSet& dsInput = hreader.GetFile().GetDataSet();
181 gdcm::Attribute<0x54,0x11> series_number_att;
182 series_number_att.SetFromDataSet(dsInput);
183 int nbEnergyWindow = series_number_att.GetValue();
186 if (dsInput.FindDataElement(gdcm::Tag(0x54, 0x52))) {
187 const gdcm::DataElement &seqRotation = dsInput.GetDataElement(gdcm::Tag(0x54, 0x52));
188 gdcm::SmartPointer<gdcm::SequenceOfItems> sqiRotation = seqRotation.GetValueAsSQ();
189 gdcm::Item &itemRotation = sqiRotation->GetItem(1);
190 const gdcm::DataElement &deRotation = itemRotation.GetDataElement(gdcm::Tag(0x54, 0x53));
191 unsigned short SamplesPerPixelRotation = *((unsigned short *) deRotation.GetByteValue()->GetPointer());
192 std::stringstream s_SamplesPerPixelRotation;
193 s_SamplesPerPixelRotation << SamplesPerPixelRotation;
194 std::string p_StrRotation = s_SamplesPerPixelRotation.str();
195 nbRotation = (int)atof(p_StrRotation.c_str());
199 int nbSlicePerEnergy = nbSlice / nbEnergyWindow;
200 if (nbSlicePerEnergy*nbEnergyWindow != nbSlice)
201 std::cerr << "Error: The number of slices is not correct !!" << std::endl;
203 for (unsigned int energy = 0; energy < nbEnergyWindow; ++energy) {
206 typename InputImageType::Pointer energyImage = InputImageType::New();
207 energyImage->SetRegions(input->GetLargestPossibleRegion());
208 typename InputImageType::IndexType startIndex;
211 startIndex[2] = 0;//energy*nbSlicePerEnergy;
212 typename InputImageType::SizeType regionSize;
213 regionSize[0] = input->GetLargestPossibleRegion().GetSize()[0];
214 regionSize[1] = input->GetLargestPossibleRegion().GetSize()[1];
215 regionSize[2] = nbSlicePerEnergy;
216 typename InputImageType::RegionType region;
217 region.SetSize(regionSize);
218 region.SetIndex(startIndex);
219 energyImage->SetRegions(region);
220 energyImage->Allocate();
222 //Create the iterator on the output
223 itk::ImageRegionIterator<InputImageType> imageOutputIterator(energyImage,energyImage->GetLargestPossibleRegion());
225 //Create the iterator on the region of the input
226 typename InputImageType::IndexType startIndexIterator;
227 startIndexIterator[0] = 0;
228 startIndexIterator[1] = 0;
229 startIndexIterator[2] = energy*nbSlicePerEnergy;
230 typename InputImageType::RegionType regionIterator;
231 regionIterator.SetSize(regionSize);
232 regionIterator.SetIndex(startIndexIterator);
233 itk::ImageRegionIterator<InputImageType> imageInputIterator(input,regionIterator);
235 //Copy the requested region
236 while(!imageInputIterator.IsAtEnd())
238 imageOutputIterator.Set(imageInputIterator.Get());
240 ++imageInputIterator;
241 ++imageOutputIterator;
244 // Output directory and filenames
245 typename WriterType::Pointer writer = WriterType::New();
246 writer->SetInput( energyImage );
247 writer->SetImageIO( gdcmIO );
248 typename ReaderType::FileNamesContainer fileNamesOutput;
249 std::string extension = NumberToString(energy) + ".dcm";
250 std::string filename_out = m_ArgsInfo.output_arg + extension;
251 fileNamesOutput.push_back(filename_out);
252 writer->SetFileNames( fileNamesOutput );
253 writer->SetMetaDataDictionaryArray(&outputArray);
257 if (m_ArgsInfo.verbose_flag)
258 std::cout << writer << std::endl;
260 } catch( itk::ExceptionObject & excp ) {
261 std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
262 std::cerr << excp << std::endl;
265 //Write sequence dicom tag with gdcm
267 reader.SetFileName( filename_out.c_str() );
269 gdcm::File &fileOutput = reader.GetFile();
270 gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
271 const unsigned int ptr_len = 42;
272 char *ptr = new char[ptr_len];
273 memset(ptr,0,ptr_len);
276 gdcm::SmartPointer<gdcm::SequenceOfItems> rotationSq = new gdcm::SequenceOfItems();
277 rotationSq->SetLengthToUndefined();
278 // Create a dataelement
279 gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
280 rotationDE.SetVR( gdcm::VR::US );
281 char essai = (char)nbRotation;
282 char *p_essai = &essai;
283 rotationDE.SetByteValue(p_essai, 1);
285 gdcm::Item rotationIt;
286 rotationIt.SetVLToUndefined();
287 gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
288 rotationDS.Insert(rotationDE);
289 rotationSq->AddItem(rotationIt);
290 // Insert sequence into data set
291 gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) );
292 rotationDEParent.SetVR(gdcm::VR::SQ);
293 rotationDEParent.SetValue(*rotationSq);
294 rotationDEParent.SetVLToUndefined();
296 gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) );
297 energyDEParent.SetVR(gdcm::VR::SQ);
299 gdcm::SmartPointer<gdcm::SequenceOfItems> energySq = new gdcm::SequenceOfItems();
300 energySq->SetLengthToUndefined();
302 //Read the caracteristics of this energy window : the name and the thresholds up and down
303 const gdcm::DataElement &seqWindow = dsInput.GetDataElement(gdcm::Tag(0x54, 0x12));
304 gdcm::SmartPointer<gdcm::SequenceOfItems> sqiWindow = seqWindow.GetValueAsSQ();
305 gdcm::Item &itemWindow = sqiWindow->GetItem(energy+1);
306 const gdcm::DataElement &deNameWindow = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x18));
307 std::string nameWindow( deNameWindow.GetByteValue()->GetPointer(), deNameWindow.GetByteValue()->GetLength() );
309 const gdcm::DataElement &deThresholds = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x13));
310 gdcm::SmartPointer<gdcm::SequenceOfItems> sqiThresholds = deThresholds.GetValueAsSQ();
311 for (unsigned int nbWindow=1; nbWindow <= sqiThresholds->GetNumberOfItems() ; ++nbWindow) {
312 gdcm::Item &itemThresholds = sqiThresholds->GetItem(nbWindow);
313 const gdcm::DataElement &deThresholdDown = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x14));
314 std::string p_StrThresholdDown( deThresholdDown.GetByteValue()->GetPointer(), deThresholdDown.GetByteValue()->GetLength() );
315 double thresholdDown = (double)atof(p_StrThresholdDown.c_str());
316 const gdcm::DataElement &deThresholdUp = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x15));
317 std::string p_StrThresholdUp( deThresholdUp.GetByteValue()->GetPointer(), deThresholdUp.GetByteValue()->GetLength() );
318 double thresholdUp = (double)atof(p_StrThresholdUp.c_str());
321 gdcm::SmartPointer<gdcm::SequenceOfItems> energyThresholdSq = new gdcm::SequenceOfItems();
322 energyThresholdSq->SetLengthToUndefined();
323 // Create a dataelement
324 gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) );
325 gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) );
326 energyThresholdDE.SetVR( gdcm::VR::DS );
327 energyUpholdDE.SetVR( gdcm::VR::DS );
328 energyThresholdDE.SetByteValue(NumberToString(thresholdDown).c_str(), (uint32_t)strlen(NumberToString(thresholdDown).c_str()));
329 energyUpholdDE.SetByteValue(NumberToString(thresholdUp).c_str(), (uint32_t)strlen(NumberToString(thresholdUp).c_str()));
331 gdcm::Item energyThresholdIt;
332 energyThresholdIt.SetVLToUndefined();
333 gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
334 energyThresholdDS.Insert(energyThresholdDE);
335 energyThresholdDS.Insert(energyUpholdDE);
336 energyThresholdSq->AddItem(energyThresholdIt);
337 // Insert sequence into data set
338 gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) );
339 energyThresholdDEParent.SetVR(gdcm::VR::SQ);
340 energyThresholdDEParent.SetValue(*energyThresholdSq);
341 energyThresholdDEParent.SetVLToUndefined();
342 // Create a dataelement
343 gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) );
344 energyNameDE.SetVR( gdcm::VR::SH );
345 energyNameDE.SetByteValue(nameWindow.c_str(), (uint32_t)strlen(nameWindow.c_str()));
348 energyIt.SetVLToUndefined();
349 gdcm::DataSet &energyDS = energyIt.GetNestedDataSet();
350 energyDS.Insert(energyNameDE);
351 energyDS.Insert(energyThresholdDEParent);
352 energySq->AddItem(energyIt);
354 // Insert sequence into data set
355 energyDEParent.SetValue(*energySq);
356 energyDEParent.SetVLToUndefined();
357 dsOutput.Insert(energyDEParent);
358 dsOutput.Insert(rotationDEParent);
360 w.SetFile( fileOutput );
361 w.SetFileName( filename_out.c_str() );
365 std::cout << "Use GDCM2" << std::endl;
368 //---------------------------------------------------------------------------
371 //---------------------------------------------------------------------------
372 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
374 typedef itk::MetaDataDictionary DictionaryType;
376 DictionaryType::ConstIterator itr = fromDict.Begin();
377 DictionaryType::ConstIterator end = fromDict.End();
378 typedef itk::MetaDataObject< std::string > MetaDataStringType;
382 itk::MetaDataObjectBase::Pointer entry = itr->second;
384 MetaDataStringType::Pointer entryvalue =
385 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
388 std::string tagkey = itr->first;
389 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
390 itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
395 //---------------------------------------------------------------------------
398 //---------------------------------------------------------------------------
399 template <typename T> std::string NumberToString ( T Number )
401 std::ostringstream ss;
405 //---------------------------------------------------------------------------
409 #endif //#define clitkPartitionEnergyWindowDicomGenericFilter_txx