Skip to main content
Became Hot Network Question
Update full testing code
Source Link
JimmyHu
  • 7.6k
  • 2
  • 11
  • 48
  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Multi-dimensional Image Data Structure withThree Variadicdimensional Templatedata Functionsstructure in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }
  
        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];
        }

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}

Godbolt link is here.Godbolt link is here.

  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            constexpr std::size_t getDimensionality() const
            {
                return size.size();
            }
    
            constexpr std::size_t getWidth() const
            {
                return size[0];
            }
    
            constexpr std::size_t getHeight() const noexcept
            {
                return size[1];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Multi-dimensional Image Data Structure with Variadic Template Functions in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        constexpr std::size_t getDimensionality() const
        {
            return size.size();
        }

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

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}

Godbolt link is here.

  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Three dimensional data structure in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }
  
        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];
        }

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}

Godbolt link is here.

Update full testing code
Source Link
JimmyHu
  • 7.6k
  • 2
  • 11
  • 48
  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            constexpr std::size_t getDimensionality() const
            {
                return size.size();
            }
    
            constexpr std::size_t getWidth() const
            {
                return size[0];
            }
    
            constexpr std::size_t getHeight() const noexcept
            {
                return size[1];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Three Multi-dimensional dataImage structureData Structure with Variadic Template Functions in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        constexpr std::size_t getDimensionality() const
        {
            return size.size();
        }

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

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}
Test with 2D image:
1   1   1   
1   3   1   
1   1   1   

Test with 3D image:
4   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0    


Computation finished at Sat Dec 30 11:10:29 2023
elapsed time: 0.00150394

Godbolt link is here.Godbolt link is here.

  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            constexpr std::size_t getWidth() const
            {
                return size[0];
            }
    
            constexpr std::size_t getHeight() const noexcept
            {
                return size[1];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Three dimensional data structure in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

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

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}
Test with 2D image:
1   1   1   
1   3   1   
1   1   1   

Test with 3D image:
4   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

Godbolt link is here.

  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            constexpr std::size_t getDimensionality() const
            {
                return size.size();
            }
    
            constexpr std::size_t getWidth() const
            {
                return size[0];
            }
    
            constexpr std::size_t getHeight() const noexcept
            {
                return size[1];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    
//  Multi-dimensional Image Data Structure with Variadic Template Functions in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        constexpr std::size_t getDimensionality() const
        {
            return size.size();
        }

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

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}
Test with 2D image:
1   1   1   
1   3   1   
1   1   1   

Test with 3D image:
4   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0    


Computation finished at Sat Dec 30 11:10:29 2023
elapsed time: 0.00150394

Godbolt link is here.

Source Link
JimmyHu
  • 7.6k
  • 2
  • 11
  • 48

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

This is a follow-up question for Three dimensional data structure in C++. I am trying to implement multi-dimensional image data structure with variadic template functions. For example, image.at(4, 3) is to access the element at location x = 4 and y = 3 in two dimensional case. In three dimensional case, image.at(4, 3, 2) is to access the element at location x = 4, y = 3 and z = 2. In this way, both constructing new image and accessing elements are easy.

