Three dimensional gaussian image generator in C++

Question

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.

Answer 1

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.

Add a vector type to store coordinates

You have to pass 9 parameters to gaussianFigure3D(). That's quite a lot. The more parameters a function has, the easier it is to make mistakes and pass them in the wrong order. It would be much better if you had a type that stored 3D coordinates:

template<typename T>
struct Vec3D {
    T x;
    T y;
    T z;
};

template<typename InputT>
requires(std::floating_point<InputT> || std::integral<InputT>)
constexpr static auto guassianfigure3D(
    Vec3D<size_t> size,
    Vec3D<size_t> center,
    Vec3D<InputT> standard_deviation)
{
    …
};

Make it more generic

You already had normalDistribution1D() and guassianFigure2D(), and now you added guassianFigure3D(). What if you want to make a four-dimensional Gaussian? A good excercise is to see if you can make the number of dimensions a parameter somehow. You already demonstrate in guassianFigure3D() that the way to construct a higher-dimensional Gaussian is to take one that has one dimension less, and then just do a Cartesian product with a normalDistribution1D(). So you have your generic case and your base case, and can thus write a recursive template function to create Gaussian images of any dimension.

You'll then also want to make Image a template that has a parameter for the number of dimensions it has. Alternatively, you can choose to make the 2D Gaussian the base case, and make higher dimension images just vectors of vectors of … of Image.

C++23 introduced std::mdspan. This might be a very helpful tool here.

Add a scalar multiplication operator

Your Image class can only multiply one image by another. But in guassianFigure3D(), you actually only need to multiply guassian_image2d by a scalar. You work around this by creating a full 2D Image with a uniform value everywhere, but that is just wasting memory and hinders optimization. Consider adding an operator*(ElementT rhs) to Image, then the loop in gaussianFigure3D() can look like:

for (size_t z = 0; z < zsize; ++z)
{
    output.push_back(
        gaussian_image2d * normalDistribution1D(z - centerz, standard_deviation_z)
    );
}

Reserve space for output

You know zsize, so you can already reserve enough space for output. This prevents unnecessary copies and moves of Images.

Answer 2

I agree with everything G. Sliepen said in their answer. Especially the advice of creating an actual 3D image rather than an array of 2D images (ideally you'd create an image object that can have an arbitrary number of dimensions), and of using an array type to hold coordinates and sizes.

I would like to add a question: what is the use case of this function? I imagine you could want to add a Gaussian blob to an existing image. With your function, you’d create a new image with the Gaussian blob, and then add it to some other image. Why not take an existing image as input, and add the Gaussian to it?

Another use case could be to create a kernel for a Gaussian filter. But a Gaussian filter is separable and should be implemented as such, not as a convolution with a 3D kernel. So this is not a good use case.

You could add some efficiency by computing only the Gaussian where it is not approximately zero. Typically we crop it at 3 sigma, but depending on the application it could be a different size. Pixels outside that (hyper-)cuboid you’d set to zero, saving you some quite expensive computation. This is of course most impactful when drawing small Gaussian blobs in a large image.

I encourage you to take a look at how I implemented an arbitrary-dimensional image in DIPlib: https://github.com/DIPlib/diplib/blob/master/include/diplib/library/image.h

The DIPlib function that adds a Gaussian blob to an existing image has the following signature:

void DrawBandlimitedPoint(
      Image& out,
      FloatArray origin,
      Image::Pixel const& value = { 1 },
      FloatArray sigmas = { 1.0 },
      dfloat truncation = 3.0
);