otsu_threshold Template Function Implementation for Image in C++

Question

This is a follow-up question for An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++, histogram Template Function Implementation for Image in C++, Histogram of Image using std::map in C++, histogram_normalized and histogram_with_bins Template Functions Implementation for Image in C++ and normalize_histogram Template Function Implementation for Image in C++. In this post, a function for calculating OTSU threshold is implemented.

The experimental implementation

• otsu_threshold template function implementation

  namespace TinyDIP
  {
   // otsu_threshold template function implementation (with Execution Policy)
   template <class ExPo, std::integral ElementT>
   requires(std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
   constexpr static auto otsu_threshold(
   ExPo execution_policy,
   const Image<ElementT>& image)
   {
   auto probabilities = normalize_histogram(execution_policy, histogram(image));
   double maxVariance = 0.0;
   ElementT optimalThreshold = 0;
   if constexpr (std::same_as<ElementT, std::uint8_t> || std::same_as<ElementT, std::uint16_t>)
   {
   for (ElementT threshold = 0; threshold < std::numeric_limits<ElementT>::max(); ++threshold)
   {
   double w_background = 0.0;
   double w_foreground = 0.0;
   double m_background = 0.0;
   double m_foreground = 0.0;
   for (std::size_t i = 0; i <= threshold; ++i)
   {
   w_background += probabilities[i];
   m_background += i * probabilities[i];
   }
   if (w_background != 0)
   {
   m_background /= w_background;
   }
   for (std::size_t i = threshold + 1; i <= std::numeric_limits<ElementT>::max(); ++i)
   {
   w_foreground += probabilities[i];
   m_foreground += i * probabilities[i];
   }
   if (w_foreground != 0)
   {
   m_foreground /= w_foreground;
   }
   double variance = w_background * w_foreground * (m_background - m_foreground) * (m_background - m_foreground);
   if (variance > maxVariance)
   {
   maxVariance = variance;
   optimalThreshold = threshold;
   }
   }
   }
   else
   {
   auto probabilityVec = std::vector<std::pair<ElementT const, double>>(std::ranges::cbegin(probabilities), std::ranges::cend(probabilities));
   for (std::size_t i = 0; i < probabilityVec.size(); ++i)
   {
   auto const& [threshold, probability] = probabilityVec[i];
   double w_background = 0.0;
   double w_foreground = 0.0;
   double m_background = 0.0;
   double m_foreground = 0.0;
   for (std::size_t j = 0; j <= i; ++j)
   {
   w_background += probabilityVec[j].second;
   m_background += probabilityVec[j].first * probabilityVec[j].second;
   }
   if (w_background != 0)
   {
   m_background /= w_background;
   }
   for (std::size_t j = i + 1; j < probabilityVec.size(); ++j)
   {
   w_foreground += probabilityVec[j].second;
   m_foreground += probabilityVec[j].first * probabilityVec[j].second;
   }
   if (w_foreground != 0)
   {
   m_foreground /= w_foreground;
   }
   double variance = w_background * w_foreground * (m_background - m_foreground) * (m_background - m_foreground);
   if (variance >= maxVariance)
   {
   maxVariance = variance;
   optimalThreshold = threshold;
   }
   }
   }
   return optimalThreshold;
   }
   // otsu_threshold template function implementation
   template <std::integral ElementT>
   constexpr static auto otsu_threshold(const Image<ElementT>& image)
   {
   return otsu_threshold(std::execution::seq, image);
   }
  }
  

• sum_second_element template function implementation

  namespace TinyDIP
  {
   // sum_second_element Template Function Implementation
   template <typename FirstT, typename SecondT, class Function = std::plus<SecondT>>
   requires(std::regular_invocable<Function, SecondT, SecondT>)
   constexpr static SecondT sum_second_element(const std::vector<std::pair<FirstT, SecondT>>& pairs, const Function& f = Function{})
   {
   SecondT sum{};
   for (auto const& [first, second] : pairs)
   {
   sum = std::invoke(f, sum, second);
   }
   return sum;
   }
   // sum_second_element Template Function Implementation (with execution policy)
   template <class ExPo, typename FirstT, typename SecondT, class Function = std::plus<SecondT>>
   requires(std::is_execution_policy_v<std::remove_cvref_t<ExPo>> and
   std::regular_invocable<Function, SecondT, SecondT>)
   constexpr static SecondT sum_second_element(ExPo execution_policy, const std::vector<std::pair<FirstT, SecondT>>& pairs, const Function& f = Function{})
   {
   const std::size_t size = pairs.size();
   std::vector<SecondT> second_elements(size);
   std::transform(
   execution_policy,
   std::ranges::cbegin(pairs),
   std::ranges::cend(pairs),
   std::ranges::begin(second_elements),
   [&](auto const& pair) { return pair.second; });
   return std::reduce(execution_policy, std::ranges::cbegin(second_elements), std::ranges::cend(second_elements), SecondT{}, f);
   }
   // sum_second_element Template Function Implementation
   template <typename KeyT, typename ValueT, class Function = std::plus<ValueT>>
   requires(std::regular_invocable<Function, ValueT, ValueT>)
   constexpr static ValueT sum_second_element(const std::map<KeyT, ValueT>& map, const Function& f = Function{})
   {
   ValueT sum{};
   for (const auto& [key, value] : map)
   {
   sum = std::invoke(f, sum, value);
   }
   return sum;
   }
  }
  

• get_normalized_input template function implementation

  namespace TinyDIP
  {
   // get_normalized_input template function implementation
   // https://codereview.stackexchange.com/a/295540/231235
   template<class ElementT, std::size_t Count, class ProbabilityType = double>
   constexpr static auto get_normalized_input(
   const std::array<ElementT, Count>& input,
   const ProbabilityType& sum)
   {
   std::array<ProbabilityType, Count> output{};
   std::transform(std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
   [&](auto&& element)
   {
   return static_cast<ProbabilityType>(element) / sum;
   });
   return output;
   }
   // get_normalized_input template function implementation
   template<class ExPo, class ElementT, std::size_t Count, class ProbabilityType = double>
   requires(std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
   constexpr static auto get_normalized_input(
   ExPo&& execution_policy,
   const std::array<ElementT, Count>& input,
   const ProbabilityType& sum)
   {
   std::array<ProbabilityType, Count> output{};
   std::transform(
   execution_policy, std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
   [&](auto&& element)
   {
   return static_cast<ProbabilityType>(element) / sum;
   });
   return output;
   }
  }
  

• normalize_histogram template function implementation

  namespace TinyDIP
  {
   // normalize_histogram template function implementation for std::array
   template<class ElementT, std::size_t Count, class ProbabilityType = double>
   constexpr static auto normalize_histogram(const std::array<ElementT, Count>& input)
   {
   auto sum = std::reduce(std::ranges::cbegin(input), std::ranges::cend(input));
   return get_normalized_input(input, static_cast<ProbabilityType>(sum));
   }
   // normalize_histogram template function implementation for std::array (with Execution Policy)
   template<class ExecutionPolicy, class ElementT, std::size_t Count, class ProbabilityType = double>
   requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
   constexpr static auto normalize_histogram(ExecutionPolicy execution_policy, const std::array<ElementT, Count>& input)
   {
   auto sum = std::reduce(execution_policy, std::ranges::cbegin(input), std::ranges::cend(input));
   return get_normalized_input(execution_policy, input, static_cast<ProbabilityType>(sum));
   }
   // normalize_histogram template function implementation (with Execution Policy)
   template<class ExecutionPolicy, class ElementT, class CountT, class ProbabilityType = double>
   requires((std::floating_point<CountT> || std::integral<CountT>) and
   (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>))
   constexpr static auto normalize_histogram(ExecutionPolicy execution_policy, const std::map<ElementT, CountT>& input)
   {
   auto sum = sum_second_element(input);
   std::map<ElementT, ProbabilityType> output{};
   for (const auto& [key, value] : input)
   {
   output.emplace(key, static_cast<ProbabilityType>(value) / static_cast<ProbabilityType>(sum));
   }
   return output;
   }
   // normalize_histogram template function implementation
   template<class ElementT, class CountT, class ProbabilityType = double>
   requires(std::floating_point<CountT> || std::integral<CountT>)
   constexpr static auto normalize_histogram(const std::map<ElementT, CountT>& input)
   {
   return normalize_histogram(std::execution::seq, input);
   }
  }
  

• histogram template function implementation

  namespace TinyDIP
  {
   // histogram template function implementation
   // https://codereview.stackexchange.com/q/295419/231235
   template<std::integral ElementT = std::uint8_t>
   requires (std::same_as<ElementT, std::uint8_t> or
   std::same_as<ElementT, std::uint16_t>)
   constexpr static auto histogram(const Image<ElementT>& input)
   {
   std::array<std::size_t, std::numeric_limits<ElementT>::max() - std::numeric_limits<ElementT>::lowest() + 1> histogram_output{};
   auto image_data = input.getImageData();
   for (std::size_t i = 0; i < image_data.size(); ++i)
   {
   ++histogram_output[image_data[i]];
   }
   return histogram_output;
   }
   // histogram template function implementation
   // https://codereview.stackexchange.com/q/295448/231235
   template<class ElementT = int>
   constexpr static auto histogram(const Image<ElementT>& input)
   {
   std::map<ElementT, std::size_t> histogram_output{};
   auto image_data = input.getImageData();
   for (std::size_t i = 0; i < image_data.size(); ++i)
   {
   if (histogram_output.contains(image_data[i]))
   {
   ++histogram_output[image_data[i]];
   }
   else
   {
   histogram_output.emplace(image_data[i], std::size_t{ 1 });
   }
   }
   return histogram_output;
   }
  }
  

• Timer class implementation

  namespace TinyDIP
  {
   class Timer
   {
   private:
   std::chrono::system_clock::time_point start, end;
   std::chrono::duration<double> elapsed_seconds;
   std::time_t end_time;
   public:
   Timer()
   {
   start = std::chrono::system_clock::now();
   }
   ~Timer()
   {
   end = std::chrono::system_clock::now();
   elapsed_seconds = end - start;
   end_time = std::chrono::system_clock::to_time_t(end);
   if (elapsed_seconds.count() != 1)
   {
   std::print(std::cout, "Computation finished at {} elapsed time: {} seconds.\n", std::ctime(&end_time), elapsed_seconds.count());
   }
   else
   {
   std::print(std::cout, "Computation finished at {} elapsed time: {} second.\n", std::ctime(&end_time), elapsed_seconds.count());
   }
   }
   };
  }
  

• apply_threshold_openmp template function implementation

  namespace TinyDIP
  {
   // apply_threshold_openmp template function implementation
   template <typename ElementT>
   constexpr static auto apply_threshold_openmp(const Image<ElementT>& image, const ElementT threshold)
   {
   return apply_each_pixel_openmp(image,
   [&](const ElementT& each_pixel)
   {
   return (each_pixel > threshold) ? std::numeric_limits<ElementT>::max() : std::numeric_limits<ElementT>::lowest();
   });
   }
  }
  

• apply_each_pixel_openmp template function implementation

  namespace TinyDIP
  {
   // apply_each_pixel_openmp template function implementation
   template<class ElementT, class F, class... Args>
   constexpr static auto apply_each_pixel_openmp(const Image<ElementT>& input, F operation, Args&&... args)
   {
   std::vector<std::invoke_result_t<F, ElementT, Args...>> output_vector;
   auto input_count = input.count();
   auto input_vector = input.getImageData();
   output_vector.resize(input_count);
   #pragma omp parallel for
   for (std::size_t i = 0; i < input_count; ++i)
   {
   output_vector[i] = std::invoke(operation, input_vector[i], args...);
   }
   return Image<std::invoke_result_t<F, ElementT, Args...>>(output_vector, input.getSize());
   }
  }
  

The usage of otsu_threshold function:

/* Developed by Jimmy Hu */
#include <chrono>
#include <filesystem>
#include <ostream>
#include "../base_types.h"
#include "../basic_functions.h"
#include "../image.h"
#include "../image_io.h"
#include "../image_operations.h"
#include "../timer.h"
template<class ExPo, class ElementT>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
constexpr static auto otsuThresholdTest(
 ExPo execution_policy,
 const TinyDIP::Image<ElementT>& input,
 std::ostream& os = std::cout
)
{
 auto hsv_image = TinyDIP::rgb2hsv(execution_policy, input);
 TinyDIP::Timer timer1;
 auto unit8_image = TinyDIP::im2uint8(TinyDIP::getVplane(hsv_image));
 return TinyDIP::apply_threshold_openmp(unit8_image, TinyDIP::otsu_threshold(execution_policy, unit8_image));
}
int main(int argc, char* argv[])
{
 std::string image_filename = "../InputImages/1.bmp"; // Default file path
 if (argc == 2) // User has specified input file
 {
 image_filename = std::string(argv[1]);
 }
 if (!std::filesystem::is_regular_file(image_filename))
 {
 throw std::runtime_error("Error: File not found!");
 }
 TinyDIP::Timer timer1;
 auto image_input = TinyDIP::bmp_read(image_filename.c_str(), true);
 image_input = TinyDIP::copyResizeBicubic(image_input, 3 * image_input.getWidth(), 3 * image_input.getHeight());
 auto image_output = otsuThresholdTest(std::execution::par, image_input);
 TinyDIP::bmp_write("test_output", TinyDIP::constructRGB(image_output, image_output, image_output));
 return EXIT_SUCCESS;
}

The output of the test code above:

Width of the input image: 512
Height of the input image: 512
Size of the input image(Byte): 786432
Computation finished at Thu Mar 6 12:36:03 2025
 elapsed time: 0.0301528 seconds.
Computation finished at Thu Mar 6 12:36:03 2025
 elapsed time: 1.9569562 seconds.

The output test_output image:

OutputImage

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++,

  histogram Template Function Implementation for Image in C++,

  Histogram of Image using std::map in C++,

  histogram_normalized and histogram_with_bins Template Functions Implementation for Image in C++ and

  normalize_histogram Template Function Implementation for Image in C++

• What changes has been made in the code since last question?

  I implemented otsu_threshold template function in this post.

• Why a new review is being asked for?

  Please review the implementation of otsu_threshold template function and its tests.

• Cris Luengo: You implement a O(n^2) algorithm, Otsu can famously be computed in O(n). Take a look for example at the Wikipedia page for "Otsu’s method".
• Cris Luengo: Sorry, Wikipedia is not very clear, the MATLAB code is hard to read and the Python code is just dumb (doesn’t even use a histogram?). Other tutorials online seem to be just as bad, showing how to implement a O(n^2) algorithm. I didn’t know the state of the internet these days was so bad! You can take a look at my implementation here: github.com/DIPlib/diplib/blob/…

Answer 1

Avoid code duplication

There's duplicated code in otsu_threshold() because normalize_histogram() returns completely different types depending on ElementT. There are several ways to avoid this:

• Create a class Histogram that hides how the histogram is actually stored internally, and provides a uniform way to access that data.
• Always create a std::vector out of the returned histogram, since the line that declared probabilityVec would work just as well if probabilities was a std::array, at least if you let it autodeduce the value type:

  auto probabilityVec = std::vector(std::ranges::begin(probabilities), std::ranges::end(probabilities));
  

  Or even simpler:

  auto probabilityVec = probabilities | std::ranges::to<std::vector>();
  

  To ensure you get an index when applying this on a std::array, you could make your own adapter:

  auto probabilityVec = probabilities | add_index_if_missing | std::ranges::to<std::vector>();
  

  Where add_index_if_missing is a range adapter you have to write that applies std::views::enumerate if necessary.
• Realizing the above, note that instead of iterating over the histogram using an integer index and [], you can use actual iterators. This means you'll have to write code like this:

  for (auto i = std::begin(probabilities); i != std::end(probabilities); ++i) {
   ...
   for (auto j = std::begin(probabilities); j != i; ++j) {
   w_background += j->second;
   m_background += j->first * j->second;
   }
   ...
  }
  

What about non-integral element types?

Your function only deals with images with integer pixel values. However, what if ElementT is float? If you think that floating point images would result in huge histograms, please realize that std::uint32_t has just as many possible values as a float does. I think your code will handle floats just fine if you just remove the constraint.

But this also brings me to the following:

Consider that the histogram might be huge

Your algorithm is \$O(N^2)\$ where \$N\$ is the size of the histogram. For 8-bit images, that's totally fine. For 16-bit images, it will execute the code in the inner loops \2ドル^32\$ times, regardless of the image size.

You can save some time by realizing that you don't need to recompute all of w_background/w_foreground/etc each time, as only one element is added/removed from the sums each iteration of the outer loop.

For 32-bit images, your code might actually be more or less efficient, depending on the number of distinct pixel values. Still, you must consider that the size of the histogram might be equal to the total number of pixels in the image, if each pixel has a unique value.

There several possible approaches you can take for dealing with huge histograms:

• Just limit the number of histogram bins. This might result in a slightly optimal threshold value.
• Use a heuristic or bisection to find the right value in less amount of time.