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) {
  88     if (m_ArgsInfo.volume_flag)
  89       UpdateWithDim<2,2>(PixelType);
  90     else
  91       UpdateWithDim<2,1>(PixelType);
  92   } else if(Dimension==3) {
  93     if (m_ArgsInfo.volume_flag)
  94       UpdateWithDim<3,3>(PixelType);
  95     else
  96       UpdateWithDim<3,2>(PixelType);
  97   }
  98   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  99   else {
 100     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
 101     return;
 102   }
 103 }
 104
 105 //-------------------------------------------------------------------
 106 // Update with the number of dimensions
 107 //-------------------------------------------------------------------
 108 template<class args_info_type>
 109 template<unsigned int inputDimension, unsigned int outputDimension>
 110 void
 111 Image2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
 112 {
 113   if (m_Verbose) std::cout << "Image was detected to be "<<inputDimension<<"D and "<< PixelType<<"..."<<std::endl;
 114
 115   if(PixelType == "short") {
 116     if (m_Verbose) std::cout << "Launching filter in "<< inputDimension <<"D and signed short..." << std::endl;
 117     UpdateWithDimAndPixelType<inputDimension, outputDimension, signed short>();
 118   }
 119   else if(PixelType == "unsigned_short"){
 120     if (m_Verbose) std::cout  << "Launching filter in "<< inputDimension <<"D and unsigned_short..." << std::endl;
 121     UpdateWithDimAndPixelType<inputDimension, outputDimension, unsigned short>();
 122   }
 123
 124   else if (PixelType == "unsigned_char") {
 125     if (m_Verbose) std::cout  << "Launching filter in "<< inputDimension <<"D and unsigned_char..." << std::endl;
 126     UpdateWithDimAndPixelType<inputDimension, outputDimension, unsigned char>();
 127   }
 128
 129   //     else if (PixelType == "char"){
 130   //       if (m_Verbose) std::cout  << "Launching filter in "<< inputDimension <<"D and signed_char..." << std::endl;
 131   //       UpdateWithDimAndPixelType<inputDimension, outputDimension, signed char>();
 132   //     }
 133   else if (PixelType == "double") {
 134     if (m_Verbose) std::cout  << "Launching filter in "<< inputDimension <<"D and double..." << std::endl;
 135     UpdateWithDimAndPixelType<inputDimension, outputDimension, double>();
 136   }
 137   else {
 138     if (m_Verbose) std::cout  << "Launching filter in "<< inputDimension <<"D and float..." << std::endl;
 139     UpdateWithDimAndPixelType<inputDimension, outputDimension, float>();
 140   }
 141 }
 142
 143 //-------------------------------------------------------------------
 144 // Update with the number of dimensions and the pixeltype read from
 145 // the dicom files. The MHD files may be resampled to match the
 146 // dicom spacing (and number of slices). Rounding errors in resampling
 147 // are handled by removing files when generating the output dicom
 148 // series.
 149 //-------------------------------------------------------------------
 150 template<class args_info_type>
 151 template <unsigned int inputDimension, unsigned int outputDimension, class  PixelType>
 152 void
 153 Image2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 154 {
 155   // Validate input parameters
 156
 157   const unsigned int InputDimension = inputDimension;
 158   const unsigned int OutputDimension = outputDimension;
 159
 160   typedef itk::Image< PixelType, InputDimension > InputImageType;
 161   typedef itk::Image< PixelType, OutputDimension > OutputImageType;
 162   typedef itk::ImageSeriesReader< InputImageType > ReaderType;
 163   typedef itk::GDCMImageIO ImageIOType;
 164 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 165   typedef itk::ShiftScaleImageFilter< InputImageType, InputImageType > ShiftScaleType;
 166 #endif
 167   typedef itk::ImageSeriesWriter< InputImageType, OutputImageType > SeriesWriterType;
 168
 169 ////////////////////////////////////////////////
 170 // 1) Read the input series
 171
 172 // Read the input (MHD file)
 173   typedef typename InputImageType::RegionType RegionType;
 174   typedef typename RegionType::SizeType SizeType;
 175   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 176   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
 177   volumeReader->SetFileName(m_ArgsInfo.input_arg);
 178   volumeReader->Update();
 179   typename InputImageType::Pointer input = volumeReader->GetOutput();
 180   typedef itk::NumericSeriesFileNames OutputNamesGeneratorType;
 181
 182   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 183   gdcmIO->LoadPrivateTagsOn();
 184   typename ReaderType::FileNamesContainer filenames;
 185   filenames.push_back(m_ArgsInfo.inputDcm_arg);
 186   typename ReaderType::Pointer reader = ReaderType::New();
 187   reader->SetImageIO(gdcmIO);
 188   reader->SetFileNames(filenames);
 189   try {
 190     reader->Update();
 191   } catch (itk::ExceptionObject &excp) {
 192     std::cerr << "Exception thrown while reading the series" << std::endl;
 193     std::cerr << excp << std::endl;
 194     return;
 195   }
 196
 197   typename InputImageType::SpacingType outputSpacing;
 198   typename InputImageType::SizeType   outputSize;
 199   for (unsigned int i = 0; i < 3; i++) {
 200       outputSpacing[i] = input->GetSpacing()[i];
 201       outputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
 202   }
 203
 204 ////////////////////////////////////////////////
 205 // 2) Ensure to have value >= -1024
 206
 207   typedef itk::ThresholdImageFilter <InputImageType> ThresholdImageFilterType;
 208   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 209   thresholdFilter->SetInput(input);
 210   thresholdFilter->ThresholdBelow(-1024);
 211   thresholdFilter->SetOutsideValue(-1024);
 212   thresholdFilter->Update();
 213
 214   input=thresholdFilter->GetOutput();
 215
 216 ////////////////////////////////////////////////
 217 // 3) Create a MetaDataDictionary for each slice.
 218
 219   // Copy the dictionary from the first image and override slice
 220   // specific fields
 221   typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
 222   typename ReaderType::DictionaryArrayType outputArray;
 223
 224   // To keep the new series in the same study as the original we need
 225   // to keep the same study UID. But we need new series and frame of
 226   // reference UID's.
 227 #if ITK_VERSION_MAJOR >= 4
 228   gdcm::UIDGenerator fuid;
 229   std::string frameOfReferenceUID = fuid.Generate();
 230 #else
 231   std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 232 #endif
 233   std::string seriesUID;
 234   if (m_ArgsInfo.newUID_flag) {
 235 #if ITK_VERSION_MAJOR >= 4
 236     gdcm::UIDGenerator suid;
 237     seriesUID = suid.Generate();
 238 #else
 239     seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 240 #endif
 241   } else {
 242     itk::ExposeMetaData<std::string>(*inputDict, "0020|000e", seriesUID);
 243   }
 244   std::string studyUID;
 245   std::string sopClassUID;
 246   itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUID);
 247   itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUID);
 248   gdcmIO->KeepOriginalUIDOn();
 249   unsigned int numberOfKeysGiven=m_ArgsInfo.key_given;
 250   unsigned int fLimit = outputSize[2];
 251   if (m_ArgsInfo.volume_flag)
 252     fLimit = 1;
 253
 254   for (unsigned int f = 0; f < fLimit; f++)
 255   {
 256     // Create a new dictionary for this slice
 257     typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
 258
 259     typedef itk::MetaDataDictionary DictionaryType;
 260
 261     DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
 262     DictionaryType::ConstIterator endDic = (*inputDict).End();
 263     typedef itk::MetaDataObject< std::string > MetaDataStringType;
 264
 265     while( itrDic != endDic )
 266     {
 267       itk::MetaDataObjectBase::Pointer  entry = itrDic->second;
 268
 269       MetaDataStringType::Pointer entryvalue =
 270       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 271       if( entryvalue )
 272       {
 273         std::string tagkey   = itrDic->first;
 274         std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 275         itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
 276       }
 277       ++itrDic;
 278     }
 279
 280     // Set the UID's for the study, series, SOP  and frame of reference
 281     itk::EncapsulateMetaData<std::string>(*dict,"0020|000d", studyUID);
 282     itk::EncapsulateMetaData<std::string>(*dict,"0020|000e", seriesUID);
 283     itk::EncapsulateMetaData<std::string>(*dict,"0020|0052", frameOfReferenceUID);
 284
 285   #if ITK_VERSION_MAJOR >= 4
 286     gdcm::UIDGenerator sopuid;
 287     std::string sopInstanceUID = sopuid.Generate();
 288   #else
 289     std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 290   #endif
 291     itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
 292     itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
 293
 294     // Change fields that are slice specific
 295     std::ostringstream value;
 296     value.str("");
 297     value << f + 1;
 298     // Image Number
 299     itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
 300
 301     // Series Description - Append new description to current series
 302     // description
 303     std::string oldSeriesDesc;
 304     itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
 305
 306     value.str("");
 307     value << oldSeriesDesc
 308           << ": Resampled with pixel spacing "
 309           << outputSpacing[0] << ", "
 310           << outputSpacing[1] << ", "
 311           << outputSpacing[2];
 312     // This is an long string and there is a 64 character limit in the
 313     // standard
 314     unsigned lengthDesc = value.str().length();
 315
 316     std::string seriesDesc( value.str(), 0,
 317                             lengthDesc > 64 ? 64
 318                             : lengthDesc);
 319     itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
 320
 321     // Series Number
 322     value.str("");
 323     value << 1001;
 324     itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
 325
 326     // Derivation Description - How this image was derived
 327     value.str("");
 328
 329
 330
 331
 332
 333     value << ": " << ITK_SOURCE_VERSION;
 334
 335     lengthDesc = value.str().length();
 336     std::string derivationDesc( value.str(), 0,
 337                                 lengthDesc > 1024 ? 1024
 338                                 : lengthDesc);
 339     itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
 340
 341     // Image Position Patient: This is calculated by computing the
 342     // physical coordinate of the first pixel in each slice.
 343     typename InputImageType::PointType position;
 344     typename InputImageType::IndexType index;
 345     index[0] = 0;
 346     index[1] = 0;
 347     index[2] = f;
 348     input->TransformIndexToPhysicalPoint(index, position);
 349
 350     value.str("");
 351     value << position[0] << "\\" << position[1] << "\\" << position[2];
 352     itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());
 353     // Slice Location: For now, we store the z component of the Image
 354     // Position Patient.
 355     value.str("");
 356     value << position[2];
 357     itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
 358
 359     // Slice Thickness: For now, we store the z spacing
 360     value.str("");
 361     value << outputSpacing[2];
 362     itk::EncapsulateMetaData<std::string>(*dict,"0018|0050", value.str());
 363     // Spacing Between Slices
 364     itk::EncapsulateMetaData<std::string>(*dict,"0018|0088", value.str());
 365
 366     //Add ggo key
 367     for (unsigned int i = 0; i < numberOfKeysGiven; i++) {
 368       std::string entryId(m_ArgsInfo.key_arg[i]  );
 369       std::string value( m_ArgsInfo.tag_arg[i] );
 370
 371       itk::EncapsulateMetaData<std::string>(*dict, entryId, value);
 372     }
 373
 374     // Save the dictionary
 375     outputArray.push_back(dict);
 376   }
 377
 378 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 379 ////////////////////////////////////////////////
 380 // 4) Shift data to undo the effect of a rescale intercept by the
 381 //    DICOM reader
 382   std::string interceptTag("0028|1052");
 383   typedef itk::MetaDataObject<std::string> MetaDataStringType;
 384   itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
 385
 386   MetaDataStringType::ConstPointer interceptValue =
 387     dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
 388
 389   int interceptShift = 0;
 390   if(interceptValue ) {
 391     std::string tagValue = interceptValue->GetMetaDataObjectValue();
 392     interceptShift = -atoi (tagValue.c_str());
 393   }
 394
 395   ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
 396   shiftScale->SetInput(resampler->GetOutput());
 397   shiftScale->SetShift(interceptShift );
 398 #endif
 399
 400 ////////////////////////////////////////////////
 401 // 5) Write the new DICOM series
 402   // Generate the file names
 403   typename ReaderType::FileNamesContainer fileNamesOutput;
 404
 405   // Generate the file names
 406   OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
 407   std::string seriesFormat(m_ArgsInfo.outputDcm_arg);
 408   seriesFormat = seriesFormat + "/";
 409   if (m_ArgsInfo.nameDicom_given)
 410     seriesFormat = seriesFormat + m_ArgsInfo.nameDicom_arg + "_";
 411   seriesFormat = seriesFormat + "IM%d.dcm";
 412   outputNames->SetSeriesFormat(seriesFormat.c_str());
 413   outputNames->SetStartIndex(1);
 414   if (m_ArgsInfo.volume_flag) {
 415     outputNames->SetEndIndex(1);
 416   } else {
 417     outputNames->SetEndIndex(outputSize[2]);
 418   }
 419
 420   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 421 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 422   seriesWriter->SetInput(input);
 423 #else
 424   seriesWriter->SetInput(input);
 425 #endif
 426   seriesWriter->SetImageIO(gdcmIO);
 427   seriesWriter->SetFileNames(outputNames->GetFileNames());
 428   seriesWriter->SetMetaDataDictionaryArray(&outputArray);
 429   try {
 430     seriesWriter->Update();
 431   } catch(itk::ExceptionObject & excp) {
 432     std::cerr << "Exception thrown while writing the series " << std::endl;
 433     std::cerr << excp << std::endl;
 434     return;
 435   }
 436 /*
 437 ////////////////////////////////////////////////
 438 // 5) Read the new dicom data tag and copy it in the model data tag to have all dicom tags
 439   gdcm::Reader readerModel, readerOutput;
 440   readerModel.SetFileName(filenames[0].c_str());
 441   readerOutput.SetFileName(fileNamesOutput[0].c_str());
 442   readerModel.Read();
 443   readerOutput.Read();
 444   gdcm::File &fileModel = readerModel.GetFile();
 445   gdcm::File &fileOutput = readerOutput.GetFile();
 446   gdcm::DataSet &dsModel = fileModel.GetDataSet();
 447   gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
 448   const unsigned int ptr_len = 42;
 449   char *ptr = new char[ptr_len];
 450   memset(ptr,0,ptr_len);
 451
 452   const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
 453   dsModel.Replace(dataOutput);
 454   gdcm::Writer w;
 455   w.SetFile(fileModel);
 456   w.SetFileName(fileNamesOutput[0].c_str());
 457   w.Write();
 458   return;
 459 */
 460 }
 461 }//end clitk
 462
 463 #endif //#define clitkImage2DicomGenericFilter_txx