2
\$\begingroup\$

This is a follow-up question for apply_each_single_output Template Function Implementation for Image in C++. Considering the alternative approach proposed by G. Sliepen, I am trying to implement pow and sqrt for multi-channel values:

//  pow Template Function Implementation
template<std::size_t channel_count = 3, typename T, typename ExpT = double>
constexpr static auto pow(const T& input, const ExpT exp)
{
    if constexpr (Multichannel<T>)
    {
        MultiChannel<decltype(std::pow(input.channels[0], exp)), channel_count> output;
        for (std::size_t i = 0; i < channel_count; ++i)
        {
            output.channels[i] = std::pow(input.channels[i], exp);
        }
        return output;
    }
    else
    {
        return std::pow(input);
    }
}

//  sqrt Template Function Implementation
template<std::size_t channel_count = 3, typename T>
constexpr static auto sqrt(const T& input)
{
    if constexpr (Multichannel<T>)
    {
        MultiChannel<decltype(std::sqrt(input.channels[0])), channel_count> output;
        for (std::size_t i = 0; i < channel_count; ++i)
        {
            output.channels[i] = std::sqrt(input.channels[i]);
        }
        return output;
    }
    else
    {
        return std::sqrt(input);
    }
}

It seems to be working fine. However, the structure in if constexpr block is duplicated. Another function apply_multichannel comes up in my mind:

//  apply_multichannel Template Function Implementation
template<std::size_t channel_count = 3, class ElementT, class Lambda, typename... Args>
[[nodiscard]] constexpr static auto apply_multichannel(const MultiChannel<ElementT, channel_count>& input, Lambda f, Args... args)
{
    MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
    std::transform(std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
        [&](auto&& input) { return std::invoke(f, input, args...); });
    return output;
}

//  apply_multichannel Template Function Implementation
template<std::size_t channel_count = 3, class T, class Lambda, typename... Args>
requires((std::same_as<T, RGB>) || (std::same_as<T, RGB_DOUBLE>) || (std::same_as<T, HSV>))
[[nodiscard]] constexpr static auto apply_multichannel(const T& input, Lambda f, Args... args)
{
    MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
    std::transform(std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
        [&](auto&& input) { return std::invoke(f, input, args...); });
    return output;
}

The version with execution policy:

//  apply_multichannel Template Function Implementation (the version with execution policy)
template<std::size_t channel_count = 3, class ExecutionPolicy, class ElementT, class Lambda, typename... Args>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
[[nodiscard]] constexpr static auto apply_multichannel(ExecutionPolicy&& execution_policy, const MultiChannel<ElementT, channel_count>& input, Lambda f, Args... args)
{
    MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
    std::transform(execution_policy, std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
        [&](auto&& input) { return std::invoke(f, input, args...); });
    return output;
}

//  apply_multichannel Template Function Implementation (the version with execution policy)
template<std::size_t channel_count = 3, class ExecutionPolicy, class T, class Lambda, typename... Args>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)&&
         ((std::same_as<T, RGB>) || (std::same_as<T, RGB_DOUBLE>) || (std::same_as<T, HSV>))
[[nodiscard]] constexpr static auto apply_multichannel(ExecutionPolicy&& execution_policy, const T& input, Lambda f, Args... args)
{
    MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
    std::transform(execution_policy, std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
        [&](auto&& input) { return std::invoke(f, input, args...); });
    return output;
}

Then, the implementation of pow and sqrt template functions can be simpler:

//  pow Template Function Implementation
template<typename T, typename ExpT = double>
[[nodiscard]] constexpr static auto pow(const T& input, const ExpT exp)
{
    if constexpr (Multichannel<T>)
    {
        return apply_multichannel(input, [&](auto&& _input, auto&& input_exp) {return std::pow(_input, input_exp); }, exp);
    }
    else
    {
        return std::pow(input);
    }
}
//  sqrt Template Function Implementation
template<typename T>
[[nodiscard]] constexpr static auto sqrt(const T& input)
{
    if constexpr (Multichannel<T>)
    {
        return apply_multichannel(input, [&](auto&& _input) {return std::sqrt(_input); });
    }
    else
    {
        return std::sqrt(input);
    }
}

pow and sqrt template functions with execution policy parameter:

//  pow Template Function Implementation (the version with execution policy)
template<class ExecutionPolicy, typename T, typename ExpT = double>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
[[nodiscard]] constexpr static auto pow(ExecutionPolicy&& execution_policy, const T& input, const ExpT exp)
{
    if constexpr (Multichannel<T>)
    {
        return apply_multichannel(execution_policy, input, [&](auto&& _input, auto&& input_exp) {return std::pow(_input, input_exp); }, exp);
    }
    else
    {
        return std::pow(input);
    }
}

//  sqrt Template Function Implementation (the version with execution policy)
template<class ExecutionPolicy, typename T>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
[[nodiscard]] constexpr static auto sqrt(ExecutionPolicy&& execution_policy, const T& input)
{
    if constexpr (Multichannel<T>)
    {
        return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::sqrt(_input); });
    }
    else
    {
        return std::sqrt(input);
    }
}

Full Testing Code

The full testing code:

