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 clitkImage2DicomGenericFilter_txx
19 #define clitkImage2DicomGenericFilter_txx
21 /* =================================================
22 * @file clitkImage2DicomGenericFilter.txx
28 ===================================================*/
30 #include "itkVersion.h"
32 #include "itkGDCMImageIO.h"
33 #include "itkGDCMSeriesFileNames.h"
34 #include "itkNumericSeriesFileNames.h"
35 #include "itkImageSeriesReader.h"
36 #include "itkImageSeriesWriter.h"
38 #include "itkThresholdImageFilter.h"
40 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
41 #include "itkShiftScaleImageFilter.h"
46 #if GDCM_MAJOR_VERSION >= 2
47 #include <gdcmUIDGenerator.h>
48 #include <gdcmImageHelper.h>
49 #include <gdcmAttribute.h>
50 #include <gdcmReader.h>
51 #include <gdcmWriter.h>
52 #include <gdcmDataElement.h>
63 //-----------------------------------------------------------
65 //-----------------------------------------------------------
66 template<class args_info_type>
67 Image2DicomGenericFilter<args_info_type>::Image2DicomGenericFilter()
74 //-----------------------------------------------------------
76 //-----------------------------------------------------------
77 template<class args_info_type>
78 void Image2DicomGenericFilter<args_info_type>::Update()
80 // Read the Dimension and PixelType
82 std::string PixelType;
83 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
87 if(Dimension==2) UpdateWithDim<2>(PixelType);
88 else if(Dimension==3) UpdateWithDim<3>(PixelType);
89 // else if (Dimension==4)UpdateWithDim<4>(PixelType);
91 std::cout<<"Error, Only for 2 or 3 Dimensions!!!"<<std::endl ;
96 //-------------------------------------------------------------------
97 // Update with the number of dimensions
98 //-------------------------------------------------------------------
99 template<class args_info_type>
100 template<unsigned int Dimension>
102 Image2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
104 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
106 if(PixelType == "short") {
107 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
108 UpdateWithDimAndPixelType<Dimension, signed short>();
110 else if(PixelType == "unsigned_short"){
111 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
112 UpdateWithDimAndPixelType<Dimension, unsigned short>();
115 else if (PixelType == "unsigned_char") {
116 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
117 UpdateWithDimAndPixelType<Dimension, unsigned char>();
120 // else if (PixelType == "char"){
121 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
122 // UpdateWithDimAndPixelType<Dimension, signed char>();
124 else if (PixelType == "double") {
125 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
126 UpdateWithDimAndPixelType<Dimension, double>();
129 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
130 UpdateWithDimAndPixelType<Dimension, float>();
134 //-------------------------------------------------------------------
135 // Update with the number of dimensions and the pixeltype read from
136 // the dicom files. The MHD files may be resampled to match the
137 // dicom spacing (and number of slices). Rounding errors in resampling
138 // are handled by removing files when generating the output dicom
140 //-------------------------------------------------------------------
141 template<class args_info_type>
142 template <unsigned int Dimension, class PixelType>
144 Image2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
146 // Validate input parameters
148 const unsigned int InputDimension = Dimension;
149 const unsigned int OutputDimension = Dimension-1;
151 typedef itk::Image< PixelType, InputDimension > InputImageType;
152 typedef itk::Image< PixelType, OutputDimension > OutputImageType;
153 typedef itk::ImageSeriesReader< InputImageType > ReaderType;
154 typedef itk::GDCMImageIO ImageIOType;
155 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
156 typedef itk::ShiftScaleImageFilter< InputImageType, InputImageType > ShiftScaleType;
158 typedef itk::ImageSeriesWriter< InputImageType, OutputImageType > SeriesWriterType;
160 ////////////////////////////////////////////////
161 // 1) Read the input series
163 // Read the input (MHD file)
164 typedef typename InputImageType::RegionType RegionType;
165 typedef typename RegionType::SizeType SizeType;
166 typedef itk::ImageFileReader<InputImageType> InputReaderType;
167 typename InputReaderType::Pointer volumeReader = InputReaderType::New();
168 volumeReader->SetFileName(m_ArgsInfo.input_arg);
169 volumeReader->Update();
170 typename InputImageType::Pointer input = volumeReader->GetOutput();
171 typedef itk::NumericSeriesFileNames OutputNamesGeneratorType;
173 ImageIOType::Pointer gdcmIO = ImageIOType::New();
174 gdcmIO->LoadPrivateTagsOn();
175 typename ReaderType::FileNamesContainer filenames;
176 filenames.push_back(m_ArgsInfo.inputDcm_arg);
177 typename ReaderType::Pointer reader = ReaderType::New();
178 reader->SetImageIO(gdcmIO);
179 reader->SetFileNames(filenames);
182 } catch (itk::ExceptionObject &excp) {
183 std::cerr << "Exception thrown while reading the series" << std::endl;
184 std::cerr << excp << std::endl;
188 typename InputImageType::SpacingType outputSpacing;
189 typename InputImageType::SizeType outputSize;
190 for (unsigned int i = 0; i < 3; i++) {
191 outputSpacing[i] = input->GetSpacing()[i];
192 outputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
195 ////////////////////////////////////////////////
196 // 2) Ensure to have value >= -1024
198 typedef itk::ThresholdImageFilter <InputImageType> ThresholdImageFilterType;
199 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
200 thresholdFilter->SetInput(input);
201 thresholdFilter->ThresholdBelow(-1024);
202 thresholdFilter->SetOutsideValue(-1024);
203 thresholdFilter->Update();
205 input=thresholdFilter->GetOutput();
207 ////////////////////////////////////////////////
208 // 3) Create a MetaDataDictionary for each slice.
210 // Copy the dictionary from the first image and override slice
212 typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
213 typename ReaderType::DictionaryArrayType outputArray;
215 // To keep the new series in the same study as the original we need
216 // to keep the same study UID. But we need new series and frame of
218 #if ITK_VERSION_MAJOR >= 4
219 gdcm::UIDGenerator suid;
220 std::string seriesUID = suid.Generate();
221 gdcm::UIDGenerator fuid;
222 std::string frameOfReferenceUID = fuid.Generate();
224 std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
225 std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
227 std::string studyUID;
228 std::string sopClassUID;
229 itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUID);
230 itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUID);
231 gdcmIO->KeepOriginalUIDOn();
232 for (unsigned int f = 0; f < outputSize[2]; f++)
234 // Create a new dictionary for this slice
235 typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
237 typedef itk::MetaDataDictionary DictionaryType;
239 DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
240 DictionaryType::ConstIterator endDic = (*inputDict).End();
241 typedef itk::MetaDataObject< std::string > MetaDataStringType;
243 while( itrDic != endDic )
245 itk::MetaDataObjectBase::Pointer entry = itrDic->second;
247 MetaDataStringType::Pointer entryvalue =
248 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
251 std::string tagkey = itrDic->first;
252 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
253 itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
258 // Set the UID's for the study, series, SOP and frame of reference
259 itk::EncapsulateMetaData<std::string>(*dict,"0020|000d", studyUID);
260 itk::EncapsulateMetaData<std::string>(*dict,"0020|000e", seriesUID);
261 itk::EncapsulateMetaData<std::string>(*dict,"0020|0052", frameOfReferenceUID);
263 #if ITK_VERSION_MAJOR >= 4
264 gdcm::UIDGenerator sopuid;
265 std::string sopInstanceUID = sopuid.Generate();
267 std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
269 itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
270 itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
272 // Change fields that are slice specific
273 std::ostringstream value;
277 itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
279 // Series Description - Append new description to current series
281 std::string oldSeriesDesc;
282 itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
285 value << oldSeriesDesc
286 << ": Resampled with pixel spacing "
287 << outputSpacing[0] << ", "
288 << outputSpacing[1] << ", "
290 // This is an long string and there is a 64 character limit in the
292 unsigned lengthDesc = value.str().length();
294 std::string seriesDesc( value.str(), 0,
297 itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
302 itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
304 // Derivation Description - How this image was derived
311 value << ": " << ITK_SOURCE_VERSION;
313 lengthDesc = value.str().length();
314 std::string derivationDesc( value.str(), 0,
315 lengthDesc > 1024 ? 1024
317 itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
319 // Image Position Patient: This is calculated by computing the
320 // physical coordinate of the first pixel in each slice.
321 typename InputImageType::PointType position;
322 typename InputImageType::IndexType index;
326 input->TransformIndexToPhysicalPoint(index, position);
329 value << position[0] << "\\" << position[1] << "\\" << position[2];
330 itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());
331 // Slice Location: For now, we store the z component of the Image
334 value << position[2];
335 itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
337 // Slice Thickness: For now, we store the z spacing
339 value << outputSpacing[2];
340 itk::EncapsulateMetaData<std::string>(*dict,"0018|0050", value.str());
341 // Spacing Between Slices
342 itk::EncapsulateMetaData<std::string>(*dict,"0018|0088", value.str());
343 // Save the dictionary
344 outputArray.push_back(dict);
347 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
348 ////////////////////////////////////////////////
349 // 4) Shift data to undo the effect of a rescale intercept by the
351 std::string interceptTag("0028|1052");
352 typedef itk::MetaDataObject<std::string> MetaDataStringType;
353 itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
355 MetaDataStringType::ConstPointer interceptValue =
356 dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
358 int interceptShift = 0;
359 if(interceptValue ) {
360 std::string tagValue = interceptValue->GetMetaDataObjectValue();
361 interceptShift = -atoi (tagValue.c_str());
364 ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
365 shiftScale->SetInput(resampler->GetOutput());
366 shiftScale->SetShift(interceptShift );
369 ////////////////////////////////////////////////
370 // 5) Write the new DICOM series
371 // Generate the file names
372 typename ReaderType::FileNamesContainer fileNamesOutput;
374 // Generate the file names
375 OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
376 std::string seriesFormat(m_ArgsInfo.outputDcm_arg);
377 seriesFormat = seriesFormat + "/" + "IM%d.dcm";
378 outputNames->SetSeriesFormat(seriesFormat.c_str());
379 outputNames->SetStartIndex(1);
380 outputNames->SetEndIndex(outputSize[2]);
382 typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
383 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
384 seriesWriter->SetInput(input);
386 seriesWriter->SetInput(input);
388 seriesWriter->SetImageIO(gdcmIO);
389 seriesWriter->SetFileNames(outputNames->GetFileNames());
390 seriesWriter->SetMetaDataDictionaryArray(&outputArray);
392 seriesWriter->Update();
393 } catch(itk::ExceptionObject & excp) {
394 std::cerr << "Exception thrown while writing the series " << std::endl;
395 std::cerr << excp << std::endl;
399 ////////////////////////////////////////////////
400 // 5) Read the new dicom data tag and copy it in the model data tag to have all dicom tags
401 gdcm::Reader readerModel, readerOutput;
402 readerModel.SetFileName(filenames[0].c_str());
403 readerOutput.SetFileName(fileNamesOutput[0].c_str());
406 gdcm::File &fileModel = readerModel.GetFile();
407 gdcm::File &fileOutput = readerOutput.GetFile();
408 gdcm::DataSet &dsModel = fileModel.GetDataSet();
409 gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
410 const unsigned int ptr_len = 42;
411 char *ptr = new char[ptr_len];
412 memset(ptr,0,ptr_len);
414 const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
415 dsModel.Replace(dataOutput);
417 w.SetFile(fileModel);
418 w.SetFileName(fileNamesOutput[0].c_str());
425 #endif //#define clitkImage2DicomGenericFilter_txx