Return to Revisions

4 of 5

Update code

edited Mar 6 at 6:50

7.6k
2
11
48

otsu_threshold Template Function Implementation for Image in C++

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:

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.

c++performance algorithm image c++23

asked Mar 6 at 4:41

JimmyHu

7.6k
2
11
48