apply_multichannel Template Function Implementation in C++

Question

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.

• G. Sliepen: 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.