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