X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkGateSimulation2DicomGenericFilter.txx;h=9ad2857f90bd8bb263463b47b932f33cabaaed93;hb=b8e5890d37dfd64409b9694f73c0be164a089e64;hp=ce1c1e3f4fcb726b64f5897938834318244fe94c;hpb=f8ad682c6c916a31d082591d744007411639794c;p=clitk.git

diff --git a/tools/clitkGateSimulation2DicomGenericFilter.txx b/tools/clitkGateSimulation2DicomGenericFilter.txx
index ce1c1e3..9ad2857 100644
--- a/tools/clitkGateSimulation2DicomGenericFilter.txx
+++ b/tools/clitkGateSimulation2DicomGenericFilter.txx
@@ -27,15 +27,24 @@
  *
  ===================================================*/
 
+#include <sstream>
 // clitk
 #include "clitkResampleImageWithOptionsFilter.h"
 #if GDCM_MAJOR_VERSION >= 2
 #include "gdcmUIDGenerator.h"
+#include <gdcmImageHelper.h>
+#include <gdcmAttribute.h>
+#include <gdcmReader.h>
+#include <gdcmWriter.h>
 #else
 #include "gdcmFile.h"
 #include "gdcmUtil.h"
 #endif
 
+#include "itkImageRegionIterator.h"
+#include "itkMetaImageIO.h"
+#include "itkMetaDataDictionary.h"
+
 
 namespace clitk
 {
@@ -125,49 +134,33 @@ void
 GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 {
 
+#if GDCM_MAJOR_VERSION >= 2
   // ImageTypes
   typedef itk::Image<PixelType, Dimension> InputImageType;
   typedef itk::Image<PixelType, Dimension> OutputImageType;
-
-  // Read the dicom directory
-  typedef itk::ImageSeriesReader< InputImageType >     ReaderType;
+  typedef itk::ImageFileReader< InputImageType >     ReaderType;
+  typedef itk::ImageSeriesReader< InputImageType >     ReaderSeriesType;
   typedef itk::GDCMImageIO ImageIOType;
-  typedef itk::GDCMSeriesFileNames NamesGeneratorType;
 
-  ImageIOType::Pointer gdcmIO = ImageIOType::New();
-  NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
-  namesGenerator->SetInputDirectory( m_ArgsInfo.inputDir_arg );
-  namesGenerator->SetOutputDirectory( m_ArgsInfo.outputDir_arg  );
-  typename   ReaderType::FileNamesContainer filenames_in = namesGenerator->GetInputFileNames();
-  typename   ReaderType::FileNamesContainer filenames_out = namesGenerator->GetOutputFileNames();
-
-  // Output the dicom files
-  unsigned int numberOfFilenames =  filenames_in.size();
-  if (m_Verbose) {
-    std::cout << numberOfFilenames <<" were read in the directory "<<m_ArgsInfo.inputDir_arg<<"..."<<std::endl<<std::endl;
-    for(unsigned int fni = 0; fni<numberOfFilenames; fni++) {
-      std::cout << "filename # " << fni << " = ";
-      std::cout << filenames_in[fni] << std::endl;
-    }
-  }
 
-  // Read the series
-  typename ReaderType::Pointer reader = ReaderType::New();
-  if (m_ArgsInfo.preserve_flag) {
-#if GDCM_MAJOR_VERSION >= 2
-    gdcmIO->LoadPrivateTagsOn();
-    gdcmIO->KeepOriginalUIDOn();
-#endif
-  }
-  reader->SetImageIO( gdcmIO );
-  reader->SetFileNames( filenames_in );
+  // Read Dicom model file
+  typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New();
+  ImageIOType::Pointer gdcmIO = ImageIOType::New();
+  std::string filename_out = m_ArgsInfo.outputFilename_arg;
+  gdcmIO->LoadPrivateTagsOn();
+  gdcmIO->KeepOriginalUIDOn();
+  typename ReaderSeriesType::FileNamesContainer fileNames;
+  fileNames.push_back(m_ArgsInfo.inputModelFilename_arg);
+  readerSeries->SetImageIO( gdcmIO );
+  readerSeries->SetFileNames( fileNames );
   try {
-    reader->Update();
+    readerSeries->Update();
   } catch (itk::ExceptionObject &excp) {
-    std::cerr << "Error: Exception thrown while reading the DICOM series!!" << std::endl;
+    std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
     std::cerr << excp << std::endl;
   }
 
+
   // Read the input (MHD file)
   typedef typename InputImageType::RegionType RegionType;
   typedef typename RegionType::SizeType SizeType;
@@ -175,255 +168,258 @@ GateSimulation2DicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
   volumeReader->SetFileName( m_InputFileName);
   volumeReader->Update();
-  
   typename InputImageType::Pointer input = volumeReader->GetOutput();
-  if ((!m_ArgsInfo.useSizeAsReference_flag && (input->GetSpacing() != reader->GetOutput()->GetSpacing())) || 
-      (m_ArgsInfo.useSizeAsReference_flag && (input->GetLargestPossibleRegion().GetSize() != reader->GetOutput()->GetLargestPossibleRegion().GetSize()))) {
-        
-    // resampling is carried out on the fly if resolution or size between 
-    // the input mhd and input dicom series is different
-    
-    // Filter
-    typedef clitk::ResampleImageWithOptionsFilter<InputImageType, InputImageType> ResampleImageFilterType;
-    typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
-    filter->SetInput(input);
-    filter->SetVerboseOptions(m_Verbose);
-    filter->SetGaussianFilteringEnabled(false);
-    filter->SetDefaultPixelValue(0);
-    
-    if (!m_ArgsInfo.useSizeAsReference_flag) {
-      filter->SetOutputSpacing(input->GetSpacing());
-      if (m_Verbose) {
-        std::cout << "Warning: The image spacing differs between the MHD file and the input dicom series. Performing resampling with default options using spacing as reference (for advanced options, use clitkResampleImage)." << std::endl;
-        std::cout << "MHD -> " << input->GetSpacing() << std::endl;
-        std::cout << "dicom -> " << reader->GetOutput()->GetSpacing() << std::endl;
-      }
-    }
-    else {
-      const SizeType& dicom_size = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
-      SizeType output_size;
-      for (unsigned int i = 0; i < Dimension; i++)
-        output_size[i] = dicom_size[i];
-      filter->SetOutputSize(output_size);
-      if (m_Verbose) {
-          std::cout << "Warning: The image size differs between the MHD file and the input dicom series. Performing resampling with default options using size as reference (for advanced options, use clitkResampleImage)." << std::endl;
-          std::cout << "MHD -> " << input->GetLargestPossibleRegion().GetSize() << std::endl;
-          std::cout << "dicom -> " << reader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
-      }
-    }
-
-    try {
-      filter->Update();
-      input = filter->GetOutput();
-    } catch( itk::ExceptionObject & excp ) {
-    std::cerr << "Error: Exception thrown while resampling!!" << std::endl;
-    std::cerr << excp << std::endl;
-    }
-  }
-    std::cout << "RESAMPING FILTER DONE" << std::endl;
-  
-  //	In some cases, due to resampling approximation issues, 
-  //	the number of slices in the MHD file may be different (smaller)
-  //	from the number of files in the template dicom directory. 
-  //	To avoid ITK generating an exception, we reduce the number 
-  //	of DCM files to be considered, and a warning is printed
-  //	in verbose mode
-  const RegionType volumeRegion = input->GetLargestPossibleRegion();
-  const SizeType& volumeSize = volumeRegion.GetSize();
-  if (m_Verbose) {
-    std::cout << volumeRegion << volumeSize << std::endl;
-  }
-    std::cout << Dimension << volumeSize[2] << numberOfFilenames << std::endl;
-  if (Dimension == 3 && volumeSize[2] < numberOfFilenames) {
-    if (m_Verbose)
-      std::cout << "Warning: The number of files in " << m_ArgsInfo.inputDir_arg << " (" << filenames_in.size() << " files) is greater than the number of slices in MHD (" << volumeSize[2] << " slices). Using only " << volumeSize[2] << " files." << std::endl;
-    
-    filenames_in.resize(volumeSize[2]);
-    filenames_out.resize(filenames_in.size());
-    numberOfFilenames =  filenames_in.size();
-  }
 
-  // Modify the meta dictionary
-  typedef itk::MetaDataDictionary   DictionaryType;
-  const std::vector<DictionaryType*>* dictionary = reader->GetMetaDataDictionaryArray();
-
-  // Get keys
-  unsigned int numberOfKeysGiven=m_ArgsInfo.key_given;
-    if (m_ArgsInfo.verbose_flag) 
-      DD(numberOfKeysGiven);
-
-  std::string seriesUID;
-  std::string frameOfReferenceUID;
-  std::string studyUID;
-  
-  // one pass through the keys given on the cmd-line, to check what will be recreated
-  std::string seriesUIDkey = "0020|000e";
-  std::string seriesNumberKey = "0020|0011";
-  std::string seriesDescriptionKey = "0008|103e";
-  std::string frameOfReferenceUIDKey = "0020|0052";
-  std::string studyUIDKey = "0020|000d";
-  std::string studyIDKey = "0020|0010";
-  std::string studyDescriptionKey = "0008|1030";
-  bool seriesUIDGiven = false;
-  bool seriesNumberGiven = false;
-  bool seriesDescriptionGiven = false;
-  bool studyUIDGiven = false;
-  bool studyIDGiven = false;
-  bool studyDescriptionGiven = false;
-  for (unsigned int i = 0; i < numberOfKeysGiven; i++) {
-    std::string entryId( m_ArgsInfo.key_arg[i] );
-    if (m_ArgsInfo.verbose_flag) 
-      DD(entryId);
-    
-    seriesUIDGiven |= (entryId ==  seriesUIDkey || entryId ==  frameOfReferenceUIDKey);
-    seriesNumberGiven |= (entryId == seriesNumberKey);
-    seriesDescriptionGiven |= (entryId == seriesDescriptionKey);
-    studyUIDGiven |= (entryId == studyUIDKey);
-    studyIDGiven |= (entryId == studyIDKey);
-    studyDescriptionGiven |= (entryId == studyDescriptionKey);
-  }
 
-  // force the creation of a new series if a new study was specified
-  if (!studyUIDGiven && m_ArgsInfo.newStudyUID_flag) {
-    m_ArgsInfo.newSeriesUID_flag = true;
-#if GDCM_MAJOR_VERSION >= 2
-    gdcm::UIDGenerator suid;
-    studyUID = suid.Generate();
-#else
-    studyUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
-#endif
+  //Read the metadata informations in the mhd
+  MetaObject tObj(3);
+  tObj.AddUserField("NumberOfEnergyWindows", MET_STRING);
+  tObj.AddUserField("NbOfHeads", MET_STRING);
+  tObj.AddUserField("NbOfProjections", MET_STRING);
+  tObj.Read(m_InputFileName.c_str());
+  char *p_energyNumber, *p_headNumber, *p_rotationNumber;
+  int energyNumber, headNumber, rotationNumber;
+  p_energyNumber = static_cast<char*>(tObj.GetUserField("NumberOfEnergyWindows"));
+  p_headNumber = static_cast<char*>(tObj.GetUserField("NbOfHeads"));
+  p_rotationNumber = static_cast<char*>(tObj.GetUserField("NbOfProjections"));
+  energyNumber = atoi(p_energyNumber);
+  headNumber = atoi(p_headNumber);
+  rotationNumber = atoi(p_rotationNumber);
+  for (unsigned int i=0; i<energyNumber; ++i) {
+    std::string tagEnergyName("EnergyWindow_");
+    tagEnergyName += NumberToString(i).c_str();
+    std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
+    tagEnergyThreshold += "_threshold";
+    tagEnergyUphold += "_uphold";
+    tObj.AddUserField(tagEnergyName.c_str(), MET_STRING);
+    tObj.AddUserField(tagEnergyThreshold.c_str(), MET_STRING);
+    tObj.AddUserField(tagEnergyUphold.c_str(), MET_STRING);
   }
-    
-  if (!seriesUIDGiven && m_ArgsInfo.newSeriesUID_flag) {
-#if GDCM_MAJOR_VERSION >= 2
-    gdcm::UIDGenerator suid;
-    seriesUID = suid.Generate();
-    gdcm::UIDGenerator fuid;
-    frameOfReferenceUID = fuid.Generate();
-#else
-    seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
-    frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
-#endif
-  }
-
-  if (m_ArgsInfo.verbose_flag) {
-    DD(seriesUID);
-    DD(frameOfReferenceUID);
-    DD(studyUID);
+  tObj.Read(m_InputFileName.c_str());
+  std::vector<char*> p_EnergyWindowName(energyNumber), p_EnergyWindowThreshold(energyNumber), p_EnergyWindowUphold(energyNumber);
+  std::vector<int> energyWindowThreshold(energyNumber), energyWindowUphold(energyNumber);
+  for (unsigned int i=0; i<energyNumber; ++i) {
+    std::string tagEnergyName("EnergyWindow_");
+    tagEnergyName += NumberToString(i).c_str();
+    std::string tagEnergyThreshold(tagEnergyName), tagEnergyUphold(tagEnergyName);
+    tagEnergyThreshold += "_threshold";
+    tagEnergyUphold += "_uphold";
+    p_EnergyWindowName[i] = static_cast<char*>(tObj.GetUserField(tagEnergyName.c_str()));
+    p_EnergyWindowThreshold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyThreshold.c_str()));
+    p_EnergyWindowUphold[i] = static_cast<char*>(tObj.GetUserField(tagEnergyUphold.c_str()));
+    energyWindowThreshold[i] = atoi(p_EnergyWindowThreshold[i]);
+    energyWindowUphold[i] = atoi(p_EnergyWindowUphold[i]);
   }
 
-  // check if file UIDs will be be preserved
-  bool useInputFileUID = true;
-  if (m_ArgsInfo.newSeriesUID_flag || m_ArgsInfo.newStudyUID_flag || seriesUIDGiven || studyUIDGiven) {
-    useInputFileUID = false;
-  }
-  else {
-#if GDCM_MAJOR_VERSION < 2
-    gdcmIO->SetKeepOriginalUID(true);
-#endif
-    namesGenerator->SetOutputDirectory( m_ArgsInfo.outputDir_arg  );
+//TODO if the input size/spacing and dicom model ones are different
+
+  // Create a new mhd image with the correct dicom order slices
+  typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
+  mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
+  mhdCorrectOrder->Allocate();
+  unsigned int zAxis(0); //z value for the input mhd image
+  for (unsigned int energy = 0; energy < energyNumber; ++energy) {
+    for (unsigned int head = 0; head < headNumber; ++head) {
+      for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
+        std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
+
+        typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
+        startIteratorIndexCorrectOrder[0] = 0;
+        startIteratorIndexCorrectOrder[1] = 0;
+        startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
+
+        typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
+        startIteratorIndexOriginalOrder[0] = 0;
+        startIteratorIndexOriginalOrder[1] = 0;
+        startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
+
+        typename InputImageType::SizeType regionSizeIterator;
+        regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
+        regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
+        regionSizeIterator[2] = 1;
+
+        typename InputImageType::RegionType regionIteratorCorrectOrder;
+        regionIteratorCorrectOrder.SetSize(regionSizeIterator);
+        regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
+
+        typename InputImageType::RegionType regionIteratorOriginalOrder;
+        regionIteratorOriginalOrder.SetSize(regionSizeIterator);
+        regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
+
+        itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
+        itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
+        while(!CorrectOrderIterator.IsAtEnd()) {
+          CorrectOrderIterator.Set(OriginalOrderIterator.Get());
+          ++CorrectOrderIterator;
+          ++OriginalOrderIterator;
+        }
+
+        ++zAxis;
+      }
+    }
   }
-  
-  filenames_out.resize(numberOfFilenames);
-  
-  time_t t;
-  t = time(&t);
-  struct tm* instanceDateTimeTm = localtime(&t);
-  char datetime[16];
-  strftime(datetime, 16, "%Y%m%d", instanceDateTimeTm);
-  std::ostringstream instanceDate;
-  instanceDate << datetime;
-  std::ostringstream instanceTime;
-  strftime(datetime, 16, "%H%M%S", instanceDateTimeTm);
-  instanceTime << datetime;
-  
-  // update output dicom keys/tags
-  for(unsigned int fni = 0; fni<numberOfFilenames; fni++) {
-    for (unsigned int i = 0; i < numberOfKeysGiven; i++) {
-      std::string entryId(m_ArgsInfo.key_arg[i]  );
-      std::string value( m_ArgsInfo.tag_arg[i] );
 
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), entryId, value );
-    }
 
-    // series UID
-    if (!seriesUIDGiven) {
-      if (m_ArgsInfo.newSeriesUID_flag) {
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), seriesUIDkey, seriesUID );
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), frameOfReferenceUIDKey, frameOfReferenceUID );
-      }
-    }
-    
-    // study UID
-    if (!studyUIDGiven) {
-      if (m_ArgsInfo.newStudyUID_flag) 
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), studyUIDKey, studyUID );
-    }
-    
-    // study description
-    if (studyUIDGiven || m_ArgsInfo.newStudyUID_flag) {
-      if (!studyIDGiven)
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), studyIDKey,itksys::SystemTools::GetFilenameName( m_ArgsInfo.outputDir_arg ));
-      if (!studyDescriptionGiven)
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), studyDescriptionKey,itksys::SystemTools::GetFilenameName( m_ArgsInfo.outputDir_arg ));
-      
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0008|0020", instanceDate.str() );
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0008|0030", instanceTime.str() );
-    }
-    
-    // series description/number
-    if (seriesUIDGiven || m_ArgsInfo.newSeriesUID_flag) {
-      if (!seriesDescriptionGiven)
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), seriesDescriptionKey, itksys::SystemTools::GetFilenameName(m_ArgsInfo.outputDir_arg) );
-      if (!seriesNumberGiven)
-        itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), seriesNumberKey, itksys::SystemTools::GetFilenameName(m_ArgsInfo.outputDir_arg) );
-
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0008|0012", instanceDate.str() );
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0008|0013", instanceTime.str() );
-    }
+  // 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, "0054|0011", NumberToString(energyNumber));
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
+  outputArray.push_back(outputDict);
 
