tools/clitkMaskOfIntegratedIntensityGenericFilter.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
  19 #ifndef clitkMaskOfIntegratedIntensityGenericFilter_txx
  20 #define clitkMaskOfIntegratedIntensityGenericFilter_txx
  21
  22 // itk include
  23 #include "itkIntensityWindowingImageFilter.h"
  24 #include "itkLabelStatisticsImageFilter.h"
  25 #include "itkMaskImageFilter.h"
  26 #include "itkMaskNegatedImageFilter.h"
  27 #include <clitkCommon.h>
  28 #include <numeric>
  29
  30 namespace clitk
  31 {
  32
  33   //--------------------------------------------------------------------
  34   template<class args_info_type>
  35   MaskOfIntegratedIntensityGenericFilter<args_info_type>::MaskOfIntegratedIntensityGenericFilter():
  36     ImageToImageGenericFilter<Self>("MaskOfIntegratedIntensity")
  37   {
  38     InitializeImageType<2>();
  39     InitializeImageType<3>();
  40   }
  41   //--------------------------------------------------------------------
  42
  43
  44   //--------------------------------------------------------------------
  45   template<class args_info_type>
  46   template<unsigned int Dim>
  47   void MaskOfIntegratedIntensityGenericFilter<args_info_type>::InitializeImageType()
  48   {
  49     ADD_DEFAULT_IMAGE_TYPES(Dim);
  50   }
  51   //--------------------------------------------------------------------
  52
  53
  54   //--------------------------------------------------------------------
  55   template<class args_info_type>
  56   void MaskOfIntegratedIntensityGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
  57   {
  58     mArgsInfo=a;
  59     this->SetIOVerbose(mArgsInfo.verbose_flag);
  60     if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
  61
  62     this->SetInputFilename(mArgsInfo.input_arg);
  63     this->SetOutputFilename(mArgsInfo.output_arg);
  64     this->SetPercentage(mArgsInfo.percentage_arg);
  65
  66   }
  67   //--------------------------------------------------------------------
  68
  69
  70   //--------------------------------------------------------------------
  71   //  https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes
  72   template <typename T>
  73   std::vector<size_t> sort_indexes(const std::vector<T> &v) {
  74
  75     // initialize original index locations
  76     std::vector<size_t> idx(v.size());
  77     std::vector<std::pair<T, size_t> > compVector(v.size());
  78     for (size_t i = 0; i < v.size(); ++i) {
  79       compVector[i].first = v[i];
  80       compVector[i].second = i;
  81     }
  82
  83     // sort indexes based on comparing values in v
  84     std::sort(compVector.begin(), compVector.end(), comparator<T>);
  85     for (size_t i = 0; i < v.size(); ++i) {
  86       idx[i] = compVector[i].second;
  87     }
  88
  89     return idx;
  90   }
  91   //--------------------------------------------------------------------
  92
  93
  94   //--------------------------------------------------------------------
  95   // Update with the number of dimensions and the pixeltype
  96   //--------------------------------------------------------------------
  97   template<class args_info_type>
  98   template<class InputImageType>
  99   void
 100   MaskOfIntegratedIntensityGenericFilter<args_info_type>::UpdateWithInputImageType()
 101   {
 102     // Main filter
 103     typedef typename InputImageType::PixelType InputPixelType;
 104     typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
 105
 106     // Reading input
 107     typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 108
 109     typename MaskImageType::Pointer mask;
 110     mask = MaskImageType::New();
 111     mask->SetRegions(input->GetLargestPossibleRegion());
 112     mask->SetOrigin(input->GetOrigin());
 113     mask->SetSpacing(input->GetSpacing());
 114     mask->Allocate();
 115     mask->FillBuffer(0);
 116
 117     // Get a vector of all values (will be easier to sort)
 118     // And compute total sum of values
 119     std::vector<double> values;
 120     typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
 121     IteratorInputType iter(input, input->GetLargestPossibleRegion());
 122     iter.GoToBegin();
 123     double total = 0.0;
 124     while (!iter.IsAtEnd()) {
 125       values.push_back(iter.Get());
 126       total += iter.Get();
 127       ++iter;
 128     }
 129
 130     // Sort (reverse)
 131     std::vector<size_t> indices = sort_indexes(values);
 132
 133     // Get max index of pixel to reach xx percent
 134     double current = 0.0;
 135     double max = GetPercentage()/100.0*total;
 136     int i=0;
 137     int n = input->GetLargestPossibleRegion().GetNumberOfPixels();
 138     std::vector<int> should_keep(values.size());;
 139     std::fill(should_keep.begin(), should_keep.end(), 0);
 140     while (current<max && i<n) { // loop by decreasing pixel values
 141       current += values[indices[i]];
 142       should_keep[indices[i]] = 1.0;
 143       ++i;
 144     }
 145     int nb = i;
 146
 147     // Set mask values
 148     typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
 149     IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
 150     iter.GoToBegin();
 151     itm.GoToBegin();
 152     i = 0;
 153     while (!iter.IsAtEnd()) {
 154       if (should_keep[i]) itm.Set(1);
 155       ++iter;
 156       ++itm;
 157       ++i;
 158     }
 159
 160     // Verbose option
 161     if (this->m_IOVerbose)
 162       std::cout << "Sum of pixel values : " << total << std::endl
 163                 << "Percentage          : " << GetPercentage() << "%" << std::endl
 164                 << "Number of pixels    : " << nb << "/" << n << std::endl
 165                 << "Number of pixels    : " << nb/n*100.0 << "%" << std::endl;
 166
 167     // Write/Save results
 168     this->template SetNextOutput<MaskImageType>(mask);
 169   }
 170   //--------------------------------------------------------------------
 171
 172 //--------------------------------------------------------------------
 173 template <typename T>
 174 bool comparator ( const std::pair<T, size_t>& l, const std::pair<T, size_t>& r)
 175    { return l.first > r.first; }
 176 //--------------------------------------------------------------------
 177
 178 }//end clitk
 179
 180 #endif //#define clitkMaskOfIntegratedIntensityGenericFilter_txx