//  apply_multichannel Template Function Implementation in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <concepts>
#include <execution>
#include <filesystem>
#include <fstream>
#include <functional>
#include <future>
#include <iostream>
#include <iterator>
#include <numeric>
#include <ranges>
#include <random>
#include <string>
#include <type_traits>
#include <variant>
#include <vector>
#include <utility>

namespace TinyDIP
{
    struct RGB
    {
        std::uint8_t channels[3];

        inline RGB operator+(const RGB& input) const
        {
            return RGB{
                static_cast<std::uint8_t>(input.channels[0] + channels[0]),
                static_cast<std::uint8_t>(input.channels[1] + channels[1]),
                static_cast<std::uint8_t>(input.channels[2] + channels[2]) };
        }

        inline RGB operator-(const RGB& input) const
        {
            return RGB{
                static_cast<std::uint8_t>(channels[0] - input.channels[0]),
                static_cast<std::uint8_t>(channels[1] - input.channels[1]),
                static_cast<std::uint8_t>(channels[2] - input.channels[2]) };
        }

        friend std::ostream& operator<<(std::ostream& out, const RGB& _myStruct)
        {
            out << '{' << +_myStruct.channels[0] << ", " << +_myStruct.channels[1] << ", " << +_myStruct.channels[2] << '}';
            return out;
        }
    };

    struct RGB_DOUBLE
    {
        double channels[3];

        inline RGB_DOUBLE operator+(const RGB_DOUBLE& input) const
        {
            return RGB_DOUBLE{
                input.channels[0] + channels[0],
                input.channels[1] + channels[1],
                input.channels[2] + channels[2] };
        }

        inline RGB_DOUBLE operator-(const RGB_DOUBLE& input) const
        {
            return RGB_DOUBLE{
                channels[0] - input.channels[0],
                channels[1] - input.channels[1],
                channels[2] - input.channels[2] };
        }

        friend std::ostream& operator<<(std::ostream& out, const RGB_DOUBLE& _myStruct)
        {
            out << '{' << +_myStruct.channels[0] << ", " << +_myStruct.channels[1] << ", " << +_myStruct.channels[2] << '}';
            return out;
        }
    };

    using GrayScale = std::uint8_t;

    struct HSV
    {
        double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255

        inline HSV operator+(const HSV& input) const
        {
            return HSV{
                input.channels[0] + channels[0],
                input.channels[1] + channels[1],
                input.channels[2] + channels[2] };
        }

        inline HSV operator-(const HSV& input) const
        {
            return HSV{
                channels[0] - input.channels[0],
                channels[1] - input.channels[1],
                channels[2] - input.channels[2] };
        }

        friend std::ostream& operator<<(std::ostream& out, const HSV& _myStruct)
        {
            out << '{' << +_myStruct.channels[0] << ", " << +_myStruct.channels[1] << ", " << +_myStruct.channels[2] << '}';
            return out;
        }
    };

    //  MultiChannel struct implementation
    template<class ElementT, std::size_t channel_count = 3>
    struct MultiChannel
    {
        std::array<ElementT, channel_count> channels;

        inline MultiChannel operator+(const MultiChannel& input) const
        {
            std::array<ElementT, channel_count> channels_output;
            for(std::size_t i = 0; i < channels.size(); ++i)
            {
                channels_output[i] = channels[i] + input.channels[i];
            }
            return MultiChannel{channels_output};
        }

        inline MultiChannel operator-(const MultiChannel& input) const
        {
            std::array<ElementT, channel_count> channels_output;
            for(std::size_t i = 0; i < channels.size(); ++i)
            {
                channels_output[i] = channels[i] - input.channels[i];
            }
            return MultiChannel{channels_output};
        }

        friend std::ostream& operator<<(std::ostream& out, const MultiChannel& _myStruct)
        {
            out << '{';
            for(std::size_t i = 0; i < channel_count; ++i)
            {
                out << +_myStruct.channels[i] << ", ";
            }
            out << '}';
            return out;
        }
    };

    struct BMPIMAGE
    {
        std::filesystem::path FILENAME;

        unsigned int XSIZE;
        unsigned int YSIZE;
        std::uint8_t FILLINGBYTE;
        std::uint8_t* IMAGE_DATA;
    };
    
    //  Reference: https://stackoverflow.com/a/48458312/6667035
    template <typename>
    struct is_tuple : std::false_type {};

    template <typename ...T>
    struct is_tuple<std::tuple<T...>> : std::true_type {};

    //  is_MultiChannel struct implementation
    template <typename>
    struct is_MultiChannel : std::false_type {};

    template <std::size_t N, typename T>
    struct is_MultiChannel<MultiChannel<T, N>> : std::true_type {};

    template <typename, typename>
    struct check_tuple_element_type {};

    template <typename TargetType, typename ...ElementT>
    struct check_tuple_element_type<TargetType, std::tuple<ElementT...>>
        : std::bool_constant<(std::is_same_v<ElementT, TargetType> || ...)>
    {
    };

    template<typename T>
    concept image_element_standard_floating_point_type =
        std::same_as<double, T>
        or std::same_as<float, T>
        or std::same_as<long double, T>
    ;

    //  Reference: https://stackoverflow.com/a/64287611/6667035
    template <typename T>
    struct is_complex : std::false_type {};

    template <typename T>
    struct is_complex<std::complex<T>> : std::true_type {};