-    // file UIDs are recreated for new studies or series
-    if (!useInputFileUID)
-    {
-      if (m_ArgsInfo.verbose_flag)
-        std::cout << "Recreating file UIDs" << std::endl;
 
-      std::string fileUID;
-#if GDCM_MAJOR_VERSION >= 2
-      gdcm::UIDGenerator fid;
-      fileUID = fid.Generate();
-#else
-      fileUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
-#endif
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0008|0018", fileUID );
-      itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), "0002|0003", fileUID );
-      
-      filenames_out[fni] = itksys::SystemTools::CollapseFullPath(fileUID.c_str(), m_ArgsInfo.outputDir_arg) + std::string(".dcm"); 
-    }
-  }
-  
   // Output directory and filenames
-  itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputDir_arg ); // create if it doesn't exist
-  typedef itk::ImageFileWriter<OutputImageType>  SeriesWriterType;
-  typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
+  typedef itk::ImageSeriesWriter<OutputImageType, OutputImageType>  WriterSerieType;
+  typename WriterSerieType::Pointer writerSerie = WriterSerieType::New();
+  writerSerie->SetInput( mhdCorrectOrder );
+  writerSerie->SetImageIO( gdcmIO );
+  typename ReaderSeriesType::FileNamesContainer fileNamesOutput;
+  fileNamesOutput.push_back(filename_out);
+  writerSerie->SetFileNames( fileNamesOutput );
+  writerSerie->SetMetaDataDictionaryArray(&outputArray);
 
