Integrate MaskOfIntegratedIntensity into clitkBinarizeImage

diff --git a/tools/clitkBinarizeImage.ggo b/tools/clitkBinarizeImage.ggo
index dfc37b8021555f843eb454afaa9c29e01b1f8ee2..3808d33f30967b4e51ba60caeac87ca71771bf2b 100644 (file)
--- a/tools/clitkBinarizeImage.ggo
+++ b/tools/clitkBinarizeImage.ggo
@@ -17,3 +17,4 @@ option "bg"                   -       "Background (BG) or 'ouside' value"               double        no      default
 option "mode"          -       "Use FG and/or BG values (if FG, the BG is replaced by the input image values)" values="FG","BG","both"  no default="both"
 
 
+option "percentage"                p   "Percentage of total pixels values (in %)" double no

diff --git a/tools/clitkBinarizeImageGenericFilter.cxx b/tools/clitkBinarizeImageGenericFilter.cxx
index 46f846b0e772cc2ed3aad24ac683b445f5225da7..6affce00b4fb59bf535e06288b0cfd005e906d81 100644 (file)
--- a/tools/clitkBinarizeImageGenericFilter.cxx
+++ b/tools/clitkBinarizeImageGenericFilter.cxx
@@ -74,6 +74,9 @@ void BinarizeImageGenericFilter::SetArgsInfo(const args_info_type & a)
   if (mArgsInfo.output_given) {
     SetOutputFilename(mArgsInfo.output_arg);
   }
+  if (mArgsInfo.percentage_given) {
+    SetPercentage(mArgsInfo.percentage_arg);
+  }
 }
 //--------------------------------------------------------------------
 
