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