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