Return to Question

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      this->at(x, y) += rhs.at(x, y);
                  }
              }
              return *this;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          for (size_t y = 0; y < ysize; y++)
          {
              for (size_t x = 0; x < xsize; x++)
              {
                  output.at(x, y) = 
                      normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x) * 
                      normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
              }
          }
          return output;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• base_types.h: The base types

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <cmath>
  #include <cstdbool>
  #include <cstdio>
  #include <cstdlib>
  #include <string>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  typedef unsigned char BYTE;
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef BYTE GrayScale;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  #endif
  

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      this->at(x, y) += rhs.at(x, y);
                  }
              }
              return *this;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          auto row_vector_x = TinyDIP::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 = TinyDIP::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;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• base_types.h: The base types

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <cmath>
  #include <cstdbool>
  #include <cstdio>
  #include <cstdlib>
  #include <string>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  typedef unsigned char BYTE;
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef BYTE GrayScale;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  #endif
  

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          for (size_t y = 0; y < ysize; y++)
          {
              for (size_t x = 0; x < xsize; x++)
              {
                  output.at(x, y) = 
                      normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x) * 
                      normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
              }
          }
          return output;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• base_types.h: The base types

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <cmath>
  #include <cstdbool>
  #include <cstdio>
  #include <cstdlib>
  #include <string>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  typedef unsigned char BYTE;
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef BYTE GrayScale;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  #endif
  

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      this->at(x, y) += rhs.at(x, y);
                  }
              }
              return *this;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          for (size_t y = 0; y < ysize; y++)
          {
              for (size_t x = 0; x < xsize; x++)
              {
                  output.at(x, y) = 
                      normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x) * 
                      normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
              }
          }
          return output;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• base_types.h: The base types

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <cmath>
  #include <cstdbool>
  #include <cstdio>
  #include <cstdlib>
  #include <string>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  typedef unsigned char BYTE;
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef BYTE GrayScale;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  #endif
  

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          for (size_t y = 0; y < ysize; y++)
          {
              for (size_t x = 0; x < xsize; x++)
              {
                  output.at(x, y) = 
                      normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x) * 
                      normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
              }
          }
          return output;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• Which question it is a follow-up to?

  Two dimensional bicubic interpolation implementation in C++ and

  Two dimensional gaussian image generator in C

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

  Based on user673679's answer, another file image_operations.h is created for those non-member helper functions for image operations implementation. Moreover, I am attempting to build the two dimensional gaussian image generator functions in C++.

• Why a new review is being asked for?

  If there is any possible improvement, please let me know.

