Two dimensional gaussian image generator in C++

Question

This is a follow-up question for Two dimensional bicubic interpolation implementation in C++ and Two dimensional gaussian image generator in C. Based on user673679's answer, another file image_operations.h is created for those non-member helper functions for image operations implementation. Moreover, the two dimensional gaussian image generator gaussianFigure2D and gaussianFigure2D2 is implemented as below.

• 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
  

The full testing code

/* Develop by Jimmy Hu */

#include "base_types.h"
#include "basic_functions.h"
#include "image.h"
#include "image_operations.h"

void gaussianFigure2DTest();

int main()
{
    gaussianFigure2DTest();
    return 0;
}

void gaussianFigure2DTest()
{
    auto image1 = TinyDIP::gaussianFigure2D2(10, 10, 5, 5, 3.0, 1.0);
    image1.print();

    auto image2 = TinyDIP::gaussianFigure2D(10, 10, 5, 5, 3.0);
    image2.print();
    return;
}

The output of the testing code above:

9.29249e-07     1.53207e-06     2.26033e-06     2.98407e-06     3.52526e-06     3.72665e-06     3.52526e-06     2.98407e-06     2.26033e-06     1.53207e-06
8.36483e-05     0.000137913     0.000203468     0.000268617     0.000317334     0.000335463     0.000317334     0.000268617     0.000203468     0.000137913
0.00277005      0.00456705      0.00673795      0.00889539      0.0105087       0.011109        0.0105087       0.00889539      0.00673795      0.00456705
0.0337462       0.055638        0.082085        0.108368        0.128022        0.135335        0.128022        0.108368        0.082085        0.055638
0.15124 0.249352        0.367879        0.485672        0.573753        0.606531        0.573753        0.485672        0.367879        0.249352
0.249352        0.411112        0.606531        0.800737        0.945959        1       0.945959        0.800737        0.606531        0.411112
0.15124 0.249352        0.367879        0.485672        0.573753        0.606531        0.573753        0.485672        0.367879        0.249352
0.0337462       0.055638        0.082085        0.108368        0.128022        0.135335        0.128022        0.108368        0.082085        0.055638
0.00277005      0.00456705      0.00673795      0.00889539      0.0105087       0.011109        0.0105087       0.00889539      0.00673795      0.00456705
8.36483e-05     0.000137913     0.000203468     0.000268617     0.000317334     0.000335463     0.000317334     0.000268617     0.000203468     0.000137913

0.0621765       0.102512        0.15124 0.199666        0.235877        0.249352        0.235877        0.199666        0.15124 0.102512
0.102512        0.169013        0.249352        0.329193        0.388896        0.411112        0.388896        0.329193        0.249352        0.169013
0.15124 0.249352        0.367879        0.485672        0.573753        0.606531        0.573753        0.485672        0.367879        0.249352
0.199666        0.329193        0.485672        0.64118 0.757465        0.800737        0.757465        0.64118 0.485672        0.329193
0.235877        0.388896        0.573753        0.757465        0.894839        0.945959        0.894839        0.757465        0.573753        0.388896
0.249352        0.411112        0.606531        0.800737        0.945959        1       0.945959        0.800737        0.606531        0.411112
0.235877        0.388896        0.573753        0.757465        0.894839        0.945959        0.894839        0.757465        0.573753        0.388896
0.199666        0.329193        0.485672        0.64118 0.757465        0.800737        0.757465        0.64118 0.485672        0.329193
0.15124 0.249352        0.367879        0.485672        0.573753        0.606531        0.573753        0.485672        0.367879        0.249352
0.102512        0.169013        0.249352        0.329193        0.388896        0.411112        0.388896        0.329193        0.249352        0.169013

All suggestions are welcome.

The summary information:

• 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

Answer 1

You are obviously on a fantastic quest to understand how to implement image processing operations in all sorts of languages. The tag reinventing-the-wheel tells me that you know that you should use an established image processing library as the basis of any serious work.

Here are some comments on your C++ implementation.

image.h

Proper English spelling would be "Developed by Jimmy Hu".

This file includes headers that you don't use. You should only include headers that you actually need.

        Image()
        {
        }

This should be

        Image() = default;

