Which cubic spline method does scipy ndimage use for affine_transform?

Question

Background/Context

I was trying to reproduce the values output by scipy's ndimage.affine_transform function, but it seems I am using a different "cubic" interpolation scheme compared to the scipy implementation.

Example

Let's take a look at a very simple example (not data you would want to use cubic interpolation for, but easy to understand). To check the values I implemented uniform Catmull-Rom splines. My small implementation example:

import numpy as np
from scipy.ndimage import affine_transform


def catmull_rom_interp(p0, p1, p2, p3, x):
    return (
        (-0.5 * p0 + 1.5 * p1 - 1.5 * p2 + 0.5 * p3) * (x ** 3)
        + (p0 - 2.5 * p1 + 2 * p2 - 0.5 * p3) * (x ** 2)
        + (-0.5 * p0 + 0.5 * p2) * x
        + p1
    )


image = np.zeros((9,))
image[3] = 13.3

scipy_result_filtered = affine_transform(
    image, np.eye(1), offset=-1.7, order=3, prefilter=True
)
scipy_result = affine_transform(image, np.eye(1), offset=-1.7, order=3, prefilter=False)

image_padded = np.pad(image, 3, mode="constant", constant_values=0)
result_manual = np.zeros((9,))

for i in range(9):
    result_manual[i] = catmull_rom_interp(*image_padded[i : i + 4], 0.3)

print(scipy_result)
print(scipy_result_filtered)
print(result_manual)

# yields
# [0. 0. 0.          0.05985    4.63061667  7.84921667  0.76031667 0.          0.        ]
# [0. 0. 0.1675183  -1.06094923 4.43537861 11.10313479 -1.75261778 0.46923634 -0.12432758]
# [0. 0. 0.         -0.41895    3.85035    10.84615    -0.97755    0.          0.        ]


#
#    PLOTTING
#

import matplotlib.pyplot as plt

plt.gca().grid()

plots = []
for i in range(9):
    plots.append(lambda x: catmull_rom_interp(*image_padded[i : i + 4], x))

plt.plot(scipy_result, "--", label="scipy", alpha=0.5)
plt.plot(scipy_result, "o", color=plt.get_cmap("tab10")(0))

plt.plot(scipy_result_filtered, "--", label="scipy filtered", alpha=0.5)
plt.plot(scipy_result_filtered, "o", color=plt.get_cmap("tab10")(1))
plt.plot(result_manual, "o")
for i in range(9):
    plt.plot(
        np.linspace(i - 0.3, i + 1 - 0.3, 100),
        plots[i](np.linspace(0, 1, 100)),
        "--",
        alpha=0.5,
        color=plt.get_cmap("tab10")(2),
        label="Catmull-Rom spline" if i == 0 else None,
    )

plt.plot(
    np.arange(-0.3, 8.8),
    [0] * 2 + list(image[:-1]),
    "o",
    label="Data to interpolate",
    color="k",
)


plt.legend(framealpha=1)
plt.show()

will yield the following plot (note that due to not knowing the true interpolation function for the scipy functions I just plotted linear connections to better highlight the different data-points):

Observations:

• The scipy method does not use Catmull-Rom splines
• The scipy method (without filtering) does not produce overshoot which is typically associated with cubic interpolation of sharp edges, but as mentioned in the scipy documentation does lead to some blurring, also this seems related to the specific shift of the image I used in the example
• The prefiltered scipy method is closer to Catmull-Rom splines but not identical (there are visible differences)

Question

• Which interpolation scheme does scipy use?
• Where to actually find it in their documentation and/or code?
• Bonus: A simple way to implement it (for checking purposes) in Python

Answer 1

Because I can't comment yet, the source code for the function is located at: https://github.com/scipy/scipy/blob/master/scipy/ndimage/interpolation.py#L355

It appears to do some basic transformation / error checking and then feed into either zoomShift or geometricTransform depending on the parameters.

I don't have enough insight to answer more to 1 or 3 unfortunately.