5
\$\begingroup\$

This is a follow-up question for Three dimensional gaussian image generator in C++. Considering the suggestion from G. Sliepen:

Structure of a 3D image

Your 2D Image class is not storing a vector of vectors, rather it stores a single vector with all the pixels in the 2D image. But instead of creating a 3D image type, your guassianFigure3D() just returns a std::vector<Image>. This is a less efficient way of storing the pixel data, and also lacks all of the properties of Image, like having an at() function that takes \$x\$, \$y\$ and \$z\$ coordinates, or an operator+() that can add two 3D images together.

I am trying to create 3D data structure Cube in this post and there are several properties:

  • at() function that takes \$x\$, \$y\$ and \$z\$ coordinates for element access in Cube class.

  • print function to print values.

  • operator+=, operator-=, operator+(), operator-() for adding / subtracting two Cubes.

The experimental implementation

  • Cube Template Class (in file "cube.h")

    namespace TinyDIP
    {
        template <typename ElementT>
        class Cube
        {
        public:
            Cube() = default;
    
            Cube(const std::size_t newWidth, const std::size_t newHeight, const std::size_t newDepth):
                width(width),
                height(height),
                depth(newDepth),
                data(width * height * depth) { }
    
            Cube(const int newWidth, const int newHeight, const int newDepth, ElementT initVal):
                width(newWidth),
                height(newHeight),
                depth(newDepth),
                data(width * height * depth, initVal) {}
    
            Cube(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight, const std::size_t newDepth):
                width(newWidth),
                height(newHeight),
                depth(newDepth)
            {
                if (input.size() != newWidth * newHeight * newDepth)
                {
                    throw std::runtime_error("Data input and the given size are mismatched!");
                }
                data = input;
            }    
    
            Cube(const std::vector<Image<ElementT>>& input)
            {
                width = input[0].getWidth();
                height = input[0].getHeight();
                depth = input.size();
                data.resize(width * height * depth);
    
                for (std::size_t z = 0; z < input.size(); ++z)
                {
                    auto image = input[z];
                    for (std::size_t y = 0; y < height; ++y)
                    {
                        for (std::size_t x = 0; x < width; ++x)
                        {
                            data[z * width * height + y * width + x] = image.at(x, y);
                        }
                    }
                }
                return;
            }
    
            constexpr ElementT& at(const unsigned int x, const unsigned int y, const unsigned int z)
            { 
                checkBoundary(x, y, z);
                return data[z * width * height + y * width + x];
            }
    
            constexpr ElementT const& at(const unsigned int x, const unsigned int y, const unsigned int z) const
            {
                checkBoundary(x, y, z);
                return data[z * width * height + y * width + x];
            }
    
            constexpr auto getSizeX() const noexcept
            {
                return width;
            }
    
            constexpr auto getSizeY() const noexcept
            {
                return height;
            }
    
            constexpr auto getSizeZ() const noexcept
            {
                return depth;
            }
    
            constexpr auto getWidth() const noexcept
            {
                return width;
            }
    
            constexpr auto getHeight() const noexcept
            {
                return height;
            }
    
            constexpr auto getDepth() const noexcept
            {
                return depth;
            }
    
            std::vector<ElementT> const& getData() const noexcept { return data; }      //  expose the internal data
    
            void print(std::string separator = "\t", std::ostream& os = std::cout) const
            {
                for(std::size_t z = 0; z < depth; ++z)
                {
                    for (std::size_t y = 0; y < height; ++y)
                    {
                        for (std::size_t x = 0; x < width; ++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";
                return;
            }
    
            friend std::ostream& operator<<(std::ostream& os, const Cube<ElementT>& rhs)
            {
                const std::string separator = "\t";
                rhs.print(separator, os);
                return os;
            }
    
            Cube<ElementT>& operator+=(const Cube<ElementT>& rhs)
            {
                check_size_same(rhs, *this);
                std::transform(std::ranges::cbegin(data), std::ranges::cend(data), std::ranges::cbegin(rhs.data),
                       std::ranges::begin(data), std::plus<>{});
                return *this;
            }
    
            Cube<ElementT>& operator-=(const Cube<ElementT>& rhs)
            {
                check_size_same(rhs, *this);
                std::transform(std::ranges::cbegin(data), std::ranges::cend(data), std::ranges::cbegin(rhs.data),
                       std::ranges::begin(data), std::minus<>{});
                return *this;
            }
    
            Cube<ElementT>& operator*=(const Cube<ElementT>& rhs)
            {
                check_size_same(rhs, *this);
                std::transform(std::ranges::cbegin(data), std::ranges::cend(data), std::ranges::cbegin(rhs.data),
                       std::ranges::begin(data), std::multiplies<>{});
                return *this;
            }
    
            Cube<ElementT>& operator/=(const Cube<ElementT>& rhs)
            {
                check_size_same(rhs, *this);
                std::transform(std::ranges::cbegin(data), std::ranges::cend(data), std::ranges::cbegin(rhs.data),
                       std::ranges::begin(data), std::divides<>{});
                return *this;
            }
    
            friend bool operator==(Cube<ElementT> const&, Cube<ElementT> const&) = default;
    
            friend bool operator!=(Cube<ElementT> const&, Cube<ElementT> const&) = default;
    
            friend Cube<ElementT> operator+(Cube<ElementT> input1, const Cube<ElementT>& input2)
            {
                return input1 += input2;
            }
    
            friend Cube<ElementT> operator-(Cube<ElementT> input1, const Cube<ElementT>& input2)
            {
                return input1 -= input2;
            }
    
            friend Cube<ElementT> operator*(Cube<ElementT> input1, ElementT input2)
            {
                return multiplies(input1, input2);
            }
    
            Cube<ElementT>& operator=(Cube<ElementT> const& input) = default;   //  Copy Assign
    
            Cube<ElementT>& operator=(Cube<ElementT>&& other) = default;        //  Move Assign
    
            Cube(const Cube<ElementT> &input) = default;                        //  Copy Constructor
    
            Cube(Cube<ElementT> &&input) = default;                             //  Move Constructor
    
        private:
            std::size_t width;
            std::size_t height;
            std::size_t depth;
            std::vector<ElementT> data;
    
            void checkBoundary(const size_t x, const size_t y, const size_t z) const
            {
                if (x >= width)
                    throw std::out_of_range("Given x out of range!");
                if (y >= height)
                    throw std::out_of_range("Given y out of range!");
                if (z >= depth)
                    throw std::out_of_range("Given z out of range!");
            }
        };
    }
    
  • Helper Functions for Cube (in file "cube_operations.h")

    namespace TinyDIP
    {
        template<typename ElementT>
        constexpr bool is_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            return x.getSizeX() == y.getSizeX();
        }
    
        template<typename ElementT>
        constexpr bool is_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            return is_width_same(x, y) && is_width_same(y, z);
        }
    
        template<typename ElementT>
        constexpr bool is_height_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            return x.getSizeY() == y.getSizeY();
        }
    
        template<typename ElementT>
        constexpr bool is_height_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            return is_height_same(x, y) && is_height_same(y, z);
        }
    
        template<typename ElementT>
        constexpr bool is_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            return x.getSizeZ() == y.getSizeZ();
        }
    
        template<typename ElementT>
        constexpr bool is_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            return is_depth_same(x, y) && is_depth_same(y, z);
        }
    
        template<typename ElementT>
        constexpr bool is_size_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            return is_width_same(x, y) && is_height_same(x, y) && is_depth_same(x, y);
        }
    
        template<typename ElementT>
        constexpr bool is_size_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            return is_size_same(x, y) && is_size_same(y, z);
        }
    
        template<typename ElementT>
        constexpr void assert_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            assert(is_width_same(x, y));
        }
    
        template<typename ElementT>
        constexpr void assert_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            assert(is_width_same(x, y, z));
        }
    
        template<typename ElementT>
        constexpr void assert_height_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            assert(is_height_same(x, y));
        }
    
        template<typename ElementT>
        constexpr void assert_height_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            assert(is_height_same(x, y, z));
        }
    
        template<typename ElementT>
        constexpr void assert_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            assert(is_depth_same(x, y));
        }
    
        template<typename ElementT>
        constexpr void assert_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            assert(is_depth_same(x, y));
            assert(is_depth_same(y, z));
        }
    
        template<typename ElementT>
        constexpr void assert_size_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            assert_width_same(x, y);
            assert_height_same(x, y);
            assert_depth_same(x, y);
        }
    
        template<typename ElementT>
        constexpr void assert_size_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
        {
            assert_size_same(x, y);
            assert_size_same(y, z);
        }
    
        template<typename ElementT>
        constexpr void check_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            if (!is_width_same(x, y))
                throw std::runtime_error("Width mismatched!");
        }
    
        template<typename ElementT>
        constexpr void check_height_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            if (!is_height_same(x, y))
                throw std::runtime_error("Height mismatched!");
        }
    
        template<typename ElementT>
        constexpr void check_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            if (!is_depth_same(x, y))
                throw std::runtime_error("Depth mismatched!");
        }
    
        template<typename ElementT>
        constexpr void check_size_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
        {
            check_width_same(x, y);
            check_height_same(x, y);
            check_depth_same(x, y);
        }
    
        //  voxelwiseOperation template function implementation
        template<std::size_t unwrap_level = 1, class... Args>
        constexpr static auto voxelwiseOperation(auto op, const Args&... inputs)
        {
            auto output = Cube(
                recursive_transform<unwrap_level>(
                    [&](auto&& element1, auto&&... elements) 
                        {
                            return op(element1, elements...);
                        },
                    inputs.getData()...),
                first_of(inputs...).getWidth(),
                first_of(inputs...).getHeight(),
                first_of(inputs...).getDepth());
            return output;
        }
    
        template<std::size_t unwrap_level = 1, class ExPo, class InputT>
        requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
        constexpr static auto voxelwiseOperation(ExPo execution_policy, auto op, const Image<InputT>& input1)
        {
            auto output = Cube(
                recursive_transform<unwrap_level>(
                    execution_policy,
                    [&](auto&& element1) 
                        {
                            return op(element1);
                        },
                    (input1.getData())),
                input1.getWidth(),
                input1.getHeight(),
                input1.getDepth());
            return output;
        }
    
        //  plus template function implementation
        template<class InputT>
        constexpr static Cube<InputT> plus(const Cube<InputT>& input1)
        {
            return input1;
        }
    
        template<class InputT, class... Args>
        constexpr static Cube<InputT> plus(const Cube<InputT>& input1, const Args&... inputs)
        {
            return voxelwiseOperation(std::plus<>{}, input1, plus(inputs...));
        }
    
        template<class InputT, class... Args>
        constexpr static auto plus(const std::vector<Cube<InputT>>& input1, const Args&... inputs)
        {
            return recursive_transform<1>(
                [](auto&& input1_element, auto&&... inputs_element)
                {
                    return plus(input1_element, inputs_element...);
                }, input1, inputs...);
        }
    
        //  subtract template function implementation
        template<class InputT>
        constexpr static Cube<InputT> subtract(const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            check_size_same(input1, input2);
            return voxelwiseOperation(std::minus<>{}, input1, input2);
        }
    
        template<class InputT>
        constexpr static auto subtract(const std::vector<Cube<InputT>>& input1, const std::vector<Cube<InputT>>& input2)
        {
            assert(input1.size() == input2.size());
            return recursive_transform<1>(
                [](auto&& input1_element, auto&& input2_element)
                {
                    return subtract(input1_element, input2_element);
                }, input1, input2);
        }
    
        //  multiplies template function implementation
        template<class InputT>
        constexpr static Cube<InputT> multiplies(const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            return voxelwiseOperation(std::multiplies<>{}, input1, input2);
        }
    
        template<class InputT, class TimesT>
        requires(std::floating_point<TimesT> || std::integral<TimesT>)
        constexpr static Cube<InputT> multiplies(const Cube<InputT>& input1, const TimesT times)
        {
            return multiplies(
                input1,
                Cube(input1.getWidth(), input1.getHeight(), input1.getDepth(), times)
            );
        }
    
        template<class InputT, class TimesT>
        requires(std::floating_point<TimesT> || std::integral<TimesT>)
        constexpr static Cube<InputT> multiplies(const TimesT times, const Cube<InputT>& input1)
        {
            return multiplies(input1, times);
        }
    
        template<class ExPo, class InputT>
        requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
        constexpr static Cube<InputT> multiplies(ExPo execution_policy, const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            return voxelwiseOperation(execution_policy, std::multiplies<>{}, input1, input2);
        }
    
        template<class InputT>
        constexpr static Cube<InputT> divides(const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            return voxelwiseOperation(std::divides<>{}, input1, input2);
        }
    
        template<class InputT>
        constexpr static auto divides(const std::vector<Cube<InputT>>& input1, const std::vector<Cube<InputT>>& input2)
        {
            assert(input1.size() == input2.size());
            return recursive_transform<1>(
                [](auto&& input1_element, auto&& input2_element)
                {
                    return divides(input1_element, input2_element);
                }, input1, input2);
        }
    
        template<class ExPo, class InputT>
        requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
        constexpr static Cube<InputT> divides(ExPo execution_policy, const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            return voxelwiseOperation(execution_policy, std::divides<>{}, input1, input2);
        }
    
        template<class InputT>
        constexpr static Cube<InputT> modulus(const Cube<InputT>& input1, const Cube<InputT>& input2)
        {
            return voxelwiseOperation(std::modulus<>{}, input1, input2);
        }
    
        template<class InputT>
        constexpr static Image<InputT> negate(const Image<InputT>& input1)
        {
            return voxelwiseOperation(std::negate<>{}, input1);
        }
    }
    

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;

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

        Image(const std::size_t width, const std::size_t height):
            width(width),
            height(height),
            image_data(width * height) { }

        Image(const std::size_t width, const std::size_t height, const ElementT initVal):
            width(width),
            height(height),
            image_data(width * height, initVal) {}

        Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight):
            width(newWidth),
            height(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):
            width(newWidth),
            height(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)
        {
            height = input.size();
            width = input[0].size();
            
            for (auto& rows : input)
            {
                image_data.insert(image_data.end(), std::ranges::begin(input), std::ranges::end(input));    //  flatten
            }
            return;
        }

        constexpr ElementT& at(const unsigned int x, const unsigned int y)
        { 
            checkBoundary(x, y);
            return image_data[y * width + x];
        }

        constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
        {
            checkBoundary(x, y);
            return image_data[y * width + x];
        }

        constexpr std::size_t getWidth() const
        {
            return width;
        }

        constexpr std::size_t getHeight() const noexcept
        {
            return height;
        }

        constexpr auto getSize() noexcept
        {
            return std::make_tuple(width, height);
        }

        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
        {
            for (std::size_t y = 0; y < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++x)
                {
                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                    os << +at(x, y) << separator;
                }
                os << "\n";
            }
            os << "\n";
            return;
        }

        //  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 < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++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;
        }

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

        Image<ElementT>& operator=(Image<ElementT> const& input) = default;  //  Copy Assign

        Image<ElementT>& operator=(Image<ElementT>&& other) = default;       //  Move Assign

        Image(const Image<ElementT> &input) = default;                       //  Copy Constructor

        Image(Image<ElementT> &&input) = default;                            //  Move Constructor
        
