Two dimensional bicubic interpolation implementation in C++

Question

This is a follow-up question for Two dimensional bicubic interpolation implementation in C and A recursive_transform Template Function with Unwrap Level for Various Type Arbitrary Nested Iterable Implementation in C++. Besides the C version code, I am attempting to make a C++ version bicubic interpolation function bicubicInterpolation which can be applied on two dimensional nested vectors std::vector<std::vector<>> structure.

Example input matrix:

1 1 1
1 100 1
1 1 1

Output matrix (bicubic interpolation result from the input matrix above):

1 1 1 1 1 1 1 1 1 1 1 1
1 9 19 27 30 27 19 9 1 0 0 0
1 19 39 56 62 56 39 19 1 0 0 0
1 27 56 79 89 79 56 27 1 0 0 0
1 30 62 89 100 89 62 30 1 0 0 0
1 27 56 79 89 79 56 27 1 0 0 0
1 19 39 56 62 56 39 19 1 0 0 0
1 9 19 27 30 27 19 9 1 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1
1 0 0 0 0 0 0 0 1 2 2 1
1 0 0 0 0 0 0 0 1 2 2 1
1 0 0 0 0 0 0 0 1 1 1 1

The experimental implementation

• namespace: TinyDIP

• bicubicInterpolation function implementation:

  constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY)
  {
   auto output = Image<ElementT>(newSizeX, newSizeY);
   auto ratiox = (float)this->getSizeX() / (float)newSizeX;
   auto ratioy = (float)this->getSizeY() / (float)newSizeY;
   for (size_t y = 0; y < newSizeY; y++)
   {
   for (size_t x = 0; x < newSizeX; 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)] = this->get(
   clip(xMappingToOriginFloor + ndatax, 0, this->getSizeX() - 1), 
   clip(yMappingToOriginFloor + ndatay, 0, this->getSizeY() - 1));
   }
   }
   output.set(x, y, bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac));
   }
   }
   return output;
  }
  

• Helper functions for bicubicInterpolation function:

  template<class InputT>
  constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
  }
  template<class InputT1, class InputT2>
  constexpr 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 InputT1, class InputT2, class InputT3>
  constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
  {
   if (input < lowerbound)
   {
   return static_cast<InputT1>(lowerbound);
   }
   if (input > upperbound)
   {
   return static_cast<InputT1>(upperbound);
   }
   return input;
  }
  