The experimental implementation

  • Image Class Implementation

    //  image.h
    namespace TinyDIP
    {
        template <typename ElementT>
        class Image
        {
        public:
            Image() = default;
    
            Image(const std::size_t width, const std::size_t height):
                image_data(width * height)
                {
                    size.reserve(2);
                    size.emplace_back(width);
                    size.emplace_back(height);
                }
    
            Image(const std::size_t width, const std::size_t height, const std::size_t depth):
                image_data(width * height * depth)
                {
                    size.reserve(3);
                    size.emplace_back(width);
                    size.emplace_back(height);
                    size.emplace_back(depth);
                }
    
            Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
                image_data(x * y * z * w)
                {
                    size.reserve(4);
                    size.emplace_back(x);
                    size.emplace_back(y);
                    size.emplace_back(z);
                    size.emplace_back(w);
                }
    
            Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
                image_data(a * b * c * d * e)
                {
                    size.reserve(5);
                    size.emplace_back(a);
                    size.emplace_back(b);
                    size.emplace_back(c);
                    size.emplace_back(d);
                    size.emplace_back(e);
                }
    
            Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
            {
                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 = input;
            }
    
            Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
            {
                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)
            {
                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(input), std::ranges::end(input));    //  flatten
                }
                return;
            }
    
            template<typename... Args>
            constexpr ElementT& at(const Args... indexInput)
            {
                checkBoundary(indexInput...);
                constexpr std::size_t n = sizeof...(Args);
                if(n != size.size())
                {
                    throw std::runtime_error("Dimensionality mismatched!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    return image_data[y * size[0] + x];
                }
                else if constexpr (n == 3)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    return image_data[(z * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 4)
                {
                    auto x = get_from_variadic_template<1>(indexInput...);
                    auto y = get_from_variadic_template<2>(indexInput...);
                    auto z = get_from_variadic_template<3>(indexInput...);
                    auto w = get_from_variadic_template<4>(indexInput...);
                    return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
                }
                else if constexpr (n == 5)
                {
                    auto a = get_from_variadic_template<1>(indexInput...);
                    auto b = get_from_variadic_template<2>(indexInput...);
                    auto c = get_from_variadic_template<3>(indexInput...);
                    auto d = get_from_variadic_template<4>(indexInput...);
                    auto e = get_from_variadic_template<5>(indexInput...);
                    return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
                }
            }
    
            constexpr std::size_t getWidth() const
            {
                return size[0];
            }
    
            constexpr std::size_t getHeight() const noexcept
            {
                return size[1];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                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";
                }
            }
    
            //  Enable this function if ElementT = RGB
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            requires(std::same_as<ElementT, RGB>)
            {
                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;
            }
    
            void setAllValue(const ElementT input)
            {
                std::fill(image_data.begin(), image_data.end(), input);
            }
    
            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!");
                }
                if constexpr (n == 2)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                }
                if constexpr (n == 3)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Given x out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Given y out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Given z out of range!");
                }
                if constexpr (n == 4)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                }
                if constexpr (n == 5)
                {
                    if (get_from_variadic_template<1>(indexInput...) >= size[0])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<2>(indexInput...) >= size[1])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<3>(indexInput...) >= size[2])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<4>(indexInput...) >= size[3])
                        throw std::out_of_range("Index out of range!");
                    if (get_from_variadic_template<5>(indexInput...) >= size[4])
                        throw std::out_of_range("Index out of range!");
                }
            }
        };
    
        template<typename T, typename ElementT>
        concept is_Image = std::is_same_v<T, Image<ElementT>>;
    }
    

Full Testing Code

The full testing code:

//  Three dimensional data structure in C++
//  Developed by Jimmy Hu

#include <algorithm>
#include <cassert>      //  for assert
#include <chrono>       //  for std::chrono::system_clock::now
#include <cmath>        //  for std::exp
#include <concepts>
#include <execution>    //  for std::is_execution_policy_v
#include <iostream>     //  for std::cout
#include <vector>

struct RGB
{
    std::uint8_t channels[3];
};

using GrayScale = std::uint8_t;

namespace TinyDIP
{
    //  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return std::size_t{0};
    }

    template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + std::size_t{1};
    }

    template<std::size_t index = 1, typename Arg, typename... Args>
    constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
    {
        if constexpr (index > 1)
            return get_from_variadic_template<index - 1>(inputs...);
        else
            return first;
    }

    template<typename... Args>
    constexpr static auto& first_of(Args&... inputs) {
        return get_from_variadic_template<1>(inputs...);
    }

    template<std::size_t, typename, typename...>
    struct get_from_variadic_template_struct { };

    template<typename T1, typename... Ts>
    struct get_from_variadic_template_struct<1, T1, Ts...>
    {
        using type = T1;
    };

    template<std::size_t index, typename T1, typename... Ts>
    requires ( requires { typename get_from_variadic_template_struct<index - 1, Ts...>::type; })
    struct get_from_variadic_template_struct<index, T1, Ts...>
    {
        using type = typename get_from_variadic_template_struct<index - 1, Ts...>::type;
    };

    template<std::size_t index, typename... Ts>
    using get_from_variadic_template_t = typename get_from_variadic_template_struct<index, Ts...>::type;
}

