X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkImage2DicomDoseGenericFilter.txx;h=df105fadfd1a75e1838d057dbd39ee6b9fdb0edb;hb=543b72e23ad001ac2a7743b9beacf48e2d0054ac;hp=4b865f60e0c0872b9adc41ea6954e13b8dd70c28;hpb=af9a6251e337d999b569f152dfc1b3b9f0da2759;p=clitk.git

diff --git a/tools/clitkImage2DicomDoseGenericFilter.txx b/tools/clitkImage2DicomDoseGenericFilter.txx
index 4b865f6..df105fa 100644
--- a/tools/clitkImage2DicomDoseGenericFilter.txx
+++ b/tools/clitkImage2DicomDoseGenericFilter.txx
@@ -31,6 +31,9 @@
 
 #include "clitkImage2DicomDoseGenericFilter.h"
 #include "clitkCommon.h"
+#include "itkMinimumMaximumImageCalculator.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
 #if GDCM_MAJOR_VERSION >= 2
 #include "gdcmUIDGenerator.h"
 #include <gdcmImageHelper.h>
@@ -135,15 +138,27 @@ void
 Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 {
 
-#if GDCM_MAJOR_VERSION == 2
+#if GDCM_MAJOR_VERSION >= 2
   // ImageTypes
   typedef itk::Image<PixelType, Dimension> InputImageType;
-  typedef itk::Image<PixelType, Dimension> OutputImageType;
+  typedef unsigned short int OutputPixelType;
+  typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
   typedef itk::ImageFileReader< InputImageType >     ReaderType;
   typedef itk::ImageSeriesReader< InputImageType >     ReaderSeriesType;
+  typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
   typedef itk::ImageRegionIterator< InputImageType > IteratorType;
+  typedef itk::MinimumMaximumImageCalculator <InputImageType> ImageCalculatorFilterType;
   typedef itk::GDCMImageIO ImageIOType;
 
+  //-----------------------------------------------------------------------------
+  // opening image input file
+  typename ReaderType::Pointer reader = ReaderType::New();
+  const char * filename = m_ArgsInfo.input_arg;
+  reader->SetFileName( filename );
+  reader->Update();
+  typename InputImageType::Pointer image = reader->GetOutput();
+  std::ostringstream value;
+
   // Read Dicom model file
   typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
   ImageIOType::Pointer gdcmIO = ImageIOType::New();
@@ -161,270 +176,141 @@ Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
     std::cerr << excp << std::endl;
   }
 
