• pixelwiseOperation template function implementation: based on recursive_transform

  template<typename Op, class InputT, class... Args>
  constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
  {
      Image<InputT> output(
          recursive_transform<1>(
              [&](auto&& element1, auto&&... elements) 
                  {
                      auto result = op(element1, elements...);
                      return static_cast<InputT>(std::clamp(
                          result,
                          static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                          static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                  },
              (input1.getImageData()),
              (inputs.getImageData())...),
          input1.getWidth(),
          input1.getHeight());
      return output;
  }
  

• image_operations.h: The file contains pixelwiseOperation template function and other image processing functions

  /* Developed by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include <string>
  #include "base_types.h"
  #include "image.h"
  
  //  Reference: https://stackoverflow.com/a/26065433/6667035
  #ifndef M_PI
      #define M_PI 3.14159265358979323846
  #endif
  
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<typename T>
      T normalDistribution1D(const T x, const T standard_deviation)
      {
          return std::exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      template<typename T>
      T normalDistribution2D(const T xlocation, const T ylocation, const T standard_deviation)
      {
          return std::exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      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));
      }
  
      template<class InputT = float, class ElementT>
      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 ), 
                      static_cast<InputT>(std::numeric_limits<ElementT>::min()),
                      static_cast<InputT>(std::numeric_limits<ElementT>::max()));
      }
  
      template<class FloatingType = float, class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT>& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          //  get used to the C++ way of casting
          auto ratiox = static_cast<FloatingType>(image.getWidth()) / static_cast<FloatingType>(width);
          auto ratioy = static_cast<FloatingType>(image.getHeight()) / static_cast<FloatingType>(height);
  
          for (size_t y = 0; y < height; ++y)
          {
              for (size_t x = 0; x < width; ++x)
              {
                  FloatingType xMappingToOrigin = static_cast<FloatingType>(x) * ratiox;
                  FloatingType yMappingToOrigin = static_cast<FloatingType>(y) * ratioy;
                  FloatingType xMappingToOriginFloor = std::floor(xMappingToOrigin);
                  FloatingType yMappingToOriginFloor = std::floor(yMappingToOrigin);
                  FloatingType xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  FloatingType 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, static_cast<FloatingType>(0), image.getWidth() - static_cast<FloatingType>(1)), 
                              std::clamp(yMappingToOriginFloor + ndatay, static_cast<FloatingType>(0), image.getHeight() - static_cast<FloatingType>(1)));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      //  multiple standard deviations
      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_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;
      }
  
      //  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 gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      template<typename Op, class InputT, class... Args>
      constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
      {
          Image<InputT> output(
              recursive_transform<1>(
                  [&](auto&& element1, auto&&... elements) 
                      {
                          auto result = op(element1, elements...);
                          return static_cast<InputT>(std::clamp(
                              result,
                              static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                              static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                      },
                  (input1.getImageData()),
                  (inputs.getImageData())...),
              input1.getWidth(),
              input1.getHeight());
          return output;
      }
  
      template<class InputT>
      constexpr static Image<InputT> plus(const Image<InputT>& input1)
      {
          return input1;
      }
  
      template<class InputT, class... Args>
      constexpr static Image<InputT> plus(const Image<InputT>& input1, const Args&... inputs)
      {
          return TinyDIP::pixelwiseOperation(std::plus<>{}, input1, plus(inputs...));
      }
  
      template<class InputT>
      constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
      {
          return TinyDIP::pixelwiseOperation(std::minus<>{}, input1, input2);
      }
  
      template<class InputT = RGB>
      requires (std::same_as<InputT, RGB>)
      constexpr static Image<InputT> subtract(Image<InputT>& input1, Image<InputT>& input2)
      {
          assert(input1.getWidth() == input2.getWidth());
          assert(input1.getHeight() == input2.getHeight());
          Image<InputT> output(input1.getWidth(), input1.getHeight());
          for (std::size_t y = 0; y < input1.getHeight(); ++y)
          {
              for (std::size_t x = 0; x < input1.getWidth(); ++x)
              {
                  for(std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                  {
                      output.at(x, y).channels[channel_index] = 
                      std::clamp(
                          input1.at(x, y).channels[channel_index] - 
                          input2.at(x, y).channels[channel_index],
                          0,
                          255);
                  }
              }
          }
          return output;
      }
  }
  
  #endif
  

• image.h: The file contains the definition of Image class.

  /* Developed by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #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 "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image() = default;
  
          Image(const std::size_t width, const std::size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const std::size_t width, const std::size_t height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height, initVal) {}
  
          Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight):
              width(newWidth),
              height(newHeight)
          {
              assert(input.size() == newWidth * newHeight);
              this->image_data = input;   //  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));    //  flatten
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y)
          { 
              checkBoundary(x, y);
              return this->image_data[y * width + x];
          }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
          {
              checkBoundary(x, y);
              return this->image_data[y * width + x];
          }
  
          constexpr std::size_t getWidth() const
          {
              return this->width;
          }
  
          constexpr std::size_t getHeight() const
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (std::size_t y = 0; y < this->height; ++y)
              {
                  for (std::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;
          }
  
          //  Enable this function if ElementT = RGB
          void print() requires(std::same_as<ElementT, RGB>)
          {
              for (std::size_t y = 0; y < this->height; ++y)
              {
                  for (std::size_t x = 0; x < this->width; ++x)
                  {
                      std::cout << "( ";
                      for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                      {
                          //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                          std::cout << +this->at(x, y).channels[channel_index] << "\t";
                      }
                      std::cout << ")\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::plus<>{});
              return *this;
          }
  
          Image<ElementT>& operator-=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::minus<>{});
              return *this;
          }
  
          Image<ElementT>& operator*=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::multiplies<>{});
              return *this;
          }
  
          Image<ElementT>& operator/=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::divides<>{});
              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;
  
          void checkBoundary(const size_t x, const size_t y)
          {
              assert(x < width);
              assert(y < height);
          }
      };
  }
  #endif
  

• base_types.h: The base types declaration

  /* Developed by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #define _USE_MATH_DEFINES
  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string>
  #include <utility>
  
  constexpr int MAX_PATH = 256;
  #define FILE_ROOT_PATH "./"
  
  using BYTE = unsigned char;
  
  struct RGB
  {
      unsigned char channels[3];
  };
  
  using GrayScale = BYTE;
  
  struct HSV
  {
      double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  };
  
  struct BMPIMAGE
  {
      char FILENAME[MAX_PATH];
  
      unsigned int XSIZE;
      unsigned int YSIZE;
      unsigned char FILLINGBYTE;
      unsigned char *IMAGE_DATA;
  };
  #endif
  

• basic_functions.h: The file contains the definition of recursive_transform

  /* Developed 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
  {
      //  recursive_variadic_invoke_result_t implementation
      template<std::size_t, typename, typename, typename...>
      struct recursive_variadic_invoke_result { };
  
      template<typename F, class...Ts1, template<class...>class Container1, typename... Ts>
      struct recursive_variadic_invoke_result<1, F, Container1<Ts1...>, Ts...>
      {
          using type = Container1<std::invoke_result_t<F,
              std::ranges::range_value_t<Container1<Ts1...>>,
              std::ranges::range_value_t<Ts>...>>;
      };
  
      template<std::size_t unwrap_level, typename F, class...Ts1, template<class...>class Container1, typename... Ts>
      requires (  std::ranges::input_range<Container1<Ts1...>> &&
                  requires { typename recursive_variadic_invoke_result<
                                          unwrap_level - 1,
                                          F,
                                          std::ranges::range_value_t<Container1<Ts1...>>,
                                          std::ranges::range_value_t<Ts>...>::type; })                //  The rest arguments are ranges
      struct recursive_variadic_invoke_result<unwrap_level, F, Container1<Ts1...>, Ts...>
      {
          using type = Container1<
              typename recursive_variadic_invoke_result<
              unwrap_level - 1,
              F,
              std::ranges::range_value_t<Container1<Ts1...>>,
              std::ranges::range_value_t<Ts>...
              >::type>;
      };
  
      template<std::size_t unwrap_level, typename F, typename T1, typename... Ts>
      using recursive_variadic_invoke_result_t = typename recursive_variadic_invoke_result<unwrap_level, F, T1, Ts...>::type;
  
      template<typename OutputIt, typename NAryOperation, typename InputIt, typename... InputIts>
      OutputIt transform(OutputIt d_first, NAryOperation op, InputIt first, InputIt last, InputIts... rest) {
          while (first != last) {
              *d_first++ = op(*first++, (*rest++)...);
          }
          return d_first;
      }
  
      //  recursive_transform for the multiple parameters cases (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class F, class Arg1, class... Args>
      constexpr auto recursive_transform(const F& f, const Arg1& arg1, const Args&... args)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_variadic_invoke_result_t<unwrap_level, F, Arg1, Args...> output{};
              transform(
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element1, auto&&... elements) { return recursive_transform<unwrap_level - 1>(f, element1, elements...); },
                  std::ranges::cbegin(arg1),
                  std::ranges::cend(arg1),
                  std::ranges::cbegin(args)...
              );
              return output;
          }
          else
          {
              return f(arg1, args...);
          }
      }
  }
  
  #endif
  

• pixelwiseOperation template function implementation: based on recursive_transform

  template<typename Op, class InputT, class... Args>
  constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
  {
      Image<InputT> output(
          recursive_transform<1>(
              [&](auto&& element1, auto&&... elements) 
                  {
                      auto result = op(element1, elements...);
                      return static_cast<InputT>(std::clamp(
                          result,
                          static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                          static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                  },
              (input1.getImageData()),
              (inputs.getImageData())...),
          input1.getWidth(),
          input1.getHeight());
      return output;
  }
  

• image_operations.h: The file contains pixelwiseOperation template function and other image processing functions

  /* Developed by Jimmy Hu */
  
  #ifndef ImageOperations_H
  #define ImageOperations_H
  
  #include <string>
  #include "base_types.h"
  #include "image.h"
  
  //  Reference: https://stackoverflow.com/a/26065433/6667035
  #ifndef M_PI
      #define M_PI 3.14159265358979323846
  #endif
  
  
  namespace TinyDIP
  {
      // Forward Declaration class Image
      template <typename ElementT>
      class Image;
  
      template<typename T>
      T normalDistribution1D(const T x, const T standard_deviation)
      {
          return std::exp(-x * x / (2 * standard_deviation * standard_deviation));
      }
  
      template<typename T>
      T normalDistribution2D(const T xlocation, const T ylocation, const T standard_deviation)
      {
          return std::exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
      }
  
      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));
      }
  
      template<class InputT = float, class ElementT>
      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 ), 
                      static_cast<InputT>(std::numeric_limits<ElementT>::min()),
                      static_cast<InputT>(std::numeric_limits<ElementT>::max()));
      }
  
      template<class FloatingType = float, class ElementT>
      Image<ElementT> copyResizeBicubic(Image<ElementT>& image, size_t width, size_t height)
      {
          auto output = Image<ElementT>(width, height);
          //  get used to the C++ way of casting
          auto ratiox = static_cast<FloatingType>(image.getWidth()) / static_cast<FloatingType>(width);
          auto ratioy = static_cast<FloatingType>(image.getHeight()) / static_cast<FloatingType>(height);
  
          for (size_t y = 0; y < height; ++y)
          {
              for (size_t x = 0; x < width; ++x)
              {
                  FloatingType xMappingToOrigin = static_cast<FloatingType>(x) * ratiox;
                  FloatingType yMappingToOrigin = static_cast<FloatingType>(y) * ratioy;
                  FloatingType xMappingToOriginFloor = std::floor(xMappingToOrigin);
                  FloatingType yMappingToOriginFloor = std::floor(yMappingToOrigin);
                  FloatingType xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  FloatingType 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, static_cast<FloatingType>(0), image.getWidth() - static_cast<FloatingType>(1)), 
                              std::clamp(yMappingToOriginFloor + ndatay, static_cast<FloatingType>(0), image.getHeight() - static_cast<FloatingType>(1)));
                      }
  
                  }
                  output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
              }
          }
          return output;
      }
  
      //  multiple standard deviations
      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_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;
      }
  
      //  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 gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
      }
  
      template<typename Op, class InputT, class... Args>
      constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
      {
          Image<InputT> output(
              recursive_transform<1>(
                  [&](auto&& element1, auto&&... elements) 
                      {
                          auto result = op(element1, elements...);
                          return static_cast<InputT>(std::clamp(
                              result,
                              static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                              static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                      },
                  (input1.getImageData()),
                  (inputs.getImageData())...),
              input1.getWidth(),
              input1.getHeight());
          return output;
      }
  
      template<class InputT>
      constexpr static Image<InputT> plus(const Image<InputT>& input1)
      {
          return input1;
      }
  
      template<class InputT, class... Args>
      constexpr static Image<InputT> plus(const Image<InputT>& input1, const Args&... inputs)
      {
          return TinyDIP::pixelwiseOperation(std::plus<>{}, input1, plus(inputs...));
      }
  
      template<class InputT>
      constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
      {
          assert(input1.getWidth() == input2.getWidth());
          assert(input1.getHeight() == input2.getHeight());
          return TinyDIP::pixelwiseOperation(std::minus<>{}, input1, input2);
      }
  
      template<class InputT = RGB>
      requires (std::same_as<InputT, RGB>)
      constexpr static Image<InputT> subtract(Image<InputT>& input1, Image<InputT>& input2)
      {
          assert(input1.getWidth() == input2.getWidth());
          assert(input1.getHeight() == input2.getHeight());
          Image<InputT> output(input1.getWidth(), input1.getHeight());
          for (std::size_t y = 0; y < input1.getHeight(); ++y)
          {
              for (std::size_t x = 0; x < input1.getWidth(); ++x)
              {
                  for(std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                  {
                      output.at(x, y).channels[channel_index] = 
                      std::clamp(
                          input1.at(x, y).channels[channel_index] - 
                          input2.at(x, y).channels[channel_index],
                          0,
                          255);
                  }
              }
          }
          return output;
      }
  }
  
  #endif
  

• image.h: The file contains the definition of Image class.

  /* Developed by Jimmy Hu */
  
  #ifndef Image_H
  #define Image_H
  
  #include <algorithm>
  #include <array>
  #include <cassert>
  #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 "image_operations.h"
  
  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image() = default;
  
          Image(const std::size_t width, const std::size_t height):
              width(width),
              height(height),
              image_data(width * height) { }
  
          Image(const std::size_t width, const std::size_t height, const ElementT initVal):
              width(width),
              height(height),
              image_data(width * height, initVal) {}
  
          Image(const std::vector<ElementT>& input, std::size_t newWidth, std::size_t newHeight):
              width(newWidth),
              height(newHeight)
          {
              assert(input.size() == newWidth * newHeight);
              this->image_data = input;   //  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));    //  flatten
              }
              return;
          }
  
          constexpr ElementT& at(const unsigned int x, const unsigned int y)
          { 
              checkBoundary(x, y);
              return this->image_data[y * width + x];
          }
  
          constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
          {
              checkBoundary(x, y);
              return this->image_data[y * width + x];
          }
  
          constexpr std::size_t getWidth() const
          {
              return this->width;
          }
  
          constexpr std::size_t getHeight() const
          {
              return this->height;
          }
  
          std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
  
          void print()
          {
              for (std::size_t y = 0; y < this->height; ++y)
              {
                  for (std::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;
          }
  
          //  Enable this function if ElementT = RGB
          void print() requires(std::same_as<ElementT, RGB>)
          {
              for (std::size_t y = 0; y < this->height; ++y)
              {
                  for (std::size_t x = 0; x < this->width; ++x)
                  {
                      std::cout << "( ";
                      for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                      {
                          //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                          std::cout << +this->at(x, y).channels[channel_index] << "\t";
                      }
                      std::cout << ")\t";
                  }
                  std::cout << "\n";
              }
              std::cout << "\n";
              return;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::plus<>{});
              return *this;
          }
  
          Image<ElementT>& operator-=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::minus<>{});
              return *this;
          }
  
          Image<ElementT>& operator*=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::multiplies<>{});
              return *this;
          }
  
          Image<ElementT>& operator/=(const Image<ElementT>& rhs)
          {
              assert(rhs.width == this->width);
              assert(rhs.height == this->height);
              std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                     image_data.begin(), std::divides<>{});
              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;
  
          void checkBoundary(const size_t x, const size_t y)
          {
              assert(x < width);
              assert(y < height);
          }
      };
  }
  #endif
  

