tools/clitkWriteDicomSeriesGenericFilter.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://oncora1.lyon.fnclcc.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 clitkWriteDicomSeriesGenericFilter_txx
  19 #define clitkWriteDicomSeriesGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkWriteDicomSeriesGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 // clitk
  31 #include "clitkResampleImageWithOptionsFilter.h"
  32
  33
  34 namespace clitk
  35 {
  36
  37
  38 //-----------------------------------------------------------
  39 // Constructor
  40 //-----------------------------------------------------------
  41 template<class args_info_type>
  42 WriteDicomSeriesGenericFilter<args_info_type>::WriteDicomSeriesGenericFilter()
  43 {
  44   m_Verbose=false;
  45   m_InputFileName="";
  46 }
  47
  48
  49 //-----------------------------------------------------------
  50 // Update
  51 //-----------------------------------------------------------
  52 template<class args_info_type>
  53 void WriteDicomSeriesGenericFilter<args_info_type>::Update()
  54 {
  55   // Read the Dimension and PixelType
  56   int Dimension;
  57   std::string PixelType;
  58   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  59
  60
  61   // Call UpdateWithDim
  62   if(Dimension==2) UpdateWithDim<2>(PixelType);
  63   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  64   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  65   else {
  66     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  67     return;
  68   }
  69 }
  70
  71 //-------------------------------------------------------------------
  72 // Update with the number of dimensions
  73 //-------------------------------------------------------------------
  74 template<class args_info_type>
  75 template<unsigned int Dimension>
  76 void
  77 WriteDicomSeriesGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
  78 {
  79   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  80
  81   if(PixelType == "short") {
  82     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  83     UpdateWithDimAndPixelType<Dimension, signed short>();
  84   }
  85   //    else if(PixelType == "unsigned_short"){
  86   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  87   //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  88   //     }
  89
  90   else if (PixelType == "unsigned_char") {
  91     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  92     UpdateWithDimAndPixelType<Dimension, unsigned char>();
  93   }
  94
  95   //     else if (PixelType == "char"){
  96   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  97   //       UpdateWithDimAndPixelType<Dimension, signed char>();
  98   //     }
  99   else {
 100     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 101     UpdateWithDimAndPixelType<Dimension, float>();
 102   }
 103 }
 104
 105 //-------------------------------------------------------------------
 106 // Update with the number of dimensions and the pixeltype read from
 107 // the dicom files. The MHD files may be resampled to match the
 108 // dicom spacing (and number of slices). Rounding errors in resampling
 109 // are handled by removing files when generating the output dicom
 110 // series.
 111 //-------------------------------------------------------------------
 112 template<class args_info_type>
 113 template <unsigned int Dimension, class  PixelType>
 114 void
 115 WriteDicomSeriesGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 116 {
 117
 118   // ImageTypes
 119   typedef itk::Image<PixelType, Dimension> InputImageType;
 120   typedef itk::Image<PixelType, 2> OutputImageType;
 121
 122   // Read the dicom directory
 123   typedef itk::ImageSeriesReader< InputImageType >     ReaderType;
 124   typedef itk::GDCMImageIO ImageIOType;
 125   typedef itk::GDCMSeriesFileNames NamesGeneratorType;
 126
 127   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 128   NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
 129   namesGenerator->SetInputDirectory( m_ArgsInfo.inputDir_arg );
 130   namesGenerator->SetOutputDirectory( m_ArgsInfo.outputDir_arg  );
 131   typename   ReaderType::FileNamesContainer filenames_in = namesGenerator->GetInputFileNames();
 132   typename   ReaderType::FileNamesContainer filenames_out = namesGenerator->GetOutputFileNames();
 133
 134   // Output the dicom files
 135   unsigned int numberOfFilenames =  filenames_in.size();
 136   if (m_Verbose) {
 137     std::cout << numberOfFilenames <<" were read in the directory "<<m_ArgsInfo.inputDir_arg<<"..."<<std::endl<<std::endl;
 138     for(unsigned int fni = 0; fni<numberOfFilenames; fni++) {
 139       std::cout << "filename # " << fni << " = ";
 140       std::cout << filenames_in[fni] << std::endl;
 141     }
 142   }
 143
 144   // Read the series
 145   typename ReaderType::Pointer reader = ReaderType::New();
 146   reader->SetImageIO( gdcmIO );
 147   reader->SetFileNames( filenames_in );
 148   try {
 149     reader->Update();
 150   } catch (itk::ExceptionObject &excp) {
 151     std::cerr << "Error: Exception thrown while reading the DICOM series!!" << std::endl;
 152     std::cerr << excp << std::endl;
 153   }
 154
 155   // Read the input (MHD file)
 156   typedef typename InputImageType::RegionType RegionType;
 157   typedef typename RegionType::SizeType SizeType;
 158   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 159   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
 160   volumeReader->SetFileName( m_InputFileName);
 161   volumeReader->Update();
 162
 163   typename InputImageType::Pointer input = volumeReader->GetOutput();
 164   if (input->GetSpacing() != reader->GetOutput()->GetSpacing())
 165   {
 166     // resampling is carried out on the fly if resolution between
 167     // the input mhd and input dicom series is different
 168     if (m_Verbose)
 169       std::cout << "Warning: The image spacing differs between the MHD file and the input dicom series. Performing resampling with default options (for advanced options, use clitkResampleImage)." << std::endl;
 170
 171     // Filter
 172     typedef clitk::ResampleImageWithOptionsFilter<InputImageType, InputImageType> ResampleImageFilterType;
 173     typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
 174     filter->SetInput(input);
 175     filter->SetVerboseOptions(m_Verbose);
 176     filter->SetGaussianFilteringEnabled(false);
 177     filter->SetDefaultPixelValue(0);
 178     filter->SetOutputSpacing(reader->GetOutput()->GetSpacing());
 179     const SizeType& input_size = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
 180     SizeType output_size;
 181     for (unsigned int i = 0; i < Dimension; i++)
 182       output_size[i] = input_size[i];
 183     filter->SetOutputSize(output_size);
 184     filter->Update();
 185     input = filter->GetOutput();
 186   }
 187
 188   //    In some cases, due to resampling approximation issues,
 189   //    the number of slices in the MHD file may be different (smaller)
 190   //    from the number of files in the template dicom directory.
 191   //    To avoid ITK generating an exception, we reduce the number
 192   //    of DCM files to be considered, and a warning is printed
 193   //    in verbose mode
 194   const RegionType volumeRegion = input->GetLargestPossibleRegion();
 195   const SizeType& volumeSize = volumeRegion.GetSize();
 196   if (Dimension == 3 && volumeSize[2] < numberOfFilenames)
 197   {
 198     if (m_Verbose)
 199       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;
 200
 201     filenames_in.resize(volumeSize[2]);
 202     filenames_out.resize(filenames_in.size());
 203     numberOfFilenames =  filenames_in.size();
 204   }
 205
 206   // Modify the meta dictionary
 207   typedef itk::MetaDataDictionary   DictionaryType;
 208   const std::vector<DictionaryType*>* dictionary = reader->GetMetaDataDictionaryArray();
 209
 210   // Get keys
 211   unsigned int numberOfKeysGiven=0;
 212   if(m_ArgsInfo.midP_flag && m_ArgsInfo.key_given)
 213     std::cerr<<"Error: both keys and midP option are given"<<std::endl;
 214   else if (m_ArgsInfo.midP_flag)
 215     numberOfKeysGiven=1;
 216   else
 217     numberOfKeysGiven=m_ArgsInfo.key_given;
 218
 219   for (unsigned int i = 0; i < numberOfKeysGiven; i++) {
 220     std::string entryId(m_ArgsInfo.key_arg[i]  );
 221     std::string value( m_ArgsInfo.tag_arg[i] );
 222     for(unsigned int fni = 0; fni<numberOfFilenames; fni++)
 223       itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), entryId, value );
 224   }
 225
 226   // Output directory and filenames
 227   itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputDir_arg ); // create if it doesn't exist
 228   typedef itk::ImageSeriesWriter<InputImageType, OutputImageType >  SeriesWriterType;
 229   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 230
 231   seriesWriter->SetInput( input );
 232   seriesWriter->SetImageIO( gdcmIO );
 233
 234   seriesWriter->SetFileNames( filenames_out );
 235   seriesWriter->SetMetaDataDictionaryArray( dictionary );
 236
 237   // Write
 238   try {
 239     seriesWriter->Update();
 240   } catch( itk::ExceptionObject & excp ) {
 241     std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 242     std::cerr << excp << std::endl;
 243   }
 244
 245 }
 246
 247 /*
 248 //-------------------------------------------------------------------
 249 // Update with the number of dimensions and the pixeltype
 250 //-------------------------------------------------------------------
 251 template<class args_info_type>
 252 template <unsigned int Dimension, class  PixelType>
 253 void
 254 WriteDicomSeriesGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 255 {
 256
 257   // ImageTypes
 258   typedef itk::Image<PixelType, Dimension> InputImageType;
 259   typedef itk::Image<PixelType, 2> OutputImageType;
 260
 261   // Read the input (volumetric)
 262   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 263   typename InputReaderType::Pointer volumeReader = InputReaderType::New();
 264   volumeReader->SetFileName( m_InputFileName);
 265   volumeReader->Update();
 266   typename InputImageType::Pointer input= volumeReader->GetOutput();
 267
 268   // Read the dicom directory
 269   typedef itk::ImageSeriesReader< InputImageType >     ReaderType;
 270   typedef itk::GDCMImageIO ImageIOType;
 271   typedef itk::GDCMSeriesFileNames NamesGeneratorType;
 272
 273   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 274   NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
 275   namesGenerator->SetInputDirectory( m_ArgsInfo.inputDir_arg );
 276   namesGenerator->SetOutputDirectory( m_ArgsInfo.outputDir_arg  );
 277   typename   ReaderType::FileNamesContainer filenames_in = namesGenerator->GetInputFileNames();
 278   typename   ReaderType::FileNamesContainer filenames_out = namesGenerator->GetOutputFileNames();
 279
 280   // Output the dicom files
 281   unsigned int numberOfFilenames =  filenames_in.size();
 282   if (m_Verbose) {
 283     std::cout << numberOfFilenames <<" were read in the directory "<<m_ArgsInfo.inputDir_arg<<"..."<<std::endl<<std::endl;
 284     for(unsigned int fni = 0; fni<numberOfFilenames; fni++) {
 285       std::cout << "filename # " << fni << " = ";
 286       std::cout << filenames_in[fni] << std::endl;
 287     }
 288   }
 289
 290   // RP: 16/03/2011
 291   //    In some cases, due to resampling approximation issues,
 292   //    the number of slices in the MHD file may be different
 293   //    from the number of slices in the template DCM directory.
 294   //    To avoid ITK generating an exception, we reduce the number
 295   //    of DCM files to be considered, provided the --force
 296   //    option is set.
 297   typedef typename InputImageType::RegionType RegionType;
 298   typedef typename RegionType::SizeType SizeType;
 299   const RegionType volumeRegion = input->GetLargestPossibleRegion();
 300   const SizeType& volumeSize = volumeRegion.GetSize();
 301   if (m_ArgsInfo.force_given && Dimension == 3 && volumeSize[2] < numberOfFilenames)
 302   {
 303     std::cout << "Warning: Number of files in " << m_ArgsInfo.inputDir_arg << " is greater than the number of slices in MHD. Using only " << volumeSize[2] << " files." << std::endl;
 304     filenames_in.resize(volumeSize[2]);
 305     filenames_out.resize(filenames_in.size());
 306     numberOfFilenames =  filenames_in.size();
 307   }
 308
 309   // Read the series
 310   typename ReaderType::Pointer reader = ReaderType::New();
 311   reader->SetImageIO( gdcmIO );
 312   reader->SetFileNames( filenames_in );
 313   try {
 314     reader->Update();
 315   } catch (itk::ExceptionObject &excp) {
 316     std::cerr << "Error: Exception thrown while reading the DICOM series!!" << std::endl;
 317     std::cerr << excp << std::endl;
 318   }
 319
 320
 321   // Modify the meta dictionary
 322   typedef itk::MetaDataDictionary   DictionaryType;
 323   const std::vector<DictionaryType*>* dictionary = reader->GetMetaDataDictionaryArray();
 324
 325   // Get keys
 326   unsigned int numberOfKeysGiven=0;
 327   if(m_ArgsInfo.midP_flag && m_ArgsInfo.key_given)
 328     std::cerr<<"Error: both keys and midP option are given"<<std::endl;
 329   else if (m_ArgsInfo.midP_flag)
 330     numberOfKeysGiven=1;
 331   else
 332     numberOfKeysGiven=m_ArgsInfo.key_given;
 333
 334   for (unsigned int i = 0; i < numberOfKeysGiven; i++) {
 335     std::string entryId(m_ArgsInfo.key_arg[i]  );
 336     std::string value( m_ArgsInfo.tag_arg[i] );
 337     for(unsigned int fni = 0; fni<numberOfFilenames; fni++)
 338       itk::EncapsulateMetaData<std::string>( *((*dictionary)[fni]), entryId, value );
 339   }
 340
 341   // Output directory and filenames
 342   itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputDir_arg ); // create if it doesn't exist
 343   typedef itk::ImageSeriesWriter<InputImageType, OutputImageType >  SeriesWriterType;
 344   typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
 345
 346   seriesWriter->SetInput( volumeReader->GetOutput() );
 347   seriesWriter->SetImageIO( gdcmIO );
 348
 349   seriesWriter->SetFileNames( filenames_out );
 350   seriesWriter->SetMetaDataDictionaryArray( dictionary );
 351
 352   // Write
 353   try {
 354     seriesWriter->Update();
 355   } catch( itk::ExceptionObject & excp ) {
 356     std::cerr << "Error: Exception thrown while writing the series!!" << std::endl;
 357     std::cerr << excp << std::endl;
 358   }
 359
 360 }
 361 */
 362
 363 }//end clitk
 364
 365 #endif //#define clitkWriteDicomSeriesGenericFilter_txx