Image Rotation with Shear Transformation in C++

Question

This is a follow-up question for Image Rotation and Transpose Functions Implementation in C++ and An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. I am trying to implement another image rotation function with shear transformation in this post. As another approach, the function is named rotate_detail_shear_transformation and I keep rotate_detail with trigonometric function ways for learning purpose. The example output:

Image Input	Output
Input	Output

The experimental implementation

• The first version of rotate_detail_shear_transformation template functions implementation (in file image_operations.h)

  namespace TinyDIP
  {
   // rotate_detail_shear_transformation template function implementation
   // rotate_detail_shear_transformation template function performs image rotation
   // Reference: https://gautamnagrawal.medium.com/rotating-image-by-any-angle-shear-transformation-using-only-numpy-d28d16eb5076
   template<arithmetic ElementT, std::floating_point FloatingType = double>
   constexpr static auto rotate_detail_shear_transformation(const Image<ElementT>& input, FloatingType radians)
   {
   if (input.getDimensionality()!=2)
   {
   throw std::runtime_error("Unsupported dimension!");
   }
   radians = std::fmod(radians, 2 * std::numbers::pi_v<long double>);
   // if negative degrees
   if(radians < 0)
   {
   radians = radians + 2 * std::numbers::pi_v<long double>;
   }
   // if 0° rotation case
   if (radians == 0)
   {
   return input;
   }
   // if 90° rotation case
   if(radians == std::numbers::pi_v<long double> / 2.0)
   {
   Image<ElementT> output(input.getHeight(), input.getWidth());
   for (std::size_t y = 0; y < input.getHeight(); ++y)
   {
   for (std::size_t x = 0; x < input.getWidth(); ++x)
   {
   output.at(input.getHeight() - y - 1, x) = 
   input.at(x, y);
   }
   }
   return output;
   }
   auto cosine = std::cos(radians);
   auto sine = std::sin(radians);
   auto height = input.getHeight();
   auto width = input.getWidth();
   FloatingType original_centre_width = std::round((static_cast<FloatingType>(width) + 1.0) / 2.0 - 1.0);
   FloatingType original_centre_height = std::round((static_cast<FloatingType>(height) + 1.0) / 2.0 - 1.0);
   // Define the height and width of the new image that is to be formed
   auto new_height = std::round(std::abs(height*cosine) + std::abs(width*sine)) + 1;
   auto new_width = std::round(std::abs(width*cosine) + std::abs(height*sine)) + 1;
   // Define another image variable of dimensions of new_height and new _column filled with zeros
   Image<ElementT> output(static_cast<std::size_t>(new_width), static_cast<std::size_t>(new_height));
   // Find the centre of the new image that will be obtained
   FloatingType new_centre_width = std::round((static_cast<FloatingType>(new_width) + 1.0) / 2.0 - 1.0);
   FloatingType new_centre_height = std::round((static_cast<FloatingType>(new_height) + 1.0) / 2.0 - 1.0);
   for (std::size_t i = 0; i < input.getHeight(); ++i)
   {
   for (std::size_t j = 0; j < input.getWidth(); ++j)
   {
   // co-ordinates of pixel with respect to the centre of original image
   auto y = height - 1.0 - i - original_centre_height;
   auto x = width - 1.0 - j - original_centre_width;
   // co-ordinate of pixel with respect to the rotated image
   auto new_y = std::round(-x * sine + y * cosine);
   auto new_x = std::round(x * cosine + y * sine);
   /* since image will be rotated the centre will change too, 
   so to adust to that we will need to change new_x and new_y with respect to the new centre*/
   new_y = new_centre_height - new_y;
   new_x = new_centre_width - new_x;
   if((0 <= new_x) && (new_x < new_width) &&
   (0 <= new_y) && (new_y < new_height))
   {
   output.at(
   static_cast<std::size_t>(new_x),
   static_cast<std::size_t>(new_y)) = 
   input.at(j, i);
   }
   }
   }
   return output;
   }
   // rotate_detail_shear_transformation template function implementation
   template<typename ElementT, class FloatingType = double>
   requires ((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, HSV>)) // TODO: Create a base class for both RGB and HSV
   constexpr static auto rotate_detail_shear_transformation(const Image<ElementT>& input, FloatingType radians)
   {
   if (input.getDimensionality()!=2)
   {
   throw std::runtime_error("Unsupported dimension!");
   }
   return apply_each(input, [&](auto&& planes) { return rotate_detail_shear_transformation(planes, radians); });
   }
  }
  

