lanczos_resample Template Function Implementation for Image in C++

Question

This is a follow-up question for Two dimensional bicubic interpolation implementation in C++. I am trying to implement Lanczos resampling for 2D image in this post.

The experimental implementation

• lanczos_resample template function implementation

  namespace TinyDIP
  {
   // Lanczos resampling function for 2D images
   // input_image: The original image to be resampled
   // new_width: The target width of the resampled image
   // new_height: The target height of the resampled image
   // a: The Lanczos kernel size parameter (default is 3 for high quality)
   template<typename ElementT>
   constexpr Image<ElementT> lanczos_resample(
   const Image<ElementT>& input_image,
   const std::size_t new_width,
   const std::size_t new_height,
   const int a = 3
   )
   {
   if (input_image.getDimensionality() != 2)
   {
   throw std::runtime_error("Lanczos resampling is implemented for 2D images only.");
   }
   const std::size_t old_width = input_image.getWidth();
   const std::size_t old_height = input_image.getHeight();
   Image<ElementT> output_image(new_width, new_height);
   const double x_ratio = static_cast<double>(old_width) / new_width;
   const double y_ratio = static_cast<double>(old_height) / new_height;
   // Process each pixel in the new 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)
   {
   // Map the new pixel coordinates back to the old image coordinates
   double x_old = (x_new + 0.5) * x_ratio - 0.5;
   double y_old = (y_new + 0.5) * y_ratio - 0.5;
   ElementT pixel_value{}; // Initialize with zero
   double total_weight = 0.0;
   // Determine the range of pixels in the old image that contribute to the new pixel
   int x_min = static_cast<int>(std::floor(x_old)) - a + 1;
   int x_max = static_cast<int>(std::floor(x_old)) + a;
   int y_min = static_cast<int>(std::floor(y_old)) - a + 1;
   int y_max = static_cast<int>(std::floor(y_old)) + a;
   // Iterate over the contributing pixels
   for (int j = y_min; j <= y_max; ++j)
   {
   for (int i = x_min; i <= x_max; ++i)
   {
   // Clamp coordinates to be within the bounds of the old image
   int clamped_i = std::clamp(i, 0, static_cast<int>(old_width) - 1);
   int clamped_j = std::clamp(j, 0, static_cast<int>(old_height) - 1);
   // Calculate the weight using the Lanczos kernel
   double weight_x = lanczos_kernel(x_old - i, a);
   double weight_y = lanczos_kernel(y_old - j, a);
   double weight = weight_x * weight_y;
   if (weight != 0.0)
   {
   pixel_value += input_image.at(clamped_i, clamped_j) * weight;
   total_weight += weight;
   }
   }
   }
   // Normalize the pixel value
   if (total_weight != 0.0)
   {
   output_image.at(x_new, y_new) = pixel_value / total_weight;
   }
   }
   }
   return output_image;
   }
   // lanczos_resample template function implementation
   template<typename ElementT>
   requires((std::same_as<ElementT, RGB>) ||
   (std::same_as<ElementT, RGB_DOUBLE>) ||
   (std::same_as<ElementT, HSV>) ||
   is_MultiChannel<ElementT>::value)
   constexpr static auto lanczos_resample(
   const Image<ElementT>& input_image,
   const std::size_t new_width,
   const std::size_t new_height,
   const int a = 3)
   {
   return apply_each(input_image, [&](auto&& each_plane)
   {
   return lanczos_resample(
   each_plane,
   new_width,
   new_height,
   a);
   });
   }
  }
  

• lanczos_kernel template function implementation

  namespace TinyDIP
  {
   // Lanczos kernel implementation
   // a: The kernel size parameter (e.g., 2 or 3)
   // x: The input value
   template<typename T>
   constexpr T lanczos_kernel(T x, const int a)
   {
   if (x == 0.0) {
   return T(1.0);
   }
   if (std::abs(x) < a) {
   return (a * std::sin(std::numbers::pi_v<T> *x) * std::sin(std::numbers::pi_v<T> *x / a)) / (std::numbers::pi_v<T> *std::numbers::pi_v<T> *x * x);
   }
   return T(0.0);
   }
  }
  

The usage of lanczos_resample template function:

/* Developed by Jimmy Hu */
#include <iostream>
#include "../base_types.h"
#include "../basic_functions.h"
#include "../image.h"
#include "../image_io.h"
#include "../image_operations.h"
#include "../timer.h"
void lanczosResamplingTest(
 std::string_view input_image_path = "../InputImages/1",
 std::string_view output_image_path = "../OutputImages/lanczosResamplingTest")
{
 auto input_img = TinyDIP::bmp_read(std::string(input_image_path).c_str(), false);
 auto output_img =
 TinyDIP::lanczos_resample(
 input_img,
 input_img.getWidth() * 2,
 input_img.getHeight() * 2
 );
 TinyDIP::bmp_write(
 std::string(output_image_path).c_str(),
 output_img);
 
}
int main(int argc, char* argv[])
{
 TinyDIP::Timer timer;
 lanczosResamplingTest();
 return EXIT_SUCCESS;
}

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Two dimensional bicubic interpolation implementation in C++

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

  I implemented lanczos_resample template function for 2-dimensional image in this post.

• Why a new review is being asked for?

  Please review the implementation of lanczos_resample template function and its tests.

• Cris Luengo: The computations are more efficient if you apply 1D interpolation along the rows, then along columns. If the scaling is a simple fraction, you can re-use the kernel weights too, computing those is the most expensive part of this algorithm.