-
-
-//-----------------------------------------------------------------------------------
-// opening dicom input file
-  gdcm::Reader reader2;
-  reader2.SetFileName( m_ArgsInfo.DicomInputFile_arg );
-  reader2.Read();
-  gdcm::File &mDCMFile = reader2.GetFile();
-  gdcm::DataSet &ds = mDCMFile.GetDataSet();
-//mDCMFile.SetMaxSizeLoadEntry(1006384); // important size required, otherwise some data are not loaded
-//mDCMFile.AddForceLoadElement(0x7fe0,0x0010); //Load pixel data no matter its size
-
-std::cout << "File:   "<< m_ArgsInfo.DicomInputFile_arg << "   loaded !"<< std::endl;
-
-
-
-//-----------------------------------------------------------------------------
-// opening image input file
-typename ReaderType::Pointer reader = ReaderType::New();
-const char * filename = m_ArgsInfo.input_arg;
-reader->SetFileName( filename );
-reader->Update();
-typename InputImageType::Pointer image = reader->GetOutput();
-
-// origin
-typename InputImageType::PointType origin = image->GetOrigin();
-DD(origin);
-
-// size
-typename InputImageType::SizeType imageSize = image->GetLargestPossibleRegion().GetSize();
-//DD(imageSize);
-int NbCols=imageSize[0];	// col
-int NbRows=imageSize[1];	// row
-int NbFrames=imageSize[2];	// frame
-DD(NbCols);
-DD(NbRows);
-DD(NbFrames);
-
-// spacing
-typename InputImageType::SpacingType Spacing = image->GetSpacing();
-DD(Spacing);
-
-// scaling
-float highestValue=pow(10,-10);
-IteratorType out( image, image->GetRequestedRegion() );
-for (out.GoToBegin(); !out.IsAtEnd(); ++out){
-//DD(out.Get());
-  if (out.Get()>highestValue) highestValue=out.Get();
-}
-double doseScaling = highestValue/(pow(2,16)-1);
-DD(doseScaling);
-
-// image data
-std::vector<unsigned short int> ImageData;
-typename InputImageType::IndexType pixelIndex;
-int l=0;
-unsigned short int pixelValue;
-//DD(highestValue);
-for (int i=0; i<NbFrames; i++){
-  pixelIndex[2] = i;
-  for (int j=0; j<NbRows; j++){
-    pixelIndex[1] = j;
-    for (int k=0; k<NbCols; k++){
-	pixelIndex[0] = k;
-	pixelValue=image->GetPixel(pixelIndex)/doseScaling;
-if(float(image->GetPixel(pixelIndex)/doseScaling)>(pow(2,16)-1.)) {
-std::cout<<"\n!!!!! WARNING !!!!! pixel index: "<<pixelIndex<<"unsigned short int capacity ful or overfuled => Highest value may become 0"<<std::endl;
-DD(pixelIndex);
-DD(image->GetPixel(pixelIndex));
-//DD(image->GetPixel(pixelIndex)/doseScaling);
-DD(pixelValue);
-std::cout<<"Pixel Value should be equal to "<<(pow(2,16)-1)<<" but should not be 0"<<std::endl;
-std::cout<<"\n"<<std::endl;
-//assert(pixelValue<=(pow(2,16)-1));	should work, but do not...
-}
-//DD(pixelValue);
-        ImageData.push_back(pixelValue);
-        l++;
-    }
-  }
-}
-DD(ImageData.size());
-
-// Relevant parameters inserted in the new dicom file
-/*
-ImagePosition
-NbCols
-NbRows
-NbFrames
-Spacing
-ImageData
-doseScaling
-*/
-
+  // update output dicom keys/tags
   typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
   typename ReaderSeriesType::DictionaryArrayType outputArray;
   typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
   CopyDictionary (*inputDict, *outputDict);
-  outputArray.push_back(outputDict);
+
+  // origin
+  typename InputImageType::PointType origin = image->GetOrigin();
+  value.str("");
+  value<<origin[0]<<'\\'<<origin[1]<<'\\'<<origin[2];
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0032", value.str());
+  DD(origin);
+
+  // orientation
+  typename InputImageType::DirectionType direction = image->GetDirection();
+  value.str("");
+  value<<direction[0][0]<<'\\'<<direction[0][1]<<'\\'<<direction[0][2]<<'\\'<<direction[1][0]<<'\\'<<direction[1][1]<<'\\'<<direction[1][2];
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0037", value.str());
+  DD(direction);
+
+  // size
+  typename InputImageType::SizeType imageSize = image->GetLargestPossibleRegion().GetSize();
+  //DD(imageSize);
+  int NbCols=imageSize[0];	// col
+  int NbRows=imageSize[1];	// row
+  int NbFrames=imageSize[2];	// frame
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0008", NumberToString(NbFrames));
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0010", NumberToString(NbRows));
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0011", NumberToString(NbCols));
+  DD(NbCols);
+  DD(NbRows);
+  DD(NbFrames);
+
+  // spacing
+  typename InputImageType::SpacingType Spacing = image->GetSpacing();
+  value.str("");
+  value<<Spacing[0]<<'\\'<<Spacing[1];
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0030", value.str());
+  value.str("");
+  value<<Spacing[2];
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0018|0050", value.str());
+  DD(Spacing);
+
+  // offset
+  float offset = 0.;
+  value.str("");
+  value << offset;
+  for (int i=1; i<NbFrames ; i++){
+    offset+=Spacing[2];
+    value << '\\';
+    value << offset;
+  }
+  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000c", value.str());
+  DD(value.str());
+
+  // scaling
+  typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
+  imageCalculatorFilter->SetImage(image);
+  imageCalculatorFilter->ComputeMaximum();
+  float highestValue=imageCalculatorFilter->GetMaximum();
+  double doseScaling = highestValue/(pow(2,16)-1);
+  value.str("");
+  value<<doseScaling;
+  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000e", value.str());
+  DD(value.str());
+
+  // image data
+  typename OutputImageType::Pointer imageOutput = OutputImageType::New();
+  imageOutput->SetRegions(image->GetLargestPossibleRegion());
+  imageOutput->SetSpacing(image->GetSpacing());
+  imageOutput->SetOrigin(image->GetOrigin());
+  imageOutput->SetDirection(image->GetDirection());
+  imageOutput->Allocate();
+
+  typename itk::ImageRegionConstIterator<InputImageType> imageIterator(image,image->GetLargestPossibleRegion());
+  typename itk::ImageRegionIterator<OutputImageType> imageOutputIterator(imageOutput,imageOutput->GetLargestPossibleRegion());
+  while(!imageIterator.IsAtEnd()) {
+    // Set the current pixel to white
+    imageOutputIterator.Set((signed short int)(imageIterator.Get()/doseScaling));
+
+    ++imageIterator;
+    ++imageOutputIterator;
+  }
 
   // Output directory and filenames
