X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkImage2DicomDoseGenericFilter.txx;h=df105fadfd1a75e1838d057dbd39ee6b9fdb0edb;hb=5578995d9a82792833333eeb3dd5c8ecac967293;hp=667340a07e486c02c3f7961a820b76e0e4c5b8ea;hpb=103778ec7e0d03e32d9ecfa1306312d32d72ad4c;p=clitk.git

diff --git a/tools/clitkImage2DicomDoseGenericFilter.txx b/tools/clitkImage2DicomDoseGenericFilter.txx
index 667340a..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,13 +138,16 @@ 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;
 
   //-----------------------------------------------------------------------------
@@ -151,115 +157,113 @@ Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
   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();
+  std::string filename_out = m_ArgsInfo.output_arg;
+  gdcmIO->LoadPrivateTagsOn();
+  gdcmIO->KeepOriginalUIDOn();
+  typename ReaderSeriesType::FileNamesContainer fileNames;
+  fileNames.push_back(m_ArgsInfo.DicomInputFile_arg);
+  readerSeries->SetImageIO( gdcmIO );
+  readerSeries->SetFileNames( fileNames );
+  try {
+    readerSeries->Update();
+  } catch (itk::ExceptionObject &excp) {
+    std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
+    std::cerr << excp << std::endl;
+  }
+
+  // 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);
 
   // 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.;
-  std::stringstream strOffset;
-  strOffset << offset;
+  value.str("");
+  value << offset;
   for (int i=1; i<NbFrames ; i++){
     offset+=Spacing[2];
-    strOffset << "\\";
-    strOffset << offset;
+    value << '\\';
+    value << offset;
   }
-  DD(strOffset.str().c_str());
+  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000c", value.str());
+  DD(value.str());
 
   // 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();
-  }
+  typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
+  imageCalculatorFilter->SetImage(image);
+  imageCalculatorFilter->ComputeMaximum();
+  float highestValue=imageCalculatorFilter->GetMaximum();
   double doseScaling = highestValue/(pow(2,16)-1);
-  DD(doseScaling);
+  value.str("");
+  value<<doseScaling;
+  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000e", value.str());
+  DD(value.str());
 
   // 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());*/
-
-  // Read Dicom model file
-  typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
-  ImageIOType::Pointer gdcmIO = ImageIOType::New();
-  std::string filename_out = m_ArgsInfo.output_arg;
-  gdcmIO->LoadPrivateTagsOn();
-  gdcmIO->KeepOriginalUIDOn();
-  typename ReaderSeriesType::FileNamesContainer fileNames;
-  fileNames.push_back(m_ArgsInfo.DicomInputFile_arg);
-  readerSeries->SetImageIO( gdcmIO );
-  readerSeries->SetFileNames( fileNames );
-  try {
-    readerSeries->Update();
-  } catch (itk::ExceptionObject &excp) {
-    std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
-    std::cerr << excp << std::endl;
+  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;
   }
 
-  // update output dicom keys/tags
-  // string for distinguishing items inside sequence:
-  const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE");
-  std::string tempString = ITEM_ENCAPSULATE_STRING + "01";
-  typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0];
-  typename ReaderSeriesType::DictionaryArrayType outputArray;
-  typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType;
-  CopyDictionary (*inputDict, *outputDict);
-  itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0032", NumberToString(origin));
-  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));
-  itk::EncapsulateMetaData<std::string>(*outputDict, "0028|0030", NumberToString(Spacing));
-  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000e", NumberToString(doseScaling));
-  itk::EncapsulateMetaData<std::string>(*outputDict, "3004|000c", strOffset.str());
-  outputArray.push_back(outputDict);
-
   // 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);
@@ -276,6 +280,36 @@ Image2DicomDoseGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
     std::cerr << excp << 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 BIS