#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::size_t width;
        std::size_t height;
        std::vector<ElementT> image_data;

        void checkBoundary(const size_t x, const size_t y) const
        {
            if (x >= width)
                throw std::out_of_range("Given x out of range!");
            if (y >= height)
                throw std::out_of_range("Given y out of range!");
        }

    };

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

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

        Cube(const std::size_t newWidth, const std::size_t newHeight, const std::size_t newDepth):
            width(width),
            height(height),
            depth(newDepth),
            data(width * height * depth) { }

        Cube(const int newWidth, const int newHeight, const int newDepth, ElementT initVal):
            width(newWidth),
            height(newHeight),
            depth(newDepth),
            data(width * height * depth, initVal) {}

        Cube(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight, const std::size_t newDepth):
            width(newWidth),
            height(newHeight),
            depth(newDepth)
        {
            if (input.size() != newWidth * newHeight * newDepth)
            {
                throw std::runtime_error("Data input and the given size are mismatched!");
            }
            data = input;
        }    

        Cube(const std::vector<Image<ElementT>>& input)
        {
            width = input[0].getWidth();
            height = input[0].getHeight();
            depth = input.size();
            data.resize(width * height * depth);
            
            for (std::size_t z = 0; z < input.size(); ++z)
            {
                auto image = input[z];
                for (std::size_t y = 0; y < height; ++y)
                {
                    for (std::size_t x = 0; x < width; ++x)
                    {
                        data[z * width * height + y * width + x] = image.at(x, y);
                    }
                }
            }
            return;
        }

        constexpr ElementT& at(const unsigned int x, const unsigned int y, const unsigned int z)
        { 
            checkBoundary(x, y, z);
            return data[z * width * height + y * width + x];
        }

        constexpr ElementT const& at(const unsigned int x, const unsigned int y, const unsigned int z) const
        {
            checkBoundary(x, y, z);
            return data[z * width * height + y * width + x];
        }

        constexpr auto getSizeX() const noexcept
        {
            return width;
        }

        constexpr auto getSizeY() const noexcept
        {
            return height;
        }

        constexpr auto getSizeZ() const noexcept
        {
            return depth;
        }

        constexpr auto getWidth() const noexcept
        {
            return width;
        }

        constexpr auto getHeight() const noexcept
        {
            return height;
        }

        constexpr auto getDepth() const noexcept
        {
            return depth;
        }

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

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            for(std::size_t z = 0; z < depth; ++z)
            {
                for (std::size_t y = 0; y < height; ++y)
                {
                    for (std::size_t x = 0; x < width; ++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";
            return;
        }

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

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

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

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

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

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

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

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

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

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

        Cube<ElementT>& operator=(Cube<ElementT> const& input) = default;   //  Copy Assign

        Cube<ElementT>& operator=(Cube<ElementT>&& other) = default;        //  Move Assign

        Cube(const Cube<ElementT> &input) = default;                        //  Copy Constructor

        Cube(Cube<ElementT> &&input) = default;                             //  Move Constructor
        
    private:
        std::size_t width;
        std::size_t height;
        std::size_t depth;
        std::vector<ElementT> data;
        
        void checkBoundary(const size_t x, const size_t y, const size_t z) const
        {
            if (x >= width)
                throw std::out_of_range("Given x out of range!");
            if (y >= height)
                throw std::out_of_range("Given y out of range!");
            if (z >= depth)
                throw std::out_of_range("Given z out of range!");
        }
    };
}

