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 suid;
 220   std::string seriesUID = suid.Generate();
 221   gdcm::UIDGenerator fuid;
 222   std::string frameOfReferenceUID = fuid.Generate();
 223 #else
 224   std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 225   std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 226 #endif
 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++)
 233   {
 234     // Create a new dictionary for this slice
 235     typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
 236
 237     typedef itk::MetaDataDictionary DictionaryType;
 238
 239     DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
 240     DictionaryType::ConstIterator endDic = (*inputDict).End();
 241     typedef itk::MetaDataObject< std::string > MetaDataStringType;
 242
 243     while( itrDic != endDic )
 244     {
 245       itk::MetaDataObjectBase::Pointer  entry = itrDic->second;
 246
 247       MetaDataStringType::Pointer entryvalue =
 248       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 249       if( entryvalue )
 250       {
 251         std::string tagkey   = itrDic->first;
 252         std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 253         itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
 254       }
 255       ++itrDic;
 256     }
 257
 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);
 262
 263   #if ITK_VERSION_MAJOR >= 4
 264     gdcm::UIDGenerator sopuid;
 265     std::string sopInstanceUID = sopuid.Generate();
 266   #else
 267     std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 268   #endif
 269     itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
 270     itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
 271
 272     // Change fields that are slice specific
 273     std::ostringstream value;
 274     value.str("");
 275     value << f + 1;
 276     // Image Number
 277     itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
 278
 279     // Series Description - Append new description to current series
 280     // description
 281     std::string oldSeriesDesc;
 282     itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
 283
 284     value.str("");
 285     value << oldSeriesDesc
 286           << ": Resampled with pixel spacing "
 287           << outputSpacing[0] << ", "
 288           << outputSpacing[1] << ", "
 289           << outputSpacing[2];
 290     // This is an long string and there is a 64 character limit in the
 291     // standard
 292     unsigned lengthDesc = value.str().length();
 293
 294     std::string seriesDesc( value.str(), 0,
 295                             lengthDesc > 64 ? 64
 296                             : lengthDesc);
 297     itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
 298
 299     // Series Number
 300     value.str("");
 301     value << 1001;
 302     itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
 303
 304     // Derivation Description - How this image was derived
 305     value.str("");
 306
 307
 308
 309
 310
 311     value << ": " << ITK_SOURCE_VERSION;
 312
 313     lengthDesc = value.str().length();
 314     std::string derivationDesc( value.str(), 0,
 315                                 lengthDesc > 1024 ? 1024
 316                                 : lengthDesc);
 317     itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
 318
 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;
 323     index[0] = 0;
 324     index[1] = 0;
 325     index[2] = f;
 326     input->TransformIndexToPhysicalPoint(index, position);
 327
 328     value.str("");
 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
 332     // Position Patient.
 333     value.str("");
 334     value << position[2];
 335     itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
 336
 337     // Slice Thickness: For now, we store the z spacing
 338     value.str("");
 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);
 345   }
 346
 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
 350 //    DICOM reader
 351   std::string interceptTag("0028|1052");
 352   typedef itk::MetaDataObject<std::string> MetaDataStringType;
 353   itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
 354
 355   MetaDataStringType::ConstPointer interceptValue =
 356     dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
 357
 358   int interceptShift = 0;
 359   if(interceptValue ) {
 360     std::string tagValue = interceptValue->GetMetaDataObjectValue();
 361     interceptShift = -atoi (tagValue.c_str());
 362   }
 363
 364   ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
 365   shiftScale->SetInput(resampler->GetOutput());
 366   shiftScale->SetShift(interceptShift );
 367 #endif
 368
 369 ////////////////////////////////////////////////
 370 // 5) Write the new DICOM series
 371   // Generate the file names
 372   typename ReaderType::FileNamesContainer fileNamesOutput;
 373
 374   // Generate the file names
 375   OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
 376   std::string seriesFormat(m_ArgsInfo.outputDcm_arg);
 377   seriesFormat = seriesFormat + "/";
 378   if (m_ArgsInfo.nameDicom_given)
 379     seriesFormat = seriesFormat + m_ArgsInfo.nameDicom_arg + "_";
 380   seriesFormat = seriesFormat + "IM%d.dcm";
 381   outputNames->SetSeriesFormat(seriesFormat.c_str());
 382   outputNames->SetStartIndex(1);
 383   outputNames->SetEndIndex(outputSize[2]);
 384
 385   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 386 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 387   seriesWriter->SetInput(input);
 388 #else
 389   seriesWriter->SetInput(input);
 390 #endif
 391   seriesWriter->SetImageIO(gdcmIO);
 392   seriesWriter->SetFileNames(outputNames->GetFileNames());
 393   seriesWriter->SetMetaDataDictionaryArray(&outputArray);
 394   try {
 395     seriesWriter->Update();
 396   } catch(itk::ExceptionObject & excp) {
 397     std::cerr << "Exception thrown while writing the series " << std::endl;
 398     std::cerr << excp << std::endl;
 399     return;
 400   }
 401 /*
 402 ////////////////////////////////////////////////
 403 // 5) Read the new dicom data tag and copy it in the model data tag to have all dicom tags
 404   gdcm::Reader readerModel, readerOutput;
 405   readerModel.SetFileName(filenames[0].c_str());
 406   readerOutput.SetFileName(fileNamesOutput[0].c_str());
 407   readerModel.Read();
 408   readerOutput.Read();
 409   gdcm::File &fileModel = readerModel.GetFile();
 410   gdcm::File &fileOutput = readerOutput.GetFile();
 411   gdcm::DataSet &dsModel = fileModel.GetDataSet();
 412   gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
 413   const unsigned int ptr_len = 42;
 414   char *ptr = new char[ptr_len];
 415   memset(ptr,0,ptr_len);
 416
 417   const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
 418   dsModel.Replace(dataOutput);
 419   gdcm::Writer w;
 420   w.SetFile(fileModel);
 421   w.SetFileName(fileNamesOutput[0].c_str());
 422   w.Write();
 423   return;
 424 */
 425 }
 426 }//end clitk
 427
 428 #endif //#define clitkImage2DicomGenericFilter_txx