    //  Reference: https://stackoverflow.com/a/58067611/6667035
    template <typename T>
    concept arithmetic = std::is_arithmetic_v<T> or is_complex<T>::value;

    //  recursive_print template function implementation
    template<typename T>
    constexpr void recursive_print(const T& input, const std::size_t level = 0)
    {
        std::cout << std::string(level, ' ') << input << '\n';
    }

    template<std::ranges::input_range Range>
    constexpr void recursive_print(const Range& input, const std::size_t level = 0)
    {
        std::cout << std::string(level, ' ') << "Level " << level << ":" << std::endl;
        std::ranges::for_each(input, [level](auto&& element) {
            recursive_print(element, level + 1);
        });
    }

    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        template<std::same_as<std::size_t>... Sizes>
        Image(Sizes... sizes): size{sizes...}, image_data((1 * ... * sizes))
        {}

        template<std::same_as<int>... Sizes>
        Image(Sizes... sizes)
        {
            size.reserve(sizeof...(sizes));
            (size.emplace_back(sizes), ...);
            image_data.resize(
                std::reduce(
                    std::ranges::cbegin(size),
                    std::ranges::cend(size),
                    std::size_t{1},
                    std::multiplies<>()
                    )
            );
        }

        Image(const std::vector<std::size_t>& sizes)
        {
            if (sizes.empty())
            {
                throw std::runtime_error("Image size vector is empty!");
            }
            size = std::move(sizes);
            image_data.resize(
                std::reduce(
                    std::ranges::cbegin(sizes),
                    std::ranges::cend(sizes),
                    std::size_t{1},
                    std::multiplies<>()
                    ));
        }