//  cube_operations.h
namespace TinyDIP
{
    template<typename ElementT>
    constexpr bool is_width_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
    {
        return x.getSizeX() == y.getSizeX();
    }

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

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

    template<typename ElementT>
    constexpr bool is_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
    {
        return x.getSizeZ() == y.getSizeZ();
    }

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

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

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

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

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

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

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

    template<typename ElementT>
    constexpr void assert_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
    {
        assert(is_depth_same(x, y));
    }

    template<typename ElementT>
    constexpr void assert_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y, const Cube<ElementT>& z)
    {
        assert(is_depth_same(x, y));
        assert(is_depth_same(y, z));
    }

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

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

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

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

    template<typename ElementT>
    constexpr void check_depth_same(const Cube<ElementT>& x, const Cube<ElementT>& y)
    {
        if (!is_depth_same(x, y))
            throw std::runtime_error("Depth mismatched!");
    }

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

    //  voxelwiseOperation template function implementation
    template<std::size_t unwrap_level = 1, class... Args>
    constexpr static auto voxelwiseOperation(auto op, const Args&... inputs)
    {
        auto output = Cube(
            recursive_transform<unwrap_level>(
                [&](auto&& element1, auto&&... elements) 
                    {
                        return op(element1, elements...);
                    },
                inputs.getData()...),
            first_of(inputs...).getWidth(),
            first_of(inputs...).getHeight(),
            first_of(inputs...).getDepth());
        return output;
    }

    template<std::size_t unwrap_level = 1, class ExPo, class InputT>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static auto voxelwiseOperation(ExPo execution_policy, auto op, const Image<InputT>& input1)
    {
        auto output = Cube(
            recursive_transform<unwrap_level>(
                execution_policy,
                [&](auto&& element1) 
                    {
                        return op(element1);
                    },
                (input1.getData())),
            input1.getWidth(),
            input1.getHeight(),
            input1.getDepth());
        return output;
    }

    //  plus template function implementation
    template<class InputT>
    constexpr static Cube<InputT> plus(const Cube<InputT>& input1)
    {
        return input1;
    }

    template<class InputT, class... Args>
    constexpr static Cube<InputT> plus(const Cube<InputT>& input1, const Args&... inputs)
    {
        return voxelwiseOperation(std::plus<>{}, input1, plus(inputs...));
    }

    template<class InputT, class... Args>
    constexpr static auto plus(const std::vector<Cube<InputT>>& input1, const Args&... inputs)
    {
        return recursive_transform<1>(
            [](auto&& input1_element, auto&&... inputs_element)
            {
                return plus(input1_element, inputs_element...);
            }, input1, inputs...);
    }

    //  subtract template function implementation
    template<class InputT>
    constexpr static Cube<InputT> subtract(const Cube<InputT>& input1, const Cube<InputT>& input2)
    {
        check_size_same(input1, input2);
        return voxelwiseOperation(std::minus<>{}, input1, input2);
    }

    template<class InputT>
    constexpr static auto subtract(const std::vector<Cube<InputT>>& input1, const std::vector<Cube<InputT>>& input2)
    {
        assert(input1.size() == input2.size());
        return recursive_transform<1>(
            [](auto&& input1_element, auto&& input2_element)
            {
                return subtract(input1_element, input2_element);
            }, input1, input2);
    }

    //  multiplies template function implementation
    template<class InputT>
    constexpr static Cube<InputT> multiplies(const Cube<InputT>& input1, const Cube<InputT>& input2)
    {
        return voxelwiseOperation(std::multiplies<>{}, input1, input2);
    }

    template<class InputT, class TimesT>
    requires(std::floating_point<TimesT> || std::integral<TimesT>)
    constexpr static Cube<InputT> multiplies(const Cube<InputT>& input1, const TimesT times)
    {
        return multiplies(
            input1,
            Cube(input1.getWidth(), input1.getHeight(), input1.getDepth(), times)
        );
    }

    template<class InputT, class TimesT>
    requires(std::floating_point<TimesT> || std::integral<TimesT>)
    constexpr static Cube<InputT> multiplies(const TimesT times, const Cube<InputT>& input1)
    {
        return multiplies(input1, times);
    }

    template<class ExPo, class InputT>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static Cube<InputT> multiplies(ExPo execution_policy, const Cube<InputT>& input1, const Cube<InputT>& input2)
    {
        return voxelwiseOperation(execution_policy, std::multiplies<>{}, input1, input2);
    }

    template<class InputT, class... Args>
    constexpr static Cube<InputT> multiplies(const Cube<InputT>& input1, const Args&... inputs)
    {
        return voxelwiseOperation(std::multiplies<>{}, input1, multiplies(inputs...));
    }
}

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;
    
    //  recursive_invoke_result_t implementation
    template<std::size_t, typename, typename>
    struct recursive_invoke_result { };

    template<typename T, typename F>
    struct recursive_invoke_result<0, F, T> { using type = std::invoke_result_t<F, T>; };

    template<std::size_t unwrap_level, std::copy_constructible F, template<typename...> typename Container, typename... Ts>
    requires (std::ranges::input_range<Container<Ts...>> &&
            requires { typename recursive_invoke_result<unwrap_level - 1, F, std::ranges::range_value_t<Container<Ts...>>>::type; })
    struct recursive_invoke_result<unwrap_level, F, Container<Ts...>>
    {
        using type = Container<typename recursive_invoke_result<unwrap_level - 1, F, std::ranges::range_value_t<Container<Ts...>>>::type>;
    };

    template<std::size_t unwrap_level, std::copy_constructible F, typename T>
    using recursive_invoke_result_t = typename recursive_invoke_result<unwrap_level, F, T>::type;

    //  recursive_variadic_invoke_result_t implementation
    template<std::size_t, typename, typename, typename...>
    struct recursive_variadic_invoke_result { };

    template<std::copy_constructible F, class...Ts1, template<class...>class Container1, typename... Ts>
    struct recursive_variadic_invoke_result<1, F, Container1<Ts1...>, Ts...>
    {
        using type = Container1<std::invoke_result_t<F,
            std::ranges::range_value_t<Container1<Ts1...>>,
            std::ranges::range_value_t<Ts>...>>;
    };

    template<std::size_t unwrap_level, std::copy_constructible F, class...Ts1, template<class...>class Container1, typename... Ts>
    requires (  std::ranges::input_range<Container1<Ts1...>> &&
                requires { typename recursive_variadic_invoke_result<
                                        unwrap_level - 1,
                                        F,
                                        std::ranges::range_value_t<Container1<Ts1...>>,
                                        std::ranges::range_value_t<Ts>...>::type; })                //  The rest arguments are ranges
    struct recursive_variadic_invoke_result<unwrap_level, F, Container1<Ts1...>, Ts...>
    {
        using type = Container1<
            typename recursive_variadic_invoke_result<
            unwrap_level - 1,
            F,
            std::ranges::range_value_t<Container1<Ts1...>>,
            std::ranges::range_value_t<Ts>...
            >::type>;
    };

    template<std::size_t unwrap_level, std::copy_constructible F, typename T1, typename... Ts>
    using recursive_variadic_invoke_result_t = typename recursive_variadic_invoke_result<unwrap_level, F, T1, Ts...>::type;

    template<typename OutputIt, std::copy_constructible NAryOperation, typename InputIt, typename... InputIts>
    OutputIt transform(OutputIt d_first, NAryOperation op, InputIt first, InputIt last, InputIts... rest) {
        while (first != last) {
            *d_first++ = op(*first++, (*rest++)...);
        }
        return d_first;
    }

    //  recursive_transform for the multiple parameters cases (the version with unwrap_level)
    template<std::size_t unwrap_level = 1, std::copy_constructible F, class Arg1, class... Args>
    requires(unwrap_level <= recursive_depth<Arg1>())
    constexpr auto recursive_transform(const F& f, const Arg1& arg1, const Args&... args)
    {
        if constexpr (unwrap_level > 0)
        {
            recursive_variadic_invoke_result_t<unwrap_level, F, Arg1, Args...> output{};
            transform(
                std::inserter(output, std::ranges::end(output)),
                [&f](auto&& element1, auto&&... elements) { return recursive_transform<unwrap_level - 1>(f, element1, elements...); },
                std::ranges::cbegin(arg1),
                std::ranges::cend(arg1),
                std::ranges::cbegin(args)...
            );
            return output;
        }
        else
        {
            return std::invoke(f, arg1, args...);
        }
    }

    //  recursive_transform implementation (the version with unwrap_level, with execution policy)
    template<std::size_t unwrap_level = 1, class ExPo, class T, class F>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>> &&
              unwrap_level <= recursive_depth<T>())
    constexpr auto recursive_transform(ExPo execution_policy, const F& f, const T& input)
    {
        if constexpr (unwrap_level > 0)
        {
            recursive_invoke_result_t<unwrap_level, F, T> output{};
            output.resize(input.size());
            std::mutex mutex;
            std::transform(execution_policy, std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
                [&](auto&& element)
                {
                    std::lock_guard lock(mutex);
                    return recursive_transform<unwrap_level - 1>(execution_policy, f, element);
                });
            return output;
        }
        else
        {
            return std::invoke(f, input);
        }
    }
}

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<std::size_t unwrap_level = 1, class... Args>
    constexpr static auto pixelwiseOperation(auto op, const Args&... inputs)
    {
        auto output = Image(
            recursive_transform<unwrap_level>(
                [&](auto&& element1, auto&&... elements) 
                    {
                        return op(element1, elements...);
                    },
                inputs.getImageData()...),
            first_of(inputs...).getWidth(),
            first_of(inputs...).getHeight());
        return output;
    }

    template<std::size_t unwrap_level = 1, class ExPo, class InputT>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static auto pixelwiseOperation(ExPo execution_policy, auto op, const Image<InputT>& input1)
    {
        auto output = Image(
            recursive_transform<unwrap_level>(
                execution_policy,
                [&](auto&& element1) 
                    {
                        return op(element1);
                    },
                (input1.getImageData())),
            input1.getWidth(),
            input1.getHeight());
        return output;
    }

    template<typename T>
    requires(std::floating_point<T> || std::integral<T>)
    T normalDistribution1D(const T x, const T standard_deviation)
    {
        return std::exp(-x * x / (2 * standard_deviation * standard_deviation));
    }

    template<typename T>
    requires(std::floating_point<T> || std::integral<T>)
    T normalDistribution2D(const T xlocation, const T ylocation, const T standard_deviation)
    {
        return std::exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * std::numbers::pi * standard_deviation * standard_deviation);
    }

    //  multiple standard deviations
    template<class InputT>
    requires(std::floating_point<InputT> || std::integral<InputT>)
    constexpr static Image<InputT> gaussianFigure2D(
        const size_t xsize, const size_t ysize, 
        const size_t centerx, const size_t centery,
        const InputT standard_deviation_x, const InputT standard_deviation_y)
    {
        auto output = Image<InputT>(xsize, ysize);
        auto row_vector_x = Image<InputT>(xsize, 1);
        for (size_t x = 0; x < xsize; ++x)
        {
            row_vector_x.at(x, 0) = normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x);
        }

        auto row_vector_y = Image<InputT>(ysize, 1);
        for (size_t y = 0; y < ysize; ++y)
        {
            row_vector_y.at(y, 0) = normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
        }
        
        for (size_t y = 0; y < ysize; ++y)
        {
            for (size_t x = 0; x < xsize; ++x)
            {
                output.at(x, y) = row_vector_x.at(x, 0) * row_vector_y.at(y, 0);
            }
        }
        return output;
    }

    //  single standard deviation
    template<class InputT>
    requires(std::floating_point<InputT> || std::integral<InputT>)
    constexpr static Image<InputT> gaussianFigure2D(
        const size_t xsize, const size_t ysize,
        const size_t centerx, const size_t centery,
        const InputT standard_deviation)
    {
        return gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
    }

    //  multiple standard deviations
    template<class InputT>
    requires(std::floating_point<InputT> || std::integral<InputT>)
    constexpr static auto gaussianFigure3D(
        const size_t xsize, const size_t ysize, const size_t zsize,
        const size_t centerx, const size_t centery, const size_t centerz,
        const InputT standard_deviation_x, const InputT standard_deviation_y, const InputT standard_deviation_z)
    {
        auto output = std::vector<Image<InputT>>();
        output.reserve(zsize);
        auto gaussian_image2d = gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation_x, standard_deviation_y);
        for (size_t z = 0; z < zsize; ++z)
        {
            output.emplace_back(
                normalDistribution1D(static_cast<InputT>(z) - static_cast<InputT>(centerz), standard_deviation_z) * 
                gaussian_image2d
            );
        }
        return Cube(output);
    }

    template<class InputT>
    constexpr static Image<InputT> plus(const Image<InputT>& input1)
    {
        return input1;
    }

    template<class InputT, class... Args>
    constexpr static Image<InputT> plus(const Image<InputT>& input1, const Args&... inputs)
    {
        return pixelwiseOperation(std::plus<>{}, input1, plus(inputs...));
    }

    template<class InputT, class... Args>
    constexpr static auto plus(const std::vector<Image<InputT>>& input1, const Args&... inputs)
    {
        return recursive_transform<1>(
            [](auto&& input1_element, auto&&... inputs_element)
            {
                return plus(input1_element, inputs_element...);
            }, input1, inputs...);
    }

    template<class InputT>
    constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
    {
        check_size_same(input1, input2);
        return pixelwiseOperation(std::minus<>{}, input1, input2);
    }

    template<class InputT>
    constexpr static auto subtract(const std::vector<Image<InputT>>& input1, const std::vector<Image<InputT>>& input2)
    {
        assert(input1.size() == input2.size());
        return recursive_transform<1>(
            [](auto&& input1_element, auto&& input2_element)
            {
                return subtract(input1_element, input2_element);
            }, input1, input2);
    }

    template<class InputT = RGB>
    requires (std::same_as<InputT, RGB>)
    constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
    {
        check_size_same(input1, input2);
        Image<InputT> output(input1.getWidth(), input1.getHeight());
        for (std::size_t y = 0; y < input1.getHeight(); ++y)
        {
            for (std::size_t x = 0; x < input1.getWidth(); ++x)
            {
                for(std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                {
                    output.at(x, y).channels[channel_index] = 
                    std::clamp(
                        input1.at(x, y).channels[channel_index] - 
                        input2.at(x, y).channels[channel_index],
                        0,
                        255);
                }
            }
        }
        return output;
    }

    template<class InputT>
    constexpr static Image<InputT> multiplies(const Image<InputT>& input1, const Image<InputT>& input2)
    {
        return pixelwiseOperation(std::multiplies<>{}, input1, input2);
    }

    template<class InputT, class TimesT>
    requires(std::floating_point<TimesT> || std::integral<TimesT>)
    constexpr static Image<InputT> multiplies(const Image<InputT>& input1, const TimesT times)
    {
        return multiplies(
            input1,
            Image(input1.getWidth(), input1.getHeight(), times)
        );
    }

    template<class ExPo, class InputT>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static Image<InputT> multiplies(ExPo execution_policy, const Image<InputT>& input1, const Image<InputT>& input2)
    {
        return pixelwiseOperation(execution_policy, std::multiplies<>{}, input1, input2);
    }

    template<class InputT, class... Args>
    constexpr static Image<InputT> multiplies(const Image<InputT>& input1, const Args&... inputs)
    {
        return pixelwiseOperation(std::multiplies<>{}, input1, multiplies(inputs...));
    }

    template<class ExPo, class InputT, class... Args>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static Image<InputT> multiplies(ExPo execution_policy, const Image<InputT>& input1, const Args&... inputs)
    {
        return pixelwiseOperation(execution_policy, std::multiplies<>{}, input1, multiplies(inputs...));
    }

    template<class InputT, class... Args>
    constexpr static auto multiplies(const std::vector<Image<InputT>>& input1, const Args&... inputs)
    {
        return recursive_transform<1>(
            [](auto&& input1_element, auto&&... inputs_element)
            {
                return multiplies(input1_element, inputs_element...);
            }, input1, inputs...);
    }

    template<class InputT>
    constexpr static Image<InputT> divides(const Image<InputT>& input1, const Image<InputT>& input2)
    {
        return pixelwiseOperation(std::divides<>{}, input1, input2);
    }

    template<class InputT>
    constexpr static auto divides(const std::vector<Image<InputT>>& input1, const std::vector<Image<InputT>>& input2)
    {
        assert(input1.size() == input2.size());
        return recursive_transform<1>(
            [](auto&& input1_element, auto&& input2_element)
            {
                return divides(input1_element, input2_element);
            }, input1, input2);
    }

    template<class ExPo, class InputT>
    requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
    constexpr static Image<InputT> divides(ExPo execution_policy, const Image<InputT>& input1, const Image<InputT>& input2)
    {
        return pixelwiseOperation(execution_policy, std::divides<>{}, input1, input2);
    }
    
    template<class InputT>
    constexpr static Image<InputT> modulus(const Image<InputT>& input1, const Image<InputT>& input2)
    {
        return pixelwiseOperation(std::modulus<>{}, input1, input2);
    }

    template<class InputT>
    constexpr static Image<InputT> negate(const Image<InputT>& input1)
    {
        return pixelwiseOperation(std::negate<>{}, input1);
    }
}

