Compiler Explorer

Source code

//  apply_each_single_output Template Function Implementation for Image in C++ (Rev.2)
//  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;
        }
    };

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
    template <typename>
    struct is_MultiChannel : std::false_type {};

template <typename ...T>
    struct is_MultiChannel<MultiChannel<T...>> : 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, ' ') << std::format("{}", 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(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!");
    }

//  getPlane template function implementation
    template<class OutputT = unsigned char>
    constexpr static auto getPlane(const Image<RGB>& input, std::size_t index)
    {
        auto input_data = input.getImageData();
        std::vector<OutputT> output_data;
        output_data.resize(input.count());
        #pragma omp parallel for
        for (std::size_t i = 0; i < input.count(); ++i)
        {
            output_data[i] = input_data[i].channels[index];
        }
        auto output = Image<OutputT>(output_data, input.getSize());
        return output;
    }

//  getPlane template function implementation
    template<class T = HSV, class OutputT = double>
    requires (std::same_as<T, HSV> || std::same_as<T, RGB_DOUBLE>)
    constexpr static auto getPlane(const Image<T>& input, std::size_t index)
    {
        auto input_data = input.getImageData();
        std::vector<OutputT> output_data;
        output_data.resize(input.count());
        #pragma omp parallel for
        for (std::size_t i = 0; i < input.count(); ++i)
        {
            output_data[i] = input_data[i].channels[index];
        }
        auto output = Image<OutputT>(output_data, input.getSize());
        return output;
    }

//  getRplane function implementation
    constexpr static auto getRplane(const Image<RGB>& input)
    {
        return getPlane(input, 0);
    }

//  getRplane function implementation
    constexpr static auto getRplane(const Image<RGB_DOUBLE>& input)
    {
        return getPlane(input, 0);
    }

//  getGplane function implementation
    constexpr static auto getGplane(const Image<RGB>& input)
    {
        return getPlane(input, 1);
    }

//  getGplane function implementation
    constexpr static auto getGplane(const Image<RGB_DOUBLE>& input)
    {
        return getPlane(input, 1);
    }

//  getBplane function implementation
    constexpr static auto getBplane(const Image<RGB>& input)
    {
        return getPlane(input, 2);
    }

//  getBplane function implementation
    constexpr static auto getBplane(const Image<RGB_DOUBLE>& input)
    {
        return getPlane(input, 2);
    }

constexpr static auto getHplane(const Image<HSV>& input)
    {
        return getPlane(input, 0);
    }

constexpr static auto getSplane(const Image<HSV>& input)
    {
        return getPlane(input, 1);
    }

constexpr static auto getVplane(const Image<HSV>& input)
    {
        return getPlane(input, 2);
    }

//  constructRGB template function implementation
    template<typename OutputT = RGB>
    constexpr static auto constructRGB(const Image<GrayScale>& r, const Image<GrayScale>& g, const Image<GrayScale>& b)
    {
        check_size_same(r, g);
        check_size_same(g, b);
        auto image_data_r = r.getImageData();
        auto image_data_g = g.getImageData();
        auto image_data_b = b.getImageData();
        std::vector<OutputT> new_data;
        new_data.resize(r.count());
        #pragma omp parallel for
        for (std::size_t index = 0; index < r.count(); ++index)
        {
            OutputT rgb {   image_data_r[index],
                        image_data_g[index],
                        image_data_b[index]};
            new_data[index] = rgb;
        }
        Image<OutputT> output(new_data, r.getSize());
        return output;
    }

//  constructRGBDOUBLE template function implementation
    template<typename OutputT = RGB_DOUBLE>
    constexpr static auto constructRGBDOUBLE(const Image<double>& r, const Image<double>& g, const Image<double>& b)
    {
        check_size_same(r, g);
        check_size_same(g, b);
        auto image_data_r = r.getImageData();
        auto image_data_g = g.getImageData();
        auto image_data_b = b.getImageData();
        std::vector<OutputT> new_data;
        new_data.resize(r.count());
        #pragma omp parallel for
        for (std::size_t index = 0; index < r.count(); ++index)
        {
            OutputT rgb_double { image_data_r[index],
                                    image_data_g[index],
                                    image_data_b[index]};
            new_data[index] = rgb_double;
        }
        Image<OutputT> output(new_data, r.getSize());
        return output;
    }

//  constructHSV template function implementation
    template<typename OutputT = HSV>
    constexpr static auto constructHSV(const Image<double>& h, const Image<double>& s, const Image<double>& v)
    {
        check_size_same(h, s);
        check_size_same(s, v);
        auto image_data_h = h.getImageData();
        auto image_data_s = s.getImageData();
        auto image_data_v = v.getImageData();
        std::vector<OutputT> new_data;
        new_data.resize(h.count());
        #pragma omp parallel for
        for (std::size_t index = 0; index < h.count(); ++index)
        {
            OutputT hsv {   image_data_h[index],
                        image_data_s[index],
                        image_data_v[index]};
            new_data[index] = hsv;
        }
        Image<OutputT> output(new_data, h.getSize());
        return output;
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<RGB>& input, F operation, Args&&... args)
    {
        auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input), args...); });
        auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input), args...); });
        auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input), args...); });
        return constructRGB(Rplane.get(), Gplane.get(), Bplane.get());
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<RGB_DOUBLE>& input, F operation, Args&&... args)
    {
        auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input), args...); });
        auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input), args...); });
        auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input), args...); });
        return constructRGBDOUBLE(Rplane.get(), Gplane.get(), Bplane.get());
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<HSV>& input, F operation, Args&&... args)
    {
        auto Hplane = std::async(std::launch::async, [&] { return std::invoke(operation, getHplane(input), args...); });
        auto Splane = std::async(std::launch::async, [&] { return std::invoke(operation, getSplane(input), args...); });
        auto Vplane = std::async(std::launch::async, [&] { return std::invoke(operation, getVplane(input), args...); });
        return constructHSV(Hplane.get(), Splane.get(), Vplane.get());
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<RGB>& input1, const Image<RGB>& input2, F operation, Args&&... args)
    {
        auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input1), getRplane(input2), args...); });
        auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input1), getGplane(input2), args...); });
        auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input1), getBplane(input2), args...); });
        return constructRGB(Rplane.get(), Gplane.get(), Bplane.get());
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<RGB_DOUBLE>& input1, const Image<RGB_DOUBLE>& input2, F operation, Args&&... args)
    {
        auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input1), getRplane(input2), args...); });
        auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input1), getGplane(input2), args...); });
        auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input1), getBplane(input2), args...); });
        return constructRGBDOUBLE(Rplane.get(), Gplane.get(), Bplane.get());
    }

