tools/clitkSpect2DicomGenericFilter.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 clitkSpect2DicomGenericFilter_txx
  19 #define clitkSpect2DicomGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkSpect2DicomGenericFilter.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 Spect2DicomGenericFilter<args_info_type>::Spect2DicomGenericFilter()
  68 {
  69   m_Verbose=false;
  70   m_InputFileName="";
  71 }
  72
  73
  74 //-----------------------------------------------------------
  75 // Update
  76 //-----------------------------------------------------------
  77 template<class args_info_type>
  78 void Spect2DicomGenericFilter<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 Spect2DicomGenericFilter<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 Spect2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 145 {
 146   // Validate input parameters
 147
 148   const unsigned int InputDimension = Dimension;
 149   const unsigned int OutputDimension = Dimension;
 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
 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);
 179   try {
 180     reader->Update();
 181   } catch (itk::ExceptionObject &excp) {
 182     std::cerr << "Exception thrown while reading the series" << std::endl;
 183     std::cerr << excp << std::endl;
 184     return;
 185   }
 186
 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];
 192   }
 193
 194 ////////////////////////////////////////////////
 195 // 2) Ensure to have value >= -1024
 196
 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();
 203
 204   input=thresholdFilter->GetOutput();
 205
 206 ////////////////////////////////////////////////
 207 // 3) Create a MetaDataDictionary for each slice.
 208
 209   // Copy the dictionary from the first image and override slice
 210   // specific fields
 211   typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
 212   typename ReaderType::DictionaryArrayType outputArray;
 213
 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
 216   // reference UID's.
 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();
 222 #else
 223   std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 224   std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 225 #endif
 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();
 231
 232   // Create a new dictionary for this slice
 233   typename ReaderType::DictionaryRawPointer dict = new typename ReaderType::DictionaryType;
 234
 235   typedef itk::MetaDataDictionary DictionaryType;
 236
 237   DictionaryType::ConstIterator itrDic = (*inputDict).Begin();
 238   DictionaryType::ConstIterator endDic = (*inputDict).End();
 239   typedef itk::MetaDataObject< std::string > MetaDataStringType;
 240
 241   while( itrDic != endDic )
 242   {
 243     itk::MetaDataObjectBase::Pointer  entry = itrDic->second;
 244
 245     MetaDataStringType::Pointer entryvalue =
 246     dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 247     if( entryvalue )
 248     {
 249       std::string tagkey   = itrDic->first;
 250       std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 251       itk::EncapsulateMetaData<std::string>(*dict, tagkey, tagvalue);
 252     }
 253     ++itrDic;
 254   }
 255
 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);
 260
 261 #if ITK_VERSION_MAJOR >= 4
 262   gdcm::UIDGenerator sopuid;
 263   std::string sopInstanceUID = sopuid.Generate();
 264 #else
 265   std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 266 #endif
 267   itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
 268   itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);
 269
 270   // Change fields that are slice specific
 271   std::ostringstream value;
 272   value.str("1");
 273
 274   // Image Number
 275   itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());
 276
 277   // Series Description - Append new description to current series
 278   // description
 279   std::string oldSeriesDesc;
 280   itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);
 281
 282   value.str("");
 283   value << oldSeriesDesc
 284         << ": Resampled with pixel spacing "
 285         << outputSpacing[0] << ", "
 286         << outputSpacing[1] << ", "
 287         << outputSpacing[2];
 288   // This is an long string and there is a 64 character limit in the
 289   // standard
 290   unsigned lengthDesc = value.str().length();
 291
 292   std::string seriesDesc( value.str(), 0,
 293                           lengthDesc > 64 ? 64
 294                           : lengthDesc);
 295   itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);
 296
 297   // Series Number
 298   value.str("");
 299   value << 1001;
 300   itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());
 301
 302   // Derivation Description - How this image was derived
 303   value.str("");
 304
 305   value << ": " << ITK_SOURCE_VERSION;
 306
 307   lengthDesc = value.str().length();
 308   std::string derivationDesc( value.str(), 0,
 309                               lengthDesc > 1024 ? 1024
 310                               : lengthDesc);
 311   itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
 312
 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;
 317   index[0] = 0;
 318   index[1] = 0;
 319   index[2] = 0;
 320   input->TransformIndexToPhysicalPoint(index, position);
 321
 322   value.str("");
 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
 326   // Position Patient.
 327   value.str("");
 328   value << position[2];
 329   itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
 330
 331   // Slice Thickness: For now, we store the z spacing
 332   value.str("");
 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);
 339
 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
 343 //    DICOM reader
 344   std::string interceptTag("0028|1052");
 345   typedef itk::MetaDataObject<std::string> MetaDataStringType;
 346   itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
 347
 348   MetaDataStringType::ConstPointer interceptValue =
 349     dynamic_cast<const MetaDataStringType *>(entry.GetPointer());
 350
 351   int interceptShift = 0;
 352   if(interceptValue ) {
 353     std::string tagValue = interceptValue->GetMetaDataObjectValue();
 354     interceptShift = -atoi (tagValue.c_str());
 355   }
 356
 357   ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
 358   shiftScale->SetInput(resampler->GetOutput());
 359   shiftScale->SetShift(interceptShift );
 360 #endif
 361
 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);
 367
 368   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 369 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR < 6 ) )
 370   seriesWriter->SetInput(input);
 371 #else
 372   seriesWriter->SetInput(input);
 373 #endif
 374   seriesWriter->SetImageIO(gdcmIO);
 375   seriesWriter->SetFileNames(fileNamesOutput);
 376   seriesWriter->SetMetaDataDictionaryArray(&outputArray);
 377   try {
 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;
 382     return;
 383   }
 384
 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());
 390   readerModel.Read();
 391   readerOutput.Read();
 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);
 399
 400   const gdcm::DataElement &dataOutput = dsOutput.GetDataElement(gdcm::Tag(0x7fe0, 0x10));
 401   dsModel.Replace(dataOutput);
 402   gdcm::Writer w;
 403   w.SetFile(fileModel);
 404   w.SetFileName(fileNamesOutput[0].c_str());
 405   w.Write();
 406   return;
 407 }
 408 }//end clitk
 409
 410 #endif //#define clitkSpect2DicomGenericFilter_txx