tools/clitkPartitionEnergyWindowDicomGenericFilter.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 clitkPartitionEnergyWindowDicomGenericFilter_txx
  19 #define clitkPartitionEnergyWindowDicomGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkPartitionEnergyWindowDicomGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include <sstream>
  31 // clitk
  32 #include "clitkResampleImageWithOptionsFilter.h"
  33 #if GDCM_MAJOR_VERSION >= 2
  34 #include "gdcmUIDGenerator.h"
  35 #include <gdcmImageHelper.h>
  36 #include <gdcmAttribute.h>
  37 #include <gdcmReader.h>
  38 #include <gdcmWriter.h>
  39 #include <gdcmDataElement.h>
  40 #include <gdcmTag.h>
  41 #else
  42 #include "gdcmFile.h"
  43 #include "gdcmUtil.h"
  44 #endif
  45
  46 #include "itkImageRegionIterator.h"
  47 #include "itkMetaImageIO.h"
  48 #include "itkMetaDataDictionary.h"
  49
  50
  51 namespace clitk
  52 {
  53
  54
  55 //-----------------------------------------------------------
  56 // Constructor
  57 //-----------------------------------------------------------
  58 template<class args_info_type>
  59 PartitionEnergyWindowDicomGenericFilter<args_info_type>::PartitionEnergyWindowDicomGenericFilter()
  60 {
  61   m_Verbose=false;
  62   m_InputFileName="";
  63 }
  64
  65
  66 //-----------------------------------------------------------
  67 // Update
  68 //-----------------------------------------------------------
  69 template<class args_info_type>
  70 void PartitionEnergyWindowDicomGenericFilter<args_info_type>::Update()
  71 {
  72   // Read the Dimension and PixelType
  73   int Dimension;
  74   std::string PixelType;
  75   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  76
  77
  78   // Call UpdateWithDim
  79   if(Dimension==2) UpdateWithDim<2>(PixelType);
  80   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  81   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  82   else {
  83     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  84     return;
  85   }
  86 }
  87
  88 //-------------------------------------------------------------------
  89 // Update with the number of dimensions
  90 //-------------------------------------------------------------------
  91 template<class args_info_type>
  92 template<unsigned int Dimension>
  93 void
  94 PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
  95 {
  96   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  97
  98   if(PixelType == "short") {
  99     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
 100     UpdateWithDimAndPixelType<Dimension, signed short>();
 101   }
 102   else if(PixelType == "unsigned_short"){
 103     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
 104     UpdateWithDimAndPixelType<Dimension, unsigned short>();
 105   }
 106
 107   else if (PixelType == "unsigned_char") {
 108     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
 109     UpdateWithDimAndPixelType<Dimension, unsigned char>();
 110   }
 111
 112   //     else if (PixelType == "char"){
 113   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
 114   //       UpdateWithDimAndPixelType<Dimension, signed char>();
 115   //     }
 116   else if (PixelType == "double") {
 117     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
 118     UpdateWithDimAndPixelType<Dimension, double>();
 119   }
 120   else {
 121     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 122     UpdateWithDimAndPixelType<Dimension, float>();
 123   }
 124 }
 125
 126 //-------------------------------------------------------------------
 127 // Update with the number of dimensions and the pixeltype read from
 128 // the dicom files. The MHD files may be resampled to match the
 129 // dicom spacing (and number of slices). Rounding errors in resampling
 130 // are handled by removing files when generating the output dicom
 131 // series.
 132 //-------------------------------------------------------------------
 133 template<class args_info_type>
 134 template <unsigned int Dimension, class  PixelType>
 135 void
 136 PartitionEnergyWindowDicomGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 137 {
 138
 139 #if GDCM_MAJOR_VERSION == 2
 140   // ImageTypes
 141   typedef itk::Image<PixelType, Dimension> InputImageType;
 142   typedef itk::Image<PixelType, Dimension> OutputImageType;
 143   typedef itk::ImageSeriesReader< InputImageType >     ReaderType;
 144   typedef itk::ImageSeriesWriter< InputImageType, InputImageType >     WriterType;
 145   typedef itk::GDCMImageIO ImageIOType;
 146   typedef typename InputImageType::RegionType RegionType;
 147   typedef typename RegionType::SizeType SizeType;
 148
 149
 150   // Read Dicom model file
 151   typename ReaderType::Pointer reader = ReaderType::New();
 152   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 153   gdcmIO->LoadPrivateTagsOn();
 154   gdcmIO->KeepOriginalUIDOn();
 155   reader->SetImageIO( gdcmIO );
 156   typename ReaderType::FileNamesContainer fileNamesInput;
 157   fileNamesInput.push_back(m_ArgsInfo.input_arg);
 158   reader->SetFileNames( fileNamesInput );
 159   try {
 160     reader->Update();
 161   } catch (itk::ExceptionObject &excp) {
 162     std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl;
 163     std::cerr << excp << std::endl;
 164   }
 165   typename InputImageType::Pointer input = reader->GetOutput();
 166
 167   typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
 168   typename ReaderType::DictionaryArrayType outputArray;
 169   typename ReaderType::DictionaryRawPointer outputDict = new typename ReaderType::DictionaryType;
 170   CopyDictionary (*inputDict, *outputDict);
 171   itk::EncapsulateMetaData<std::string>(*outputDict, "0054|0011", NumberToString(1));
 172   outputArray.push_back(outputDict);
 173
 174
 175   //Read the number of slices and energy window, rotation
 176   int nbSlice = input->GetLargestPossibleRegion().GetSize()[2];
 177   gdcm::Reader hreader;
 178   hreader.SetFileName(m_ArgsInfo.input_arg);
 179   hreader.Read();
 180   gdcm::DataSet& dsInput = hreader.GetFile().GetDataSet();
 181   gdcm::Attribute<0x54,0x11> series_number_att;
 182   series_number_att.SetFromDataSet(dsInput);
 183   int nbEnergyWindow = series_number_att.GetValue();
 184
 185   int nbRotation;
 186   if (dsInput.FindDataElement(gdcm::Tag(0x54, 0x52))) {
 187     const gdcm::DataElement &seqRotation = dsInput.GetDataElement(gdcm::Tag(0x54, 0x52));
 188     gdcm::SmartPointer<gdcm::SequenceOfItems> sqiRotation = seqRotation.GetValueAsSQ();
 189     gdcm::Item &itemRotation = sqiRotation->GetItem(1);
 190     const gdcm::DataElement &deRotation = itemRotation.GetDataElement(gdcm::Tag(0x54, 0x53));
 191     unsigned short SamplesPerPixelRotation = *((unsigned short *) deRotation.GetByteValue()->GetPointer());
 192     std::stringstream s_SamplesPerPixelRotation;
 193     s_SamplesPerPixelRotation << SamplesPerPixelRotation;
 194     std::string p_StrRotation = s_SamplesPerPixelRotation.str();
 195     nbRotation = (int)atof(p_StrRotation.c_str());
 196   } else
 197     nbRotation = 1;
 198
 199   int nbSlicePerEnergy = nbSlice / nbEnergyWindow;
 200   if (nbSlicePerEnergy*nbEnergyWindow != nbSlice)
 201     std::cerr << "Error: The number of slices is not correct !!" << std::endl;
 202
 203   for (unsigned int energy = 0; energy < nbEnergyWindow; ++energy) {
 204
 205     //Create the output
 206     typename InputImageType::Pointer energyImage = InputImageType::New();
 207     energyImage->SetRegions(input->GetLargestPossibleRegion());
 208     typename InputImageType::IndexType startIndex;
 209     startIndex[0] = 0;
 210     startIndex[1] = 0;
 211     startIndex[2] = 0;//energy*nbSlicePerEnergy;
 212     typename InputImageType::SizeType regionSize;
 213     regionSize[0] = input->GetLargestPossibleRegion().GetSize()[0];
 214     regionSize[1] = input->GetLargestPossibleRegion().GetSize()[1];
 215     regionSize[2] = nbSlicePerEnergy;
 216     typename InputImageType::RegionType region;
 217     region.SetSize(regionSize);
 218     region.SetIndex(startIndex);
 219     energyImage->SetRegions(region);
 220     energyImage->Allocate();
 221
 222     //Create the iterator on the output
 223     itk::ImageRegionIterator<InputImageType> imageOutputIterator(energyImage,energyImage->GetLargestPossibleRegion());
 224
 225     //Create the iterator on the region of the input
 226     typename InputImageType::IndexType startIndexIterator;
 227     startIndexIterator[0] = 0;
 228     startIndexIterator[1] = 0;
 229     startIndexIterator[2] = energy*nbSlicePerEnergy;
 230     typename InputImageType::RegionType regionIterator;
 231     regionIterator.SetSize(regionSize);
 232     regionIterator.SetIndex(startIndexIterator);
 233     itk::ImageRegionIterator<InputImageType> imageInputIterator(input,regionIterator);
 234
 235     //Copy the requested region
 236     while(!imageInputIterator.IsAtEnd())
 237     {
 238       imageOutputIterator.Set(imageInputIterator.Get());
 239
 240       ++imageInputIterator;
 241       ++imageOutputIterator;
 242     }
 243
 244     // Output directory and filenames
 245     typename WriterType::Pointer writer = WriterType::New();
 246     writer->SetInput( energyImage );
 247     writer->SetImageIO( gdcmIO );
 248     typename ReaderType::FileNamesContainer fileNamesOutput;
 249     std::string extension = NumberToString(energy) + ".dcm";
 250     std::string filename_out = m_ArgsInfo.output_arg + extension;
 251     fileNamesOutput.push_back(filename_out);
 252     writer->SetFileNames( fileNamesOutput );
 253     writer->SetMetaDataDictionaryArray(&outputArray);
 254
 255     // Write
 256     try {
 257       if (m_ArgsInfo.verbose_flag)
 258         std::cout << writer << std::endl;
 259       writer->Update();
 260     } catch( itk::ExceptionObject & excp ) {
 261       std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 262       std::cerr << excp << std::endl;
 263     }
 264
 265     //Write sequence dicom tag with gdcm
 266     gdcm::Reader reader;
 267     reader.SetFileName( filename_out.c_str() );
 268     reader.Read();
 269     gdcm::File &fileOutput = reader.GetFile();
 270     gdcm::DataSet &dsOutput = fileOutput.GetDataSet();
 271     const unsigned int ptr_len = 42;
 272     char *ptr = new char[ptr_len];
 273     memset(ptr,0,ptr_len);
 274     //Write rotation tag
 275     // Create a Sequence
 276     gdcm::SmartPointer<gdcm::SequenceOfItems> rotationSq = new gdcm::SequenceOfItems();
 277     rotationSq->SetLengthToUndefined();
 278     // Create a dataelement
 279     gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) );
 280     rotationDE.SetVR( gdcm::VR::US );
 281     char essai = (char)nbRotation;
 282     char *p_essai = &essai;
 283     rotationDE.SetByteValue(p_essai, 1);
 284     // Create an item
 285     gdcm::Item rotationIt;
 286     rotationIt.SetVLToUndefined();
 287     gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet();
 288     rotationDS.Insert(rotationDE);
 289     rotationSq->AddItem(rotationIt);
 290     // Insert sequence into data set
 291     gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) );
 292     rotationDEParent.SetVR(gdcm::VR::SQ);
 293     rotationDEParent.SetValue(*rotationSq);
 294     rotationDEParent.SetVLToUndefined();
 295     //Write energy
 296     gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) );
 297     energyDEParent.SetVR(gdcm::VR::SQ);
 298     // Create a Sequence
 299     gdcm::SmartPointer<gdcm::SequenceOfItems> energySq = new gdcm::SequenceOfItems();
 300     energySq->SetLengthToUndefined();
 301
 302     //Read the caracteristics of this energy window : the name and the thresholds up and down
 303     const gdcm::DataElement &seqWindow = dsInput.GetDataElement(gdcm::Tag(0x54, 0x12));
 304     gdcm::SmartPointer<gdcm::SequenceOfItems> sqiWindow = seqWindow.GetValueAsSQ();
 305     gdcm::Item &itemWindow = sqiWindow->GetItem(energy+1);
 306     const gdcm::DataElement &deNameWindow = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x18));
 307     std::string nameWindow( deNameWindow.GetByteValue()->GetPointer(), deNameWindow.GetByteValue()->GetLength() );
 308
 309     const gdcm::DataElement &deThresholds = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x13));
 310     gdcm::SmartPointer<gdcm::SequenceOfItems> sqiThresholds = deThresholds.GetValueAsSQ();
 311     for (unsigned int nbWindow=1; nbWindow <= sqiThresholds->GetNumberOfItems() ; ++nbWindow) {
 312       gdcm::Item &itemThresholds = sqiThresholds->GetItem(nbWindow);
 313       const gdcm::DataElement &deThresholdDown = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x14));
 314       std::string p_StrThresholdDown( deThresholdDown.GetByteValue()->GetPointer(), deThresholdDown.GetByteValue()->GetLength() );
 315       double thresholdDown = (double)atof(p_StrThresholdDown.c_str());
 316       const gdcm::DataElement &deThresholdUp = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x15));
 317       std::string p_StrThresholdUp( deThresholdUp.GetByteValue()->GetPointer(), deThresholdUp.GetByteValue()->GetLength() );
 318       double thresholdUp = (double)atof(p_StrThresholdUp.c_str());
 319
 320       //Now write it !
 321       gdcm::SmartPointer<gdcm::SequenceOfItems> energyThresholdSq = new gdcm::SequenceOfItems();
 322       energyThresholdSq->SetLengthToUndefined();
 323       // Create a dataelement
 324       gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) );
 325       gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) );
 326       energyThresholdDE.SetVR( gdcm::VR::DS );
 327       energyUpholdDE.SetVR( gdcm::VR::DS );
 328       energyThresholdDE.SetByteValue(NumberToString(thresholdDown).c_str(), (uint32_t)strlen(NumberToString(thresholdDown).c_str()));
 329       energyUpholdDE.SetByteValue(NumberToString(thresholdUp).c_str(), (uint32_t)strlen(NumberToString(thresholdUp).c_str()));
 330       // Create an item
 331       gdcm::Item energyThresholdIt;
 332       energyThresholdIt.SetVLToUndefined();
 333       gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet();
 334       energyThresholdDS.Insert(energyThresholdDE);
 335       energyThresholdDS.Insert(energyUpholdDE);
 336       energyThresholdSq->AddItem(energyThresholdIt);
 337       // Insert sequence into data set
 338       gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) );
 339       energyThresholdDEParent.SetVR(gdcm::VR::SQ);
 340       energyThresholdDEParent.SetValue(*energyThresholdSq);
 341       energyThresholdDEParent.SetVLToUndefined();
 342       // Create a dataelement
 343       gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) );
 344       energyNameDE.SetVR( gdcm::VR::SH );
 345       energyNameDE.SetByteValue(nameWindow.c_str(), (uint32_t)strlen(nameWindow.c_str()));
 346       // Create an item
 347       gdcm::Item energyIt;
 348       energyIt.SetVLToUndefined();
 349       gdcm::DataSet &energyDS = energyIt.GetNestedDataSet();
 350       energyDS.Insert(energyNameDE);
 351       energyDS.Insert(energyThresholdDEParent);
 352       energySq->AddItem(energyIt);
 353     }
 354     // Insert sequence into data set
 355     energyDEParent.SetValue(*energySq);
 356     energyDEParent.SetVLToUndefined();
 357     dsOutput.Insert(energyDEParent);
 358     dsOutput.Insert(rotationDEParent);
 359     gdcm::Writer w;
 360     w.SetFile( fileOutput );
 361     w.SetFileName( filename_out.c_str() );
 362     w.Write();
 363   }
 364 #else
 365   std::cout << "Use GDCM2" << std::endl;
 366 #endif
 367 }
 368 //---------------------------------------------------------------------------
 369
 370
 371 //---------------------------------------------------------------------------
 372 void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict)
 373 {
 374   typedef itk::MetaDataDictionary DictionaryType;
 375
 376   DictionaryType::ConstIterator itr = fromDict.Begin();
 377   DictionaryType::ConstIterator end = fromDict.End();
 378   typedef itk::MetaDataObject< std::string > MetaDataStringType;
 379
 380   while( itr != end )
 381     {
 382     itk::MetaDataObjectBase::Pointer  entry = itr->second;
 383
 384     MetaDataStringType::Pointer entryvalue =
 385       dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
 386     if( entryvalue )
 387       {
 388       std::string tagkey   = itr->first;
 389       std::string tagvalue = entryvalue->GetMetaDataObjectValue();
 390       itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue);
 391       }
 392     ++itr;
 393     }
 394 }
 395 //---------------------------------------------------------------------------
 396
 397
 398 //---------------------------------------------------------------------------
 399 template <typename T> std::string NumberToString ( T Number )
 400 {
 401    std::ostringstream ss;
 402    ss << Number;
 403    return ss.str();
 404 }
 405 //---------------------------------------------------------------------------
 406
 407 }//end clitk
 408
 409 #endif //#define clitkPartitionEnergyWindowDicomGenericFilter_txx