        template<std::ranges::input_range Range,
                 std::same_as<std::size_t>... Sizes>
        Image(const Range& input, Sizes... sizes):
            size{sizes...}, image_data(begin(input), end(input))
        {
            if (image_data.size() != (1 * ... * sizes)) {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
        }

        template<std::same_as<std::size_t>... Sizes>
        Image(std::vector<ElementT>&& input, Sizes... sizes):
            size{sizes...}, image_data(begin(input), end(input))
        {
            if (input.empty())
            {
                throw std::runtime_error("Input vector is empty!");
            }
            if (image_data.size() != (1 * ... * sizes)) {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
        }

        Image(const std::vector<ElementT>& input, const std::vector<std::size_t>& sizes)
        {
            if (input.empty())
            {
                throw std::runtime_error("Input vector is empty!");
            }
            size = std::move(sizes);
            image_data = std::move(input);
            auto count = std::reduce(std::ranges::cbegin(sizes), std::ranges::cend(sizes), 1, std::multiplies());
            if (image_data.size() != count) {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
        }


        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            if (input.empty())
            {
                throw std::runtime_error("Input vector is empty!");
            }
            size.reserve(2);
            size.emplace_back(newWidth);
            size.emplace_back(newHeight);
            if (input.size() != newWidth * newHeight)
            {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
            image_data = std::move(input);              //  Reference: https://stackoverflow.com/a/51706522/6667035
        }

        Image(const std::vector<std::vector<ElementT>>& input)
        {
            if (input.empty())
            {
                throw std::runtime_error("Input vector is empty!");
            }
            size.reserve(2);
            size.emplace_back(input[0].size());
            size.emplace_back(input.size());
            for (auto& rows : input)
            {
                image_data.insert(image_data.end(), std::ranges::begin(rows), std::ranges::end(rows));    //  flatten
            }
            return;
        }

        //  at template function implementation
        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            return const_cast<ElementT&>(static_cast<const Image &>(*this).at(indexInput...));
        }

        //  at template function implementation
        //  Reference: https://codereview.stackexchange.com/a/288736/231235
        template<typename... Args>
        constexpr ElementT const& at(const Args... indexInput) const
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            std::size_t i = 0;
            std::size_t stride = 1;
            std::size_t position = 0;

            auto update_position = [&](auto index) {
                position += index * stride;
                stride *= size[i++];
            };
            (update_position(indexInput), ...);

            return image_data[position];
        }

        //  at_without_boundary_check template function implementation
        template<typename... Args>
        constexpr ElementT& at_without_boundary_check(const Args... indexInput)
        {
            return const_cast<ElementT&>(static_cast<const Image &>(*this).at_without_boundary_check(indexInput...));
        }

        template<typename... Args>
        constexpr ElementT const& at_without_boundary_check(const Args... indexInput) const
        {
            std::size_t i = 0;
            std::size_t stride = 1;
            std::size_t position = 0;

            auto update_position = [&](auto index) {
                position += index * stride;
                stride *= size[i++];
            };
            (update_position(indexInput), ...);

            return image_data[position];
        }

        //  get function implementation
        constexpr ElementT get(std::size_t index) const noexcept
        {
            return image_data[index];
        }

        //  set function implementation
        constexpr ElementT& set(const std::size_t index) noexcept
        {
            return image_data[index];
        }

        //  set template function implementation
        template<class TupleT>
        requires(is_tuple<TupleT>::value and
                 check_tuple_element_type<std::size_t, TupleT>::value)
        constexpr bool set(const TupleT location, const ElementT draw_value)
        {
            if (checkBoundaryTuple(location))
            {
                image_data[tuple_location_to_index(location)] = draw_value;
                return true;
            }
            return false;
        }

        //  cast template function implementation
        template<typename TargetT>
        constexpr Image<TargetT> cast()
        {
            std::vector<TargetT> output_data;
            output_data.resize(image_data.size());
            std::transform(
                std::ranges::cbegin(image_data),
                std::ranges::cend(image_data),
                std::ranges::begin(output_data),
                [](auto& input){ return static_cast<TargetT>(input); }
                );
            Image<TargetT> output(output_data, size);
            return output;
        }

        constexpr std::size_t count() const noexcept
        {
            return std::reduce(std::ranges::cbegin(size), std::ranges::cend(size), 1, std::multiplies());
        }
  
        constexpr std::size_t getDimensionality() const noexcept
        {
            return size.size();
        }

        constexpr std::size_t getWidth() const noexcept
        {
            return size[0];
        }

        constexpr std::size_t getHeight() const noexcept
        {
            return size[1];
        }

        //  getSize function implementation
        constexpr auto getSize() const noexcept
        {
            return size;
        }

        //  getSize function implementation
        constexpr auto getSize(std::size_t index) const noexcept
        {
            return size[index];
        }

        //  getStride function implementation
        constexpr std::size_t getStride(std::size_t index) const noexcept
        {
            if(index == 0)
            {
                return std::size_t{1};
            }
            std::size_t output = std::size_t{1};
            for(std::size_t i = 0; i < index; ++i)
            {
                output *= size[i];
            }
            return output;
        }

        std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data

        //  print function implementation
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            if constexpr (is_MultiChannel<ElementT>::value)
            {
                if (size.size() == 1)
                {
                    for (std::size_t x = 0; x < size[0]; ++x)
                    {
                        os << at(x) << separator;
                    }
                    os << "\n";
                }
                else if (size.size() == 2)
                {
                    for (std::size_t y = 0; y < size[1]; ++y)
                    {
                        for (std::size_t x = 0; x < size[0]; ++x)
                        {
                            os << at(x, y) << separator;
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
                else if (size.size() == 3)
                {
                    for (std::size_t z = 0; z < size[2]; ++z)
                    {
                        for (std::size_t y = 0; y < size[1]; ++y)
                        {
                            for (std::size_t x = 0; x < size[0]; ++x)
                            {
                                os << at(x, y, z) << separator;
                            }
                            os << "\n";
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
                else if (size.size() == 4)
                {
                    for (std::size_t a = 0; a < size[3]; ++a)
                    {
                        os << "group = " << a << "\n";
                        for (std::size_t z = 0; z < size[2]; ++z)
                        {
                            for (std::size_t y = 0; y < size[1]; ++y)
                            {
                                for (std::size_t x = 0; x < size[0]; ++x)
                                {
                                    os << at(x, y, z, a) << separator;
                                }
                                os << "\n";
                            }
                            os << "\n";
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
            }
            else
            {
                if (size.size() == 1)
                {
                    for (std::size_t x = 0; x < size[0]; ++x)
                    {
                        //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                        os << +at(x) << separator;
                    }
                    os << "\n";
                }
                else if (size.size() == 2)
                {
                    for (std::size_t y = 0; y < size[1]; ++y)
                    {
                        for (std::size_t x = 0; x < size[0]; ++x)
                        {
                            //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                            os << +at(x, y) << separator;
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
                else if (size.size() == 3)
                {
                    for (std::size_t z = 0; z < size[2]; ++z)
                    {
                        for (std::size_t y = 0; y < size[1]; ++y)
                        {
                            for (std::size_t x = 0; x < size[0]; ++x)
                            {
                                //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                                os << +at(x, y, z) << separator;
                            }
                            os << "\n";
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
                else if (size.size() == 4)
                {
                    for (std::size_t a = 0; a < size[3]; ++a)
                    {
                        os << "group = " << a << "\n";
                        for (std::size_t z = 0; z < size[2]; ++z)
                        {
                            for (std::size_t y = 0; y < size[1]; ++y)
                            {
                                for (std::size_t x = 0; x < size[0]; ++x)
                                {
                                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                                    os << +at(x, y, z, a) << separator;
                                }
                                os << "\n";
                            }
                            os << "\n";
                        }
                        os << "\n";
                    }
                    os << "\n";
                }
            }
        }

        //  Enable this function if ElementT = RGB or RGB_DOUBLE or HSV
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB> or std::same_as<ElementT, RGB_DOUBLE> or std::same_as<ElementT, HSV> or is_MultiChannel<ElementT>::value)
        {
            if (size.size() == 1)
            {
                for (std::size_t x = 0; x < size[0]; ++x)
                {
                    os << "( ";
                    for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                    {
                        //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                        os << +at(x).channels[channel_index] << separator;
                    }
                    os << ")" << separator;
                }
                os << "\n";
            }
            else if (size.size() == 2)
            {
                for (std::size_t y = 0; y < size[1]; ++y)
                {
                    for (std::size_t x = 0; x < size[0]; ++x)
                    {
                        os << "( ";
                        for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                        {
                            //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                            os << +at(x, y).channels[channel_index] << separator;
                        }
                        os << ")" << separator;
                    }
                    os << "\n";
                }
                os << "\n";
            }
            return;
        }

        Image<ElementT>& setAllValue(const ElementT input)
        {
            std::fill(std::ranges::begin(image_data), std::ranges::end(image_data), input);
            return *this;
        }

        friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
        {
            const std::string separator = "\t";
            rhs.print(separator, os);
            return os;
        }

        Image<ElementT>& operator+=(const Image<ElementT>& rhs)
        {
            check_size_same(rhs, *this);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::plus<>{});
            return *this;
        }

        Image<ElementT>& operator-=(const Image<ElementT>& rhs)
        {
            check_size_same(rhs, *this);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::minus<>{});
            return *this;
        }

        Image<ElementT>& operator*=(const Image<ElementT>& rhs)
        {
            check_size_same(rhs, *this);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::multiplies<>{});
            return *this;
        }

        Image<ElementT>& operator/=(const Image<ElementT>& rhs)
        {
            check_size_same(rhs, *this);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::divides<>{});
            return *this;
        }

        friend bool operator==(Image<ElementT> const&, Image<ElementT> const&) = default;

        friend bool operator!=(Image<ElementT> const&, Image<ElementT> const&) = default;

        friend Image<ElementT> operator+(Image<ElementT> input1, const Image<ElementT>& input2)
        {
            return input1 += input2;
        }

        friend Image<ElementT> operator-(Image<ElementT> input1, const Image<ElementT>& input2)
        {
            return input1 -= input2;
        }

        friend Image<ElementT> operator*(Image<ElementT> input1, ElementT input2)
        {
            return multiplies(input1, input2);
        }

        friend Image<ElementT> operator*(ElementT input1, Image<ElementT> input2)
        {
            return multiplies(input2, input1);
        }
        
#ifdef USE_BOOST_SERIALIZATION

        void Save(std::string filename)
        {
            const std::string filename_with_extension = filename + ".dat";
            //  Reference: https://stackoverflow.com/questions/523872/how-do-you-serialize-an-object-in-c
            std::ofstream ofs(filename_with_extension, std::ios::binary);
            boost::archive::binary_oarchive ArchiveOut(ofs);
            //  write class instance to archive
            ArchiveOut << *this;
            //  archive and stream closed when destructors are called
            ofs.close();
        }
        
#endif
    private:
        std::vector<std::size_t> size;
        std::vector<ElementT> image_data;

        template<typename... Args>
        void checkBoundary(const Args... indexInput) const
        {
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            std::size_t parameter_pack_index = 0;
            auto function = [&](auto index) {
                if (index >= size[parameter_pack_index])
                    throw std::out_of_range("Given index out of range!");
                parameter_pack_index = parameter_pack_index + 1;
            };

            (function(indexInput), ...);
        }

        //  checkBoundaryTuple template function implementation
        template<class TupleT>
        requires(TinyDIP::is_tuple<TupleT>::value)
        constexpr bool checkBoundaryTuple(const TupleT location)
        {
            constexpr std::size_t n = std::tuple_size<TupleT>{};
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            std::size_t parameter_pack_index = 0;
            auto function = [&](auto index) {
                if (std::cmp_greater_equal(index, size[parameter_pack_index]))
                    return false;
                parameter_pack_index = parameter_pack_index + 1;
                return true;
            };
            return std::apply([&](auto&&... args) { return ((function(args))&& ...);}, location);
        }

        //  tuple_location_to_index template function implementation
        template<class TupleT>
        requires(TinyDIP::is_tuple<TupleT>::value)
        constexpr std::size_t tuple_location_to_index(TupleT location)
        {
            std::size_t i = 0;
            std::size_t stride = 1;
            std::size_t position = 0;
            auto update_position = [&](auto index) {
                    position += index * stride;
                    stride *= size[i++];
                };
            std::apply([&](auto&&... args) {((update_position(args)), ...);}, location);
            return position;
        }

#ifdef USE_BOOST_SERIALIZATION
        friend class boost::serialization::access;
        template<class Archive>
        void serialize(Archive& ar, const unsigned int version)
        {
            ar& size;
            ar& image_data;
        }
#endif

    };

    template<typename T, typename ElementT>
    concept is_Image = std::is_same_v<T, Image<ElementT>>;

    //  zeros template function implementation
    template<typename ElementT, std::same_as<std::size_t>... Sizes>
    constexpr static auto zeros(Sizes... sizes)
    {
        auto output = Image<ElementT>(sizes...);
        return output;
    }

    //  ones template function implementation
    template<typename ElementT, std::same_as<std::size_t>... Sizes>
    constexpr static auto ones(Sizes... sizes)
    {
        auto output = zeros<ElementT>(sizes...);
        output.setAllValue(1);
        return output;
    }

    //  rand template function implementation
    template<image_element_standard_floating_point_type ElementT = double, typename Urbg, std::same_as<std::size_t>... Sizes>
    requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>>
    constexpr static auto rand(Urbg&& urbg, Sizes... sizes)
    {
        if constexpr (sizeof...(Sizes) == 1)
        {
            return rand(std::forward<Urbg>(urbg), sizes..., sizes...);
        }
        else
        {
            std::vector<ElementT> image_data((... * sizes));
            //  Reference: https://stackoverflow.com/a/23143753/6667035
            //  Reference: https://codereview.stackexchange.com/a/294739/231235
            auto dist = std::uniform_real_distribution<ElementT>{};
            std::ranges::generate(image_data, [&dist, &urbg]() { return dist(urbg); });

            return Image<ElementT>{std::move(image_data), sizes...};
        }
    }

    //  rand template function implementation
    template<image_element_standard_floating_point_type ElementT = double, std::same_as<std::size_t>... Size>
    inline auto rand(Size... size)
    {
        return rand<ElementT>(std::mt19937{std::random_device{}()}, size...);
    }

    //  rand template function implementation
    template<image_element_standard_floating_point_type ElementT = double, typename Urbg>
    requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>>
    constexpr auto rand(Urbg&& urbg) -> ElementT
    {
        auto dist = std::uniform_real_distribution<ElementT>{};
        return dist(urbg);
    }

    //  rand template function implementation
    template<image_element_standard_floating_point_type ElementT = double>
    inline auto rand()
    {
        return rand<ElementT>(std::mt19937{std::random_device{}()});
    }

    template<typename ElementT>
    constexpr void check_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        if (!is_width_same(x, y))
            throw std::runtime_error("Width mismatched!");
    }

    template<typename ElementT>
    constexpr void check_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        if (!is_height_same(x, y))
            throw std::runtime_error("Height mismatched!");
    }

    //  check_size_same template function implementation
    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        if(x.getSize() != y.getSize())
            throw std::runtime_error("Size mismatched!");
    }

    //  Multichannel Concept Implementation
    template<typename T>
    concept Multichannel = requires(T a)
    {
        { a.channels }; // or whatever is best to check for multiple channels
    };

    //  apply_multichannel Template Function Implementation
    template<std::size_t channel_count = 3, class ElementT, class Lambda, typename... Args>
    [[nodiscard]] constexpr static auto apply_multichannel(const MultiChannel<ElementT, channel_count>& input, Lambda f, Args... args)
    {
        MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
        std::transform(std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
            [&](auto&& input) { return std::invoke(f, input, args...); });
        return output;
    }

    //  apply_multichannel Template Function Implementation (the version with execution policy)
    template<std::size_t channel_count = 3, class ExecutionPolicy, class ElementT, class Lambda, typename... Args>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto apply_multichannel(ExecutionPolicy&& execution_policy, const MultiChannel<ElementT, channel_count>& input, Lambda f, Args... args)
    {
        MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
        std::transform(execution_policy, std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
            [&](auto&& input) { return std::invoke(f, input, args...); });
        return output;
    }

    //  apply_multichannel Template Function Implementation
    template<std::size_t channel_count = 3, class T, class Lambda, typename... Args>
    requires((std::same_as<T, RGB>) || (std::same_as<T, RGB_DOUBLE>) || (std::same_as<T, HSV>))
    [[nodiscard]] constexpr static auto apply_multichannel(const T& input, Lambda f, Args... args)
    {
        MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
        std::transform(std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
            [&](auto&& input) { return std::invoke(f, input, args...); });
        return output;
    }

    //  apply_multichannel Template Function Implementation (the version with execution policy)
    template<std::size_t channel_count = 3, class ExecutionPolicy, class T, class Lambda, typename... Args>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)&&
            ((std::same_as<T, RGB>) || (std::same_as<T, RGB_DOUBLE>) || (std::same_as<T, HSV>))
    [[nodiscard]] constexpr static auto apply_multichannel(ExecutionPolicy&& execution_policy, const T& input, Lambda f, Args... args)
    {
        MultiChannel<decltype(std::invoke(f, input.channels[0], args...)), channel_count> output;
        std::transform(execution_policy, std::ranges::cbegin(input.channels), std::ranges::cend(input.channels), std::ranges::begin(output.channels),
            [&](auto&& input) { return std::invoke(f, input, args...); });
        return output;
    }

    //  abs Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto abs(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::abs(_input); });
        }
        else
        {
            return std::abs(input);
        }
    }

    //  abs Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto abs(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::abs(_input); });
        }
        else
        {
            return std::abs(input);
        }
    }

    //  pow Template Function Implementation
    template<typename T, typename ExpT = double>
    [[nodiscard]] constexpr static auto pow(const T& input, const ExpT exp)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input, auto&& input_exp) {return std::pow(_input, input_exp); }, exp);
        }
        else
        {
            return std::pow(input);
        }
    }

    //  pow Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T, typename ExpT = double>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto pow(ExecutionPolicy&& execution_policy, const T& input, const ExpT exp)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input, auto&& input_exp) {return std::pow(_input, input_exp); }, exp);
        }
        else
        {
            return std::pow(input);
        }
    }

    //  sqrt Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto sqrt(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::sqrt(_input); });
        }
        else
        {
            return std::sqrt(input);
        }
    }

    //  sqrt Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto sqrt(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::sqrt(_input); });
        }
        else
        {
            return std::sqrt(input);
        }
    }

    //  cbrt Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto cbrt(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::cbrt(_input); });
        }
        else
        {
            return std::cbrt(input);
        }
    }

    //  cbrt Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto cbrt(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::cbrt(_input); });
        }
        else
        {
            return std::cbrt(input);
        }
    }

    //  sin Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto sin(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::sin(_input); });
        }
        else
        {
            return std::sin(input);
        }
    }

    //  sin Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto sin(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::sin(_input); });
        }
        else
        {
            return std::sin(input);
        }
    }

    //  cos Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto cos(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::cos(_input); });
        }
        else
        {
            return std::cos(input);
        }
    }

    //  cos Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto cos(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::cos(_input); });
        }
        else
        {
            return std::cos(input);
        }
    }

    //  tan Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto tan(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::tan(_input); });
        }
        else
        {
            return std::tan(input);
        }
    }

    //  tan Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto tan(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::tan(_input); });
        }
        else
        {
            return std::tan(input);
        }
    }

    //  cot Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto cot(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return 1 / std::tan(_input); });
        }
        else
        {
            return 1 / std::tan(input);
        }
    }

    //  cot Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto cot(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return 1 / std::tan(_input); });
        }
        else
        {
            return 1 / std::tan(input);
        }
    }

    //  sec Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto sec(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return 1 / std::cos(_input); });
        }
        else
        {
            return 1 / std::cos(input);
        }
    }

    //  sec Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto sec(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return 1 / std::cos(_input); });
        }
        else
        {
            return 1 / std::cos(input);
        }
    }

    //  csc Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto csc(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return 1 / std::sin(_input); });
        }
        else
        {
            return 1 / std::sin(input);
        }
    }

    //  csc Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto csc(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return 1 / std::sin(_input); });
        }
        else
        {
            return 1 / std::sin(input);
        }
    }

    //  asin Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto asin(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::asin(_input); });
        }
        else
        {
            return std::asin(input);
        }
    }

    //  asin Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto asin(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::asin(_input); });
        }
        else
        {
            return std::asin(input);
        }
    }

    //  acos Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto acos(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::acos(_input); });
        }
        else
        {
            return std::acos(input);
        }
    }

    //  acos Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto acos(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::acos(_input); });
        }
        else
        {
            return std::acos(input);
        }
    }

    //  atan Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto atan(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::atan(_input); });
        }
        else
        {
            return std::atan(input);
        }
    }

    //  atan Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto atan(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::atan(_input); });
        }
        else
        {
            return std::atan(input);
        }
    }

    //  sinh Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto sinh(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::sinh(_input); });
        }
        else
        {
            return std::sinh(input);
        }
    }

    //  sinh Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto sinh(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::sinh(_input); });
        }
        else
        {
            return std::sinh(input);
        }
    }

    //  cosh Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto cosh(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::cosh(_input); });
        }
        else
        {
            return std::cosh(input);
        }
    }

    //  cosh Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto cosh(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::cosh(_input); });
        }
        else
        {
            return std::cosh(input);
        }
    }

    //  tanh Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto tanh(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::tanh(_input); });
        }
        else
        {
            return std::tanh(input);
        }
    }

    //  tanh Template Function Implementation (the version with execution policy)
    template<class ExecutionPolicy, typename T>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    [[nodiscard]] constexpr static auto tanh(ExecutionPolicy&& execution_policy, const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(execution_policy, input, [&](auto&& _input) {return std::tanh(_input); });
        }
        else
        {
            return std::tanh(input);
        }
    }

    //  asinh Template Function Implementation
    template<typename T>
    [[nodiscard]] constexpr static auto asinh(const T& input)
    {
        if constexpr (Multichannel<T>)
        {
            return apply_multichannel(input, [&](auto&& _input) {return std::asinh(_input); });
        }
        else
        {
            return std::asinh(input);
        }
    }

    //  two_input_map_reduce Template Function Implementation
    template<
        class ExecutionPolicy,
        std::ranges::input_range Input1,
        std::ranges::input_range Input2,
        class T,
        class BinaryOp1 = std::minus<T>,
        class BinaryOp2 = std::plus<T>
    >
    requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
    constexpr auto two_input_map_reduce(
        ExecutionPolicy execution_policy,
        const Input1& input1,
        const Input2& input2,
        const T init = {},
        const BinaryOp1& binop1 = std::minus<T>(),
        const BinaryOp2& binop2 = std::plus<T>())
    {
        if (input1.size() != input2.size())
        {
            throw std::runtime_error("Size mismatched!");
        }
        auto transformed = std::views::zip(input1, input2)
                     | std::views::transform([&](auto input) {
                           return std::invoke(binop1, std::get<0>(input), std::get<1>(input));
                       });
        return std::reduce(
            execution_policy,
            transformed.begin(), transformed.end(),
            init,
            binop2
        );
    }

    //  euclidean_distance Template Function Implementation
    template<
        arithmetic OutputT = double,
        class ExPo,
        arithmetic ElementT1 = double,
        arithmetic ElementT2 = double
        >
    requires(std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static auto euclidean_distance(
        ExPo execution_policy,
        const Image<ElementT1>& input1,
        const Image<ElementT2>& input2,
        const OutputT output = 0.0
        )
    {
        if (input1.getSize() != input2.getSize())
        {
            throw std::runtime_error("Size mismatched!");
        }
        return std::sqrt(two_input_map_reduce(execution_policy, input1.getImageData(), input2.getImageData(), OutputT{},
            [&](auto&& element1, auto&& element2) {
                return std::pow(element1 - element2, 2.0);
            }));
    }

    //  euclidean_distance Template Function Implementation
    template<
        arithmetic OutputT = double,
        arithmetic ElementT1 = double,
        arithmetic ElementT2 = double
        >
    constexpr static auto euclidean_distance(
        const Image<ElementT1>& input1,
        const Image<ElementT2>& input2,
        const OutputT output = 0.0
        )
    {
        return euclidean_distance(std::execution::seq, input1, input2, output);
    }

    //  euclidean_distance Template Function Implementation for multiple channel image
    template<
        class ExPo,
        class ElementT1,
        class ElementT2
        >
    requires((std::same_as<ElementT1, RGB>) || (std::same_as<ElementT1, RGB_DOUBLE>) || (std::same_as<ElementT1, HSV>)) and
            ((std::same_as<ElementT2, RGB>) || (std::same_as<ElementT2, RGB_DOUBLE>) || (std::same_as<ElementT2, HSV>)) and
            (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static auto euclidean_distance(
        ExPo execution_policy,
        const Image<ElementT1>& input1,
        const Image<ElementT2>& input2
        )
    {
        return sqrt(execution_policy, two_input_map_reduce(execution_policy, input1.getImageData(), input2.getImageData(), MultiChannel<double>{},
            [&](auto&& element1, auto&& element2) {
                return pow(execution_policy, element1 - element2, 2.0);
            }));
    }

    //  euclidean_distance Template Function Implementation for multiple channel image
    template<
        class ExPo,
        class ElementT1,
        class ElementT2,
        std::size_t Size
        >
    constexpr static auto euclidean_distance(
        ExPo execution_policy,
        const Image<MultiChannel<ElementT1, Size>>& input1,
        const Image<MultiChannel<ElementT2, Size>>& input2
        )
    {
        return sqrt(execution_policy, two_input_map_reduce(execution_policy, input1.getImageData(), input2.getImageData(), MultiChannel<double, Size>{},
            [&](auto&& element1, auto&& element2) {
                return pow(execution_policy, element1 - element2, 2.0);
            }));
    }

    //  euclidean_distance Template Function Implementation for multiple channel image
    template<
        class ElementT1,
        class ElementT2
        >
    requires((std::same_as<ElementT1, RGB>) || (std::same_as<ElementT1, RGB_DOUBLE>) || (std::same_as<ElementT1, HSV>) || (is_MultiChannel<ElementT1>::value)) and
            ((std::same_as<ElementT2, RGB>) || (std::same_as<ElementT2, RGB_DOUBLE>) || (std::same_as<ElementT2, HSV>) || (is_MultiChannel<ElementT2>::value))
    constexpr static auto euclidean_distance(
        const Image<ElementT1>& input1,
        const Image<ElementT2>& input2
        )
    {
        return euclidean_distance(std::execution::seq, input1, input2);
    }

    //  hypot Template Function Implementation
    template<typename... Args>
    constexpr auto hypot(Args... args)
    {
        return std::sqrt((std::pow(args, 2.0) + ...));
    }

    //  Copy from https://stackoverflow.com/a/41171552
    template<class TupType, std::size_t... I>
    void print_tuple(const TupType& _tup, std::index_sequence<I...>)
    {
        std::cout << "(";
        (..., (std::cout << (I == 0? "" : ", ") << std::get<I>(_tup)));
        std::cout << ")\n";
    }

    template<class... T>
    void print_tuple(const std::tuple<T...>& _tup)
    {
        print_tuple(_tup, std::make_index_sequence<sizeof...(T)>());
    }
}

void euclidean_distanceRGBTest(
    const std::size_t sizex = 3,
    const std::size_t sizey = 2
    )
{
    TinyDIP::Image<TinyDIP::RGB> image1(sizex, sizey);
    image1.setAllValue(TinyDIP::RGB{1, 2, 3});
    image1.print(" ");
    TinyDIP::Image<TinyDIP::RGB> image2(sizex, sizey);
    image2.print(" ");
    std::cout << "euclidean_distance of image1 and image2: " 
    << TinyDIP::euclidean_distance(image1, image2) << '\n';
    return;
}

int main()
{
  auto start = std::chrono::system_clock::now();
  euclidean_distanceRGBTest();
  auto end = std::chrono::system_clock::now();
  std::chrono::duration<double> elapsed_seconds = end - start;
  std::time_t end_time = std::chrono::system_clock::to_time_t(end);
  if (elapsed_seconds.count() != 1)
  {
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << " seconds.\n";
  }
  else
  {
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << " second.\n";
  }
  return EXIT_SUCCESS;
}

The output of the test code above:

( 1 2 3 ) ( 1 2 3 ) ( 1 2 3 ) 
( 1 2 3 ) ( 1 2 3 ) ( 1 2 3 ) 

( 0 0 0 ) ( 0 0 0 ) ( 0 0 0 ) 
( 0 0 0 ) ( 0 0 0 ) ( 0 0 0 ) 

euclidean_distance of image1 and image2: {2.44949, 4.89898, 7.34847, }
Computation finished at Mon Jan 13 15:40:09 2025
elapsed time: 9.2951e-05 seconds.

Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

  • Which question it is a follow-up to?

    apply_each_single_output Template Function Implementation for Image in C++.

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

    apply_multichannel template functions implementation is proposed in this post.

  • Why a new review is being asked for?

    Please check the implementation of apply_multichannel template functions and the related code.

\$\endgroup\$
1
  • \$\begingroup\$ If T has 2 channels and you forget to explicitly specify that, the result will have 3 channels and it will do an out-of-bounds read. \$\endgroup\$ Commented Jan 10 at 19:39

0

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.