]> Creatis software - clitk.git/blobdiff - tools/clitkPartitionEnergyWindowDicomGenericFilter.txx
Debug RTStruct conversion with empty struc
[clitk.git] / tools / clitkPartitionEnergyWindowDicomGenericFilter.txx
index a32dc96d81a1076b8a0ab5cf42537216cf7fe41a..8e8fa33f37f23359d785f9d6531b8e4aab316382 100644 (file)
@@ -36,6 +36,8 @@
 #include <gdcmAttribute.h>
 #include <gdcmReader.h>
 #include <gdcmWriter.h>
+#include <gdcmDataElement.h>
+#include <gdcmTag.h>
 #else
 #include "gdcmFile.h"
 #include "gdcmUtil.h"
@@ -134,7 +136,7 @@ void
 PartitionEnergyWindowDicomGenericFilter<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;
@@ -148,7 +150,6 @@ PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelTy
   // Read Dicom model file
   typename ReaderType::Pointer reader = ReaderType::New();
   ImageIOType::Pointer gdcmIO = ImageIOType::New();
-  std::string filename_out = m_ArgsInfo.output_arg;
   gdcmIO->LoadPrivateTagsOn();
   gdcmIO->KeepOriginalUIDOn();
   reader->SetImageIO( gdcmIO );
@@ -167,22 +168,33 @@ PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelTy
   typename ReaderType::DictionaryArrayType outputArray;
   typename ReaderType::DictionaryRawPointer outputDict = new typename ReaderType::DictionaryType;
   CopyDictionary (*inputDict, *outputDict);
-  //itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(energyNumber));
-  //itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0021", NumberToString(headNumber));
+  itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(1));
   outputArray.push_back(outputDict);
 
 
-  //Read the number of slices and energy window
+  //Read the number of slices and energy window, rotation
   int nbSlice = input->GetLargestPossibleRegion().GetSize()[2];
   gdcm::Reader hreader;
   hreader.SetFileName(m_ArgsInfo.input_arg);
   hreader.Read();
-  gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
+  gdcm::DataSet& dsInput = hreader.GetFile().GetDataSet();
   gdcm::Attribute<0x54,0x11> series_number_att;
-  series_number_att.SetFromDataSet(ds);
+  series_number_att.SetFromDataSet(dsInput);
   int nbEnergyWindow = series_number_att.GetValue();
 
-  std::cout << nbSlice << " " << nbEnergyWindow << std::endl;
+  int nbRotation;
+  if (dsInput.FindDataElement(gdcm::Tag(0x54, 0x52))) {
+    const gdcm::DataElement &seqRotation = dsInput.GetDataElement(gdcm::Tag(0x54, 0x52));
+    gdcm::SmartPointer<gdcm::SequenceOfItems> sqiRotation = seqRotation.GetValueAsSQ();
+    gdcm::Item &itemRotation = sqiRotation->GetItem(1);
+    const gdcm::DataElement &deRotation = itemRotation.GetDataElement(gdcm::Tag(0x54, 0x53));
+    unsigned short SamplesPerPixelRotation = *((unsigned short *) deRotation.GetByteValue()->GetPointer());
+    std::stringstream s_SamplesPerPixelRotation;
+    s_SamplesPerPixelRotation << SamplesPerPixelRotation;
+    std::string p_StrRotation = s_SamplesPerPixelRotation.str();
+    nbRotation = (int)atof(p_StrRotation.c_str());
+  } else
+    nbRotation = 1;
 
   int nbSlicePerEnergy = nbSlice / nbEnergyWindow;
   if (nbSlicePerEnergy*nbEnergyWindow != nbSlice)
@@ -234,6 +246,8 @@ PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelTy
     writer->SetInput( energyImage );
     writer->SetImageIO( gdcmIO );
     typename ReaderType::FileNamesContainer fileNamesOutput;