-  seriesWriter->SetInput( input );
-  seriesWriter->SetImageIO( gdcmIO );
-  
-  seriesWriter->SetFileName( filenames_out[0] );
-  seriesWriter->SetMetaDataDictionary( *((*dictionary)[0]) );
 
   // Write
   try {
     if (m_ArgsInfo.verbose_flag)
-      std::cout << seriesWriter << std::endl;
-    seriesWriter->Update();
+      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;
   }
 
+
+  //Write sequence dicom tag with gdcm
+  gdcm::Reader reader;
+  reader.SetFileName( fileNamesOutput[0].c_str() );
+  reader.Read();
+  gdcm::File &file = reader.GetFile();
+  gdcm::DataSet &ds2 = file.GetDataSet();
+  const unsigned int ptr_len = 42;
+  char *ptr = new char[ptr_len];
+  memset(ptr,0,ptr_len);
+
+  //Write rotation tag
+  // Create a Sequence
+  gdcm::SmartPointer<gdcm::SequenceOfItems> rotationSq = new gdcm::SequenceOfItems();
+  rotationSq->SetLengthToUndefined();
+  // Create a dataelement
+  gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
+  rotationDE.SetVR( gdcm::VR::US );
+  char essai = (char)rotationNumber;
+  char *p_essai = &essai;
+  rotationDE.SetByteValue(p_essai, 1);
+  // Create an item
+  gdcm::Item rotationIt;
+  rotationIt.SetVLToUndefined();
+  gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
+  rotationDS.Insert(rotationDE);
+  rotationSq->AddItem(rotationIt);
+  // Insert sequence into data set
+  gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) );
+  rotationDEParent.SetVR(gdcm::VR::SQ);
+  rotationDEParent.SetValue(*rotationSq);
+  rotationDEParent.SetVLToUndefined();
+
+  //Write energy
+  gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) );
+  energyDEParent.SetVR(gdcm::VR::SQ);
+  // Create a Sequence
+  gdcm::SmartPointer<gdcm::SequenceOfItems> energySq = new gdcm::SequenceOfItems();
+  energySq->SetLengthToUndefined();
+  for (unsigned int i=0; i<energyNumber; ++i) {
+    gdcm::SmartPointer<gdcm::SequenceOfItems> energyThresholdSq = new gdcm::SequenceOfItems();
+    energyThresholdSq->SetLengthToUndefined();
+    // Create a dataelement
+    gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) );
+    gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) );
+    energyThresholdDE.SetVR( gdcm::VR::DS );
+    energyUpholdDE.SetVR( gdcm::VR::DS );
+    energyThresholdDE.SetByteValue(p_EnergyWindowThreshold[i], (uint32_t)strlen(p_EnergyWindowThreshold[i]));
+    energyUpholdDE.SetByteValue(p_EnergyWindowUphold[i], (uint32_t)strlen(p_EnergyWindowUphold[i]));
+    // Create an item
+    gdcm::Item energyThresholdIt;
+    energyThresholdIt.SetVLToUndefined();
+    gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
+    energyThresholdDS.Insert(energyThresholdDE);
+    energyThresholdDS.Insert(energyUpholdDE);
+    energyThresholdSq->AddItem(energyThresholdIt);
+    // Insert sequence into data set
+    gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) );
+    energyThresholdDEParent.SetVR(gdcm::VR::SQ);
+    energyThresholdDEParent.SetValue(*energyThresholdSq);
+    energyThresholdDEParent.SetVLToUndefined();
+    // Create a dataelement
+    gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) );
+    energyNameDE.SetVR( gdcm::VR::SH );
+    energyNameDE.SetByteValue(p_EnergyWindowName[i], (uint32_t)strlen(p_EnergyWindowName[i]));
+    // Create an item
+    gdcm::Item energyIt;
+    energyIt.SetVLToUndefined();
+    gdcm::DataSet &energyDS = energyIt.GetNestedDataSet();
+    energyDS.Insert(energyNameDE);
+    energyDS.Insert(energyThresholdDEParent);
+    energySq->AddItem(energyIt);
+  }
+  // Insert sequence into data set
+  energyDEParent.SetValue(*energySq);
+  energyDEParent.SetVLToUndefined();
+  ds2.Insert(energyDEParent);
+  ds2.Insert(rotationDEParent);
+
+  gdcm::Writer w;
+  w.SetFile( file );
+  w.SetFileName( fileNamesOutput[0].c_str() );
+  w.Write();
+
+#else
+  std::cout << "Use GDCM2" << std::endl;
+#endif
+}
+//---------------------------------------------------------------------------
+
+
+//---------------------------------------------------------------------------
+void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
+{
+  typedef itk::MetaDataDictionary DictionaryType;
+
+  DictionaryType::ConstIterator itr = fromDict.Begin();
+  DictionaryType::ConstIterator end = fromDict.End();
+  typedef itk::MetaDataObject< std::string > MetaDataStringType;
+
+  while( itr != end )
+    {
+    itk::MetaDataObjectBase::Pointer  entry = itr->second;
+
+    MetaDataStringType::Pointer entryvalue =
+      dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
+    if( entryvalue )
+      {
+      std::string tagkey   = itr->first;
+      std::string tagvalue = entryvalue->GetMetaDataObjectValue();
+      itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
+      }
+    ++itr;
+    }
+}
+//---------------------------------------------------------------------------
+
+
+//---------------------------------------------------------------------------
+template <typename T> std::string NumberToString ( T Number )
+{
+   std::ostringstream ss;
+   ss << Number;
+   return ss.str();
 }
+//---------------------------------------------------------------------------
 
 }//end clitk