Gaussian Fisheye Image Generator Implementation in C++

Question

This is a follow-up question for An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++ and Three dimensional gaussian image generator in C++. I am trying to make a Gaussian fisheye image generator in this post. The example output:

Image Input	Output

The experimental implementation

• gaussian_fisheye template function implementation (in file image_operations.h)

  namespace TinyDIP
  {
      //  gaussian_fisheye template function implementation
      template<arithmetic ElementT, std::floating_point FloatingType = double>
      constexpr static auto gaussian_fisheye(const Image<ElementT>& input, FloatingType D0)
      {
          if (input.getDimensionality()!=2)
          {
              throw std::runtime_error("Unsupported dimension!");
          }
  
          Image<ElementT> output(input.getWidth(), input.getHeight());
          for (std::size_t y = 0; y < input.getHeight(); ++y)
          {
              for (std::size_t x = 0; x < input.getWidth(); ++x)
              {
                  FloatingType distance_x = x - static_cast<FloatingType>(input.getWidth()) / 2.0;
                  FloatingType distance_y = y - static_cast<FloatingType>(input.getHeight()) / 2.0;
                  FloatingType distance = std::hypot(distance_x, distance_y);
                  FloatingType angle = std::atan2(distance_y, distance_x);
                  FloatingType weight = normalDistribution2D(std::fabs(distance_x), std::fabs(distance_y), D0) / normalDistribution2D(0.0, 0.0, D0);
                  FloatingType new_distance = distance * weight;
                  FloatingType new_distance_x = new_distance * std::cos(angle);
                  FloatingType new_distance_y = new_distance * std::sin(angle);
                  output.at(
                      static_cast<std::size_t>(new_distance_x + static_cast<FloatingType>(input.getWidth()) / 2.0),
                      static_cast<std::size_t>(new_distance_y + static_cast<FloatingType>(input.getHeight()) / 2.0)) = 
                      input.at(x, y);
              }
          }
          return output;
      }
  
      //  gaussian_fisheye template function implementation
      template<typename ElementT, class FloatingType = double>
      requires ((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, HSV>))
      constexpr static auto gaussian_fisheye(const Image<ElementT>& input, FloatingType D0)
      {
          return apply_each(input, [&](auto&& planes) { return gaussian_fisheye(planes, D0); });
      }
  }
  

• normalDistribution2D function implementation (in file image_operations.h)

  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 * std::numbers::pi * standard_deviation * standard_deviation);
  }
  

• apply_each function implementation (in file image_operations.h)

  template<class T = RGB, class F, class... Args>
  requires (std::same_as<T, RGB>)
  constexpr static auto apply_each(Image<T> input, F operation, Args&&... args)
  {
      return constructRGB(operation(getRplane(input), args...), operation(getGplane(input), args...), operation(getBplane(input), args...));
  }
  