//  apply_each template function implementation
    template<class F, class... Args>
    constexpr static auto apply_each(const Image<HSV> input1, const Image<HSV> input2, F operation, Args&&... args)
    {
        auto Hplane = std::async(std::launch::async, [&] { return std::invoke(operation, getHplane(input1), getHplane(input2), args...); });
        auto Splane = std::async(std::launch::async, [&] { return std::invoke(operation, getSplane(input1), getSplane(input2), args...); });
        auto Vplane = std::async(std::launch::async, [&] { return std::invoke(operation, getVplane(input1), getVplane(input2), args...); });
        return constructHSV(Hplane.get(), Splane.get(), Vplane.get());
    }

//  apply_each_single_output template function implementation
    template<class ElementT, class F, class... Args>
    constexpr static auto apply_each_single_output(const std::size_t channel_count, const Image<ElementT>& input1, const Image<ElementT>& input2, F operation, Args&&... args)
    {
        std::vector<decltype(std::invoke(operation, getPlane(input1, 0), getPlane(input2, 0), args...))> output;
        output.reserve(channel_count);
        for (std::size_t channel_index = 0; channel_index < channel_count; ++channel_index)
        {
            output.emplace_back(std::invoke(operation, getPlane(input1, channel_index), getPlane(input2, channel_index), args...));
        }
        return output;
    }

//  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<
        arithmetic OutputT = double,
        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>))
    constexpr static auto euclidean_distance(
        const Image<ElementT1>& input1,
        const Image<ElementT2>& input2,
        const OutputT output = 0.0
        )
    {
        return apply_each_single_output(3, input1, input2, [&](auto&& planes1, auto&& planes2) { return euclidean_distance(planes1, planes2, output); });
    }

//  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: " << '\n';
    TinyDIP::recursive_print(TinyDIP::euclidean_distance(image1, image2));
    return;
}

void euclidean_distanceRGB_DOUBLETest(
    const std::size_t sizex = 3,
    const std::size_t sizey = 2
    )
{
    TinyDIP::Image<TinyDIP::RGB_DOUBLE> image1(sizex, sizey);
    image1.setAllValue(TinyDIP::RGB_DOUBLE{1.5, 2.5, 3.5});
    image1.print();
    TinyDIP::Image<TinyDIP::RGB_DOUBLE> image2(sizex, sizey);
    image2.print();
    std::cout << "euclidean_distance of image1 and image2: " << '\n';
    TinyDIP::recursive_print(TinyDIP::euclidean_distance(image1, image2));
    return;
}

int main()
{
  auto start = std::chrono::system_clock::now();
  euclidean_distanceRGBTest();
  euclidean_distanceRGB_DOUBLETest();
  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;
}