//  image.h
namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            image_data(width * height)
            {
                size.reserve(2);
                size.emplace_back(width);
                size.emplace_back(height);
            }

        Image(const std::size_t width, const std::size_t height, const std::size_t depth):
            image_data(width * height * depth)
            {
                size.reserve(3);
                size.emplace_back(width);
                size.emplace_back(height);
                size.emplace_back(depth);
            }

        Image(const std::size_t x, const std::size_t y, const std::size_t z, const std::size_t w):
            image_data(x * y * z * w)
            {
                size.reserve(4);
                size.emplace_back(x);
                size.emplace_back(y);
                size.emplace_back(z);
                size.emplace_back(w);
            }

        Image(const std::size_t a, const std::size_t b, const std::size_t c, const std::size_t d, const std::size_t e):
            image_data(a * b * c * d * e)
            {
                size.reserve(5);
                size.emplace_back(a);
                size.emplace_back(b);
                size.emplace_back(c);
                size.emplace_back(d);
                size.emplace_back(e);
            }

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight)
        {
            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 = input;
        }

        Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
        {
            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)
        {
            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(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        template<typename... Args>
        constexpr ElementT& at(const Args... indexInput)
        {
            checkBoundary(indexInput...);
            constexpr std::size_t n = sizeof...(Args);
            if(n != size.size())
            {
                throw std::runtime_error("Dimensionality mismatched!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

        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!");
            }
            if constexpr (n == 2)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                return image_data[y * size[0] + x];
            }
            else if constexpr (n == 3)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                return image_data[(z * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 4)
            {
                auto x = get_from_variadic_template<1>(indexInput...);
                auto y = get_from_variadic_template<2>(indexInput...);
                auto z = get_from_variadic_template<3>(indexInput...);
                auto w = get_from_variadic_template<4>(indexInput...);
                return image_data[((w * size[2] + z) * size[1] + y) * size[0] + x];
            }
            else if constexpr (n == 5)
            {
                auto a = get_from_variadic_template<1>(indexInput...);
                auto b = get_from_variadic_template<2>(indexInput...);
                auto c = get_from_variadic_template<3>(indexInput...);
                auto d = get_from_variadic_template<4>(indexInput...);
                auto e = get_from_variadic_template<5>(indexInput...);
                return image_data[(((e * size[3] + d) * size[2] + c) * size[1] + b) * size[0] + a];
            }
        }

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

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            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";
            }
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            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;
        }

        void setAllValue(const ElementT input)
        {
            std::fill(image_data.begin(), image_data.end(), input);
        }

        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!");
            }
            if constexpr (n == 2)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
            }
            if constexpr (n == 3)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Given x out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Given y out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Given z out of range!");
            }
            if constexpr (n == 4)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
            }
            if constexpr (n == 5)
            {
                if (get_from_variadic_template<1>(indexInput...) >= size[0])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<2>(indexInput...) >= size[1])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<3>(indexInput...) >= size[2])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<4>(indexInput...) >= size[3])
                    throw std::out_of_range("Index out of range!");
                if (get_from_variadic_template<5>(indexInput...) >= size[4])
                    throw std::out_of_range("Index out of range!");
            }
        }
    };

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

namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getWidth() == y.getWidth();
    }

    template<typename ElementT>
    constexpr bool is_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_width_same(x, y) && is_width_same(y, z) && is_width_same(x, z);
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return x.getHeight() == y.getHeight();
    }

    template<typename ElementT>
    constexpr bool is_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_height_same(x, y) && is_height_same(y, z) && is_height_same(x, z);
    }
    
    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        return is_width_same(x, y) && is_height_same(x, y);
    }

    template<typename ElementT>
    constexpr bool is_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        return is_size_same(x, y) && is_size_same(y, z) && is_size_same(x, z);
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_width_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_width_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_width_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert(is_height_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_height_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert(is_height_same(x, y, z));
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        assert_width_same(x, y);
        assert_height_same(x, y);
    }

    template<typename ElementT>
    constexpr void assert_size_same(const Image<ElementT>& x, const Image<ElementT>& y, const Image<ElementT>& z)
    {
        assert_size_same(x, y);
        assert_size_same(y, z);
        assert_size_same(x, z);
    }

    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!");
    }

    template<typename ElementT>
    constexpr void check_size_same(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        check_width_same(x, y);
        check_height_same(x, y);
    }
}

template<typename ElementT>
void multidimensionalImageTest(const size_t size = 3)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(1, 1) = 3;
    image2d.print();
    std::cout << "Test with 3D image:\n";
    auto image3d = TinyDIP::Image<double>(size, size, size);
    image3d.setAllValue(0);
    image3d.at(0, 0, 0) = 4;
    image3d.print();
    return;
}

int main()
{
    auto start = std::chrono::system_clock::now();
    multidimensionalImageTest<double>();
    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);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}

The output of the test code above:

Test with 2D image:
1   1   1   
1   3   1   
1   1   1   

Test with 3D image:
4   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

0   0   0   
0   0   0   
0   0   0   

Godbolt link is here.

All suggestions are welcome.

The summary information:

  • Which question it is a follow-up to?

    Three dimensional data structure in C++

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

    I am trying to create multi-dimensional image data structure with variadic template functions in this post.

  • Why a new review is being asked for?

    Despite the usage of Image class is easy and the feasible (dimensionality of an Image object can be changed in run-time), there are several if constexpr blocks to deal with various dimension case. Is there any better way to improve the implementation?