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 clitkSpect2DicomGenericFilter_txx
19 #define clitkSpect2DicomGenericFilter_txx
21 /* =================================================
22 * @file clitkSpect2DicomGenericFilter.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 Spect2DicomGenericFilter<args_info_type>::Spect2DicomGenericFilter()
74 //-----------------------------------------------------------
76 //-----------------------------------------------------------
77 template<class args_info_type>
78 void Spect2DicomGenericFilter<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 Spect2DicomGenericFilter<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 Spect2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
146 // Validate input parameters
148 const unsigned int InputDimension = Dimension;
149 const unsigned int OutputDimension = Dimension;
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();
172 ImageIOType::Pointer gdcmIO = ImageIOType::New();
173 gdcmIO->LoadPrivateTagsOn();
174 typename ReaderType::FileNamesContainer filenames;
175 filenames.push_back(m_ArgsInfo.inputDcm_arg);
176 typename ReaderType::Pointer reader = ReaderType::New();
177 reader->SetImageIO(gdcmIO);
178 reader->SetFileNames(filenames);
181 } catch (itk::ExceptionObject &excp) {
182 std::cerr << "Exception thrown while reading the series" << std::endl;
183 std::cerr << excp << std::endl;
187 typename InputImageType::SpacingType outputSpacing;
188 typename InputImageType::SizeType outputSize;
189 for (unsigned int i = 0; i < 3; i++) {
190 outputSpacing[i] = input->GetSpacing()[i];
191 outputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
194 ////////////////////////////////////////////////
195 // 2) Ensure to have value >= -1024
197 typedef itk::ThresholdImageFilter <InputImageType> ThresholdImageFilterType;
198 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
199 thresholdFilter->SetInput(input);
200 thresholdFilter->ThresholdBelow(-1024);
201 thresholdFilter->SetOutsideValue(-1024);
202 thresholdFilter->Update();
204 input=thresholdFilter->GetOutput();
206 ////////////////////////////////////////////////
207 // 3) Create a MetaDataDictionary for each slice.
209 // Copy the dictionary from the first image and override slice
211 typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
212 typename ReaderType::DictionaryArrayType outputArray;
214 // To keep the new series in the same study as the original we need
215 // to keep the same study UID. But we need new series and frame of
217 #if ITK_VERSION_MAJOR >= 4
218 gdcm::UIDGenerator suid;
219 std::string seriesUID = suid.Generate();
220 gdcm::UIDGenerator fuid;
221 std::string frameOfReferenceUID = fuid.Generate();
223 std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
224 std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
226 std::string studyUID;
227 std::string sopClassUID;
228 itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUID);
229 itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUID);
230 gdcmIO->KeepOriginalUIDOn();
232 // Create a new dictionary for this slice
233 typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
235 typedef itk::MetaDataDictionary DictionaryType;
237 DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
238 DictionaryType::ConstIterator endDic = (*inputDict).End();
239 typedef itk::MetaDataObject< std::string > MetaDataStringType;
241 while( itrDic != endDic )
243 itk::MetaDataObjectBase::Pointer entry = itrDic->second;
245 MetaDataStringType::Pointer entryvalue =
246 dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
249 std::string tagkey = itrDic->first;
250 std::string tagvalue = entryvalue->GetMetaDataObjectValue();
251 itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
256 // Set the UID's for the study, series, SOP and frame of reference
257 itk::EncapsulateMetaData<std::string>(*dict,"0020|000d", studyUID);
258 itk::EncapsulateMetaData<std::string>(*dict,"0020|000e", seriesUID);
259 itk::EncapsulateMetaData<std::string>(*dict,"0020|0052", frameOfReferenceUID);
261 #if ITK_VERSION_MAJOR >= 4
262 gdcm::UIDGenerator sopuid;
263 std::string sopInstanceUID = sopuid.Generate();
265 std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
267 itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
268 itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
270 // Change fields that are slice specific
271 std::ostringstream value;
275 itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
277 // Series Description - Append new description to current series
279 std::string oldSeriesDesc;
280 itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
283 value << oldSeriesDesc
284 << ": Resampled with pixel spacing "
285 << outputSpacing[0] << ", "
286 << outputSpacing[1] << ", "
288 // This is an long string and there is a 64 character limit in the
290 unsigned lengthDesc = value.str().length();
292 std::string seriesDesc( value.str(), 0,
295 itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
300 itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
302 // Derivation Description - How this image was derived
305 value << ": " << ITK_SOURCE_VERSION;
307 lengthDesc = value.str().length();
308 std::string derivationDesc( value.str(), 0,
309 lengthDesc > 1024 ? 1024
311 itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
313 // Image Position Patient: This is calculated by computing the
314 // physical coordinate of the first pixel in each slice.
315 typename InputImageType::PointType position;
316 typename InputImageType::IndexType index;
320 input->TransformIndexToPhysicalPoint(index, position);
323 value << position[0] << "\\" << position[1] << "\\" << position[2];
324 itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());
325 // Slice Location: For now, we store the z component of the Image
328 value << position[2];
329 itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
331 // Slice Thickness: For now, we store the z spacing
333 value << outputSpacing[2];
334 itk::EncapsulateMetaData<std::string>(*dict,"0018|0050", value.str());
335 // Spacing Between Slices
336 itk::EncapsulateMetaData<std::string>(*dict,"0018|0088", value.str());
337 // Save the dictionary
338 outputArray.push_back(inputDict);
340 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
341 ////////////////////////////////////////////////
342 // 4) Shift data to undo the effect of a rescale intercept by the
344 std::string interceptTag("0028|1052");
345 typedef itk::MetaDataObject<std::string> MetaDataStringType;
346 itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
348 MetaDataStringType::ConstPointer interceptValue =
349 dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
351 int interceptShift = 0;
352 if(interceptValue ) {
353 std::string tagValue = interceptValue->GetMetaDataObjectValue();
354 interceptShift = -atoi (tagValue.c_str());
357 ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
358 shiftScale->SetInput(resampler->GetOutput());
359 shiftScale->SetShift(interceptShift );
362 ////////////////////////////////////////////////
363 // 5) Write the new DICOM series
364 // Generate the file names
365 typename ReaderType::FileNamesContainer fileNamesOutput;
366 fileNamesOutput.push_back(m_ArgsInfo.outputDcm_arg);
368 typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
369 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
370 seriesWriter->SetInput(input);
372 seriesWriter->SetInput(input);
374 seriesWriter->SetImageIO(gdcmIO);
375 seriesWriter->SetFileNames(fileNamesOutput);
376 seriesWriter->SetMetaDataDictionaryArray(&outputArray);
378 seriesWriter->Update();
379 } catch(itk::ExceptionObject & excp) {
380 std::cerr << "Exception thrown while writing the series " << std::endl;
381 std::cerr << excp << std::endl;
385 ////////////////////////////////////////////////
386 // 5) Read the new dicom data tag and copy it in the model data tag to have all dicom tags
387 gdcm::Reader readerModel, readerOutput;
388 readerModel.SetFileName(filenames[0].c_str());
389 readerOutput.SetFileName(fileNamesOutput[0].c_str());
392 gdcm::File &fileModel = readerModel.GetFile();
393 gdcm::File &fileOutput = readerOutput.GetFile();
394 gdcm::DataSet &dsModel = fileModel.GetDataSet();
395 gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
396 const unsigned int ptr_len = 42;
397 char *ptr = new char[ptr_len];
398 memset(ptr,0,ptr_len);
400 const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
401 dsModel.Replace(dataOutput);
403 w.SetFile(fileModel);
404 w.SetFileName(fileNamesOutput[0].c_str());
410 #endif //#define clitkSpect2DicomGenericFilter_txx