• Image class implementation (in file image.h)

  namespace TinyDIP
  {
      template <typename ElementT>
      class Image
      {
      public:
          Image() = default;
  
          template<std::same_as<std::size_t>... Sizes>
          Image(Sizes... sizes): size{sizes...}, image_data((1 * ... * sizes))
          {}
  
          template<std::same_as<int>... Sizes>
          Image(Sizes... sizes)
          {
              size.reserve(sizeof...(sizes));
              (size.push_back(sizes), ...);
              image_data.resize(
                  std::reduce(
                      std::ranges::cbegin(size),
                      std::ranges::cend(size),
                      std::size_t{1},
                      std::multiplies<>()
                      )
              );
          }
  
          template<std::ranges::input_range Range,
                   std::same_as<std::size_t>... Sizes>
          Image(const Range& input, Sizes... sizes):
              size{sizes...}, image_data(begin(input), end(input))
          {
              if (image_data.size() != (1 * ... * sizes)) {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
          }
  
          Image(std::vector<ElementT>&& input, std::size_t newWidth, std::size_t newHeight)
          {
              size.reserve(2);
              size.emplace_back(newWidth);
              size.emplace_back(newHeight);
              if (input.size() != newWidth * newHeight)
              {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
              image_data = std::move(input);              //  Reference: https://stackoverflow.com/a/51706522/6667035
          }
  
          Image(const std::vector<std::vector<ElementT>>& input)
          {
              size.reserve(2);
              size.emplace_back(input[0].size());
              size.emplace_back(input.size());
              for (auto& rows : input)
              {
                  image_data.insert(image_data.end(), std::ranges::begin(input), std::ranges::end(input));    //  flatten
              }
              return;
          }
  
          //  at template function implementation
          template<typename... Args>
          constexpr ElementT& at(const Args... indexInput)
          {
              return const_cast<ElementT&>(static_cast<const Image &>(*this).at(indexInput...));
          }
  
          //  at template function implementation
          //  Reference: https://codereview.stackexchange.com/a/288736/231235
          template<typename... Args>
          constexpr ElementT const& at(const Args... indexInput) const
          {
              checkBoundary(indexInput...);
              constexpr std::size_t n = sizeof...(Args);
              if(n != size.size())
              {
                  throw std::runtime_error("Dimensionality mismatched!");
              }
              std::size_t i = 0;
              std::size_t stride = 1;
              std::size_t position = 0;
  
              auto update_position = [&](auto index) {
                  position += index * stride;
                  stride *= size[i++];
              };
              (update_position(indexInput), ...);
  
              return image_data[position];
          }
  
          constexpr std::size_t count() const noexcept
          {
              return std::reduce(std::ranges::cbegin(size), std::ranges::cend(size), 1, std::multiplies());
          }
  
          constexpr std::size_t getDimensionality() const noexcept
          {
              return size.size();
          }
  
          constexpr std::size_t getWidth() const noexcept
          {
              return size[0];
          }
  
          constexpr std::size_t getHeight() const noexcept
          {
              return size[1];
          }
  
          constexpr auto getSize() noexcept
          {
              return size;
          }
  
          std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
  
          void print(std::string separator = "\t", std::ostream& os = std::cout) const
          {
              if(size.size() == 1)
              {
                  for(std::size_t x = 0; x < size[0]; ++x)
                  {
                      //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                      os << +at(x) << separator;
                  }
                  os << "\n";
              }
              else if(size.size() == 2)
              {
                  for (std::size_t y = 0; y < size[1]; ++y)
                  {
                      for (std::size_t x = 0; x < size[0]; ++x)
                      {
                          //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                          os << +at(x, y) << separator;
                      }
                      os << "\n";
                  }
                  os << "\n";
              }
              else if (size.size() == 3)
              {
                  for(std::size_t z = 0; z < size[2]; ++z)
                  {
                      for (std::size_t y = 0; y < size[1]; ++y)
                      {
                          for (std::size_t x = 0; x < size[0]; ++x)
                          {
                              //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                              os << +at(x, y, z) << separator;
                          }
                          os << "\n";
                      }
                      os << "\n";
                  }
                  os << "\n";
              }
          }
  
          //  Enable this function if ElementT = RGB
          void print(std::string separator = "\t", std::ostream& os = std::cout) const
          requires(std::same_as<ElementT, RGB>)
          {
              for (std::size_t y = 0; y < size[1]; ++y)
              {
                  for (std::size_t x = 0; x < size[0]; ++x)
                  {
                      os << "( ";
                      for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                      {
                          //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                          os << +at(x, y).channels[channel_index] << separator;
                      }
                      os << ")" << separator;
                  }
                  os << "\n";
              }
              os << "\n";
              return;
          }
  
          Image<ElementT>& setAllValue(const ElementT input)
          {
              std::fill(std::ranges::begin(image_data), std::ranges::end(image_data), input);
              return *this;
          }
  
          friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
          {
              const std::string separator = "\t";
              rhs.print(separator, os);
              return os;
          }
  
          Image<ElementT>& operator+=(const Image<ElementT>& rhs)
          {
              check_size_same(rhs, *this);
              std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                     std::ranges::begin(image_data), std::plus<>{});
              return *this;
          }
  
          Image<ElementT>& operator-=(const Image<ElementT>& rhs)
          {
              check_size_same(rhs, *this);
              std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                     std::ranges::begin(image_data), std::minus<>{});
              return *this;
          }
  
          Image<ElementT>& operator*=(const Image<ElementT>& rhs)
          {
              check_size_same(rhs, *this);
              std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                     std::ranges::begin(image_data), std::multiplies<>{});
              return *this;
          }
  
          Image<ElementT>& operator/=(const Image<ElementT>& rhs)
          {
              check_size_same(rhs, *this);
              std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                     std::ranges::begin(image_data), std::divides<>{});
              return *this;
          }
  
          friend bool operator==(Image<ElementT> const&, Image<ElementT> const&) = default;
  
          friend bool operator!=(Image<ElementT> const&, Image<ElementT> const&) = default;
  
          friend Image<ElementT> operator+(Image<ElementT> input1, const Image<ElementT>& input2)
          {
              return input1 += input2;
          }
  
          friend Image<ElementT> operator-(Image<ElementT> input1, const Image<ElementT>& input2)
          {
              return input1 -= input2;
          }
  
          friend Image<ElementT> operator*(Image<ElementT> input1, ElementT input2)
          {
              return multiplies(input1, input2);
          }
  
          friend Image<ElementT> operator*(ElementT input1, Image<ElementT> input2)
          {
              return multiplies(input2, input1);
          }
  
  #ifdef USE_BOOST_SERIALIZATION
  
          void Save(std::string filename)
          {
              const std::string filename_with_extension = filename + ".dat";
              //    Reference: https://stackoverflow.com/questions/523872/how-do-you-serialize-an-object-in-c
              std::ofstream ofs(filename_with_extension, std::ios::binary);
              boost::archive::binary_oarchive ArchiveOut(ofs);
              //    write class instance to archive
              ArchiveOut << *this;
              //    archive and stream closed when destructors are called
              ofs.close();
          }
  
  #endif
      private:
          std::vector<std::size_t> size;
          std::vector<ElementT> image_data;
  
          template<typename... Args>
          void checkBoundary(const Args... indexInput) const
          {
              constexpr std::size_t n = sizeof...(Args);
              if(n != size.size())
              {
                  throw std::runtime_error("Dimensionality mismatched!");
              }
              std::size_t parameter_pack_index = 0;
              auto function = [&](auto index) {
                  if (index >= size[parameter_pack_index])
                      throw std::out_of_range("Given index out of range!");
                  parameter_pack_index = parameter_pack_index + 1;
              };
  
              (function(indexInput), ...);
          }
  
  #ifdef USE_BOOST_SERIALIZATION
          friend class boost::serialization::access;
          template<class Archive>
          void serialize(Archive& ar, const unsigned int version)
          {
              ar& size;
              ar& image_data;
          }
  #endif
  
      };
  
      template<typename T, typename ElementT>
      concept is_Image = std::is_same_v<T, Image<ElementT>>;
  }
  
  
  #endif
  