template<typename T>
void gaussianFigure3DTest(const size_t size = 3)
{
    auto gaussian_image3d = TinyDIP::gaussianFigure3D(
        size,
        size,
        size,
        1, 1, 1, 
        static_cast<T>(3), static_cast<T>(3), static_cast<T>(3));
    gaussian_image3d.print();
    return;
}

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

0.846482  0.894839  0.846482  
0.894839  0.945959  0.894839  
0.846482  0.894839  0.846482  

0.894839  0.945959  0.894839  
0.945959  1  0.945959  
0.894839  0.945959  0.894839  

0.846482  0.894839  0.846482  
0.894839  0.945959  0.894839  
0.846482  0.894839  0.846482  


Computation finished at Tue Dec 26 06:37:57 2023
elapsed time: 6.6216e-05

Godbolt link is here.

All suggestions are welcome.

The summary information:

  • Which question it is a follow-up to?

    Three dimensional gaussian image generator in C++

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

    I am trying to create 3D data structure Cube in this post.

  • Why a new review is being asked for?

    Please review Cube template class implementation and all suggestions are welcome.

\$\endgroup\$

2 Answers 2

4
\$\begingroup\$

Lots of small issues

  • Missing checks for overflow of width * height * depth.
  • Inconsistent types for coordinates: one constructor uses std::size_t for the dimensions, another int, and at() takes unsigned ints.
  • Sometimes you use size_t instead of std::size_t.
  • Unnecessary copies, for example in:
    auto image = input[z];
    
  • Unnecessary return statement in the last constructor.
  • Why have both getSizeX() and getWidth()?
  • Why not use std::ranges::transform(), like so:
    std::ranges::transform(data, rhs.data, std::ranges::begin(data), std::plus<>{});
    
  • Remove the copy/move constructors/assignment operators; if they are all default you don't need to specify them at all (also see the rule of zero).
  • Do you really need all of is_foo_same(), assert_foo_same() and check_foo_same()?
  • Why do you even need assert_size_same() that takes three parameters? And even if you do, you only need two calls to assert_size_same() in it.
  • Multiplying a Cube by an ElementT still generates another Cube unnecessarily.