+    std::string extension = NumberToString(energy) + ".dcm";
+    std::string filename_out = m_ArgsInfo.output_arg + extension;
     fileNamesOutput.push_back(filename_out);
     writer->SetFileNames( fileNamesOutput );
     writer->SetMetaDataDictionaryArray(&outputArray);
@@ -247,217 +261,106 @@ PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelTy
       std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
       std::cerr << excp << std::endl;
     }
-  }
-
-/*
-
-  //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);
-  }
-  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]);
-  }
-
-//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;
-      }
-    }
-  }
-
-
-  // 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);
 
-
-  // Output directory and filenames
-  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);
-
-
-  // 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;
-  }
-
-
-  //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();
+    //Write sequence dicom tag with gdcm
+    gdcm::Reader reader;
+    reader.SetFileName( filename_out.c_str() );
+    reader.Read();
+    gdcm::File &fileOutput = reader.GetFile();
+    gdcm::DataSet &dsOutput = fileOutput.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 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]));
+    gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
+    rotationDE.SetVR( gdcm::VR::US );
+    char essai = (char)nbRotation;
+    char *p_essai = &essai;
+    rotationDE.SetByteValue(p_essai, 1);
     // Create an item
-    gdcm::Item energyThresholdIt;
-    energyThresholdIt.SetVLToUndefined();
-    gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
-    energyThresholdDS.Insert(energyThresholdDE);
-    energyThresholdDS.Insert(energyUpholdDE);
-    energyThresholdSq->AddItem(energyThresholdIt);
+    gdcm::Item rotationIt;
+    rotationIt.SetVLToUndefined();
+    gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
+    rotationDS.Insert(rotationDE);
+    rotationSq->AddItem(rotationIt);
     // 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);
+    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();
+
+    //Read the caracteristics of this energy window : the name and the thresholds up and down
+    const gdcm::DataElement &seqWindow = dsInput.GetDataElement(gdcm::Tag(0x54, 0x12));
+    gdcm::SmartPointer<gdcm::SequenceOfItems> sqiWindow = seqWindow.GetValueAsSQ();
+    gdcm::Item &itemWindow = sqiWindow->GetItem(energy+1);
+    const gdcm::DataElement &deNameWindow = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x18));
+    std::string nameWindow( deNameWindow.GetByteValue()->GetPointer(), deNameWindow.GetByteValue()->GetLength() );
+
+    const gdcm::DataElement &deThresholds = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x13));
+    gdcm::SmartPointer<gdcm::SequenceOfItems> sqiThresholds = deThresholds.GetValueAsSQ();
+    for (unsigned int nbWindow=1; nbWindow <= sqiThresholds->GetNumberOfItems() ; ++nbWindow) {
+      gdcm::Item &itemThresholds = sqiThresholds->GetItem(nbWindow);
+      const gdcm::DataElement &deThresholdDown = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x14));
+      std::string p_StrThresholdDown( deThresholdDown.GetByteValue()->GetPointer(), deThresholdDown.GetByteValue()->GetLength() );
+      double thresholdDown = (double)atof(p_StrThresholdDown.c_str());
+      const gdcm::DataElement &deThresholdUp = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x15));
+      std::string p_StrThresholdUp( deThresholdUp.GetByteValue()->GetPointer(), deThresholdUp.GetByteValue()->GetLength() );
+      double thresholdUp = (double)atof(p_StrThresholdUp.c_str());
+
+      //Now write it !
+      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(NumberToString(thresholdDown).c_str(), (uint32_t)strlen(NumberToString(thresholdDown).c_str()));
+      energyUpholdDE.SetByteValue(NumberToString(thresholdUp).c_str(), (uint32_t)strlen(NumberToString(thresholdUp).c_str()));
+      // 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(nameWindow.c_str(), (uint32_t)strlen(nameWindow.c_str()));
+      // 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();
+    dsOutput.Insert(energyDEParent);
+    dsOutput.Insert(rotationDEParent);
+    gdcm::Writer w;
+    w.SetFile( fileOutput );
+    w.SetFileName( filename_out.c_str() );
+    w.Write();
   }
-  // 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