@@ -85,60 +88,170 @@ template<class InputImageType>
 void
 BinarizeImageGenericFilter::UpdateWithInputImageType()
 {
+  if (mArgsInfo.percentage_given)
+    MaskOfIntegratedIntensity<InputImageType>();
+  else {
+    // Reading input
+    typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 
-  // Reading input
-  typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
+    // Main filter
+    typedef typename InputImageType::PixelType PixelType;
+    typedef itk::Image<uchar, InputImageType::ImageDimension> OutputImageType;
 
-  // Main filter
-  typedef typename InputImageType::PixelType PixelType;
-  typedef itk::Image<uchar, InputImageType::ImageDimension> OutputImageType;
-
-  // Filter
-  typedef itk::BinaryThresholdImageFilter<InputImageType, OutputImageType> BinaryThresholdImageFilterType;
-  typename BinaryThresholdImageFilterType::Pointer thresholdFilter=BinaryThresholdImageFilterType::New();
-  thresholdFilter->SetInput(input);
-  thresholdFilter->SetInsideValue(mArgsInfo.fg_arg);
-
-  if (mArgsInfo.lower_given) thresholdFilter->SetLowerThreshold(static_cast<PixelType>(mArgsInfo.lower_arg));
-  if (mArgsInfo.upper_given) thresholdFilter->SetUpperThreshold(static_cast<PixelType>(mArgsInfo.upper_arg));
-
-  /* Three modes :
-     - FG -> only use FG value for pixel in the Foreground (or Inside), keep input values for outside
-     - BG -> only use BG value for pixel in the Background (or Outside), keep input values for inside
-     - both -> use FG and BG (real binary image)
-  */
-  if (mArgsInfo.mode_arg == std::string("both")) {
-    thresholdFilter->SetOutsideValue(mArgsInfo.bg_arg);
-    thresholdFilter->Update();
-    typename OutputImageType::Pointer outputImage = thresholdFilter->GetOutput();
-    this->template SetNextOutput<OutputImageType>(outputImage);
-  } else {
-    typename InputImageType::Pointer outputImage;
-    thresholdFilter->SetOutsideValue(0);
-    if (mArgsInfo.mode_arg == std::string("BG")) {
-      typedef itk::MaskImageFilter<InputImageType,OutputImageType> maskFilterType;
-      typename maskFilterType::Pointer maskFilter = maskFilterType::New();
-      maskFilter->SetInput1(input);
-      maskFilter->SetInput2(thresholdFilter->GetOutput());
-      maskFilter->SetOutsideValue(mArgsInfo.bg_arg);
-      maskFilter->Update();
-      outputImage = maskFilter->GetOutput();
+    // Filter
+    typedef itk::BinaryThresholdImageFilter<InputImageType, OutputImageType> BinaryThresholdImageFilterType;
+    typename BinaryThresholdImageFilterType::Pointer thresholdFilter=BinaryThresholdImageFilterType::New();
+    thresholdFilter->SetInput(input);
+    thresholdFilter->SetInsideValue(mArgsInfo.fg_arg);
+
+    if (mArgsInfo.lower_given) thresholdFilter->SetLowerThreshold(static_cast<PixelType>(mArgsInfo.lower_arg));
+    if (mArgsInfo.upper_given) thresholdFilter->SetUpperThreshold(static_cast<PixelType>(mArgsInfo.upper_arg));
+
+    /* Three modes :
+       - FG -> only use FG value for pixel in the Foreground (or Inside), keep input values for outside
+       - BG -> only use BG value for pixel in the Background (or Outside), keep input values for inside
+       - both -> use FG and BG (real binary image)
+    */
+    if (mArgsInfo.mode_arg == std::string("both")) {
+      thresholdFilter->SetOutsideValue(mArgsInfo.bg_arg);
+      thresholdFilter->Update();
+      typename OutputImageType::Pointer outputImage = thresholdFilter->GetOutput();
+      this->template SetNextOutput<OutputImageType>(outputImage);
     } else {
-      typedef itk::MaskNegatedImageFilter<InputImageType,OutputImageType> maskFilterType;
-      typename maskFilterType::Pointer maskFilter = maskFilterType::New();
-      maskFilter->SetInput1(input);
-      maskFilter->SetInput2(thresholdFilter->GetOutput());
-      maskFilter->SetOutsideValue(mArgsInfo.fg_arg);
-      maskFilter->Update();
-      outputImage = maskFilter->GetOutput();
+      typename InputImageType::Pointer outputImage;
+      thresholdFilter->SetOutsideValue(0);
+      if (mArgsInfo.mode_arg == std::string("BG")) {
+        typedef itk::MaskImageFilter<InputImageType,OutputImageType> maskFilterType;
+        typename maskFilterType::Pointer maskFilter = maskFilterType::New();
+        maskFilter->SetInput1(input);
+        maskFilter->SetInput2(thresholdFilter->GetOutput());
+        maskFilter->SetOutsideValue(mArgsInfo.bg_arg);
+        maskFilter->Update();
+        outputImage = maskFilter->GetOutput();
+      } else {
+        typedef itk::MaskNegatedImageFilter<InputImageType,OutputImageType> maskFilterType;
+        typename maskFilterType::Pointer maskFilter = maskFilterType::New();
+        maskFilter->SetInput1(input);
+        maskFilter->SetInput2(thresholdFilter->GetOutput());
+        maskFilter->SetOutsideValue(mArgsInfo.fg_arg);
+        maskFilter->Update();
+        outputImage = maskFilter->GetOutput();
+      }
+      // Write/Save results
+      this->template SetNextOutput<InputImageType>(outputImage);
     }
-    // Write/Save results
-    this->template SetNextOutput<InputImageType>(outputImage);
   }
 }
 //--------------------------------------------------------------------
 
 
+//--------------------------------------------------------------------
+//  https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes
+template <typename T>
+std::vector<size_t> sort_indexes(const std::vector<T> &v) {
+
+  // initialize original index locations
+  std::vector<size_t> idx(v.size());
+  std::vector<std::pair<T, size_t> > compVector(v.size());
+  for (size_t i = 0; i < v.size(); ++i) {
+    compVector[i].first = v[i];
+    compVector[i].second = i;
+  }
+
+  // sort indexes based on comparing values in v
+  std::sort(compVector.begin(), compVector.end(), comparator<T>);
+  for (size_t i = 0; i < v.size(); ++i) {
+    idx[i] = compVector[i].second;
+  }
+
+  return idx;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+// Update with the number of dimensions and the pixeltype
+//--------------------------------------------------------------------
+template<class InputImageType>
+void
+BinarizeImageGenericFilter::MaskOfIntegratedIntensity()
+{
+  // Main filter
+  typedef typename InputImageType::PixelType InputPixelType;
+  typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
+
+  // Reading input
+  typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
+
+  typename MaskImageType::Pointer mask;
+  mask = MaskImageType::New();
+  mask->SetRegions(input->GetLargestPossibleRegion());
+  mask->SetOrigin(input->GetOrigin());
+  mask->SetSpacing(input->GetSpacing());
+  mask->Allocate();
+  mask->FillBuffer(0);
+
+  // Get a vector of all values (will be easier to sort)
+  // And compute total sum of values
+  std::vector<double> values;
+  typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
+  IteratorInputType iter(input, input->GetLargestPossibleRegion());
+  iter.GoToBegin();
+  double total = 0.0;
+  while (!iter.IsAtEnd()) {
+    values.push_back(iter.Get());
+    total += iter.Get();
+    ++iter;
+  }
+
+  // Sort (reverse)
+  std::vector<size_t> indices = sort_indexes(values);
+
+  // Get max index of pixel to reach xx percent
+  double current = 0.0;
+  double max = GetPercentage()/100.0*total;
+  int i=0;
+  int n = input->GetLargestPossibleRegion().GetNumberOfPixels();
+  std::vector<int> should_keep(values.size());;
+  std::fill(should_keep.begin(), should_keep.end(), 0);
+  while (current<max and i<n) { // loop by decreasing pixel values
+    current += values[indices[i]];
+    should_keep[indices[i]] = 1.0;
+    ++i;
+  }
+  int nb = i;
+
+  // Set mask values
+  typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
+  IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
+  iter.GoToBegin();
+  itm.GoToBegin();
+  i = 0;
+  while (!iter.IsAtEnd()) {
+    if (should_keep[i]) itm.Set(1);
+    ++iter;
+    ++itm;
+    ++i;
+  }
+
+  // Verbose option
+  if (this->m_IOVerbose)
+    std::cout << "Sum of pixel values : " << total << std::endl
+              << "Percentage          : " << GetPercentage() << "%" << std::endl
+              << "Number of pixels    : " << nb << "/" << n << std::endl
+              << "Number of pixels    : " << nb/n*100.0 << "%" << std::endl;
+
+  // Write/Save results
+  this->template SetNextOutput<MaskImageType>(mask);
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <typename T>
+bool comparator ( const std::pair<T, size_t>& l, const std::pair<T, size_t>& r)
+ { return l.first > r.first; }
+//--------------------------------------------------------------------
+
 }//end clitk
 
 #endif  //#define clitkBinarizeImageGenericFilter_cxx

diff --git a/tools/clitkBinarizeImageGenericFilter.h b/tools/clitkBinarizeImageGenericFilter.h
index c21e31e2b4de1c7548f3301522cb5498b3c0f54e..4a919f5e7f80f51018358d74430159a551806210 100644 (file)
--- a/tools/clitkBinarizeImageGenericFilter.h
+++ b/tools/clitkBinarizeImageGenericFilter.h
@@ -46,20 +46,28 @@ namespace clitk
 
     //--------------------------------------------------------------------
     void SetArgsInfo(const args_info_type & a);
+    void SetPercentage(double p) { mPercentage = p; }
+    double GetPercentage() const { return mPercentage; }
 
     //--------------------------------------------------------------------
     // Main function called each time the filter is updated
     template<class InputImageType>  
     void UpdateWithInputImageType();
 
+    template<class InputImageType>
+    void MaskOfIntegratedIntensity();
+
   protected:
     BinarizeImageGenericFilter();
     template<unsigned int Dim> void InitializeImageType();
     args_info_type mArgsInfo;
+    double mPercentage;
     
   }; // end class
   //--------------------------------------------------------------------
     
+//Implementation of the pair comparative function
+template <typename T> bool comparator ( const std::pair<T, size_t>& l, const std::pair<T, size_t>& r);
 } // end namespace clitk
 //--------------------------------------------------------------------