voxelwiseOperation() seems unnecessarily complex

It's weird to see recursive_transform() used here, but I guess if you want to support operations with more than two arguments and you couldn't use C++23's zip() yet, it might make sense. I would still have just written two simpler versions, one for unary operations and one for binary operations. However, even your version can be simplified a tiny bit this way:

template<std::size_t unwrap_level = 1, class First, class... Rest>
constexpr static auto voxelwiseOperation(auto op, const First& first, const Rest&... rest)
{
    return Cube(
        recursive_transform<unwrap_level>(
            [&](auto&&... elements) 
            {
                return op(elements...);
            },
            first.getData(),
            rest.getData()...
        ),
        first.getWidth(),
        first.getHeight(),
        first.getDepth()
    );
}

There is also the issue that recursive_transform() will generate a new std::vector<ElementT>, and since you don't have a constructor that can move from a std::vector<ElementT>, that whole vector will be copied, but that could be avoided.

What about hypercubes?

What if you need even higher-dimensional images? Adding a four-dimensional version would be a lot of work and lines of code. You could instead try to make Image take a template parameter that tells it the number of dimensions it has.

\$\endgroup\$
2
\$\begingroup\$

The word “cube” implies that the three dimensions are equal (x=y=z). I would have chosen a name like Image3D or VolumetricImage for this class.