• 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"
  namespace TinyDIP
  {
   template <typename ElementT>
   class Image
   {
   public:
   Image()
   {
   }
   Image(const int newWidth, const int newHeight)
   {
   this->image_data.resize(newHeight);
   for (size_t i = 0; i < newHeight; ++i) {
   this->image_data[i].resize(newWidth);
   }
   this->image_data = recursive_transform<2>(this->image_data, [](ElementT element) { return ElementT{}; });
   return;
   }
   Image(const int newWidth, const int newHeight, const ElementT initVal)
   {
   this->image_data.resize(newHeight);
   for (size_t i = 0; i < newHeight; ++i) {
   this->image_data[i].resize(newWidth);
   }
   this->image_data = recursive_transform<2>(this->image_data, [initVal](ElementT element) { return initVal; });
   return;
   }
   Image(const std::vector<std::vector<ElementT>>& input)
   {
   this->image_data = recursive_transform<2>(input, [](ElementT element) {return element; } ); // Deep copy
   return;
   }
   template<class OutputT>
   constexpr auto cast()
   {
   return this->transform([](ElementT element) { return static_cast<OutputT>(element); });
   }
   constexpr auto get(const unsigned int locationx, const unsigned int locationy)
   {
   return this->image_data[locationy][locationx];
   }
   constexpr auto set(const unsigned int locationx, const unsigned int locationy, const ElementT& element)
   {
   this->image_data[locationy][locationx] = element;
   return *this;
   }
   template<class InputT>
   constexpr auto set(const unsigned int locationx, const unsigned int locationy, const InputT& element)
   {
   this->image_data[locationy][locationx] = static_cast<ElementT>(element);
   return *this;
   }
   constexpr auto getSizeX()
   {
   return this->image_data[0].size();
   }
   constexpr auto getSizeY()
   {
   return this->image_data.size();
   }
   constexpr auto getData()
   {
   return this->transform([](ElementT element) { return element; }); // Deep copy
   }
   void print()
   {
   for (auto& row_element : this->toString())
   {
   for (auto& element : row_element)
   {
   std::cout << element << "\t";
   }
   std::cout << "\n";
   }
   std::cout << "\n";
   return;
   }
   constexpr auto toString()
   {
   return this->transform([](ElementT element) { return std::to_string(element); });
   }
   constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY)
   {
   auto output = Image<ElementT>(newSizeX, newSizeY);
   auto ratiox = (float)this->getSizeX() / (float)newSizeX;
   auto ratioy = (float)this->getSizeY() / (float)newSizeY;
   for (size_t y = 0; y < newSizeY; y++)
   {
   for (size_t x = 0; x < newSizeX; 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)] = this->get(
   clip(xMappingToOriginFloor + ndatax, 0, this->getSizeX() - 1), 
   clip(yMappingToOriginFloor + ndatay, 0, this->getSizeY() - 1));
   }
   }
   output.set(x, y, bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac));
   }
   }
   return output;
   }
   Image<ElementT>& operator=(Image<ElementT> const& input) // Copy Assign
   {
   this->image_data = input.getData();
   return *this;
   }
   Image<ElementT>& operator=(Image<ElementT>&& other) // Move Assign
   {
   this->image_data = std::move(other.image_data);
   std::cout << "move assigned\n";
   return *this;
   }
   Image(const Image<ElementT> &input) // Copy Constructor
   {
   this->image_data = input.getData();
   }
   /* Move Constructor
   */
   Image(Image<ElementT> &&input) : image_data(std::move(input.image_data))
   {
   }
   private:
   std::vector<std::vector<ElementT>> image_data;
   template<class F>
   constexpr auto transform(const F& f)
   {
   return recursive_transform<2>(this->image_data, f);
   }
   template<class InputT>
   constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
   }
   template<class InputT1, class InputT2>
   constexpr 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 InputT1, class InputT2, class InputT3>
   constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
   {
   if (input < lowerbound)
   {
   return static_cast<InputT1>(lowerbound);
   }
   if (input > upperbound)
   {
   return static_cast<InputT1>(upperbound);
   }
   return input;
   }
   };
  }
  #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 "./"
  #define True true
  #define False false
  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
  

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

The full testing code

The grayscale type data has been tested here.

/* Develop by Jimmy Hu */
#include "base_types.h"
#include "basic_functions.h"
#include "image.h"
void bicubicInterpolationTest();
int main()
{
 bicubicInterpolationTest();
 return 0;
}
void bicubicInterpolationTest()
{
 TinyDIP::Image<GrayScale> image1(3, 3, 1);
 std::cout << "Width: " + std::to_string(image1.getSizeX()) + "\n";
 std::cout << "Height: " + std::to_string(image1.getSizeY()) + "\n";
 image1 = image1.set(1, 1, 100);
 image1.print();
 auto image2 = image1.bicubicInterpolation(12, 12);
 image2.print();
}

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Two dimensional bicubic interpolation implementation in C and

  A recursive_transform Template Function with Unwrap Level for Various Type Arbitrary Nested Iterable Implementation in C++.

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

  I am attempting to make a C++ version bicubic interpolation function bicubicInterpolation which can be applied on two dimensional nested vectors std::vector<std::vector<>> structure.

• Why a new review is being asked for?

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

• user555045: Can't help but notice that the bug from last time (filling 9 out of 16 entries in ndata) is still in here
• JimmyHu: @harold I apologize for the missing point you mentioned. The question has been updated.

Answer 1

I like auto for lots of things, but not as a return type unless it's actually necessary.

For example, to understand what cast returns it's necessary to dig through several template functions, layers of meta-templates, and then understand what std::invoke_result_t does in this particular case.

Please just write out the return type.

---

 template<class InputT1, class InputT2, class InputT3>
 constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
 {
 if (input < lowerbound)
 {
 return static_cast<InputT1>(lowerbound);
 }
 if (input > upperbound)
 {
 return static_cast<InputT1>(upperbound);
 }
 return input;
 }

The standard library already provides std::clamp so this isn't necessary.

Note that providing a single template argument is safer than separate arguments for all the inputs. Conversions (and thus comparisons) between different types may not be well-defined since types have different ranges.

---

#define True true
#define False false

Why? C++ uses true and false.

---

 std::vector<std::vector<ElementT>> image_data;

