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     iota(idx.begin(), idx.end(), 0);
  78
  79     // sort indexes based on comparing values in v
  80     std::sort(idx.begin(), idx.end(),
  81               [&v](size_t i1, size_t i2) {return v[i1] > v[i2];});
  82
  83     return idx;
  84   }
  85   //--------------------------------------------------------------------
  86
  87
  88   //--------------------------------------------------------------------
  89   // Update with the number of dimensions and the pixeltype
  90   //--------------------------------------------------------------------
  91   template<class args_info_type>
  92   template<class InputImageType>
  93   void
  94   MaskOfIntegratedIntensityGenericFilter<args_info_type>::UpdateWithInputImageType()
  95   {
  96     // Main filter
  97     typedef typename InputImageType::PixelType InputPixelType;
  98     typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
  99
 100     // Reading input
 101     typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 102
 103     typename MaskImageType::Pointer mask;
 104     mask = MaskImageType::New();
 105     mask->SetRegions(input->GetLargestPossibleRegion());
 106     mask->SetOrigin(input->GetOrigin());
 107     mask->SetSpacing(input->GetSpacing());
 108     mask->Allocate();
 109     mask->FillBuffer(0);
 110
 111     // Get a vector of all values (will be easier to sort)
 112     // And compute total sum of values
 113     std::vector<double> values;
 114     typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
 115     IteratorInputType iter(input, input->GetLargestPossibleRegion());
 116     iter.GoToBegin();
 117     double total = 0.0;
 118     while (!iter.IsAtEnd()) {
 119       values.push_back(iter.Get());
 120       total += iter.Get();
 121       ++iter;
 122     }
 123
 124     // Sort (reverse)
 125     auto indices = sort_indexes(values);
 126
 127     // Get max index of pixel to reach xx percent
 128     double current = 0.0;
 129     double max = GetPercentage()/100.0*total;
 130     int i=0;
 131     auto n = input->GetLargestPossibleRegion().GetNumberOfPixels();
 132     std::vector<int> should_keep(values.size());;
 133     std::fill(should_keep.begin(), should_keep.end(), 0);
 134     while (current<max and i<n) { // loop by decreasing pixel values
 135       current += values[indices[i]];
 136       should_keep[indices[i]] = 1.0;
 137       ++i;
 138     }
 139     int nb = i;
 140
 141     // Set mask values
 142     typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
 143     IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
 144     iter.GoToBegin();
 145     itm.GoToBegin();
 146     i = 0;
 147     while (!iter.IsAtEnd()) {
 148       if (should_keep[i]) itm.Set(1);
 149       ++iter;
 150       ++itm;
 151       ++i;
 152     }
 153
 154     // Verbose option
 155     if (this->m_IOVerbose)
 156       std::cout << "Sum of pixel values : " << total << std::endl
 157                 << "Percentage          : " << GetPercentage() << "%" << std::endl
 158                 << "Number of pixels    : " << nb << "/" << n << std::endl
 159                 << "Number of pixels    : " << nb/n*100.0 << "%" << std::endl;
 160
 161     // Write/Save results
 162     this->template SetNextOutput<MaskImageType>(mask);
 163   }
 164   //--------------------------------------------------------------------
 165
 166
 167 }//end clitk
 168
 169 #endif //#define clitkMaskOfIntegratedIntensityGenericFilter_txx