But note how much code you now have duplicated from class Image. Code duplication is bad for many reasons, most importantly it increases the cost of maintenance. A 2D image is just a 3D image where the z dimension is 1. You could now rewrite the old Image class to a specialization of this one where you fill in z=1 for the sizes and z=0 when indexing. Most code won’t change…

But it would be even easier if you do as I suggested in the previous review: create an image class that has an arbitrary number of dimensions. Make sizes and coordinates an array type instead of passing each value separately, then it is easier to define the API.

You could look at how libraries like OpenCV, CImg, Vigra, or ITK implement the image class (besides my own DIPlib that I already pointed to in the previous review).

For example,

  • CImg defines an image to be 4D, where the 3rd and 4th dimensions default to size 1.
  • OpenCV and DIPlib define an image to have an arbitrary number of dimensions, and this dimensionality can be chosen at run time (using a vector-like type to specify the dimensions). But OpenCV adds functionality specific to 2D images (ie some API works differently if the image has 2 dimensions or if it has more).
  • Vigra and ITK define the dimensionality as a template parameter.

Any of these choices is fine, depending on your use cases.


Code like this:

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            for(std::size_t z = 0; z < depth; ++z)
            {
                for (std::size_t y = 0; y < height; ++y)
                {
                    for (std::size_t x = 0; x < width; ++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";
            return;
        }

looks fine on first read, but at() hides a lot of duplicate operations. Sure, printing is not efficient anyway, it doesn’t matter in this case. But I wanted to point out this being deceptive.

This is true for any image processing library: calling an at() function in a loop is suboptimal. Instead, add an iterator to your image class that visits all pixels. Or do like OpenCV, which has a function that returns a pointer to a specific line (e.g. line_at(y,z)). This function does some repeated work when called in a loop, but at least the inner loop is just incrementing a pointer.

The at() function itself can be a bit more efficient if you rewrite

data[z * width * height + y * width + x]

into

data[(z * height + y) * width + x]

Also, to reduce the code duplication between the const and non-const versions of the at() function, write a separate function that computes the linear index from the coordinates. This could be a private member function.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.