Image Bicubic Interpolation does not match OpenCV and Scikit-image implementations

Question

I am trying to implement bicubic convolution interpolation for images from the paper "Cubic Convolution Interpolation for Digital Image Processing" in Python. However, my implementation which looks like a proper scale still differs from reference implementations and I do not understand why. This is especially noticeable in smaller images, like this one:

Here's an image generated by the MWE with original unscaled image, my bad bicubic, opencv/skimage bicubic scales, and their differences from my scaled image.

Here's the code I have so far turned into a MWE without multiprocessing:

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()

It can be used with e.g. python scaling.py to scale test_sharpening.png to 2x.

I made the implementation so far and it seems to work okay, but still differs. I also tried changing the value of a , but that's not the problem.

Answer 1

It seems my method of scaling the coordinates was wrong. For example with ratio 2 the new points on the y-axis were 2.0, 2.5, 3.0 and so on.

This is wrong as the coordinates are supposed to be in-between the old points and not directly on top of them. I changed the scaling to be:

# scale new image's coordinate to be in old image based on its midpoint
y = ((j + 0.5) / ratio) - 0.5 + 2
x = ((i + 0.5) / ratio) - 0.5 + 2

and now the new point coordinates are 1.75, 2.25, 2.75 and so on. My intuition would tell me that the points should start from 2.25 instead, but in that case the image seems to be shifted slightly compared to reference.

This implementation now matches up near-perfectly with cv2's implementation when a=-0.75 along with other implementations, with the only exception being my reflected border which seems to instead be copied on other implementations.

I put the final code on Github along with a Rust version that is about 200 times faster to enable testing on larger images.