• The first version of rotate_detail_shear_transformation template function implementation using the technique called forward mapping. It iterates through each pixel of the source image, calculates its new position in the destination image. There is an issue of aliasing in this way. To improve it, another better approach (second version) is to use reverse mapping combined with bilinear interpolation.

  namespace TinyDIP
  {
   // rotate_detail_shear_transformation template function implementation
   // rotate_detail_shear_transformation template function performs image rotation
   // Reference: https://gautamnagrawal.medium.com/rotating-image-by-any-angle-shear-transformation-using-only-numpy-d28d16eb5076
   template<arithmetic ElementT, arithmetic FloatingType = double>
   constexpr static auto rotate_detail_shear_transformation(const Image<ElementT>& input, FloatingType radians)
   {
   if (input.getDimensionality() != 2)
   {
   throw std::runtime_error("Unsupported dimension!");
   }
   radians = std::fmod(radians, 2 * std::numbers::pi_v<long double>);
   if (radians < 0)
   {
   radians = radians + 2 * std::numbers::pi_v<long double>;
   }
   // Handle 0-degree rotation
   if (radians == 0)
   {
   return input;
   }
   // Handle 90-degree rotation (can still be optimized for speed)
   if (radians == std::numbers::pi_v<long double> / 2.0)
   {
   Image<ElementT> output(input.getHeight(), input.getWidth());
   for (std::size_t y = 0; y < input.getHeight(); ++y)
   {
   for (std::size_t x = 0; x < input.getWidth(); ++x)
   {
   output.at(input.getHeight() - y - 1, x) = input.at(x, y);
   }
   }
   return output;
   }
   auto cosine = std::cos(radians);
   auto sine = std::sin(radians);
   auto height = static_cast<FloatingType>(input.getHeight());
   auto width = static_cast<FloatingType>(input.getWidth());
   // Calculate original image center
   FloatingType original_centre_width = (width - 1.0) / 2.0;
   FloatingType original_centre_height = (height - 1.0) / 2.0;
   // Calculate the dimensions of the new image
   auto new_height_f = std::abs(height * cosine) + std::abs(width * sine);
   auto new_width_f = std::abs(width * cosine) + std::abs(height * sine);
   auto new_height = static_cast<std::size_t>(std::ceil(new_height_f));
   auto new_width = static_cast<std::size_t>(std::ceil(new_width_f));
   Image<ElementT> output(new_width, new_height);
   output.setAllValue(ElementT{ 0 }); // Initialize with a background color
   // Calculate new image center
   FloatingType new_centre_width = (static_cast<FloatingType>(new_width) - 1.0) / 2.0;
   FloatingType new_centre_height = (static_cast<FloatingType>(new_height) - 1.0) / 2.0;
   // Reverse Mapping: Iterate over the destination image
   for (std::size_t y_new = 0; y_new < new_height; ++y_new)
   {
   for (std::size_t x_new = 0; x_new < new_width; ++x_new)
   {
   // Translate coordinates to be relative to the new image's center
   auto x_prime = static_cast<FloatingType>(x_new) - new_centre_width;
   auto y_prime = static_cast<FloatingType>(y_new) - new_centre_height;
   // Apply the inverse rotation to find the corresponding source coordinate
   // x = x'cos(θ) + y'sin(θ)
   // y = -x'sin(θ) + y'cos(θ)
   auto x_original_centered = x_prime * cosine + y_prime * sine;
   auto y_original_centered = -x_prime * sine + y_prime * cosine;
   // Translate back to the original image's coordinate system
   auto x_source = x_original_centered + original_centre_width;
   auto y_source = y_original_centered + original_centre_height;
   // Use bilinear interpolation to get the pixel value
   output.at(x_new, y_new) = bilinear_interpolate(input, x_source, y_source);
   }
   }
   return output;
   }
  }
  

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

  namespace TinyDIP
  {
   template<arithmetic ElementT, arithmetic FloatingType>
   constexpr ElementT bilinear_interpolate(const Image<ElementT>& image, const FloatingType x, const FloatingType y)
   {
   const auto width = static_cast<FloatingType>(image.getWidth());
   const auto height = static_cast<FloatingType>(image.getHeight());
   // Boundary check
   if (x < 0 || x > width - 1 || y < 0 || y > height - 1)
   {
   return ElementT{0}; // Return black/zero for out-of-bounds pixels
   }
   // Get the integer coordinates of the four surrounding pixels
   auto x1 = static_cast<std::size_t>(x);
   auto y1 = static_cast<std::size_t>(y);
   auto x2 = x1 + 1;
   auto y2 = y1 + 1;
   // Ensure the second set of coordinates are within bounds
   if (x2 >= image.getWidth()) x2 = x1;
   if (y2 >= image.getHeight()) y2 = y1;
   // Get the values of the four corner pixels
   const ElementT& q11 = image.at(x1, y1);
   const ElementT& q21 = image.at(x2, y1);
   const ElementT& q12 = image.at(x1, y2);
   const ElementT& q22 = image.at(x2, y2);
   // Calculate fractional parts (weights)
   FloatingType dx = x - x1;
   FloatingType dy = y - y1;
   // Interpolate in the x-direction
   auto r1 = q11 * (1.0 - dx) + q21 * dx;
   auto r2 = q12 * (1.0 - dx) + q22 * dx;
   // Interpolate in the y-direction and cast back to the element type
   return static_cast<ElementT>(r1 * (1.0 - dy) + r2 * dy);
   }
  }
  

• 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 rotate function:

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

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Image Rotation and Transpose Functions Implementation in C++ and

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

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

  I am trying to implement image rotation with shear transformation in this post.

• Why a new review is being asked for?

  Please review the implementation of rotate_detail_shear_transformation template functions.

• Loki Astari: Rotation without rotating. by Matt Parker YouTube
• Loki Astari: You have obvious color averaging issues. There are obvious artifacts in the rotated image. If you are using the obvious color averaging technique this is probably the cause. See here for a better technique. youtube.com/watch?v=LKnqECcg6Gw&t=2s
• Loki Astari: If memory servies me: bad => RGB = {(R1+R2)/2, (G1+G2)/2, (B1+B2)/2} Good => RGB = { SQRT(R1^2 + R2^2)/2, SQRT(G1^2 + G2^2)/2, SQRT(B1^2 + B2^2)/2}
• ALX23z: Normally, people fill in the new image by iterating over each element of the target image and figuring out its source point. Otherwise, you end up with various problems of skipping points and having multiple source points match the same target point. Of course, for efficiency purposes, you should add fast filters to sift out-bounds points.