自研图像双三次插值与OpenCV、Scikit-image实现结果不符求助

我正尝试基于论文《Cubic Convolution Interpolation for Digital Image Processing》在Python中实现图像双三次卷积插值。尽管实现能正常缩放，但结果与OpenCV和Scikit-image的参考存在差异，小图像上尤为明显：

以下是最小可复现示例生成的对比图，包含原始图、我的插值结果、OpenCV/Scikit-image结果及差异：

最小可复现代码：

import math
import time
from functools import cache

import cv2 as cv
import matplotlib.pyplot as plt
import numpy as np
import skimage


def u(s: float):
    # bicubic convolution kernel aka catmull-rom spline
    # the value of a here is -0.5 as that was used in Keys' version
    a: float = -0.5
    s = abs(s)
    if 0 <= s < 1:
        return (a + 2) * s**3 - (a + 3) * s**2 + 1
    elif 1 <= s < 2:
        return a * s**3 - 5 * a * s**2 + 8 * a * s - 4 * a
    return 0


in_file = "test_sharpen.png"
ratio = 2.0

im_data = cv.imread(str(in_file))

# because plt uses rgb
im_data = cv.cvtColor(im_data, cv.COLOR_RGB2BGR)

start = time.perf_counter()
print("Scaling image...")

H, W, C = im_data.shape

# pad by 2 px
image = cv.copyMakeBorder(im_data, 2, 2, 2, 2, cv.BORDER_REFLECT)

image = image.astype(np.float64) / 255

# create new image
new_H = math.floor(H * ratio)
new_W = math.floor(W * ratio)
big_image = np.zeros((new_H, new_W, C))
for c in range(C):
    for j in range(new_H):
        # scale new image's coordinate to be in old image
        y = j * (1 / ratio) + 2
        # we separate x and y to integer and fractional parts
        iy = int(y)
        # ix and iy are essentially the closest original pixels
        # as all the old pixels are in integer positions
        # decx and decy as the fractional parts are then the distances
        # to the original pixels on the left and above
        decy = iy - y
        for i in range(new_W):
            x = i * (1 / ratio) + 2
            ix = int(x)
            decx = ix - x

            pix = sum(
                sum(
                    image[iy + M, ix + L, c] * u(decx + L) * u(decy + M)
                    for L in range(-1, 2 + 1)
                )
                for M in range(-1, 2 + 1)
            )

            # we limit results to [0, 1] because bicubic interpolation
            # can produce pixel values outside the original range
            big_image[j, i, c] = max(min(1, pix), 0)

big_image = (big_image * 255).astype(np.uint8)

print(f"Finished scaling in {time.perf_counter() - start} seconds")


# generate proper bicubic scales with opencv and skimage
# and compare them to my scale with plt
proper_cv = cv.resize(im_data, None, None, ratio, ratio, cv.INTER_CUBIC)
proper_skimage = skimage.util.img_as_ubyte(
    skimage.transform.rescale(im_data, ratio, channel_axis=-1, order=3)
)


fig, ax = plt.subplots(nrows=4, ncols=2)
ax[0, 0].imshow(im_data)
ax[0, 0].set_title("Original")
ax[0, 1].imshow(big_image)
ax[0, 1].set_title("My scale")

ax[1, 0].set_title("Proper OpenCV")
ax[1, 0].imshow(proper_cv)
ax[1, 1].set_title("Proper Skimage")
ax[1, 1].imshow(proper_cv)

print("my scale vs proper_cv psnr:", cv.PSNR(big_image, proper_cv))

ax[2, 0].set_title("Absdiff OpenCV vs My")
diffy_cv = cv.absdiff(big_image, proper_cv)
ax[2, 0].imshow(diffy_cv)
ax[2, 1].set_title("Absdiff Skimage vs My")
diffy_skimage = cv.absdiff(big_image, proper_skimage)
ax[2, 1].imshow(diffy_skimage)

ax[3, 1].set_title("Absdiff CV vs Skimage")
ax[3, 1].imshow(cv.absdiff(proper_cv, proper_skimage))
ax[3, 0].set_title("Absdiff CV vs Skimage")
ax[3, 0].imshow(cv.absdiff(proper_cv, proper_skimage))

print("diffy_cv", diffy_cv.min(), diffy_cv.max(), diffy_cv.dtype, diffy_cv.shape)
print(
    "diffy_skimage",
    diffy_skimage.min(),
    diffy_skimage.max(),
    diffy_skimage.dtype,
    diffy_skimage.shape,
)
print(
    "proper_skimage vs proper_opencv psnr:",
    cv.PSNR(big_image, proper_cv),
    cv.absdiff(proper_cv, proper_skimage).max(),
)
plt.show()

可通过python scaling.py将test_sharpen.png缩放至2倍。调整参数a无法解决差异，希望排查原因。

核心问题排查与修正 #

1. 像素中心坐标映射错误 #

这是导致差异的最关键原因。你的代码中，新图像坐标到原图像的映射没有考虑像素中心对齐：

• 原图像每个像素的中心位置是k + 0.5（k为0-based像素索引），新图像第j个像素的中心对应原图像的连续坐标应为(j + 0.5) / ratio（而非j / ratio）。
• 加上padding的偏移后，正确的坐标计算应为：

  y = (j + 0.5) / ratio + 2
  x = (i + 0.5) / ratio + 2
  

2. 整数部分取整逻辑问题 #

使用int(y)取整会在y为负数时出现错误（尽管padding后不会触发，但逻辑上不准确），应改用math.floor(y)来获取双三次窗口的左上角起始索引：

iy = math.floor(y)
ix = math.floor(x)

3. 目标尺寸计算不一致 #

OpenCV的cv.resize使用四舍五入计算目标尺寸，而你用math.floor会导致尺寸偏差，应改为：

new_H = round(H * ratio)
new_W = round(W * ratio)

修正后的核心循环代码 #

for c in range(C):
    for j in range(new_H):
        # 修正坐标映射：对齐像素中心
        y = (j + 0.5) / ratio + 2
        iy = math.floor(y)
        decy = y - iy  # 改为正的小数部分，不影响u函数结果，但逻辑更清晰
        for i in range(new_W):
            x = (i + 0.5) / ratio + 2
            ix = math.floor(x)
            decx = x - ix

            pix = sum(
                sum(
                    image[iy + M, ix + L, c] * u(decx - L) * u(decy - M)
                    for L in range(0, 4)  # 从0到3，对应ix到ix+3
                )
                for M in range(0, 4)  # 从0到3，对应iy到iy+3
            )

            big_image[j, i, c] = max(min(1, pix), 0)

注：将L/M从-1到2改为0到3，配合math.floor(y)，窗口范围更直观，计算结果一致。

其他细节 #

• OpenCV的INTER_CUBIC默认使用BORDER_REPLICATE处理边界，而你用BORDER_REFLECT，若需完全对齐，可将padding方式改为cv.BORDER_REPLICATE。
• Scikit-image的rescale默认会对图像进行归一化，你的代码中已经做了相同处理，修正坐标后结果会基本对齐。

内容的提问来源于stack exchange，提问作者tzv