Skip to main content
Became Hot Network Question
deprettify program output
Source Link
Toby Speight
  • 88.3k
  • 14
  • 104
  • 327
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 Sat Dec 23 13:00:22 2023
elapsed time: 0.00157651
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 Sat Dec 23 13:00:22 2023
elapsed time: 0.00157651
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 Sat Dec 23 13:00:22 2023
elapsed time: 0.00157651
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 Sat Dec 23 13:00:22 2023
elapsed time: 0.00157651
Source Link
JimmyHu
  • 7.5k
  • 2
  • 11
  • 48

Three dimensional gaussian image generator in C++

This is a follow-up question for Two dimensional gaussian image generator in C++. Besides the two dimensional case, I am trying to implement three dimensional gaussian image generator which with additional zsize, centerz and standard_deviation_z parameters.

The experimental implementation

  • gaussianFigure3D Template Function Implementation

    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>>();
        auto gaussian_image2d = gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation_x, standard_deviation_y);
        for (size_t z = 0; z < zsize; ++z)
        {
            output.push_back(
                multiplies(gaussian_image2d,
                Image(xsize, ysize, normalDistribution1D(static_cast<InputT>(z) - static_cast<InputT>(centerz), standard_deviation_z)))
            );
        }
        return output;
    }
    
  • Image class

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

Full Testing Code

The full testing code:

//  Three dimensional gaussian image generator in C++
//  Developed by Jimmy Hu

#include <algorithm>
#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;
        }

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

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<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>>();
        auto gaussian_image2d = gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation_x, standard_deviation_y);
        for (size_t z = 0; z < zsize; ++z)
        {
            output.push_back(
                multiplies(gaussian_image2d,
                Image(xsize, ysize, normalDistribution1D(static_cast<InputT>(z) - static_cast<InputT>(centerz), standard_deviation_z)))
            );
        }
        return 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 Image<InputT> 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 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));
    for(auto each_plane : gaussian_image3d)
    {
        each_plane.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 Sat Dec 23 13:00:22 2023
elapsed time: 0.00157651

Godbolt link is here.

All suggestions are welcome.

The summary information:

  • Which question it is a follow-up to?

    Two dimensional gaussian image generator in C++

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

    I am trying to implement gaussianFigure3D template function in this post.

  • Why a new review is being asked for?

    Please review gaussianFigure3D template function implementation and all suggestions are welcome.