Lanczos Resampling грешка

Написах програма за преоразмеряване на изображения, използвайки повторна семплиране на Lanczos. Взех изпълнението направо от указанията в wikipedia. Резултатите изглеждат добре визуално, но по някаква причина не съвпадат много добре с резултата от преоразмеряването на Matlab с Lanczos (в пикселна грешка).

Някой вижда ли грешки? Това изобщо не е моята област на компетентност...

Ето моят филтър (използвам Lanczos3 по подразбиране):

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

И моят код за преоразмеряване на изображението:

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