tools/clitkBinarizeImageGenericFilter.cxx

   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 clitkBinarizeImageGenericFilter_cxx
  19 #define clitkBinarizeImageGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkBinarizeImageGenericFilter.cxx
  23  * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
  24  * @date   29 june 2009
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkBinarizeImageGenericFilter.h"
  31
  32 // itk include
  33 #include "itkBinaryThresholdImageFilter.h"
  34 #include "itkMaskImageFilter.h"
  35 #include "itkMaskNegatedImageFilter.h"
  36
  37 #include <clitkCommon.h>
  38
  39 namespace clitk
  40 {
  41
  42 //--------------------------------------------------------------------
  43 BinarizeImageGenericFilter::BinarizeImageGenericFilter():
  44   ImageToImageGenericFilter<Self>("BinarizeImage")
  45 {
  46   InitializeImageType<2>();
  47   InitializeImageType<3>();
  48   InitializeImageType<4>();
  49 }
  50 //--------------------------------------------------------------------
  51
  52
  53 //--------------------------------------------------------------------
  54 template<unsigned int Dim>
  55 void BinarizeImageGenericFilter::InitializeImageType()
  56 {
  57   ADD_DEFAULT_IMAGE_TYPES(Dim);
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 void BinarizeImageGenericFilter::SetArgsInfo(const args_info_type & a)
  64 {
  65   mArgsInfo=a;
  66   if (mArgsInfo.verbose_given)
  67     SetIOVerbose(mArgsInfo.verbose_flag);
  68   if (mArgsInfo.imagetypes_given && mArgsInfo.imagetypes_flag)
  69     this->PrintAvailableImageTypes();
  70
  71   if (mArgsInfo.input_given) {
  72     SetInputFilename(mArgsInfo.input_arg);
  73   }
  74   if (mArgsInfo.output_given) {
  75     SetOutputFilename(mArgsInfo.output_arg);
  76   }
  77   if (mArgsInfo.percentage_given) {
  78     SetPercentage(mArgsInfo.percentage_arg);
  79   }
  80 }
  81 //--------------------------------------------------------------------
  82
  83
  84 //--------------------------------------------------------------------
  85 // Update with the number of dimensions and the pixeltype
  86 //--------------------------------------------------------------------
  87 template<class InputImageType>
  88 void
  89 BinarizeImageGenericFilter::UpdateWithInputImageType()
  90 {
  91   if (mArgsInfo.percentage_given)
  92     MaskOfIntegratedIntensity<InputImageType>();
  93   else {
  94     // Reading input
  95     typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
  96
  97     // Main filter
  98     typedef typename InputImageType::PixelType PixelType;
  99     typedef itk::Image<uchar, InputImageType::ImageDimension> OutputImageType;
 100
 101     // Filter
 102     typedef itk::BinaryThresholdImageFilter<InputImageType, OutputImageType> BinaryThresholdImageFilterType;
 103     typename BinaryThresholdImageFilterType::Pointer thresholdFilter=BinaryThresholdImageFilterType::New();
 104     thresholdFilter->SetInput(input);
 105     thresholdFilter->SetInsideValue(mArgsInfo.fg_arg);
 106
 107     if (mArgsInfo.lower_given) thresholdFilter->SetLowerThreshold(static_cast<PixelType>(mArgsInfo.lower_arg));
 108     if (mArgsInfo.upper_given) thresholdFilter->SetUpperThreshold(static_cast<PixelType>(mArgsInfo.upper_arg));
 109
 110     /* Three modes :
 111        - FG -> only use FG value for pixel in the Foreground (or Inside), keep input values for outside
 112        - BG -> only use BG value for pixel in the Background (or Outside), keep input values for inside
 113        - both -> use FG and BG (real binary image)
 114     */
 115     if (mArgsInfo.mode_arg == std::string("both")) {
 116       thresholdFilter->SetOutsideValue(mArgsInfo.bg_arg);
 117       thresholdFilter->Update();
 118       typename OutputImageType::Pointer outputImage = thresholdFilter->GetOutput();
 119       this->template SetNextOutput<OutputImageType>(outputImage);
 120     } else {
 121       typename InputImageType::Pointer outputImage;
 122       thresholdFilter->SetOutsideValue(0);
 123       if (mArgsInfo.mode_arg == std::string("BG")) {
 124         typedef itk::MaskImageFilter<InputImageType,OutputImageType> maskFilterType;
 125         typename maskFilterType::Pointer maskFilter = maskFilterType::New();
 126         maskFilter->SetInput1(input);
 127         maskFilter->SetInput2(thresholdFilter->GetOutput());
 128         maskFilter->SetOutsideValue(mArgsInfo.bg_arg);
 129         maskFilter->Update();
 130         outputImage = maskFilter->GetOutput();
 131       } else {
 132         typedef itk::MaskNegatedImageFilter<InputImageType,OutputImageType> maskFilterType;
 133         typename maskFilterType::Pointer maskFilter = maskFilterType::New();
 134         maskFilter->SetInput1(input);
 135         maskFilter->SetInput2(thresholdFilter->GetOutput());
 136         maskFilter->SetOutsideValue(mArgsInfo.fg_arg);
 137         maskFilter->Update();
 138         outputImage = maskFilter->GetOutput();
 139       }
 140       // Write/Save results
 141       this->template SetNextOutput<InputImageType>(outputImage);
 142     }
 143   }
 144 }
 145 //--------------------------------------------------------------------
 146
 147
 148 //--------------------------------------------------------------------
 149 //  https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes
 150 template <typename T>
 151 std::vector<size_t> sort_indexes(const std::vector<T> &v) {
 152
 153   // initialize original index locations
 154   std::vector<size_t> idx(v.size());
 155   std::vector<std::pair<T, size_t> > compVector(v.size());
 156   for (size_t i = 0; i < v.size(); ++i) {
 157     compVector[i].first = v[i];
 158     compVector[i].second = i;
 159   }
 160
 161   // sort indexes based on comparing values in v
 162   std::sort(compVector.begin(), compVector.end(), comparator<T>);
 163   for (size_t i = 0; i < v.size(); ++i) {
 164     idx[i] = compVector[i].second;
 165   }
 166
 167   return idx;
 168 }
 169 //--------------------------------------------------------------------
 170
 171
 172 //--------------------------------------------------------------------
 173 // Update with the number of dimensions and the pixeltype
 174 //--------------------------------------------------------------------
 175 template<class InputImageType>
 176 void
 177 BinarizeImageGenericFilter::MaskOfIntegratedIntensity()
 178 {
 179   // Main filter
 180   typedef typename InputImageType::PixelType InputPixelType;
 181   typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
 182
 183   // Reading input
 184   typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 185
 186   typename MaskImageType::Pointer mask;
 187   mask = MaskImageType::New();
 188   mask->SetRegions(input->GetLargestPossibleRegion());
 189   mask->SetOrigin(input->GetOrigin());
 190   mask->SetSpacing(input->GetSpacing());
 191   mask->Allocate();
 192   mask->FillBuffer(0);
 193
 194   // Get a vector of all values (will be easier to sort)
 195   // And compute total sum of values
 196   std::vector<double> values;
 197   typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
 198   IteratorInputType iter(input, input->GetLargestPossibleRegion());
 199   iter.GoToBegin();
 200   double total = 0.0;
 201   while (!iter.IsAtEnd()) {
 202     values.push_back(iter.Get());
 203     total += iter.Get();
 204     ++iter;
 205   }
 206
 207   // Sort (reverse)
 208   std::vector<size_t> indices = sort_indexes(values);
 209
 210   // Get max index of pixel to reach xx percent
 211   double current = 0.0;
 212   double max = GetPercentage()/100.0*total;
 213   int i=0;
 214   int n = input->GetLargestPossibleRegion().GetNumberOfPixels();
 215   std::vector<int> should_keep(values.size());;
 216   std::fill(should_keep.begin(), should_keep.end(), 0);
 217   while (current<max and i<n) { // loop by decreasing pixel values
 218     current += values[indices[i]];
 219     should_keep[indices[i]] = 1.0;
 220     ++i;
 221   }
 222   int nb = i;
 223
 224   // Set mask values
 225   typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
 226   IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
 227   iter.GoToBegin();
 228   itm.GoToBegin();
 229   i = 0;
 230   while (!iter.IsAtEnd()) {
 231     if (should_keep[i]) itm.Set(1);
 232     ++iter;
 233     ++itm;
 234     ++i;
 235   }
 236
 237   // Verbose option
 238   if (this->m_IOVerbose)
 239     std::cout << "Sum of pixel values : " << total << std::endl
 240               << "Percentage          : " << GetPercentage() << "%" << std::endl
 241               << "Number of pixels    : " << nb << "/" << n << std::endl
 242               << "Number of pixels    : " << nb/n*100.0 << "%" << std::endl;
 243
 244   // Write/Save results
 245   this->template SetNextOutput<MaskImageType>(mask);
 246 }
 247 //--------------------------------------------------------------------
 248
 249 //--------------------------------------------------------------------
 250 template <typename T>
 251 bool comparator ( const std::pair<T, size_t>& l, const std::pair<T, size_t>& r)
 252  { return l.first > r.first; }
 253 //--------------------------------------------------------------------
 254
 255 }//end clitk
 256
 257 #endif  //#define clitkBinarizeImageGenericFilter_cxx