tools/clitkImage2DicomDoseGenericFilter.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 clitkImage2DicomDoseGenericFilter_txx
  19 #define clitkImage2DicomDoseGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkImage2DicomDosemGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "math.h"
  31
  32 #include "clitkImage2DicomDoseGenericFilter.h"
  33 #include "clitkCommon.h"
  34 #include "itkMinimumMaximumImageCalculator.h"
  35 #include "itkImageRegionIterator.h"
  36 #include "itkImageRegionConstIterator.h"
  37 #if GDCM_MAJOR_VERSION >= 2
  38 #include "gdcmUIDGenerator.h"
  39 #include <gdcmImageHelper.h>
  40 #include <gdcmAttribute.h>
  41 #include <gdcmReader.h>
  42 #include <gdcmWriter.h>
  43 #else
  44 #include "gdcmFile.h"
  45 #include "gdcmUtil.h"
  46 #endif
  47
  48 //#include "gdcmBase.h"
  49 //#include "gdcmDocEntry.h"
  50
  51
  52
  53 namespace clitk
  54 {
  55
  56
  57 //-----------------------------------------------------------
  58 // Constructor
  59 //-----------------------------------------------------------
  60 template<class args_info_type>
  61 Image2DicomDoseGenericFilter<args_info_type>::Image2DicomDoseGenericFilter()
  62 {
  63   m_Verbose=false;
  64   m_InputFileName="";
  65 }
  66
  67
  68 //-----------------------------------------------------------
  69 // Update
  70 //-----------------------------------------------------------
  71 template<class args_info_type>
  72 void Image2DicomDoseGenericFilter<args_info_type>::Update()
  73 {
  74   // Read the Dimension and PixelType
  75   int Dimension;
  76   std::string PixelType;
  77   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  78
  79
  80   // Call UpdateWithDim
  81   if(Dimension==2) UpdateWithDim<2>(PixelType);
  82   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  83   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  84   else {
  85     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  86     return;
  87   }
  88 }
  89
  90 //-------------------------------------------------------------------
  91 // Update with the number of dimensions
  92 //-------------------------------------------------------------------
  93 template<class args_info_type>
  94 template<unsigned int Dimension>
  95 void
  96 Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
  97 {
  98   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  99
 100   if(PixelType == "short") {
 101     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
 102     UpdateWithDimAndPixelType<Dimension, signed short>();
 103   }
 104   else if(PixelType == "unsigned_short"){
 105     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
 106     UpdateWithDimAndPixelType<Dimension, unsigned short>();
 107   }
 108
 109   else if (PixelType == "unsigned_char") {
 110     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
 111     UpdateWithDimAndPixelType<Dimension, unsigned char>();
 112   }
 113
 114   //     else if (PixelType == "char"){
 115   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
 116   //       UpdateWithDimAndPixelType<Dimension, signed char>();
 117   //     }
 118   else if (PixelType == "double") {
 119     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
 120     UpdateWithDimAndPixelType<Dimension, double>();
 121   }
 122   else {
 123     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 124     UpdateWithDimAndPixelType<Dimension, float>();
 125   }
 126 }
 127
 128 //-------------------------------------------------------------------
 129 // Update with the number of dimensions and the pixeltype read from
 130 // the dicom files. The MHD files may be resampled to match the
 131 // dicom spacing (and number of slices). Rounding errors in resampling
 132 // are handled by removing files when generating the output dicom
 133 // series.
 134 //-------------------------------------------------------------------
 135 template<class args_info_type>
 136 template <unsigned int Dimension, class  PixelType>
 137 void
 138 Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 139 {
 140
 141 #if GDCM_MAJOR_VERSION >= 2
 142   // ImageTypes
 143   typedef itk::Image<PixelType, Dimension> InputImageType;
 144   typedef unsigned short int OutputPixelType;
 145   typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
 146   typedef itk::ImageFileReader< InputImageType >     ReaderType;
 147   typedef itk::ImageSeriesReader< InputImageType >     ReaderSeriesType;
 148   typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
 149   typedef itk::ImageRegionIterator< InputImageType > IteratorType;
 150   typedef itk::MinimumMaximumImageCalculator <InputImageType> ImageCalculatorFilterType;
 151   typedef itk::GDCMImageIO ImageIOType;
 152
 153   //-----------------------------------------------------------------------------
 154   // opening image input file
 155   typename ReaderType::Pointer reader = ReaderType::New();
 156   const char * filename = m_ArgsInfo.input_arg;
 157   reader->SetFileName( filename );
 158   reader->Update();
 159   typename InputImageType::Pointer image = reader->GetOutput();
 160   std::ostringstream value;
 161
 162   // Read Dicom model file
 163   typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
 164   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 165   std::string filename_out = m_ArgsInfo.output_arg;
 166   gdcmIO->LoadPrivateTagsOn();
 167   gdcmIO->KeepOriginalUIDOn();
 168   typename ReaderSeriesType::FileNamesContainer fileNames;
 169   fileNames.push_back(m_ArgsInfo.DicomInputFile_arg);
 170   readerSeries->SetImageIO( gdcmIO );
 171   readerSeries->SetFileNames( fileNames );
 172   try {
 173     readerSeries->Update();
 174   } catch (itk::ExceptionObject &excp) {
 175     std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
 176     std::cerr << excp << std::endl;
 177   }
 178
 179   // update output dicom keys/tags
 180   typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
 181   typename ReaderSeriesType::DictionaryArrayType outputArray;
 182   typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
 183   CopyDictionary (*inputDict, *outputDict);
 184
 185   // origin
 186   typename InputImageType::PointType origin = image->GetOrigin();
 187   value.str("");
 188   value<<origin[0]<<'\\'<<origin[1]<<'\\'<<origin[2];
 189   itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0032", value.str());
 190   DD(origin);
 191
 192   // orientation
 193   typename InputImageType::DirectionType direction = image->GetDirection();
 194   value.str("");
 195   value<<direction[0][0]<<'\\'<<direction[0][1]<<'\\'<<direction[0][2]<<'\\'<<direction[1][0]<<'\\'<<direction[1][1]<<'\\'<<direction[1][2];
 196   itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0037", value.str());
 197   DD(direction);
 198
 199   // size
 200   typename InputImageType::SizeType imageSize = image->GetLargestPossibleRegion().GetSize();
 201   //DD(imageSize);
 202   int NbCols=imageSize[0];      // col
 203   int NbRows=imageSize[1];      // row
 204   int NbFrames=imageSize[2];    // frame
 205   itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0008", NumberToString(NbFrames));
 206   itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0010", NumberToString(NbRows));
 207   itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0011", NumberToString(NbCols));
 208   DD(NbCols);
 209   DD(NbRows);
 210   DD(NbFrames);
 211
 212   // spacing
 213   typename InputImageType::SpacingType Spacing = image->GetSpacing();
 214   value.str("");
 215   value<<Spacing[0]<<'\\'<<Spacing[1];
 216   itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0030", value.str());
 217   value.str("");
 218   value<<Spacing[2];
 219   itk::EncapsulateMetaData<std::string>(*outputDict, "0018|0050", value.str());
 220   DD(Spacing);
 221
 222   // offset
 223   float offset = 0.;
 224   value.str("");
 225   value << offset;
 226   for (int i=1; i<NbFrames ; i++){
 227     offset+=Spacing[2];
 228     value << '\\';
 229     value << offset;
 230   }
 231   itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000c", value.str());
 232   DD(value.str());
 233
 234   // scaling
 235   typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
 236   imageCalculatorFilter->SetImage(image);
 237   imageCalculatorFilter->ComputeMaximum();
 238   float highestValue=imageCalculatorFilter->GetMaximum();
 239   double doseScaling = highestValue/(pow(2,16)-1);
 240   value.str("");
 241   value<<doseScaling;
 242   itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000e", value.str());
 243   DD(value.str());
 244
 245   // image data
 246   typename OutputImageType::Pointer imageOutput = OutputImageType::New();
 247   imageOutput->SetRegions(image->GetLargestPossibleRegion());
 248   imageOutput->SetSpacing(image->GetSpacing());
 249   imageOutput->SetOrigin(image->GetOrigin());
 250   imageOutput->SetDirection(image->GetDirection());
 251   imageOutput->Allocate();
 252
 253   typename itk::ImageRegionConstIterator<InputImageType> imageIterator(image,image->GetLargestPossibleRegion());
 254   typename itk::ImageRegionIterator<OutputImageType> imageOutputIterator(imageOutput,imageOutput->GetLargestPossibleRegion());
 255   while(!imageIterator.IsAtEnd()) {
 256     // Set the current pixel to white
 257     imageOutputIterator.Set((signed short int)(imageIterator.Get()/doseScaling));
 258
 259     ++imageIterator;
 260     ++imageOutputIterator;
 261   }
 262
 263   // Output directory and filenames
 264   typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
 265   outputArray.push_back(outputDict);
 266   writerSerie->SetInput( imageOutput );
 267   writerSerie->SetImageIO( gdcmIO );
 268   typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
 269   fileNamesOutput.push_back(filename_out);
 270   writerSerie->SetFileNames( fileNamesOutput );
 271   writerSerie->SetMetaDataDictionaryArray(&outputArray);
 272
 273   // Write
 274   try {
 275     if (m_ArgsInfo.verbose_flag)
 276       std::cout << writerSerie << std::endl;
 277     writerSerie->Update();
 278   } catch( itk::ExceptionObject & excp ) {
 279     std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 280     std::cerr << excp << std::endl;
 281   }
 282
 283   //Read sequence dicom tag with gdcm
 284   gdcm::Reader readerTemplateGDCM;
 285   readerTemplateGDCM.SetFileName( fileNames[0].c_str() );
 286   readerTemplateGDCM.Read();
 287   gdcm::File &fileTemplate = readerTemplateGDCM.GetFile();
 288   gdcm::DataSet &dsTemplate = fileTemplate.GetDataSet();
 289   const unsigned int ptr_len = 42;
 290   char *ptrTemplate = new char[ptr_len];
 291   memset(ptrTemplate,0,ptr_len);
 292
 293   const gdcm::DataElement &referenceRTPlanSq = dsTemplate.GetDataElement(gdcm::Tag(0x300c, 0x02));
 294
 295   //Create the Dose Grid Scaling data element (ITK 4.13 do not take into account - it works well with ITK 4.5.1)
 296   gdcm::DataElement deDoseGridScaling( gdcm::Tag(0x3004,0x0e) );
 297   deDoseGridScaling.SetVR( gdcm::VR::DS );
 298   deDoseGridScaling.SetByteValue(NumberToString(doseScaling).c_str(), (uint32_t)strlen(NumberToString(doseScaling).c_str()));
 299
 300   //Copy/Write sequence dicom tag with gdcm
 301   gdcm::Reader readerOutputGDCM;
 302   readerOutputGDCM.SetFileName( fileNamesOutput[0].c_str() );
 303   readerOutputGDCM.Read();
 304   gdcm::File &file = readerOutputGDCM.GetFile();
 305   gdcm::DataSet &dsOutput = file.GetDataSet();
 306
 307   dsOutput.Insert(referenceRTPlanSq);
 308   dsOutput.Replace(deDoseGridScaling);
 309   gdcm::Writer w;
 310   w.SetFile( file );
 311   w.SetFileName( fileNamesOutput[0].c_str() );
 312   w.Write();
 313
 314 //---------------------------------------------------------------------------------------
 315 //WRITE DICOM BIS
 316 // The previous way of writting DICOM-RT-DOSE works only for ITK
 317 // and resulting RT-DOSE files are not readable by commercial systems.
 318 // The nex step is to copy again the output file with another syntax, allowing to make RT-DOSE files readable by commercial systems.
 319 // see Jean-Pierre Roux comment below.
 320
 321 /*gdcm::FileHelper *fh = new gdcm::FileHelper(args_info.OutputFile_arg);
 322    void *imageData;
 323    int dataSize;
 324
 325    dataSize  = fh->GetImageDataRawSize();
 326    imageData = fh->GetImageDataRaw();// somewhat important : Loads the Pixels in memory !
 327
 328    fh->SetWriteModeToRaw();
 329    fh->SetWriteTypeToDcmExplVR();
 330
 331    bool res = fh->Write(args_info.OutputFile_arg);
 332
 333    if(!res)
 334       std::cout <<"Fail to write [" << args_info.OutputFile_arg << "]" <<std::endl;
 335    else std::cout<<"\n DICOM File re-written, using the FileHelper syntax, in order to be processed by commercial systems !"<<std::endl;
 336
 337 delete fh; */
 338 /*  gdcm::Writer w;
 339   w.SetFile(mDCMFile);
 340   w.SetFileName(m_ArgsInfo.output_arg);
 341   w.Write();*/
 342 //   fh->Delete();
 343
 344 //---------------------------------------------------------------------------------------
 345 //---------------------------------------------------------------------------------------
 346 // Jean-Pierre Roux help for DICOM-RT-DOSE writting
 347 //---------------------------------------------------------------------------------------
 348 //---------------------------------------------------------------------------------------
 349
 350 //      de      Jean-Pierre Roux <jpr@creatis.insa-lyon.fr>
 351 // répondre à jpr@creatis.insa-lyon.fr
 352 // à   Loic Grevillot <loic.grevillot@gmail.com>
 353 // cc   jpr@creatis.insa-lyon.fr,
 354 // Joël Schaerer <joel.schaerer@gmail.com>,
 355 // David Sarrut <David.Sarrut@creatis.insa-lyon.fr>,
 356 // Joel Schaerer <joel.schaerer@creatis.insa-lyon.fr>
 357 // date 12 juillet 2010 12:23
 358 // objet        Re: DICOM RT DOSE
 359 //
 360 // masquer les détails 12:23 (Il y a 21 heures)
 361 //
 362 // Bonjour,
 363 //
 364 // J'aurais écrit à peut prèt la même chose.
 365 // (Ci après un extrait de mon code -Example/ReWrite.cxx-)
 366 //
 367 // L'utilisation d'un FileHelper (qui ne changera rien dans ce cas précis) est une mesure de précaution, car,  l'élément 7FE0|0010 peut être compressé (ce qui n'est pas le cas pour tes images), que la manière de stocker les pixels ainsi compressés était parfois un peu ... curieuse.
 368 // Je décompresse, et réécrit non compressé.
 369 //
 370 // ============================
 371 //    gdcm::File *f = new gdcm::File();
 372 //    f->SetMaxSizeLoadEntry(0x7fffffff);
 373 //    f->SetFileName( fileName );
 374 //    bool res = f->Load();
 375 //    if ( !res )
 376 //    {
 377 //       delete f;
 378 //       return 0;
 379 //    }
 380 //
 381 //    if (!f->IsReadable())
 382 //    {
 383 //        std::cerr << "Sorry, not a Readable DICOM / ACR File"  <<std::endl;
 384 //        delete f;
 385 //        return 0;
 386 //    }
 387 //
 388 //    gdcm::FileHelper *fh = new gdcm::FileHelper(f);
 389 //    void *imageData;
 390 //    int dataSize;
 391 //
 392 //    dataSize  = fh->GetImageDataRawSize();
 393 //    imageData = fh->GetImageDataRaw();// somewhat important : Loads the Pixels in memory !
 394 //
 395 //    fh->SetWriteModeToRaw();
 396 //    fh->SetWriteTypeToDcmExplVR();
 397 //    res = fh->Write(outputFileName);
 398 //
 399 //    if(!res)
 400 //       std::cout <<"Fail to write [" << outputFileName << "]" <<std::endl;
 401 //
 402 //    f->Delete();
 403 //    fh->Delete();
 404
 405 //---------------------------------------------------------------------------------------
 406 //---------------------------------------------------------------------------------------
 407 // pour supprimer des tags:
 408 //---------------------------------------------------------------------------------------
 409 //---------------------------------------------------------------------------------------
 410 //      SOLUCE 1 QUI MARCHE POUR UN SQItem SIMPLE
 411 /*
 412 gdcm::DocEntry *a;
 413 //a= ((gdcm::SQItem*)mDCMFile)->GetDocEntry(0xfffe,0xe00d);
 414 a= ((gdcm::SQItem*)mDCMFile)->GetDocEntry(0x3004,0x000e);
 415 ((gdcm::SQItem*)mDCMFile)->RemoveEntry(a);
 416 mDCMFile->Write (args_info.OutputFile_arg, type);
 417 */
 418
 419 //      SOLUCE 2 QUI MARCHE POUR UNE SeqEntry->SQItem
 420 /*
 421 std::cout<<"\ntest correction fichier apres ecriture\n"<<std::endl;
 422 gdcm::SeqEntry *seqEntry = mDCMFile->GetSeqEntry(0x300c,0x0002);
 423 gdcm::SQItem* currentItem = seqEntry->GetFirstSQItem();
 424 gdcm::DocEntry *a;
 425 //a= currentItem->GetDocEntry(0x0008,0x1155);
 426 a= currentItem->GetDocEntry(0xfffe,0xe00d);
 427 currentItem->RemoveEntry(a);
 428 mDCMFile->Write (args_info.OutputFile_arg, type);
 429 */
 430
 431 //gdcm::DocEntry *a;
 432 //a=GetDocEntry(0x7fe0,0x0000);
 433 //((gdcm::SQItem*)mDCMFile)->RemoveEntry(a);
 434
 435 //-----------------------------------------------------------------------------------
 436   std::cout <<"\n## DICOM Image to RT DOSE application is ended..."<<std::endl;
 437   std::cout <<"#########################################################\n" << std::endl;
 438
 439 #else
 440   std::cout << "Use GDCM2" << std::endl;
 441 #endif
 442 }
 443 //---------------------------------------------------------------------------
 444
 445
 446 //---------------------------------------------------------------------------
 447 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
 448 {
 449   typedef itk::MetaDataDictionary DictionaryType;
 450
 451   DictionaryType::ConstIterator itr = fromDict.Begin();
 452   DictionaryType::ConstIterator end = fromDict.End();
 453   typedef itk::MetaDataObject< std::string > MetaDataStringType;
 454
 455   while( itr != end )
 456     {
 457     itk::MetaDataObjectBase::Pointer  entry = itr->second;
 458
 459     MetaDataStringType::Pointer entryvalue =
 460       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 461     if( entryvalue )
 462       {
 463       std::string tagkey   = itr->first;
 464       std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 465       itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
 466       }
 467     ++itr;
 468     }
 469 }
 470 //---------------------------------------------------------------------------
 471
 472 //---------------------------------------------------------------------------
 473 template <typename T> std::string NumberToString ( T Number )
 474 {
 475    std::ostringstream ss;
 476    ss << Number;
 477    return ss.str();
 478 }
 479 //---------------------------------------------------------------------------
 480
 481
 482 }//end clitk
 483
 484 #endif //#define clitkImage2DicomDoseGenericFilter_txx