• Image template class implementation (image.h):

  /* Develop by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <numeric>
  #include <string>
  #include <type_traits>
  #include <variant>
  #include <vector>
  #include "basic_functions.h"
  #include "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image()
          {
          }
  
          Image(const size_t width, const size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const int width, const int height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height)
          {
              this->image_data = recursive_transform<1>(this->image_data, [initVal](ElementT element) { return initVal; });
              return;
          }
  
          Image(const std::vector<ElementT>& input, size_t newWidth, size_t newHeight)
          {
              this->width = newWidth;
              this->height = newHeight;
              this->image_data = recursive_transform<1>(input, [](ElementT element) { return element; });   //  Deep copy
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              this->height = input.size();
              this->width = input[0].size();
  
              for (auto& rows : input)
              {
                  this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y) { return this->image_data[y * width + x]; }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return this->image_data[y * width + x]; }
  
          constexpr size_t getWidth()
          {
              return this->width;
          }
  
          constexpr size_t getHeight()
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (size_t y = 0; y < this->height; y++)
              {
                  for (size_t x = 0; x < this->width; x++)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      std::cout << +this->at(x, y) << "\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          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
  
      private:
          size_t width;
          size_t height;
          std::vector<ElementT> image_data;
      };
  }
  #endif
  

• image_operations.h: non-member helper functions for image operations.

  /* Develop by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include "base_types.h"
  #include "image.h"
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          auto ratiox = (float)image.getWidth() / (float)width;
          auto ratioy = (float)image.getHeight() / (float)height;
  
          for (size_t y = 0; y < height; y++)
          {
              for (size_t x = 0; x < width; x++)
              {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  ElementT ndata[4 * 4];
                  for (int ndatay = -1; ndatay <= 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax <= 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
                              std::clamp(xMappingToOriginFloor + ndatax, 0.0f, image.getWidth() - 1.0f), 
                              std::clamp(yMappingToOriginFloor + ndatay, 0.0f, image.getHeight() - 1.0f));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      template<class ElementT, class InputT>
      constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
      {
          auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
          auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
          auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
          auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
          return std::clamp(cubicPolate( x1, x2, x3, x4, fracy ), 0.0f, 255.0f);
      }
  
      template<class InputT1, class InputT2>
      constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
      {
          auto A = (v3-v2)-(v0-v1);
          auto B = (v0-v1)-A;
          auto C = v2-v0;
          auto D = v1;
          return D + frac * (C + frac * (B + frac * A));
      }
  
      //  single standard deviation
      template<class 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 gaussianFigure2D2(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      //  multiple standard deviations
      template<class InputT>
      constexpr static Image<InputT> gaussianFigure2D2(
          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 = TinyDIP::Image<InputT>(xsize, ysize);
          for (size_t y = 0; y < ysize; y++)
          {
              for (size_t x = 0; x < xsize; x++)
              {
                  output.at(x, y) = 
                      normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x) * 
                      normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
              }
          }
          return output;
      }
  
      float normalDistribution1D(const float x, const float standard_deviation)
      {
          return expf(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      double normalDistribution1D(const double x, const double standard_deviation)
      {
          return exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      long double normalDistribution1D(const long double x, const long double standard_deviation)
      {
          return expl(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      float normalDistribution2D(const float xlocation, const float ylocation, const float standard_deviation)
      {
          return expf(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      double normalDistribution2D(const double xlocation, const double ylocation, const double standard_deviation)
      {
          return exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  }
  
  #endif
  

• basic_functions.h: The basic functions

  /* Develop by Jimmy Hu */
  
  #ifndef BasicFunctions_H
  #define BasicFunctions_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #include <chrono>
  #include <complex>
  #include <concepts>
  #include <deque>
  #include <execution>
  #include <exception>
  #include <functional>
  #include <iostream>
  #include <iterator>
  #include <list>
  #include <map>
  #include <mutex>
  #include <numeric>
  #include <optional>
  #include <ranges>
  #include <stdexcept>
  #include <string>
  #include <tuple>
  #include <type_traits>
  #include <utility>
  #include <variant>
  #include <vector>
  
  namespace TinyDIP
  {
      template<typename T>
      concept is_back_inserterable = requires(T x)
      {
          std::back_inserter(x);
      };
  
      template<typename T>
      concept is_inserterable = requires(T x)
      {
          std::inserter(x, std::ranges::end(x));
      };
  
      //  recursive_invoke_result_t implementation
      template<typename, typename>
      struct recursive_invoke_result { };
  
      template<typename T, std::invocable<T> F>
      struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };
  
      template<typename F, template<typename...> typename Container, typename... Ts>
      requires (
          !std::invocable<F, Container<Ts...>>&&
          std::ranges::input_range<Container<Ts...>>&&
          requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
          struct recursive_invoke_result<F, Container<Ts...>>
      {
          using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
      };
  
      template<typename F, typename T>
      using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;
  
      //  recursive_transform implementation (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class T, class F>
      constexpr auto recursive_transform(const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::ranges::transform(
                  std::ranges::cbegin(input),
                  std::ranges::cend(input),
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
              );
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      //  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>>)
      constexpr auto recursive_transform(ExPo execution_policy, const T& input, const F& f)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_invoke_result_t<F, T> output{};
              std::mutex mutex;
  
              //  Reference: https://en.cppreference.com/w/cpp/algorithm/for_each
              std::for_each(execution_policy, input.cbegin(), input.cend(),
                  [&](auto&& element)
                  {
                      auto result = recursive_transform<unwrap_level - 1>(execution_policy, element, f);
                      std::lock_guard lock(mutex);
                      output.emplace_back(std::move(result));
                  }
              );
  
              return output;
          }
          else
          {
              return f(input);
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_vector_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_vector_generator<dim - 1>(input, times);
              std::vector<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, std::size_t times, class T>
      constexpr auto n_dim_array_generator(T input)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_array_generator<dim - 1, times>(input);
              std::array<decltype(element), times> output;
              std::fill(std::ranges::begin(output), std::ranges::end(output), element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_deque_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_deque_generator<dim - 1>(input, times);
              std::deque<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, class T>
      constexpr auto n_dim_list_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              auto element = n_dim_list_generator<dim - 1>(input, times);
              std::list<decltype(element)> output(times, element);
              return output;
          }
      }
  
      template<std::size_t dim, template<class...> class Container = std::vector, class T>
      constexpr auto n_dim_container_generator(T input, std::size_t times)
      {
          if constexpr (dim == 0)
          {
              return input;
          }
          else
          {
              return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
          }
      }
  }
  
  #endif
  

• base_types.h: The base types

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <cmath>
  #include <cstdbool>
  #include <cstdio>
  #include <cstdlib>
  #include <string>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  typedef unsigned char BYTE;
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef BYTE GrayScale;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  #endif
  

• Which question it is a follow-up to?

  Two dimensional bicubic interpolation implementation in C++ and

  Two dimensional gaussian image generator in C

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

  Based on user673679's answer, another file image_operations.h is created for those non-member helper functions for image operations implementation. Moreover, I am attempting to build the two dimensional gaussian image generator functions in C++.

• Why a new review is being asked for?

  If there is any possible improvement, please let me know.

Reference

• Multivariate normal distribution

  https://en.wikipedia.org/wiki/Multivariate_normal_distribution