-  typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
   typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
-  writerSerie->SetInput( image );
+  outputArray.push_back(outputDict);
+  writerSerie->SetInput( imageOutput );
   writerSerie->SetImageIO( gdcmIO );
   typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
   fileNamesOutput.push_back(filename_out);
   writerSerie->SetFileNames( fileNamesOutput );
   writerSerie->SetMetaDataDictionaryArray(&outputArray);
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-//--------------------------------------------------------------
-std::cout<<"\nECRITURE DU FICHIER DICOM !"<<std::endl;
-
-//gdcm::ValEntry *b;
-std::string Value("");
-
-
-gdcm::DataElement DE;
-
-DE = gdcm::Tag(0x20, 0x32);
-Value = origin[0];
-Value += "\\";
-Value += origin[1];
-Value += "\\";
-Value += origin[2];
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x28, 0x11);
-Value = NbCols;
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x28, 0x10);
-Value = NbRows;
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x28, 0x08);
-Value = NbFrames;
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x3004, 0x0e);
-Value = doseScaling;
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x28, 0x30);
-Value = Spacing[0];
-Value += "\\";
-Value += Spacing[1];
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-DE = gdcm::Tag(0x3004, 0x000c);
-float offset = 0.;
-Value = offset;
-  for (int i=1; i<NbFrames ; i++){
-    offset+=Spacing[2];
-    Value += "\\";
-    Value += offset;
+  // Write
+  try {
+    if (m_ArgsInfo.verbose_flag)
+      std::cout << writerSerie << std::endl;
+    writerSerie->Update();
+  } catch( itk::ExceptionObject & excp ) {
+    std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
+    std::cerr << excp << std::endl;
   }
 
-DE.SetVR( gdcm::VR::US );
-DE.SetByteValue(Value.c_str(), 1);
-ds.Insert(DE);
-DD(Value);
-Value = "";
-
-/*
-// NbCols
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x0028,0x0011);
-Value<<NbCols;
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");
-
-// NbRows
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x0028,0x0010);
-Value<<NbRows;
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");
-//DD(Value.str());
-
-// NbFrames
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x0028,0x0008);
-Value<<NbFrames;
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");
-
-// doseScaling
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x3004,0x000e);
-Value<<doseScaling;
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");
-
-// Spacing X Y
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x0028, 0x0030);
-Value<<Spacing[0]<<'\\'<<Spacing[1];
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");
-
-// Spacing Z ([Grid Frame Offset Vector])
-b = ((gdcm::SQItem*)mDCMFile)->GetValEntry(0x3004, 0x000c);
-float offset=0.;
-Value<<offset;
-  for (int i=1; i<NbFrames ; i++){
-    offset+=Spacing[2];
-    Value<<'\\'<<offset;
-  }
-b->SetValue(Value.str());
-DD(Value.str());
-Value.str("");	
-*/
-//ImageData
-//bool data = mDCMFile->SetBinEntry(reinterpret_cast<uint8_t*>( &(ImageData[0]) ) , (int)(sizeof(unsigned short int) * ImageData.size()) , 0x7fe0, 0x0010);
-//if (data)  std::cout<<"\n DICOM dose data written !"<<std::endl;
+  //Read sequence dicom tag with gdcm
+  gdcm::Reader readerTemplateGDCM;
+  readerTemplateGDCM.SetFileName( fileNames[0].c_str() );
+  readerTemplateGDCM.Read();
+  gdcm::File &fileTemplate = readerTemplateGDCM.GetFile();
+  gdcm::DataSet &dsTemplate = fileTemplate.GetDataSet();
+  const unsigned int ptr_len = 42;
+  char *ptrTemplate = new char[ptr_len];
+  memset(ptrTemplate,0,ptr_len);
+
+  const gdcm::DataElement &referenceRTPlanSq = dsTemplate.GetDataElement(gdcm::Tag(0x300c, 0x02));
+
+  //Create the Dose Grid Scaling data element (ITK 4.13 do not take into account - it works well with ITK 4.5.1)
+  gdcm::DataElement deDoseGridScaling( gdcm::Tag(0x3004,0x0e) );
+  deDoseGridScaling.SetVR( gdcm::VR::DS );
+  deDoseGridScaling.SetByteValue(NumberToString(doseScaling).c_str(), (uint32_t)strlen(NumberToString(doseScaling).c_str()));
+
+  //Copy/Write sequence dicom tag with gdcm
+  gdcm::Reader readerOutputGDCM;
+  readerOutputGDCM.SetFileName( fileNamesOutput[0].c_str() );
+  readerOutputGDCM.Read();
+  gdcm::File &file = readerOutputGDCM.GetFile();
+  gdcm::DataSet &dsOutput = file.GetDataSet();
+
+  dsOutput.Insert(referenceRTPlanSq);
+  dsOutput.Replace(deDoseGridScaling);
+  gdcm::Writer w;
+  w.SetFile( file );
+  w.SetFileName( fileNamesOutput[0].c_str() );
+  w.Write();
 
-//---------------------------------------------------------------------------------------
-//WRITE DICOM
-/*gdcm::FileType type = mDCMFile->GetFileType();
-//type=(gdcm::FileType)2;
-bool ecriture = mDCMFile->Write (args_info.OutputFile_arg, type);
-if (ecriture) std::cout<<"\n DICOM File written !"<<std::endl;
-*/
 //---------------------------------------------------------------------------------------
 //WRITE DICOM BIS
 // The previous way of writting DICOM-RT-DOSE works only for ITK
@@ -449,10 +335,10 @@ if (ecriture) std::cout<<"\n DICOM File written !"<<std::endl;
    else std::cout<<"\n DICOM File re-written, using the FileHelper syntax, in order to be processed by commercial systems !"<<std::endl;
 
 delete fh; */
-  gdcm::Writer w;
+/*  gdcm::Writer w;
   w.SetFile(mDCMFile);
   w.SetFileName(m_ArgsInfo.output_arg);
-  w.Write();
+  w.Write();*/
 //   fh->Delete();
 
 //---------------------------------------------------------------------------------------
@@ -547,8 +433,8 @@ mDCMFile->Write (args_info.OutputFile_arg, type);
 //((gdcm::SQItem*)mDCMFile)->RemoveEntry(a);
 
 //-----------------------------------------------------------------------------------
-std::cout <<"\n## DICOM Image to RT DOSE application is ended..."<<std::endl;
-std::cout <<"#########################################################\n" << std::endl;
+  std::cout <<"\n## DICOM Image to RT DOSE application is ended..."<<std::endl;
+  std::cout <<"#########################################################\n" << std::endl;
 
 #else
   std::cout << "Use GDCM2" << std::endl;
@@ -582,6 +468,17 @@ void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary
     }
 }
 //---------------------------------------------------------------------------
+
+//---------------------------------------------------------------------------
+template <typename T> std::string NumberToString ( T Number )
+{
+   std::ostringstream ss;
+   ss << Number;
+   return ss.str();
+}
+//---------------------------------------------------------------------------
+
+
 }//end clitk
 
 #endif //#define clitkImage2DicomDoseGenericFilter_txx