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 ReaderType::Pointer reader = ReaderType::New();
 148   typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
 149   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 150   std::string filename_out = m_ArgsInfo.outputFilename_arg;
 151   gdcmIO->LoadPrivateTagsOn();
 152   gdcmIO->KeepOriginalUIDOn();
 153   reader->SetImageIO( gdcmIO );
 154   reader->SetFileName( m_ArgsInfo.inputModelFilename_arg );
 155   typename ReaderSeriesType::FileNamesContainer fileNames;
 156   fileNames.push_back(m_ArgsInfo.inputModelFilename_arg);
 157   readerSeries->SetImageIO( gdcmIO );
 158   readerSeries->SetFileNames( fileNames );
 159   try {
 160     reader->Update();
 161     readerSeries->Update();
 162   } catch (itk::ExceptionObject &excp) {
 163     std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
 164     std::cerr << excp << std::endl;
 165   }
 166
 167
 168   // Read the input (MHD file)
 169   typedef typename InputImageType::RegionType RegionType;
 170   typedef typename RegionType::SizeType SizeType;
 171   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 172   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
 173   volumeReader->SetFileName( m_InputFileName);
 174   volumeReader->Update();
 175   typename InputImageType::Pointer input = volumeReader->GetOutput();
 176
 177
 178   //Read the metadata informations in the mhd
 179   MetaObject tObj(3);
 180   tObj.AddUserField("NumberOfEnergyWindows", MET_STRING);
 181   tObj.AddUserField("NbOfHeads", MET_STRING);
 182   tObj.AddUserField("NbOfProjections", MET_STRING);
 183   tObj.Read(m_InputFileName.c_str());
 184   char *p_energyNumber, *p_headNumber, *p_rotationNumber;
 185   int energyNumber, headNumber, rotationNumber;
 186   p_energyNumber = static_cast<char*>(tObj.GetUserField("NumberOfEnergyWindows"));
 187   p_headNumber = static_cast<char*>(tObj.GetUserField("NbOfHeads"));
 188   p_rotationNumber = static_cast<char*>(tObj.GetUserField("NbOfProjections"));
 189   energyNumber = atoi(p_energyNumber);
 190   headNumber = atoi(p_headNumber);
 191   rotationNumber = atoi(p_rotationNumber);
 192   for (unsigned int i=0; i<energyNumber; ++i) {
 193     std::string tagEnergyName("EnergyWindow_");
 194     tagEnergyName += NumberToString(i).c_str();
 195     std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
 196     tagEnergyThreshold += "_threshold";
 197     tagEnergyUphold += "_uphold";
 198     tObj.AddUserField(tagEnergyName.c_str(), MET_STRING);
 199     tObj.AddUserField(tagEnergyThreshold.c_str(), MET_STRING);
 200     tObj.AddUserField(tagEnergyUphold.c_str(), MET_STRING);
 201   }
 202   tObj.Read(m_InputFileName.c_str());
 203   std::vector<char*> p_EnergyWindowName(energyNumber), p_EnergyWindowThreshold(energyNumber), p_EnergyWindowUphold(energyNumber);
 204   std::vector<int> energyWindowThreshold(energyNumber), energyWindowUphold(energyNumber);
 205   for (unsigned int i=0; i<energyNumber; ++i) {
 206     std::string tagEnergyName("EnergyWindow_");
 207     tagEnergyName += NumberToString(i).c_str();
 208     std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
 209     tagEnergyThreshold += "_threshold";
 210     tagEnergyUphold += "_uphold";
 211     p_EnergyWindowName[i] = static_cast<char*>(tObj.GetUserField(tagEnergyName.c_str()));
 212     p_EnergyWindowThreshold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyThreshold.c_str()));
 213     p_EnergyWindowUphold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyUphold.c_str()));
 214     energyWindowThreshold[i] = atoi(p_EnergyWindowThreshold[i]);
 215     energyWindowUphold[i] = atoi(p_EnergyWindowUphold[i]);
 216   }
 217
 218 //TODO if the input size/spacing and dicom model ones are different
 219
 220   // Create a new mhd image with the correct dicom order slices
 221   typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
 222   mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
 223   mhdCorrectOrder->Allocate();
 224   unsigned int zAxis(0); //z value for the input mhd image
 225   for (unsigned int energy = 0; energy < energyNumber; ++energy) {
 226     for (unsigned int head = 0; head < headNumber; ++head) {
 227       for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
 228         std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
 229
 230         typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
 231         startIteratorIndexCorrectOrder[0] = 0;
 232         startIteratorIndexCorrectOrder[1] = 0;
 233         startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
 234
 235         typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
 236         startIteratorIndexOriginalOrder[0] = 0;
 237         startIteratorIndexOriginalOrder[1] = 0;
 238         startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
 239
 240         typename InputImageType::SizeType regionSizeIterator;
 241         regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
 242         regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
 243         regionSizeIterator[2] = 1;
 244
 245         typename InputImageType::RegionType regionIteratorCorrectOrder;
 246         regionIteratorCorrectOrder.SetSize(regionSizeIterator);
 247         regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
 248
 249         typename InputImageType::RegionType regionIteratorOriginalOrder;
 250         regionIteratorOriginalOrder.SetSize(regionSizeIterator);
 251         regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
 252
 253         itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
 254         itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
 255         while(!CorrectOrderIterator.IsAtEnd()) {
 256           CorrectOrderIterator.Set(OriginalOrderIterator.Get());
 257           ++CorrectOrderIterator;
 258           ++OriginalOrderIterator;
 259         }
 260
 261         ++zAxis;
 262       }
 263     }
 264   }
 265
 266
 267   // update output dicom keys/tags
 268   // string for distinguishing items inside sequence:
 269   const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE");
 270   std::string tempString = ITEM_ENCAPSULATE_STRING + "01";
 271   typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
 272   typename ReaderSeriesType::DictionaryArrayType outputArray;
 273   typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
 274   CopyDictionary (*inputDict, *outputDict);
 275   itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(energyNumber));
 276   itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
 277   outputArray.push_back(outputDict);
 278
 279
 280   // Output directory and filenames
 281   typedef itk::ImageFileWriter<OutputImageType>  WriterType;
 282   typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
 283   typename WriterType::Pointer writer = WriterType::New();
 284   typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
 285   writer->SetInput( mhdCorrectOrder );
 286   writer->SetImageIO( gdcmIO );
 287   writer->SetFileName( filename_out );
 288   writerSerie->SetInput( mhdCorrectOrder );
 289   writerSerie->SetImageIO( gdcmIO );
 290   typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
 291   fileNamesOutput.push_back(filename_out);
 292   writerSerie->SetFileNames( fileNamesOutput );
 293   writerSerie->SetMetaDataDictionaryArray(&outputArray);
 294
 295
 296   // Write
 297   try {
 298     if (m_ArgsInfo.verbose_flag)
 299       std::cout << writer << std::endl;
 300     //writer->Update();
 301     writerSerie->Update();
 302   } catch( itk::ExceptionObject & excp ) {
 303     std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 304     std::cerr << excp << std::endl;
 305   }
 306
 307
 308   //Write sequence dicom tag with gdcm
 309   gdcm::Reader reader2;
 310   reader2.SetFileName( fileNamesOutput[0].c_str() );
 311   reader2.Read();
 312   gdcm::File &file = reader2.GetFile();
 313   gdcm::DataSet &ds2 = file.GetDataSet();
 314   const unsigned int ptr_len = 42;
 315   char *ptr = new char[ptr_len];
 316   memset(ptr,0,ptr_len);
 317
 318   //Write rotation tag
 319   // Create a Sequence
 320   gdcm::SmartPointer<gdcm::SequenceOfItems> rotationSq = new gdcm::SequenceOfItems();
 321   rotationSq->SetLengthToUndefined();
 322   // Create a dataelement
 323   gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
 324   rotationDE.SetVR( gdcm::VR::US );
 325   char essai = (char)rotationNumber;
 326   char *p_essai = &essai;
 327   rotationDE.SetByteValue(p_essai, 1);
 328   // Create an item
 329   gdcm::Item rotationIt;
 330   rotationIt.SetVLToUndefined();
 331   gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
 332   rotationDS.Insert(rotationDE);
 333   rotationSq->AddItem(rotationIt);
 334   // Insert sequence into data set
 335   gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) );
 336   rotationDEParent.SetVR(gdcm::VR::SQ);
 337   rotationDEParent.SetValue(*rotationSq);
 338   rotationDEParent.SetVLToUndefined();
 339
 340   //Write energy
 341   gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) );
 342   energyDEParent.SetVR(gdcm::VR::SQ);
 343   // Create a Sequence
 344   gdcm::SmartPointer<gdcm::SequenceOfItems> energySq = new gdcm::SequenceOfItems();
 345   energySq->SetLengthToUndefined();
 346   for (unsigned int i=0; i<energyNumber; ++i) {
 347     gdcm::SmartPointer<gdcm::SequenceOfItems> energyThresholdSq = new gdcm::SequenceOfItems();
 348     energyThresholdSq->SetLengthToUndefined();
 349     // Create a dataelement
 350     gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) );
 351     gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) );
 352     energyThresholdDE.SetVR( gdcm::VR::DS );
 353     energyUpholdDE.SetVR( gdcm::VR::DS );
 354     energyThresholdDE.SetByteValue(p_EnergyWindowThreshold[i], (uint32_t)strlen(p_EnergyWindowThreshold[i]));
 355     energyUpholdDE.SetByteValue(p_EnergyWindowUphold[i], (uint32_t)strlen(p_EnergyWindowUphold[i]));
 356     // Create an item
 357     gdcm::Item energyThresholdIt;
 358     energyThresholdIt.SetVLToUndefined();
 359     gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
 360     energyThresholdDS.Insert(energyThresholdDE);
 361     energyThresholdDS.Insert(energyUpholdDE);
 362     energyThresholdSq->AddItem(energyThresholdIt);
 363     // Insert sequence into data set
 364     gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) );
 365     energyThresholdDEParent.SetVR(gdcm::VR::SQ);
 366     energyThresholdDEParent.SetValue(*energyThresholdSq);
 367     energyThresholdDEParent.SetVLToUndefined();
 368     // Create a dataelement
 369     gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) );
 370     energyNameDE.SetVR( gdcm::VR::SH );
 371     energyNameDE.SetByteValue(p_EnergyWindowName[i], (uint32_t)strlen(p_EnergyWindowName[i]));
 372     // Create an item
 373     gdcm::Item energyIt;
 374     energyIt.SetVLToUndefined();
 375     gdcm::DataSet &energyDS = energyIt.GetNestedDataSet();
 376     energyDS.Insert(energyNameDE);
 377     energyDS.Insert(energyThresholdDEParent);
 378     energySq->AddItem(energyIt);
 379   }
 380   // Insert sequence into data set
 381   energyDEParent.SetValue(*energySq);
 382   energyDEParent.SetVLToUndefined();
 383   ds2.Insert(energyDEParent);
 384   ds2.Insert(rotationDEParent);
 385
 386   gdcm::Writer w;
 387   w.SetFile( file );
 388   w.SetFileName( fileNamesOutput[0].c_str() );
 389   w.Write();
 390
 391 #else
 392   std::cout << "Use GDCM2" << std::endl;
 393 #endif
 394 }
 395 //---------------------------------------------------------------------------
 396
 397
 398 //---------------------------------------------------------------------------
 399 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
 400 {
 401   typedef itk::MetaDataDictionary DictionaryType;
 402
 403   DictionaryType::ConstIterator itr = fromDict.Begin();
 404   DictionaryType::ConstIterator end = fromDict.End();
 405   typedef itk::MetaDataObject< std::string > MetaDataStringType;
 406
 407   while( itr != end )
 408     {
 409     itk::MetaDataObjectBase::Pointer  entry = itr->second;
 410
 411     MetaDataStringType::Pointer entryvalue =
 412       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 413     if( entryvalue )
 414       {
 415       std::string tagkey   = itr->first;
 416       std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 417       itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
 418       }
 419     ++itr;
 420     }
 421 }
 422 //---------------------------------------------------------------------------
 423
 424
 425 //---------------------------------------------------------------------------
 426 template <typename T> std::string NumberToString ( T Number )
 427 {
 428    std::ostringstream ss;
 429    ss << Number;
 430    return ss.str();
 431 }
 432 //---------------------------------------------------------------------------
 433
 434 }//end clitk
 435
 436 #endif //#define clitkGateSimulation2DicomGenericFilter_txx