• base_types.h: The base types declaration

  /* Developed by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #define _USE_MATH_DEFINES
  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string>
  #include <utility>
  
  constexpr int MAX_PATH = 256;
  #define FILE_ROOT_PATH "./"
  
  using BYTE = unsigned char;
  
  struct RGB
  {
      unsigned char channels[3];
  };
  
  using GrayScale = BYTE;
  
  struct HSV
  {
      double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  };
  
  struct BMPIMAGE
  {
      char FILENAME[MAX_PATH];
  
      unsigned int XSIZE;
      unsigned int YSIZE;
      unsigned char FILLINGBYTE;
      unsigned char *IMAGE_DATA;
  };
  #endif
  

• basic_functions.h: The file contains the definition of recursive_transform

  /* Developed 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
  {
      //  recursive_variadic_invoke_result_t implementation
      template<std::size_t, typename, typename, typename...>
      struct recursive_variadic_invoke_result { };
  
      template<typename F, class...Ts1, template<class...>class Container1, typename... Ts>
      struct recursive_variadic_invoke_result<1, F, Container1<Ts1...>, Ts...>
      {
          using type = Container1<std::invoke_result_t<F,
              std::ranges::range_value_t<Container1<Ts1...>>,
              std::ranges::range_value_t<Ts>...>>;
      };
  
      template<std::size_t unwrap_level, typename F, class...Ts1, template<class...>class Container1, typename... Ts>
      requires (  std::ranges::input_range<Container1<Ts1...>> &&
                  requires { typename recursive_variadic_invoke_result<
                                          unwrap_level - 1,
                                          F,
                                          std::ranges::range_value_t<Container1<Ts1...>>,
                                          std::ranges::range_value_t<Ts>...>::type; })                //  The rest arguments are ranges
      struct recursive_variadic_invoke_result<unwrap_level, F, Container1<Ts1...>, Ts...>
      {
          using type = Container1<
              typename recursive_variadic_invoke_result<
              unwrap_level - 1,
              F,
              std::ranges::range_value_t<Container1<Ts1...>>,
              std::ranges::range_value_t<Ts>...
              >::type>;
      };
  
      template<std::size_t unwrap_level, typename F, typename T1, typename... Ts>
      using recursive_variadic_invoke_result_t = typename recursive_variadic_invoke_result<unwrap_level, F, T1, Ts...>::type;
  
      template<typename OutputIt, typename NAryOperation, typename InputIt, typename... InputIts>
      OutputIt transform(OutputIt d_first, NAryOperation op, InputIt first, InputIt last, InputIts... rest) {
          while (first != last) {
              *d_first++ = op(*first++, (*rest++)...);
          }
          return d_first;
      }
  
      //  recursive_transform for the multiple parameters cases (the version with unwrap_level)
      template<std::size_t unwrap_level = 1, class F, class Arg1, class... Args>
      constexpr auto recursive_transform(const F& f, const Arg1& arg1, const Args&... args)
      {
          if constexpr (unwrap_level > 0)
          {
              recursive_variadic_invoke_result_t<unwrap_level, F, Arg1, Args...> output{};
              transform(
                  std::inserter(output, std::ranges::end(output)),
                  [&f](auto&& element1, auto&&... elements) { return recursive_transform<unwrap_level - 1>(f, element1, elements...); },
                  std::ranges::cbegin(arg1),
                  std::ranges::cend(arg1),
                  std::ranges::cbegin(args)...
              );
              return output;
          }
          else
          {
              return f(arg1, args...);
          }
      }
  }
  
  #endif
  

Test platform

MacOS: g++-11 (Homebrew GCC 11.1.0_1) 11.1.0