tools/clitkImage2DicomGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   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
   8
   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.
  12
  13   It is distributed under dual licence
  14
  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
  20
  21 /* =================================================
  22  * @file   clitkImage2DicomGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "itkVersion.h"
  31 #include "itkImage.h"
  32 #include "itkGDCMImageIO.h"
  33 #include "itkGDCMSeriesFileNames.h"
  34 #include "itkNumericSeriesFileNames.h"
  35 #include "itkImageSeriesReader.h"
  36 #include "itkImageSeriesWriter.h"
  37
  38 #include "itkThresholdImageFilter.h"
  39
  40 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
  41 #include "itkShiftScaleImageFilter.h"
  42 #endif
  43
  44 #include <string>
  45 #include <sstream>
  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>
  53 #include <gdcmTag.h>
  54 #else
  55 #include "gdcmFile.h"
  56 #include "gdcmUtil.h"
  57 #endif
  58
  59 namespace clitk
  60 {
  61
  62
  63 //-----------------------------------------------------------
  64 // Constructor
  65 //-----------------------------------------------------------
  66 template<class args_info_type>
  67 Image2DicomGenericFilter<args_info_type>::Image2DicomGenericFilter()
  68 {
  69   m_Verbose=false;
  70   m_InputFileName="";
  71 }
  72
  73
  74 //-----------------------------------------------------------
  75 // Update
  76 //-----------------------------------------------------------
  77 template<class args_info_type>
  78 void Image2DicomGenericFilter<args_info_type>::Update()
  79 {
  80   // Read the Dimension and PixelType
  81   int Dimension;
  82   std::string PixelType;
  83   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  84
  85
  86   // Call UpdateWithDim
  87   if(Dimension==2) UpdateWithDim<2>(PixelType);
  88   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  89   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  90   else {
  91     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  92     return;
  93   }
  94 }
  95
  96 //-------------------------------------------------------------------
  97 // Update with the number of dimensions
  98 //-------------------------------------------------------------------
  99 template<class args_info_type>
 100 template<unsigned int Dimension>
 101 void
 102 Image2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
 103 {
 104   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
 105
 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>();
 109   }
 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>();
 113   }
 114
 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>();
 118   }
 119
 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>();
 123   //     }
 124   else if (PixelType == "double") {
 125     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
 126     UpdateWithDimAndPixelType<Dimension, double>();
 127   }
 128   else {
 129     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 130     UpdateWithDimAndPixelType<Dimension, float>();
 131   }
 132 }
 133
 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
 139 // series.
 140 //-------------------------------------------------------------------
 141 template<class args_info_type>
 142 template <unsigned int Dimension, class  PixelType>
 143 void
 144 Image2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 145 {
 146   // Validate input parameters
 147
 148   const unsigned int InputDimension = Dimension;
 149   const unsigned int OutputDimension = Dimension-1;
 150
 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;
 157 #endif
 158   typedef itk::ImageSeriesWriter< InputImageType, OutputImageType > SeriesWriterType;
 159
 160 ////////////////////////////////////////////////
 161 // 1) Read the input series
 162
 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;
 172
 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);
 180   try {
 181     reader->Update();
 182   } catch (itk::ExceptionObject &excp) {
 183     std::cerr << "Exception thrown while reading the series" << std::endl;
 184     std::cerr << excp << std::endl;
 185     return;
 186   }
 187
 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];
 193   }
 194
 195 ////////////////////////////////////////////////
 196 // 2) Ensure to have value >= -1024
 197
 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();
 204
 205   input=thresholdFilter->GetOutput();
 206
 207 ////////////////////////////////////////////////
 208 // 3) Create a MetaDataDictionary for each slice.
 209
 210   // Copy the dictionary from the first image and override slice
 211   // specific fields
 212   typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
 213   typename ReaderType::DictionaryArrayType outputArray;
 214
 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
 217   // reference UID's.
 218 #if ITK_VERSION_MAJOR >= 4
 219   gdcm::UIDGenerator fuid;
 220   std::string frameOfReferenceUID = fuid.Generate();
 221 #else
 222   std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 223 #endif
 224   std::string seriesUID;
 225   if (m_ArgsInfo.newUID_flag) {
 226 #if ITK_VERSION_MAJOR >= 4
 227     gdcm::UIDGenerator suid;
 228     seriesUID = suid.Generate();
 229 #else
 230     seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 231 #endif
 232   } else {
 233     itk::ExposeMetaData<std::string>(*inputDict, "0020|000e", seriesUID);
 234   }
 235   std::string studyUID;
 236   std::string sopClassUID;
 237   itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUID);
 238   itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUID);
 239   gdcmIO->KeepOriginalUIDOn();
 240   for (unsigned int f = 0; f < outputSize[2]; f++)
 241   {
 242     // Create a new dictionary for this slice
 243     typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
 244
 245     typedef itk::MetaDataDictionary DictionaryType;
 246
 247     DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
 248     DictionaryType::ConstIterator endDic = (*inputDict).End();
 249     typedef itk::MetaDataObject< std::string > MetaDataStringType;
 250
 251     while( itrDic != endDic )
 252     {
 253       itk::MetaDataObjectBase::Pointer  entry = itrDic->second;
 254
 255       MetaDataStringType::Pointer entryvalue =
 256       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 257       if( entryvalue )
 258       {
 259         std::string tagkey   = itrDic->first;
 260         std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 261         itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
 262       }
 263       ++itrDic;
 264     }
 265
 266     // Set the UID's for the study, series, SOP  and frame of reference
 267     itk::EncapsulateMetaData<std::string>(*dict,"0020|000d", studyUID);
 268     itk::EncapsulateMetaData<std::string>(*dict,"0020|000e", seriesUID);
 269     itk::EncapsulateMetaData<std::string>(*dict,"0020|0052", frameOfReferenceUID);
 270
 271   #if ITK_VERSION_MAJOR >= 4
 272     gdcm::UIDGenerator sopuid;
 273     std::string sopInstanceUID = sopuid.Generate();
 274   #else
 275     std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 276   #endif
 277     itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
 278     itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
 279
 280     // Change fields that are slice specific
 281     std::ostringstream value;
 282     value.str("");
 283     value << f + 1;
 284     // Image Number
 285     itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
 286
 287     // Series Description - Append new description to current series
 288     // description
 289     std::string oldSeriesDesc;
 290     itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
 291
 292     value.str("");
 293     value << oldSeriesDesc
 294           << ": Resampled with pixel spacing "
 295           << outputSpacing[0] << ", "
 296           << outputSpacing[1] << ", "
 297           << outputSpacing[2];
 298     // This is an long string and there is a 64 character limit in the
 299     // standard
 300     unsigned lengthDesc = value.str().length();
 301
 302     std::string seriesDesc( value.str(), 0,
 303                             lengthDesc > 64 ? 64
 304                             : lengthDesc);
 305     itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
 306
 307     // Series Number
 308     value.str("");
 309     value << 1001;
 310     itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
 311
 312     // Derivation Description - How this image was derived
 313     value.str("");
 314
 315
 316
 317
 318
 319     value << ": " << ITK_SOURCE_VERSION;
 320
 321     lengthDesc = value.str().length();
 322     std::string derivationDesc( value.str(), 0,
 323                                 lengthDesc > 1024 ? 1024
 324                                 : lengthDesc);
 325     itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
 326
 327     // Image Position Patient: This is calculated by computing the
 328     // physical coordinate of the first pixel in each slice.
 329     typename InputImageType::PointType position;
 330     typename InputImageType::IndexType index;
 331     index[0] = 0;
 332     index[1] = 0;
 333     index[2] = f;
 334     input->TransformIndexToPhysicalPoint(index, position);
 335
 336     value.str("");
 337     value << position[0] << "\\" << position[1] << "\\" << position[2];
 338     itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());
 339     // Slice Location: For now, we store the z component of the Image
 340     // Position Patient.
 341     value.str("");
 342     value << position[2];
 343     itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
 344
 345     // Slice Thickness: For now, we store the z spacing
 346     value.str("");
 347     value << outputSpacing[2];
 348     itk::EncapsulateMetaData<std::string>(*dict,"0018|0050", value.str());
 349     // Spacing Between Slices
 350     itk::EncapsulateMetaData<std::string>(*dict,"0018|0088", value.str());
 351     // Save the dictionary
 352     outputArray.push_back(dict);
 353   }
 354
 355 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 356 ////////////////////////////////////////////////
 357 // 4) Shift data to undo the effect of a rescale intercept by the
 358 //    DICOM reader
 359   std::string interceptTag("0028|1052");
 360   typedef itk::MetaDataObject<std::string> MetaDataStringType;
 361   itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
 362
 363   MetaDataStringType::ConstPointer interceptValue =
 364     dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
 365
 366   int interceptShift = 0;
 367   if(interceptValue ) {
 368     std::string tagValue = interceptValue->GetMetaDataObjectValue();
 369     interceptShift = -atoi (tagValue.c_str());
 370   }
 371
 372   ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
 373   shiftScale->SetInput(resampler->GetOutput());
 374   shiftScale->SetShift(interceptShift );
 375 #endif
 376
 377 ////////////////////////////////////////////////
 378 // 5) Write the new DICOM series
 379   // Generate the file names
 380   typename ReaderType::FileNamesContainer fileNamesOutput;
 381
 382   // Generate the file names
 383   OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
 384   std::string seriesFormat(m_ArgsInfo.outputDcm_arg);
 385   seriesFormat = seriesFormat + "/";
 386   if (m_ArgsInfo.nameDicom_given)
 387     seriesFormat = seriesFormat + m_ArgsInfo.nameDicom_arg + "_";
 388   seriesFormat = seriesFormat + "IM%d.dcm";
 389   outputNames->SetSeriesFormat(seriesFormat.c_str());
 390   outputNames->SetStartIndex(1);
 391   outputNames->SetEndIndex(outputSize[2]);
 392
 393   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 394 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 395   seriesWriter->SetInput(input);
 396 #else
 397   seriesWriter->SetInput(input);
 398 #endif
 399   seriesWriter->SetImageIO(gdcmIO);
 400   seriesWriter->SetFileNames(outputNames->GetFileNames());
 401   seriesWriter->SetMetaDataDictionaryArray(&outputArray);
 402   try {
 403     seriesWriter->Update();
 404   } catch(itk::ExceptionObject & excp) {
 405     std::cerr << "Exception thrown while writing the series " << std::endl;
 406     std::cerr << excp << std::endl;
 407     return;
 408   }
 409 /*
 410 ////////////////////////////////////////////////
 411 // 5) Read the new dicom data tag and copy it in the model data tag to have all dicom tags
 412   gdcm::Reader readerModel, readerOutput;
 413   readerModel.SetFileName(filenames[0].c_str());
 414   readerOutput.SetFileName(fileNamesOutput[0].c_str());
 415   readerModel.Read();
 416   readerOutput.Read();
 417   gdcm::File &fileModel = readerModel.GetFile();
 418   gdcm::File &fileOutput = readerOutput.GetFile();
 419   gdcm::DataSet &dsModel = fileModel.GetDataSet();
 420   gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
 421   const unsigned int ptr_len = 42;
 422   char *ptr = new char[ptr_len];
 423   memset(ptr,0,ptr_len);
 424
 425   const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
 426   dsModel.Replace(dataOutput);
 427   gdcm::Writer w;
 428   w.SetFile(fileModel);
 429   w.SetFileName(fileNamesOutput[0].c_str());
 430   w.Write();
 431   return;
 432 */
 433 }
 434 }//end clitk
 435
 436 #endif //#define clitkImage2DicomGenericFilter_txx