This constructor initializes the array with a constant value:

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

Use the std::vector constructor for this:

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

If you further add a default value to the last parameter (const ElementT initVal = {}) then you don't need the previous constructor (Image(const size_t width, const size_t height)).

        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
        }

This constructor can be easily implemented in terms of std::copy. You do need to add a test here to verify that the input array has the right number of elements considering newWidth and newHeight.

The implementation of at doesn't check bounds. I would suggest adding assert statements for bounds, to catch bugs in debug mode. A release build would then not test bounds.

operator+= should really test that the sizes of the second image match those of the first. You can then simply iterate over the two arrays without worrying about x and y coordinates. If sizes don't match, the operation should fail.

It is not necessary to use this-> to access members. Some people like it, but I think it is superfluous.

The standard way to increment in C++ is ++x, not x++. For an int it doesn't make any difference, but for more complex objects it could be very wasteful. ++x increments the variable, then returns its value. x++ makes a copy of the value of the variable, increments the variable, then returns the copy. If you are not using the value of the expression, use the pre-increment to avoid the unnecessary copy.

image_operations.h

You don't need to forward declare class Image. You include its header file, so the class is already declared.

auto ratiox = (float)image.getWidth() / (float)width;

Please get used to the C++ way of casting:

auto ratiox = static_cast<float>(image.getWidth()) / static_cast<float>(width);

It is a more obvious cast and makes your code more readable.

Your functions gaussianFigure2D and gaussianFigure2D2 don't need different names. I would argue the library would be easier to use if they had the same name. The compiler will pick the one or the other function depending on the arguments passed.

std::exp is implemented for all three floating-point types. Consequently, normalDistribution1D can be implemented as a template:

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

The same is true for normalDistribution2D.

Functions here are declared out of order: functions use other functions that are declared later in the file. Somehow this works for you, I think, because templates are instantiated when used. It would be better practice to declare normalDistribution1D before gaussianFigure2D, and gaussianFigure2D2 before gaussianFigure2D. Note that class methods can be declared in any order -- the compiler will first collect all member names (collect the full definition of the class) before attempting to compile each of the function bodies. The same is not true for free functions such as the ones we're talking about here.

basic_functions.h

Again you have way too many headers included here.

This is very complex code, which seems out of place with the relatively straight-forward C++ code elsewhere. With the changes suggested above, you don't need this file at all.

base_types.h

You should not be using <cstdbool>. It is deprecated in C++17 and removed in C++20. The other <c...> headers should not be needed either. Use the C++ library, not the C library.

Instead of #define MAX_PATH 256, use constexpr int MAX_PATH = 256. Preprocessor macros are nice to generate code, but they should not be used for things that you can do with native C++ constructs.

Likewise, instead of the C construct typedef unsigned char BYTE, use using BYTE = unsigned char, which is, to me, much more readable.

This is also a C thing:

typedef struct RGB
{
    unsigned char channels[3];
} RGB;

In C++ it is just

struct RGB
{
    unsigned char channels[3];
};

Answer 2

"error: no matching function for call to 'gaussianFigure2D(const size_t&, const size_t&, const size_t&, const size_t&, const double&, const double&)'"

Here is how your code defines it:

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)

and the other form just has one more parameter.

Here's where it is called: auto image1 = TinyDIP::gaussianFigure2D2(10, 10, 5, 5, 3.0, 1.0);

That looks like it should work just fine. Since there are a different number of parameters, it will pick the one with 6 parameters and not get confused as to which you meant. The template argument deduction and conversions would not be any different when this is overloaded.

Having the second form simply supply the 6th argument automatically duplicating the 5th is a perfect use case for overloading!

I think your problem is that you defined them in the wrong order. The 6-argument form goes first, so the 5-argument form can call it. It worked as-written, even though gaussianFigure2D2 is not declared yet, because this is a template body and the call involves something of type T, so it puts off resolving the symbol until phase 2, when the template is instantiated.

---

Besides the structs that you apparently copied out of a C header file, as mentioned by Cris, you also have

typedef BYTE GrayScale;

Don't use typedef anymore, at all. For type alises, use using. And we have std::byte which is not spelled in all caps, so don't declare your own.

using GrayScale = std::byte;