tools/clitkGateSimulation2DicomGenericFilter.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 clitkGateSimulation2DicomGenericFilter_txx
  19 #define clitkGateSimulation2DicomGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkGateSimulation2DicomGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include <sstream>
  31 // clitk
  32 #include "clitkResampleImageWithOptionsFilter.h"
  33 #if GDCM_MAJOR_VERSION >= 2
  34 #include "gdcmUIDGenerator.h"
  35 #include <gdcmImageHelper.h>
  36 #include <gdcmAttribute.h>
  37 #include <gdcmReader.h>
  38 #else
  39 #include "gdcmFile.h"
  40 #include "gdcmUtil.h"
  41 #endif
  42
  43 #include "itkImageRegionIterator.h"
  44
  45
  46 namespace clitk
  47 {
  48
  49
  50 //-----------------------------------------------------------
  51 // Constructor
  52 //-----------------------------------------------------------
  53 template<class args_info_type>
  54 GateSimulation2DicomGenericFilter<args_info_type>::GateSimulation2DicomGenericFilter()
  55 {
  56   m_Verbose=false;
  57   m_InputFileName="";
  58 }
  59
  60
  61 //-----------------------------------------------------------
  62 // Update
  63 //-----------------------------------------------------------
  64 template<class args_info_type>
  65 void GateSimulation2DicomGenericFilter<args_info_type>::Update()
  66 {
  67   // Read the Dimension and PixelType
  68   int Dimension;
  69   std::string PixelType;
  70   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  71
  72
  73   // Call UpdateWithDim
  74   if(Dimension==2) UpdateWithDim<2>(PixelType);
  75   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  76   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  77   else {
  78     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  79     return;
  80   }
  81 }
  82
  83 //-------------------------------------------------------------------
  84 // Update with the number of dimensions
  85 //-------------------------------------------------------------------
  86 template<class args_info_type>
  87 template<unsigned int Dimension>
  88 void
  89 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
  90 {
  91   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  92
  93   if(PixelType == "short") {
  94     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  95     UpdateWithDimAndPixelType<Dimension, signed short>();
  96   }
  97   else if(PixelType == "unsigned_short"){
  98     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  99     UpdateWithDimAndPixelType<Dimension, unsigned short>();
 100   }
 101
 102   else if (PixelType == "unsigned_char") {
 103     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
 104     UpdateWithDimAndPixelType<Dimension, unsigned char>();
 105   }
 106
 107   //     else if (PixelType == "char"){
 108   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
 109   //       UpdateWithDimAndPixelType<Dimension, signed char>();
 110   //     }
 111   else if (PixelType == "double") {
 112     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
 113     UpdateWithDimAndPixelType<Dimension, double>();
 114   }
 115   else {
 116     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 117     UpdateWithDimAndPixelType<Dimension, float>();
 118   }
 119 }
 120
 121 //-------------------------------------------------------------------
 122 // Update with the number of dimensions and the pixeltype read from
 123 // the dicom files. The MHD files may be resampled to match the
 124 // dicom spacing (and number of slices). Rounding errors in resampling
 125 // are handled by removing files when generating the output dicom
 126 // series.
 127 //-------------------------------------------------------------------
 128 template<class args_info_type>
 129 template <unsigned int Dimension, class  PixelType>
 130 void
 131 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 132 {
 133
 134 #if GDCM_MAJOR_VERSION == 2
 135   // ImageTypes
 136   typedef itk::Image<PixelType, Dimension> InputImageType;
 137   typedef itk::Image<PixelType, Dimension> OutputImageType;
 138   typedef itk::ImageFileReader< InputImageType >     ReaderType;
 139   typedef itk::ImageSeriesReader< InputImageType >     ReaderSeriesType;
 140   typedef itk::GDCMImageIO ImageIOType;
 141
 142
 143   // Read Dicom model file
 144   typename ReaderType::Pointer reader = ReaderType::New();
 145   typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
 146   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 147   std::string filename_out = m_ArgsInfo.outputFilename_arg;
 148   gdcmIO->LoadPrivateTagsOn();
 149   gdcmIO->KeepOriginalUIDOn();
 150   reader->SetImageIO( gdcmIO );
 151   reader->SetFileName( m_ArgsInfo.inputModelFilename_arg );
 152   typename ReaderSeriesType::FileNamesContainer fileNames;
 153   fileNames.push_back(m_ArgsInfo.inputModelFilename_arg);
 154   readerSeries->SetImageIO( gdcmIO );
 155   readerSeries->SetFileNames( fileNames );
 156   try {
 157     reader->Update();
 158     readerSeries->Update();
 159   } catch (itk::ExceptionObject &excp) {
 160     std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
 161     std::cerr << excp << std::endl;
 162   }
 163
 164
 165   // Read the input (MHD file)
 166   typedef typename InputImageType::RegionType RegionType;
 167   typedef typename RegionType::SizeType SizeType;
 168   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 169   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
 170   volumeReader->SetFileName( m_InputFileName);
 171   volumeReader->Update();
 172   typename InputImageType::Pointer input = volumeReader->GetOutput();
 173
 174
 175 //TODO if the input size/spacing and dicom model ones are different
 176 /*  if (input->GetLargestPossibleRegion().GetSize() != reader->GetOutput()->GetLargestPossibleRegion().GetSize()) {
 177
 178     // resampling is carried out on the fly if resolution or size between
 179     // the input mhd and input dicom series is different
 180
 181     // Filter
 182     typedef clitk::ResampleImageWithOptionsFilter<InputImageType, InputImageType> ResampleImageFilterType;
 183     typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
 184     filter->SetInput(input);
 185     filter->SetVerboseOptions(m_Verbose);
 186     filter->SetGaussianFilteringEnabled(false);
 187     filter->SetDefaultPixelValue(0);
 188
 189     const SizeType& input_size = input->GetLargestPossibleRegion().GetSize();
 190     SizeType output_size;
 191     for (unsigned int i = 0; i < Dimension; i++)
 192       output_size[i] = input_size[i];
 193     filter->SetOutputSize(output_size);
 194     if (m_Verbose) {
 195         std::cout << "Warning: The image size differs between the MHD file and the input dicom series. Performing resampling with default options using mhd size as reference (for advanced options, use clitkResampleImage)." << std::endl;
 196         std::cout << "MHD -> " << input->GetLargestPossibleRegion().GetSize() << std::endl;
 197         std::cout << "dicom -> " << reader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
 198     }
 199
 200     try {
 201       filter->Update();
 202       input = filter->GetOutput();
 203     } catch( itk::ExceptionObject & excp ) {
 204     std::cerr << "Error: Exception thrown while resampling!!" << std::endl;
 205     std::cerr << excp << std::endl;
 206     }
 207   }
 208 */
 209
 210
 211   //Read the dicom file to find the energy, the head & the steps (TODO: do it on the mhd filetext)
 212   gdcm::Reader hreader;
 213   hreader.SetFileName(m_ArgsInfo.inputModelFilename_arg);
 214   hreader.Read();
 215   gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
 216   int energyNumber, headNumber, rotationNumber; //Read the number of energy, of head and rotation steps
 217   gdcm::Attribute<0x54,0x11> energyNumber_att;
 218   energyNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
 219   energyNumber = energyNumber_att.GetValue();
 220
 221   gdcm::Attribute<0x54,0x21> headNumber_att;
 222   headNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
 223   headNumber = headNumber_att.GetValue();
 224
 225   gdcm::Tag squenceTag(0x54,0x52); //rotation, read a sequence first
 226   const gdcm::DataElement &sequence = hreader.GetFile().GetDataSet().GetDataElement( squenceTag );
 227   gdcm::DataSet &sequenceDataSet = sequence.GetValueAsSQ()->GetItem(1).GetNestedDataSet();
 228   gdcm::Tag rotationNumber_Tag(0x54,0x53);
 229   const gdcm::DataElement &rotationNumber_Element = sequenceDataSet.GetDataElement( rotationNumber_Tag );
 230   gdcm::Attribute<0x54,0x53> rotationNumber_att;
 231   rotationNumber_att.SetFromDataElement(rotationNumber_Element);
 232   rotationNumber = rotationNumber_att.GetValue();
 233
 234 energyNumber = 3;
 235 headNumber = 4;
 236 rotationNumber = 6;
 237   // Create a new mhd image with the dicom order slices
 238   typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
 239   mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
 240   mhdCorrectOrder->Allocate();
 241   unsigned int zAxis(0); //z value for the input mhd image
 242   for (unsigned int energy = 0; energy < energyNumber; ++energy) {
 243     for (unsigned int head = 0; head < headNumber; ++head) {
 244       for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
 245         std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
 246
 247         typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
 248         startIteratorIndexCorrectOrder[0] = 0;
 249         startIteratorIndexCorrectOrder[1] = 0;
 250         startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
 251
 252         typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
 253         startIteratorIndexOriginalOrder[0] = 0;
 254         startIteratorIndexOriginalOrder[1] = 0;
 255         startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
 256
 257         typename InputImageType::SizeType regionSizeIterator;
 258         regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
 259         regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
 260         regionSizeIterator[2] = 1;
 261
 262         typename InputImageType::RegionType regionIteratorCorrectOrder;
 263         regionIteratorCorrectOrder.SetSize(regionSizeIterator);
 264         regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
 265
 266         typename InputImageType::RegionType regionIteratorOriginalOrder;
 267         regionIteratorOriginalOrder.SetSize(regionSizeIterator);
 268         regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
 269
 270         itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
 271         itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
 272         while(!CorrectOrderIterator.IsAtEnd()) {
 273           CorrectOrderIterator.Set(OriginalOrderIterator.Get());
 274           ++CorrectOrderIterator;
 275           ++OriginalOrderIterator;
 276         }
 277
 278         ++zAxis;
 279       }
 280     }
 281   }
 282
 283
 284   // update output dicom keys/tags
 285   // string for distinguishing items inside sequence:
 286   const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE");
 287   std::string tempString = ITEM_ENCAPSULATE_STRING + "01";
 288   typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
 289   typename ReaderSeriesType::DictionaryArrayType outputArray;
 290   typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
 291   typename ReaderSeriesType::DictionaryRawPointer sequenceDict = new typename ReaderSeriesType::DictionaryType;
 292   // Dictionary for encapsulating the sequences.
 293   itk::GDCMImageIO::Pointer gdcmIOStructureSetRoi = itk::GDCMImageIO::New();
 294   typename ReaderSeriesType::DictionaryType &structureSetRoiSeq = gdcmIOStructureSetRoi->GetMetaDataDictionary();
 295   CopyDictionary (*inputDict, *outputDict);
 296   //CopyDictionary (*inputDict, *sequenceDict);
 297
 298   itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(energyNumber));
 299   itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
 300
 301   itk::EncapsulateMetaData<std::string>(*sequenceDict, "0054|0053", NumberToString(rotationNumber));
 302   //itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0053", NumberToString(rotationNumber));
 303   itk::EncapsulateMetaData<typename ReaderSeriesType::DictionaryType>(structureSetRoiSeq, tempString, *sequenceDict);
 304   itk::EncapsulateMetaData<typename ReaderSeriesType::DictionaryType>(*outputDict, "0054|0052",structureSetRoiSeq);
 305
 306 outputArray.push_back(outputDict);
 307   // Output directory and filenames
 308   //itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputFilename_arg ); // create if it doesn't exist
 309   typedef itk::ImageFileWriter<OutputImageType>  WriterType;
 310   typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
 311   typename WriterType::Pointer writer = WriterType::New();
 312   typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
 313
 314   writer->SetInput( mhdCorrectOrder );
 315   writer->SetImageIO( gdcmIO );
 316
 317   writer->SetFileName( filename_out );
 318   //writer->SetMetaDataDictionary(outputDict);
 319   writerSerie->SetInput( mhdCorrectOrder );
 320   writerSerie->SetImageIO( gdcmIO );
 321   typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
 322   fileNamesOutput.push_back(filename_out);
 323   writerSerie->SetFileNames( fileNamesOutput );
 324   writerSerie->SetMetaDataDictionaryArray(&outputArray);
 325
 326   // Write
 327   try {
 328     if (m_ArgsInfo.verbose_flag)
 329       std::cout << writer << std::endl;
 330     //writer->Update();
 331     writerSerie->Update();
 332   } catch( itk::ExceptionObject & excp ) {
 333     std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 334     std::cerr << excp << std::endl;
 335   }
 336 #else
 337   std::cout << "Use GDCM2" << std::endl;
 338 #endif
 339 }
 340
 341 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
 342 {
 343   typedef itk::MetaDataDictionary DictionaryType;
 344
 345   DictionaryType::ConstIterator itr = fromDict.Begin();
 346   DictionaryType::ConstIterator end = fromDict.End();
 347   typedef itk::MetaDataObject< std::string > MetaDataStringType;
 348
 349   while( itr != end )
 350     {
 351     itk::MetaDataObjectBase::Pointer  entry = itr->second;
 352
 353     MetaDataStringType::Pointer entryvalue =
 354       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 355     if( entryvalue )
 356       {
 357       std::string tagkey   = itr->first;
 358       std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 359       itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
 360       }
 361     ++itr;
 362     }
 363 }
 364
 365 template <typename T>
 366   std::string NumberToString ( T Number )
 367   {
 368      std::ostringstream ss;
 369      ss << Number;
 370      return ss.str();
 371   }
 372 }//end clitk
 373
 374 #endif //#define clitkGateSimulation2DicomGenericFilter_txx