tools/clitkMedianTemporalDimensionGenericFilter.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 clitkMedianTemporalDimensionGenericFilter_txx
  19 #define clitkMedianTemporalDimensionGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkMedianTemporalDimensionGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include <algorithm>
  31
  32 namespace clitk
  33 {
  34
  35   //-----------------------------------------------------------
  36   // Constructor
  37   //-----------------------------------------------------------
  38   template<class args_info_type>
  39     MedianTemporalDimensionGenericFilter<args_info_type>::MedianTemporalDimensionGenericFilter()
  40     {
  41       m_Verbose=false;
  42       m_InputFileName="";
  43     }
  44
  45
  46   //-----------------------------------------------------------
  47   // Update
  48   //-----------------------------------------------------------
  49   template<class args_info_type>
  50     void MedianTemporalDimensionGenericFilter<args_info_type>::Update()
  51     {
  52       // Read the Dimension and PixelType
  53       int Dimension, Components;
  54       std::string PixelType;
  55       ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
  56
  57
  58       // Call UpdateWithDim
  59       if (m_ArgsInfo.input_given>1) Dimension+=1;
  60       if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
  61       else if(Dimension==3)UpdateWithDim<3>(PixelType, Components);
  62       else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
  63       else {
  64         std::cout<<"Error, Only for 2, 3 or 4D!!!"<<std::endl ;
  65         return;
  66       }
  67     }
  68
  69   //-------------------------------------------------------------------
  70   // Update with the number of dimensions
  71   //-------------------------------------------------------------------
  72   template<class args_info_type>
  73     template<unsigned int Dimension>
  74     void
  75     MedianTemporalDimensionGenericFilter<args_info_type>::UpdateWithDim(const std::string PixelType, const int Components)
  76     {
  77       if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D with "<<Components<<" component(s) of "<<  PixelType<<"..."<<std::endl;
  78
  79       if (Components==1) {
  80         if(PixelType == "short") {
  81           if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  82           UpdateWithDimAndPixelType<Dimension, signed short>();
  83         }
  84         else if (PixelType == "unsigned_char") {
  85           if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  86           UpdateWithDimAndPixelType<Dimension, unsigned char>();
  87         }
  88         else if (PixelType == "float") {
  89           if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  90           UpdateWithDimAndPixelType<Dimension, float>();
  91         }
  92         else if (PixelType == "double") {
  93           if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
  94           UpdateWithDimAndPixelType<Dimension, double>();
  95         }
  96         else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
  97       }
  98     }
  99
 100
 101   //-------------------------------------------------------------------
 102   // Update with the number of dimensions and the pixeltype
 103   //-------------------------------------------------------------------
 104   template<class args_info_type>
 105     template <unsigned int Dimension, class  PixelType>
 106     void
 107     MedianTemporalDimensionGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 108     {
 109
 110       // ImageTypes
 111       typedef itk::Image<PixelType, Dimension> InputImageType;
 112       typedef itk::Image<PixelType, Dimension-1> OutputImageType;
 113       typename InputImageType::Pointer input;
 114
 115       if (m_ArgsInfo.input_given ==1 ) {
 116         // Read the input
 117         typedef itk::ImageFileReader<InputImageType> InputReaderType;
 118         typename InputReaderType::Pointer reader = InputReaderType::New();
 119         reader->SetFileName( m_InputFileName);
 120         reader->Update();
 121         input= reader->GetOutput();
 122       }
 123
 124       else {
 125         // Read and join multiple inputs
 126         if (m_Verbose) std::cout<<m_ArgsInfo.input_given<<" inputs given..."<<std::endl;
 127         std::vector<std::string> filenames;
 128         for(unsigned int i=0; i<m_ArgsInfo.input_given; i++) {
 129           if (m_Verbose) std::cout<<m_ArgsInfo.input_arg[i]<<std::endl;
 130           filenames.push_back(m_ArgsInfo.input_arg[i]);
 131         }
 132
 133         typedef itk::ImageSeriesReader<InputImageType> ImageReaderType;
 134         typename  ImageReaderType::Pointer reader= ImageReaderType::New();
 135         reader->SetFileNames(filenames);
 136         reader->Update();
 137         input =reader->GetOutput();
 138       }
 139
 140
 141       // Output properties
 142       typename OutputImageType::RegionType region;
 143       typename OutputImageType::RegionType::SizeType size;
 144       typename OutputImageType::IndexType index;
 145       typename OutputImageType::SpacingType spacing;
 146       typename OutputImageType::PointType origin;
 147       typename InputImageType::RegionType region4D=input->GetLargestPossibleRegion();
 148       typename InputImageType::RegionType::SizeType size4D=region4D.GetSize();
 149       typename InputImageType::IndexType index4D=region4D.GetIndex();
 150       typename InputImageType::SpacingType spacing4D=input->GetSpacing();
 151       typename InputImageType::PointType origin4D=input->GetOrigin();
 152
 153       for (unsigned int i=0; i< Dimension-1; i++) {
 154         size[i]=size4D[i];
 155         index[i]=index4D[i];
 156         spacing[i]=spacing4D[i];
 157         origin[i]=origin4D[i];
 158       }
 159       region.SetSize(size);
 160       region.SetIndex(index);
 161       typename OutputImageType::Pointer output= OutputImageType::New();
 162       output->SetRegions(region);
 163       output->SetSpacing(spacing);
 164       output->SetOrigin(origin);
 165       output->Allocate();
 166
 167
 168       // Region iterators
 169       typedef itk::ImageRegionIterator<InputImageType> IteratorType;
 170       std::vector<IteratorType> iterators(size4D[Dimension-1]);
 171       for (unsigned int i=0; i< size4D[Dimension-1]; i++) {
 172         typename InputImageType::RegionType regionIt=region4D;
 173         typename InputImageType::RegionType::SizeType sizeIt=regionIt.GetSize();
 174         sizeIt[Dimension-1]=1;
 175         regionIt.SetSize(sizeIt);
 176         typename InputImageType::IndexType indexIt=regionIt.GetIndex();
 177         indexIt[Dimension-1]=i;
 178         regionIt.SetIndex(indexIt);
 179         iterators[i]=IteratorType(input, regionIt);
 180       }
 181
 182       typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
 183       OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
 184
 185       // Take the median
 186       double value;
 187       typename InputImageType::IndexValueType temporal_dimension = size4D[Dimension-1];
 188       std::vector<PixelType> temp(temporal_dimension);
 189       while (!(iterators[0]).IsAtEnd()) {
 190         value=0.;
 191         for (unsigned int i=0; i<temporal_dimension; i++) {
 192           temp[i] = iterators[i].Get();
 193           ++(iterators[i]);
 194         }
 195         if (temporal_dimension & 1) {
 196           nth_element(temp.begin(), temp.begin() + temporal_dimension/2,temp.end());
 197           value = temp[temporal_dimension/2];
 198         } else {
 199           nth_element(temp.begin(), temp.begin() + temporal_dimension/2 - 1, temp.end());
 200           value = temp[temporal_dimension/2 - 1];
 201           nth_element(temp.begin(), temp.begin() + temporal_dimension/2, temp.end());
 202           value += temp[temporal_dimension/2];
 203           value /= 2;
 204         }
 205         avIt.Set(value);
 206         ++avIt;
 207       }
 208
 209   // Output
 210   typedef itk::ImageFileWriter<OutputImageType> WriterType;
 211   typename WriterType::Pointer writer = WriterType::New();
 212   writer->SetFileName(m_ArgsInfo.output_arg);
 213   writer->SetInput(output);
 214   writer->Update();
 215
 216 }
 217
 218
 219 }//end clitk
 220
 221 #endif //#define clitkMedianTemporalDimensionGenericFilter_txx