Bicubic interpolation Python
I have developed Bicubic interpolation for demonstration to some undergraduate students using Python Programming language.
The methodology is as explained in wikipedia, The code is working fine except the results I am getting are slightly different than what is obtained when using scipy library.
The interpolation code is shown below in the function bicubic_interpolation.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from scipy import interpolate
import sympy as syp
import pandas as pd
pd.options.display.max_colwidth = 200
%matplotlib inline
def bicubic_interpolation(xi, yi, zi, xnew, ynew):
# check sorting
if np.any(np.diff(xi) < 0) and np.any(np.diff(yi) < 0) and\
np.any(np.diff(xnew) < 0) and np.any(np.diff(ynew) < 0):
raise ValueError('data are not sorted')
if zi.shape != (xi.size, yi.size):
raise ValueError('zi is not set properly use np.meshgrid(xi, yi)')
z = np.zeros((xnew.size, ynew.size))
deltax = xi[1] - xi[0]
deltay = yi[1] - yi[0]
for n, x in enumerate(xnew):
for m, y in enumerate(ynew):
if xi.min() <= x <= xi.max() and yi.min() <= y <= yi.max():
i = np.searchsorted(xi, x) - 1
j = np.searchsorted(yi, y) - 1
x0 = xi[i-1]
x1 = xi[i]
x2 = xi[i+1]
x3 = x1+2*deltax
y0 = yi[j-1]
y1 = yi[j]
y2 = yi[j+1]
y3 = y1+2*deltay
px = (x-x1)/(x2-x1)
py = (y-y1)/(y2-y1)
f00 = zi[i-1, j-1] #row0 col0 >> x0,y0
f01 = zi[i-1, j] #row0 col1 >> x1,y0
f02 = zi[i-1, j+1] #row0 col2 >> x2,y0
f10 = zi[i, j-1] #row1 col0 >> x0,y1
f11 = p00 = zi[i, j] #row1 col1 >> x1,y1
f12 = p01 = zi[i, j+1] #row1 col2 >> x2,y1
f20 = zi[i+1,j-1] #row2 col0 >> x0,y2
f21 = p10 = zi[i+1,j] #row2 col1 >> x1,y2
f22 = p11 = zi[i+1,j+1] #row2 col2 >> x2,y2
if 0 < i < xi.size-2 and 0 < j < yi.size-2:
f03 = zi[i-1, j+2] #row0 col3 >> x3,y0
f13 = zi[i,j+2] #row1 col3 >> x3,y1
f23 = zi[i+1,j+2] #row2 col3 >> x3,y2
f30 = zi[i+2,j-1] #row3 col0 >> x0,y3
f31 = zi[i+2,j] #row3 col1 >> x1,y3
f32 = zi[i+2,j+1] #row3 col2 >> x2,y3
f33 = zi[i+2,j+2] #row3 col3 >> x3,y3
elif i<=0:
f03 = f02 #row0 col3 >> x3,y0
f13 = f12 #row1 col3 >> x3,y1
f23 = f22 #row2 col3 >> x3,y2
f30 = zi[i+2,j-1] #row3 col0 >> x0,y3
f31 = zi[i+2,j] #row3 col1 >> x1,y3
f32 = zi[i+2,j+1] #row3 col2 >> x2,y3
f33 = f32 #row3 col3 >> x3,y3
elif j<=0:
f03 = zi[i-1, j+2] #row0 col3 >> x3,y0
f13 = zi[i,j+2] #row1 col3 >> x3,y1
f23 = zi[i+1,j+2] #row2 col3 >> x3,y2
f30 = f20 #row3 col0 >> x0,y3
f31 = f21 #row3 col1 >> x1,y3
f32 = f22 #row3 col2 >> x2,y3
f33 = f23 #row3 col3 >> x3,y3
elif i == xi.size-2 or j == yi.size-2:
f03 = f02 #row0 col3 >> x3,y0
f13 = f12 #row1 col3 >> x3,y1
f23 = f22 #row2 col3 >> x3,y2
f30 = f20 #row3 col0 >> x0,y3
f31 = f21 #row3 col1 >> x1,y3
f32 = f22 #row3 col2 >> x2,y3
f33 = f23 #row3 col3 >> x3,y3
px00 = (f12 - f10)/2*deltax
px01 = (f22 - f20)/2*deltax
px10 = (f13 - f11)/2*deltax
px11 = (f23 - f21)/2*deltax
py00 = (f21 - f01)/2*deltay
py01 = (f22 - f02)/2*deltay
py10 = (f31 - f11)/2*deltay
py11 = (f32 - f12)/2*deltay
pxy00 = ((f22-f20) - (f02-f00))/4*deltax*deltay
pxy01 = ((f32-f30) - (f12-f10))/4*deltax*deltay
pxy10 = ((f23-f21) - (f03-f01))/4*deltax*deltay
pxy11 = ((f33-f31) - (f13-f11))/4*deltax*deltay
f = np.array([p00, p01, p10, p11,
px00, px01, px10, px11,
py00, py01, py10, py11,
pxy00, pxy01, pxy10, pxy11])
a = A@f
a = a.reshape(4,4).transpose()
z[n,m] = np.array([1, px, px**2, px**3]) @ a @ np.array([1, py, py**2, py**3])
return z
In the function bicubic_interpolation the inputs are xi= old x data range, yi= old y range, zi= old values at grids points (x,y), xnew, and ynew are the new horizontal data ranges. All inputs are 1D numpy arrays except zi which is 2D numpy array.
The data I am testing the function on are shown below. I am comparing the results as well with scipy and true model (function f).
def f(x,y):
return np.sin(np.sqrt(x ** 2 + y ** 2))
x = np.linspace(-6, 6, 11)
y = np.linspace(-6, 6, 11)
xx, yy = np.meshgrid(x, y)
z = f(xx, yy)
x_new = np.linspace(-6, 6, 100)
y_new = np.linspace(-6, 6, 100)
xx_new, yy_new = np.meshgrid(x_new, y_new)
z_new = bicubic_interpolation(x, y, z, x_new, y_new)
z_true = f(xx_new, yy_new)
f_scipy = interpolate.interp2d(x, y, z, kind='cubic')
z_scipy = f_scipy(x_new, y_new)
fig, ax = plt.subplots(2, 2, sharey=True, figsize=(16,12))
img0 = ax[0, 0].scatter(xx, yy, c=z, s=100)
ax[0, 0].set_title('original points')
fig.colorbar(img0, ax=ax[0, 0], orientation='vertical', shrink=1, pad=0.01)
img1 = ax[0, 1].imshow(z_new, vmin=z_new.min(), vmax=z_new.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[0, 1].set_title('bicubic our code')
fig.colorbar(img1, ax=ax[0, 1], orientation='vertical', shrink=1, pad=0.01)
img2 = ax[1, 0].imshow(z_scipy, vmin=z_scipy.min(), vmax=z_scipy.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[1, 0].set_title('bicubic scipy')
fig.colorbar(img2, ax=ax[1, 0], orientation='vertical', shrink=1, pad=0.01)
img3 = ax[1, 1].imshow(z_true, vmin=z_true.min(), vmax=z_true.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[1, 1].set_title('true model')
fig.colorbar(img3, ax=ax[1, 1], orientation='vertical', shrink=1, pad=0.01)
plt.subplots_adjust(wspace=0.05, hspace=0.15)
plt.show()
The results are shown below:
The matrix A (inside funtion bicubic_interpolation) is as explained in the Wikipedia site and simply can be obtained using the following code:
x = syp.Symbol('x')
y = syp.Symbol('y')
a00, a01, a02, a03, a10, a11, a12, a13 = syp.symbols('a00 a01 a02 a03 a10 a11 a12 a13')
a20, a21, a22, a23, a30, a31, a32, a33 = syp.symbols('a20 a21 a22 a23 a30 a31 a32 a33')
p = a00 + a01*y + a02*y**2 + a03*y**3\
+ a10*x + a11*x*y + a12*x*y**2 + a13*x*y**3\
+ a20*x**2 + a21*x**2*y + a22*x**2*y**2 + a23*x**2*y**3\
+ a30*x**3 + a31*x**3*y + a32*x**3*y**2 + a33*x**3*y**3
px = syp.diff(p, x)
py = syp.diff(p, y)
pxy = syp.diff(p, x, y)
df = pd.DataFrame(columns=['function', 'evaluation'])
for i in range(2):
for j in range(2):
function = 'p({}, {})'.format(j,i)
df.loc[len(df)] = [function, p.subs({x:j, y:i})]
for i in range(2):
for j in range(2):
function = 'px({}, {})'.format(j,i)
df.loc[len(df)] = [function, px.subs({x:j, y:i})]
for i in range(2):
for j in range(2):
function = 'py({}, {})'.format(j,i)
df.loc[len(df)] = [function, py.subs({x:j, y:i})]
for i in range(2):
for j in range(2):
function = 'pxy({}, {})'.format(j,i)
df.loc[len(df)] = [function, pxy.subs({x:j, y:i})]
eqns = df['evaluation'].tolist()
symbols = [a00,a01,a02,a03,a10,a11,a12,a13,a20,a21,a22,a23,a30,a31,a32,a33]
A = syp.linear_eq_to_matrix(eqns, *symbols)[0]
A = np.array(A.inv()).astype(np.float64)
print(df)
print(A)
I would like to know where the problem is with the bicubic_interpolation function and why it is slightly different than the result obtained by scipy?
Any help is greatly appreciated!
Answer
Not sure why Wikipedia implementation is not working as expected. Probably, the reason is that these values might be approximated in a different way than what is explained in their site.
px00 = (f12 - f10)/2*deltax
px01 = (f22 - f20)/2*deltax
px10 = (f13 - f11)/2*deltax
px11 = (f23 - f21)/2*deltax
py00 = (f21 - f01)/2*deltay
py01 = (f22 - f02)/2*deltay
py10 = (f31 - f11)/2*deltay
py11 = (f32 - f12)/2*deltay
pxy00 = ((f22-f20) - (f02-f00))/4*deltax*deltay
pxy01 = ((f32-f30) - (f12-f10))/4*deltax*deltay
pxy10 = ((f23-f21) - (f03-f01))/4*deltax*deltay
pxy11 = ((f33-f31) - (f13-f11))/4*deltax*deltay
However, I found this document which have a different implementation and it is well explained and easier understood than Wikipedia. The result I obtain using this implementation are quite similar to what obtained by SciPy.
def bicubic_interpolation2(xi, yi, zi, xnew, ynew):
# check sorting
if np.any(np.diff(xi) < 0) and np.any(np.diff(yi) < 0) and\
np.any(np.diff(xnew) < 0) and np.any(np.diff(ynew) < 0):
raise ValueError('data are not sorted')
if zi.shape != (xi.size, yi.size):
raise ValueError('zi is not set properly use np.meshgrid(xi, yi)')
z = np.zeros((xnew.size, ynew.size))
deltax = xi[1] - xi[0]
deltay = yi[1] - yi[0]
for n, x in enumerate(xnew):
for m, y in enumerate(ynew):
if xi.min() <= x <= xi.max() and yi.min() <= y <= yi.max():
i = np.searchsorted(xi, x) - 1
j = np.searchsorted(yi, y) - 1
x1 = xi[i]
x2 = xi[i+1]
y1 = yi[j]
y2 = yi[j+1]
px = (x-x1)/(x2-x1)
py = (y-y1)/(y2-y1)
f00 = zi[i-1, j-1] #row0 col0 >> x0,y0
f01 = zi[i-1, j] #row0 col1 >> x1,y0
f02 = zi[i-1, j+1] #row0 col2 >> x2,y0
f10 = zi[i, j-1] #row1 col0 >> x0,y1
f11 = p00 = zi[i, j] #row1 col1 >> x1,y1
f12 = p01 = zi[i, j+1] #row1 col2 >> x2,y1
f20 = zi[i+1,j-1] #row2 col0 >> x0,y2
f21 = p10 = zi[i+1,j] #row2 col1 >> x1,y2
f22 = p11 = zi[i+1,j+1] #row2 col2 >> x2,y2
if 0 < i < xi.size-2 and 0 < j < yi.size-2:
f03 = zi[i-1, j+2] #row0 col3 >> x3,y0
f13 = zi[i,j+2] #row1 col3 >> x3,y1
f23 = zi[i+1,j+2] #row2 col3 >> x3,y2
f30 = zi[i+2,j-1] #row3 col0 >> x0,y3
f31 = zi[i+2,j] #row3 col1 >> x1,y3
f32 = zi[i+2,j+1] #row3 col2 >> x2,y3
f33 = zi[i+2,j+2] #row3 col3 >> x3,y3
elif i<=0:
f03 = f02 #row0 col3 >> x3,y0
f13 = f12 #row1 col3 >> x3,y1
f23 = f22 #row2 col3 >> x3,y2
f30 = zi[i+2,j-1] #row3 col0 >> x0,y3
f31 = zi[i+2,j] #row3 col1 >> x1,y3
f32 = zi[i+2,j+1] #row3 col2 >> x2,y3
f33 = f32 #row3 col3 >> x3,y3
elif j<=0:
f03 = zi[i-1, j+2] #row0 col3 >> x3,y0
f13 = zi[i,j+2] #row1 col3 >> x3,y1
f23 = zi[i+1,j+2] #row2 col3 >> x3,y2
f30 = f20 #row3 col0 >> x0,y3
f31 = f21 #row3 col1 >> x1,y3
f32 = f22 #row3 col2 >> x2,y3
f33 = f23 #row3 col3 >> x3,y3
elif i == xi.size-2 or j == yi.size-2:
f03 = f02 #row0 col3 >> x3,y0
f13 = f12 #row1 col3 >> x3,y1
f23 = f22 #row2 col3 >> x3,y2
f30 = f20 #row3 col0 >> x0,y3
f31 = f21 #row3 col1 >> x1,y3
f32 = f22 #row3 col2 >> x2,y3
f33 = f23 #row3 col3 >> x3,y3
Z = np.array([f00, f01, f02, f03,
f10, f11, f12, f13,
f20, f21, f22, f23,
f30, f31, f32, f33]).reshape(4,4).transpose()
X = np.tile(np.array([-1, 0, 1, 2]), (4,1))
X[0,:] = X[0,:]**3
X[1,:] = X[1,:]**2
X[-1,:] = 1
Cr = [email protected](X)
R = [email protected]([px**3, px**2, px, 1])
Y = np.tile(np.array([-1, 0, 1, 2]), (4,1)).transpose()
Y[:,0] = Y[:,0]**3
Y[:,1] = Y[:,1]**2
Y[:,-1] = 1
Cc = np.linalg.inv(Y)@R
z[n,m]=([email protected]([py**3, py**2, py, 1]))
return z
def f(x,y):
return np.sin(np.sqrt(x ** 2 + y ** 2))
x = np.linspace(-6, 6, 11)
y = np.linspace(-6, 6, 11)
xx, yy = np.meshgrid(x, y)
z = f(xx, yy)
x_new = np.linspace(-6, 6, 100)
y_new = np.linspace(-6, 6, 100)
xx_new, yy_new = np.meshgrid(x_new, y_new)
z_new = bicubic_interpolation2(x, y, z, x_new, y_new)
z_true = f(xx_new, yy_new)
f_scipy = interpolate.interp2d(x, y, z, kind='cubic')
z_scipy = f_scipy(x_new, y_new)
fig, ax = plt.subplots(2, 2, sharey=True, figsize=(16,12))
img0 = ax[0, 0].scatter(xx, yy, c=z, s=100)
ax[0, 0].set_title('original points')
fig.colorbar(img0, ax=ax[0, 0], orientation='vertical', shrink=1, pad=0.01)
img1 = ax[0, 1].imshow(z_new, vmin=z_new.min(), vmax=z_new.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[0, 1].set_title('bicubic our code')
fig.colorbar(img1, ax=ax[0, 1], orientation='vertical', shrink=1, pad=0.01)
img2 = ax[1, 0].imshow(z_scipy, vmin=z_scipy.min(), vmax=z_scipy.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[1, 0].set_title('bicubic scipy')
fig.colorbar(img2, ax=ax[1, 0], orientation='vertical', shrink=1, pad=0.01)
img3 = ax[1, 1].imshow(z_true, vmin=z_true.min(), vmax=z_true.max(), origin='lower',
extent=[x_new.min(), x_new.max(), y_new.max(), y_new.min()])
ax[1, 1].set_title('true model')
fig.colorbar(img3, ax=ax[1, 1], orientation='vertical', shrink=1, pad=0.01)
plt.subplots_adjust(wspace=0.05, hspace=0.15)
plt.show()
fig, ax = plt.subplots(1, 2, sharey=True, figsize=(10, 6))
ax[0].plot(xx[0,:], z[5,:], 'or', label='original')
ax[0].plot(xx_new[0,:], z_true[int(100/10*5),:], label='true')
ax[0].plot(xx_new[0,:], z_new[int(100/10*5), :], label='our interpolation')
ax[1].plot(xx[0,:], z[5,:], 'or', label='original')
ax[1].plot(xx_new[0,:], z_true[int(100/10*5),:], label='true')
ax[1].plot(xx_new[0,:], z_scipy[int(100/10*5), :], label='scipy interpolation')
for axes in ax:
axes.legend()
axes.grid()
plt.show()