The usage of gaussian_fisheye function:

std::string file_path = "InputImages/1";
auto bmp1 = TinyDIP::bmp_read(file_path.c_str(), false);
bmp1 = gaussian_fisheye(bmp1, 800.0);
TinyDIP::bmp_write("test", bmp1);

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++ and

  Three dimensional gaussian image generator in C++

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

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

• Why a new review is being asked for?

  Please review the implementation of gaussian_fisheye template function.

Answer 1

Unnecessary calculations

You calculate angle using std::atan2(), but later use std::cos() and std::sin() to convert it back into x and y coordinates. This is unnecessary, you already had those to begin with, you only need to scale them.

Furthermore, you don't need to use std::fabs() before passing the x and y position to normalDistribution2D(), since it will square those values anyway.

Finally, you don't need any of the static casts to FloatingType, only the final casts back to std::size_t. So:

FloatingType distance_x = x - input.getWidth() / 2.0;
FloatingType distance_y = y - input.getHeight() / 2.0;
FloatingType weight = normalDistribution2D(distance_x, distance_y, D0)
                    / normalDistribution2D(0, 0, D0);
FloatingType new_distance_x = distance_x * weight;
FloatingType new_distance_y = distance_y * weight;
output.at(static_cast<std::size_t>(new_distance_x + input.getWidth() / 2.0),
          static_cast<std::size_t>(new_distance_y + input.getHeight() / 2.0)) =
    input.at(x, y);