Lanczos Resampling error

Question

I have written an image resizer using Lanczos re-sampling. I've taken the implementation straight from the directions on wikipedia. The results look good visually, but for some reason it does not match the result from Matlab's resize with Lanczos very well (in pixel error).

Does anybody see any errors? This is not my area of expertise at all...

Here is my filter (I'm using Lanczos3 by default):

double lanczos_size_ = 3.0;
inline double sinc(double x) {
  double pi = 3.1415926;
  x = (x * pi);
  if (x < 0.01 && x > -0.01)
    return 1.0 + x*x*(-1.0/6.0 + x*x*1.0/120.0);
  return sin(x)/x;
}

inline double LanczosFilter(double x) {
  if (std::abs(x) < lanczos_size_) {
    double pi = 3.1415926;
    return sinc(x)*sinc(x/lanczos_size_);
  } else {
    return 0.0;
  }
}

And my code to resize the image:

Image Resize(Image& image, int new_rows, int new_cols) {
  int old_cols = image.size().cols;
  int old_rows = image.size().rows;
  double col_ratio =
      static_cast<double>(old_cols)/static_cast<double>(new_cols);
  double row_ratio =
      static_cast<double>(old_rows)/static_cast<double>(new_rows);

  // Apply filter first in width, then in height.
  Image horiz_image(new_cols, old_rows);
  for (int r = 0; r < old_rows; r++) {
    for (int c = 0; c < new_cols; c++) {
      // x is the new col in terms of the old col coordinates.
      double x = static_cast<double>(c)*col_ratio;
      // The old col corresponding to the closest new col.
      int floor_x = static_cast<int>(x);

      horiz_image[r][c] = 0.0;
      double weight = 0.0;
      // Add up terms across the filter.
      for (int i = floor_x - lanczos_size_ + 1; i < floor_x + lanczos_size_; i++) {
        if (i >= 0 && i < old_cols) {
          double lanc_term = LanczosFilter(x - i);
          horiz_image[r][c] += image[r][i]*lanc_term;
          weight += lanc_term;
        }
      }
      // Normalize the filter.
      horiz_image[r][c] /= weight;
      // Strap the pixel values to valid values.
      horiz_image[r][c] = (horiz_image[r][c] > 1.0) ? 1.0 : horiz_image[r][c];
      horiz_image[r][c] = (horiz_image[r][c] < 0.0) ? 0.0 : horiz_image[r][c];
    }
  }

  // Now apply a vertical filter to the horiz image.
  Image new_image(new_cols, new_rows);
  for (int r = 0; r < new_rows; r++) {
    double x = static_cast<double>(r)*row_ratio;
    int floor_x = static_cast<int>(x);
    for (int c = 0; c < new_cols; c++) {      
      new_image[r][c] = 0.0;
      double weight = 0.0;
      for (int i = floor_x - lanczos_size_ + 1; i < floor_x + lanczos_size_; i++) {
        if (i >= 0 && i < old_rows) {
          double lanc_term = LanczosFilter(x - i);
          new_image[r][c] += horiz_image[i][c]*lanc_term;
          weight += lanc_term;
        }
      }
      new_image[r][c] /= weight;
      new_image[r][c] = (new_image[r][c] > 1.0) ? 1.0 : new_image[r][c];
      new_image[r][c] = (new_image[r][c] < 0.0) ? 0.0 : new_image[r][c];
    }
  }
  return new_image;
}

Answer 1

I think there is a mistake in your sinc function. Below the fraction bar you have to square pi and x. Additional you have to multiply the function with lanczos size L(x) = **a***sin(pi*x)*sin(pi*x/a) * (pi**²**x**²**)^-1

Edit: My mistake, there is all right.

Answer 2

Here is Lanczosh in one single loop. no errors. Uses mentioned at top procedures.

void ResizeDD(
double* const pixelsSrc, 
const int old_cols, 
const int old_rows,
double* const pixelsTarget, 
int const new_rows, int const new_cols)
{
double col_ratio =
    static_cast<double>(old_cols) / static_cast<double>(new_cols);
double row_ratio =
    static_cast<double>(old_rows) / static_cast<double>(new_rows);

// Now apply a filter to the image.
for (int r = 0; r < new_rows; ++r)
{
    const double row_within = static_cast<double>(r)* row_ratio;
    int floor_row = static_cast<int>(row_within);
    for (int c = 0; c < new_cols; ++c)
    {
        // x is the new col in terms of the old col coordinates.
        double col_within = static_cast<double>(c)* col_ratio;
        // The old col corresponding to the closest new col.
        int floor_col = static_cast<int>(col_within);

        double& v_toSet = pixelsTarget[r * new_cols + c];
        v_toSet = 0.0;
        double weight = 0.0;
        for (int i = floor_row - lanczos_size_ + 1; i <= floor_row + lanczos_size_; ++i)
        {
            for (int j = floor_col - lanczos_size_ + 1; j <= floor_col + lanczos_size_; ++j)
            {
                if (i >= 0 && i < old_rows && j >= 0 && j < old_cols)
                {
                    const double lanc_term = LanczosFilter(row_within - i + col_within - j);
                    v_toSet += pixelsSrc[i * old_rows + j] * lanc_term;
                    weight += lanc_term;
                }
            }
        }

        v_toSet /= weight;
        v_toSet = (v_toSet > 1.0) ? 1.0 : v_toSet;
        v_toSet = (v_toSet < 0.0) ? 0.0 : v_toSet;
    }
}

}

Answer 3

The line

for (int i = floor_x - lanczos_size_ + 1; i < floor_x + lanczos_size_; i++)

should be

for (int i = floor_x - lanczos_size_ + 1; i <= floor_x + lanczos_size_; i++)

Do not know but perhaps other mistakes linger too.