It's faster (and usually easier) to use a one-dimensional vector: std::vector<ElementT> image_data; containing width * height elements, and calculate indices when necessary as y * width + x.

So this:

 Image(const int newWidth, const int newHeight)
 {
 this->image_data.resize(newHeight);
 for (size_t i = 0; i < newHeight; ++i) {
 this->image_data[i].resize(newWidth);
 }
 this->image_data = recursive_transform<2>(this->image_data, [](ElementT element) { return ElementT{}; });
 return;
 }

Should simply be:

 Image(const unsigned int width, const unsigned int height):
 m_width(width),
 m_height(height),
 m_image_data(width * height) { }

Note that the std::vector constructor (and resize function) will "value-initialize" elements, so we don't need to do that separately.

We should also be using an unsigned type for the size, or at least checking that the integer isn't less than or equal to zero.

---

 template<class OutputT>
 constexpr auto cast()
 {
 return this->transform([](ElementT element) { return static_cast<OutputT>(element); });
 }

I don't think this belongs inside the class. If we want to do something with the internal data, it would be better to expose the internal data in a function like:

 std::vector<ElementT> const& getImageData() const { return m_imageData; }

The user can then use std::transform or whatever as they need to.

---

 constexpr auto get(const unsigned int locationx, const unsigned int locationy)
 {
 return this->image_data[locationy][locationx];
 }
 constexpr auto set(const unsigned int locationx, const unsigned int locationy, const ElementT& element)
 {
 this->image_data[locationy][locationx] = element;
 return *this;
 }
 template<class InputT>
 constexpr auto set(const unsigned int locationx, const unsigned int locationy, const InputT& element)
 {
 this->image_data[locationy][locationx] = static_cast<ElementT>(element);
 return *this;
 }

It's more idiomatic in C++ to provide two versions of a get function, one returning a mutable reference, and the other a const-reference.

constexpr ElementT& at(const unsigned int x, const unsigned int y) { return m_image_data[y * width + x]; }
constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return m_image_data[y * width + x]; }

---

 constexpr auto getData()
 {
 return this->transform([](ElementT element) { return element; }); // Deep copy
 }

As with cast above we don't need this function. We can just expose a const& of the image data.

---

 void print()
 {
 for (auto& row_element : this->toString())
 {
 for (auto& element : row_element)
 {
 std::cout << element << "\t";
 }
 std::cout << "\n";
 }
 std::cout << "\n";
 return;
 }

It's more flexible to provide an output stream operator, since we can specify the stream to print to.

We definitely don't want to be creating a temporary copy filled with strings to do this!!!

We can simply use ElementT's operator<< and send the elements directly to the output stream.

---

 constexpr auto toString()
 {
 return this->transform([](ElementT element) { return std::to_string(element); });
 }

This seems an unlikely thing for anyone to want to do. Again, we can expose the image data by const& and let the user use std::transform if they really want to.

---

 constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY) ...

There's no point in using references for the size parameters. Again, they should be unsigned, or we have to check that the values are above 0.

This function would probably be better as a free function, e.g.:

 template<ElementT>
 Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, int width, int height);

---

 Image<ElementT>& operator=(Image<ElementT> const& input) // Copy Assign
 {
 this->image_data = input.getData();
 return *this;
 }
 Image<ElementT>& operator=(Image<ElementT>&& other) // Move Assign
 {
 this->image_data = std::move(other.image_data);
 std::cout << "move assigned\n";
 return *this;
 }
 Image(const Image<ElementT> &input) // Copy Constructor
 {
 this->image_data = input.getData();
 }
 /* Move Constructor
 */
 Image(Image<ElementT> &&input) : image_data(std::move(input.image_data))
 {
 }

These can all be = default no?

Note that we can write Image<ElementT> as just Image inside the class.

---

 template<class InputT>
 constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
 }
 template<class InputT1, class InputT2>
 constexpr 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));
 }

These functions don't need to access any class data. They be static, or even non-member helper functions in another file.

As small POD types, the parameters should be passed by value, not const&.

• JimmyHu: Thank you for the answer. We can simply use ElementT's operator<< and send the elements directly to the output stream. How about TinyDIP::Image<GrayScale> image1 here? The GrayScale type is defined from unsigned char (8-bits per pixel), but operator<< in unsigned char can't print number correctly. The std::to_string is still needed here.
• user673679: True. I'd forgotten about that. It is possible to use the + trick to force numeric output instead, see: isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number