目录
简介
nearest插值核
linear插值核
cubic插值核
lanczos插值核
图像插值几何中心对齐
图像插值
展望
参考资料
简介 主要介绍常用的插值核(interpolation kernel)和对应的插值算法，主要是应用于图像上的，也就是二维的插值，对应的就是图像处理中经常遇到的nearest插值，bilinear插值，bicubic插值和lanzcos插值算法。
nearest插值核顾名思义，就是最靠近哪个点，就用哪个点的值，一般就是用四舍五入的方法来判断。
以缩放比例为ratio为例，基本公式如下：
srcx = dstx / ratio
srcx = (int)(srcx + 0.5)
其插值核如下图所示，可以看到，-0.5~0.5中间恒为1， 其他区域都为0，相位phase不同，其取值也不同，但存在跳变点。
此次的相位计算为
phase=1 - (srcx - (int)srcx)
在做图像nearest插值时，是二维处理的，那就有两个方向，可以理解为核的大小为2x2的矩阵，但这个矩阵很特殊，一定是某一个数值为1， 其他的都为0。
核心代码如下
def GetNearest(x, a=1):
out = np.zeros_like(x)
index = (np.abs(x) <= 0.5)
out[index] = 1
return out
def DrawNearest():
a = 2
xx = np.linspace(-a-2, a+2, 16*(a+2)*2+1)
yynearest = GetNearest(xx, a)
x1 = np.linspace(-a, a, 2 * a, endpoint=False)
phase1 = 0.25
phase2 = 0.75
y0 = GetNearest(x1 + 0, a)
y1 = GetNearest(x1 + phase1, a)
y2 = GetNearest(x1 + phase2, a)
plt.subplot(1, 3, 1)
plt.title('nearest phase=0.0')
plt.plot(xx, yynearest, 'b')
plt.scatter(x1, y0, c='r')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.title('nearest phase=0.25')
plt.plot(xx, yynearest, 'b')
plt.scatter(x1+phase1, y1, c='r')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(xx, yynearest, 'b')
plt.title('nearest phase=0.75')
plt.scatter(x1+phase2, y2, c='r')
plt.grid(True)
plt.show()
s = 1
x0 = np.linspace(-s, s, 2 * s, endpoint=False)
for p1 in range(1, 4, 2):
for p2 in range(1, 4, 2):
yh = GetNearest(x0+p1/4, s)
yw = GetNearest(x0+p2/4, s)
yh = yh.reshape((1, 2*s))
yw = yw.reshape((1, 2*s))
mat0 = yh.T @ yw
print(mat0)代码中，对X和Y方向分别以1/4,3/4为相位，计算得到对应插值核，分别如下：
phasey: 0.25 phasex: 0.25
[[0. 0.]
[0. 1.]]
phasey: 0.25 phasex: 0.75
[[0. 0.]
[1. 0.]]
phasey: 0.75 phasex: 0.25
[[0. 1.]
[0. 0.]]
phasey: 0.75 phasex: 0.75
[[1. 0.]
[0. 0.]]
Linear插值核是取邻近的两个点进行插值得到，靠近哪个点，哪个点的权重就更大，总权重为1，如图所示，p1和p2之间，插值出p，则有：
p=(p1*x2 + p2*x1)/(x1+x2)
Linear插值核如下图所示，同样的算了几个不同的phase值。
对应到图像处理时，是二维的，也就是bilinear插值，那么X和Y方向都计算出系数，然后两个方向乘起来，得到一个2x2的矩阵。
核心代码如下
def GetLinear(x, a=1):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 1 - tmp
return out
def DrawBilinear():
a = 2
xx = np.linspace(-a - 2, a + 2, 16 * (a + 2) * 2 + 1)
yy = GetLinear(xx, a)
x1 = np.linspace(-a, a, 2 * a, endpoint=False)
phase1 = 0.25
phase2 = 0.75
y0 = GetLinear(x1 + 0, a)
y1 = GetLinear(x1 + phase1, a)
y2 = GetLinear(x1 + phase2, a)
plt.subplot(1, 3, 1)
plt.title('linear phase=0.0')
plt.plot(xx, yy, 'b')
plt.scatter(x1, y0, c='r')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.title('linear phase=0.25')
plt.plot(xx, yy, 'b')
plt.scatter(x1 + phase1, y1, c='r')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(xx, yy, 'b')
plt.title('linear phase=0.75')
plt.scatter(x1 + phase2, y2, c='r')
plt.grid(True)
plt.show()
s = 1
x0 = np.linspace(-s, s, 2 * s, endpoint=False)
for p1 in range(1, 4, 2):
for p2 in range(1, 4, 2):
yh = GetLinear(x0 + p1 / 4, s)
yw = GetLinear(x0 + p2 / 4, s)
yh = yh.reshape((1, 2 * s))
yw = yw.reshape((1, 2 * s))
mat0 = yh.T @ yw
print('phasey:', p1 / 4)
print('phasex:', p2 / 4)
print(mat0)同样的，算了相位为0.25和0.75时的插值核矩阵，结果如下：
phasey: 0.25 phasex: 0.25
[[0.0625 0.1875]
[0.1875 0.5625]]
phasey: 0.25 phasex: 0.75
[[0.1875 0.0625]
[0.5625 0.1875]]
phasey: 0.75 phasex: 0.25
[[0.1875 0.5625]
[0.0625 0.1875]]
phasey: 0.75 phasex: 0.75
[[0.5625 0.1875]
[0.1875 0.0625]]
cubic插值核[1]使用了三次方运算，对应的核大小和Linear不一样，Linear是参考邻近的2个点，cubic至少参考4个点，还有一种是6个点，那么对应的二维图像处理时，对应的bicubic插值核分别为4x4和6x6，由一维变二维，只需要把两个一维向量相乘即可。
4个节点的公式如下
其中a取-0.5时，有
至于是如何得到的，可以参考下论文，主要还是根据方程求解出来的系数，最后只剩下了a。
其对应的插值核如下图所示。
核心代码如下
def GetCubic(x, a=2):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 1.5 * tmp * tmp * tmp - 2.5 * tmp * tmp + 1
index10 = (np.abs(x) <= 2)
index11 = (np.abs(x) > 1)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = -0.5 * tmp * tmp * tmp + 2.5 * tmp * tmp - 4 * tmp + 2
return out
def DrawBiCubic():
a = 3
xx = np.linspace(-a - 2, a + 2, 16 * (a + 2) * 2 + 1)
yy = GetCubic(xx, a)
x1 = np.linspace(-a, a, 2 * a, endpoint=False)
phase1 = 0.25
phase2 = 0.75
y0 = GetCubic(x1 + 0, a)
y1 = GetCubic(x1 + phase1, a)
y2 = GetCubic(x1 + phase2, a)
plt.subplot(1, 3, 1)
plt.title('cubic phase=0.0')
plt.plot(xx, yy, 'b')
plt.scatter(x1, y0, c='r')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.title('cubic phase=0.25')
plt.plot(xx, yy, 'b')
plt.scatter(x1 + phase1, y1, c='r')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(xx, yy, 'b')
plt.title('cubic phase=0.75')
plt.scatter(x1 + phase2, y2, c='r')
plt.grid(True)
plt.show()6个节点的公式如下：
其插值核如下
核心代码如下
def GetCubic3(x, a=3):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 4 * tmp * tmp * tmp / 3 - 7 * tmp * tmp / 3 + 1
index10 = (np.abs(x) <= 2)
index11 = (np.abs(x) > 1)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = -7 * tmp * tmp * tmp / 12 + 3 * tmp * tmp - 59 * tmp / 12 + 15/6
index10 = (np.abs(x) <= 3)
index11 = (np.abs(x) > 2)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = tmp * tmp * tmp / 12 - 2 * tmp * tmp / 3 + 21 * tmp / 12 - 3 / 2
return out
def DrawBiCubic3():
a = 4
xx = np.linspace(-a - 2, a + 2, 16 * (a + 2) * 2 + 1)
yy = GetCubic3(xx, a)
x1 = np.linspace(-a, a, 2 * a, endpoint=False)
phase1 = 0.25
phase2 = 0.75
y0 = GetCubic3(x1 + 0, a)
y1 = GetCubic3(x1 + phase1, a)
y2 = GetCubic3(x1 + phase2, a)
plt.subplot(1, 3, 1)
plt.title('cubic3 phase=0.0')
plt.plot(xx, yy, 'b')
plt.scatter(x1, y0, c='r')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.title('cubic3 phase=0.25')
plt.plot(xx, yy, 'b')
plt.scatter(x1 + phase1, y1, c='r')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(xx, yy, 'b')
plt.title('cubic3 phase=0.75')
plt.scatter(x1 + phase2, y2, c='r')
plt.grid(True)
plt.show() lanczos插值核[2]主要是使用了sinc函数，构造了lanzcos函数，其公式如下，其中a表示窗口的半径，不同a，得到的插值核形状不一样，参考的点的数目也不一样，参考点数目为2*a个。
其插值核如下图所示，a=2时
a=3时
a=4时
可以发现，不同的a，曲线上的花瓣的数目也不一样，为2*a-1个瓣。
而且和cubic插值核对比，a=2时，和cubic的4个节点很接近，a=3时，和cubic的6个节点很接近。
对比曲线如下
其核心代码如下：
def GetLanczosCoef(x, a):
out = np.zeros_like(x)
out[x == 0] = 1
index = (np.abs(x) <= a)
index1 = (x != 0)
index2 = np.zeros_like(index)
for i in range(len(index)):
index2[i] = index[i] and index1[i]
# print(index2)
out[index2] = a * np.sin(np.pi * x[index2]) * np.sin(np.pi * x[index2] / a) / (np.pi * np.pi * x[index2] * x[index2])
return out
def DrawLanczos():
a = 4
xx = np.linspace(-a - 2, a + 2, 16 * (a + 2) * 2 + 1)
yy = GetLanczosCoef(xx, a)
x1 = np.linspace(-a, a, 2 * a, endpoint=False)
phase1 = 0.25
phase2 = 0.75
y0 = GetLanczosCoef(x1 + 0, a)
y1 = GetLanczosCoef(x1 + phase1, a)
y2 = GetLanczosCoef(x1 + phase2, a)
plt.subplot(1, 3, 1)
plt.title('lanczos a=4 phase=0.0')
plt.plot(xx, yy, 'b')
plt.scatter(x1, y0, c='r')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.title('lanczos a=4 phase=0.25')
plt.plot(xx, yy, 'b')
plt.scatter(x1 + phase1, y1, c='r')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(xx, yy, 'b')
plt.title('lanczos a=4 phase=0.75')
plt.scatter(x1 + phase2, y2, c='r')
plt.grid(True)
plt.show() 介绍了4种插值的方法，可以发现，如果相位为0时，得到的系数都是1，但是这样就会出现图像不对齐的问题[3]。
还是以缩放ratio为例来说明，
srcx=dstx / ratio
srcy=dsty / ratio
如果源图像和目标图像的原点（0，0）均选择左上角，然后根据插值公式计算目标图像每点像素，假设你需要将一幅5x5的图像缩小成3x3，那么源图像和目标图像各个像素之间的对应关系如下：
也就是说，图像当以左上角为对齐点来看，会使偏右偏下的原像素点的价值变低，新插值点受偏左偏上的原像素点影响更大。
那么，让坐标加1或者选择右下角为原点怎么样呢？很不幸，还是一样的效果，不过这次得到的图像将偏右偏下。
最好的方法就是，两个图像的几何中心重合，并且目标图像的每个像素之间都是等间隔的，并且都和两边有一定的边距。
公式如下：
srcx=(dstx + 0.5) / ratio - 0.5
srcy=(dsty +0.5) / ratio - 0.5
图像插值 根据上面介绍的几种插值方法，和几何中心对齐，然后实现图像插值，并且和Python Image Library模块里的结果尽心对比。
主要代码如下，由于跑了for循环，整体速度还是比较慢，代码以2倍上采样为例。
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.image as img
from scipy.ndimage import filters, measurements, interpolation
import glob
from scipy.io import savemat
import ntpath
import math
import cv2
from skimage.metrics import mean_squared_error as compare_mse
from skimage.metrics import peak_signal_noise_ratio as compare_psnr
from skimage.metrics import structural_similarity as compare_ssim
from PIL import Image
def GetLanczosCoef(x, a):
out = np.zeros_like(x)
out[x == 0] = 1
index = (np.abs(x) <= a)
index1 = (x != 0)
index2 = np.zeros_like(index)
for i in range(len(index)):
index2[i] = index[i] and index1[i]
# print(index2)
out[index2] = a * np.sin(np.pi * x[index2]) * np.sin(np.pi * x[index2] / a) / (np.pi * np.pi * x[index2] * x[index2])
return out
def GetCubic(x, a=2):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 1.5 * tmp * tmp * tmp - 2.5 * tmp * tmp + 1
index10 = (np.abs(x) <= 2)
index11 = (np.abs(x) > 1)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = -0.5 * tmp * tmp * tmp + 2.5 * tmp * tmp - 4 * tmp + 2
return out
def GetCubic3(x, a=3):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 4 * tmp * tmp * tmp / 3 - 7 * tmp * tmp / 3 + 1
index10 = (np.abs(x) <= 2)
index11 = (np.abs(x) > 1)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = -7 * tmp * tmp * tmp / 12 + 3 * tmp * tmp - 59 * tmp / 12 + 15/6
index10 = (np.abs(x) <= 3)
index11 = (np.abs(x) > 2)
index1 = np.zeros_like(index10)
for i in range(len(index10)):
index1[i] = index10[i] and index11[i]
tmp = np.abs(x[index1])
out[index1] = tmp * tmp * tmp / 12 - 2 * tmp * tmp / 3 + 21 * tmp / 12 - 3 / 2
return out
def GetLinear(x, a=1):
out = np.zeros_like(x)
index = (np.abs(x) <= 1)
tmp = np.abs(x[index])
out[index] = 1 - tmp
return out
def GetNearest(x, a=1):
out = np.zeros_like(x)
index = (np.abs(x) <= 0.5)
out[index] = 1
return out
# mode 0 nearest
# 1 bilinear
# 2 bicubic 2
# 3 bicubic 3
# 4 lanczos
def MyInterpolation(img, ratio, a, mode):
h, w, c = img.shape
dsth = (int)(h*ratio)
dstw = (int)(w*ratio)
if mode == 0:
a = 1
fun = GetNearest
elif mode == 1:
a = 1
fun = GetLinear
elif mode == 2:
a = 2
fun = GetCubic
elif mode == 3:
a = 3
fun = GetCubic3
else:
fun = GetLanczosCoef
imgBorder = cv2.copyMakeBorder(img, a, a, a, a, cv2.BORDER_REPLICATE)
imgout = np.zeros((dsth, dstw, c), np.uint8)
x1 = np.linspace(-a, a, 2*a, endpoint=False)
for j in range(dsth):
srcj = (j + 0.5) / ratio - 0.5 + a
cj = (int)(srcj)
phasej = 1 - (srcj - cj)
yj = fun(x1 + phasej, a)
yj = yj.reshape((1, 2*a))
for i in range(dstw):
srci = (i + 0.5) / ratio - 0.5 + a
ci = (int)(srci)
phasei = 1 - (srci - ci)
yi = fun(x1 + phasei, a)
yi = yi.reshape((1, 2*a))
mat = yj.T @ yi
coefs = np.sum(mat)
b = imgBorder[cj - a + 1: cj + a + 1, ci - a + 1:ci + a + 1, 0] * mat
g = imgBorder[cj - a + 1: cj + a + 1, ci - a + 1:ci + a + 1, 1] * mat
r = imgBorder[cj - a + 1: cj + a + 1, ci - a + 1:ci + a + 1, 2] * mat
b = np.sum(b) / coefs
g = np.sum(g) / coefs
r = np.sum(r) / coefs
imgout[j, i, 0] = np.clip(b, 0, 255)
imgout[j, i, 1] = np.clip(g, 0, 255)
imgout[j, i, 2] = np.clip(r, 0, 255)
return imgout
def TestImage():
img = cv2.imread('lennax2.png', 1)
gt = cv2.imread('lenna.png', 1)
print(img.shape)
print(gt.shape)
h, w, c = img.shape
img1 = Image.fromarray(img)
imgoutpil = img1.resize((w * 2, h * 2), Image.LANCZOS)
imgoutpil = np.array(imgoutpil)
psnrpil = compare_psnr(gt, imgoutpil)
cv2.imwrite('out-pil-lanzcos.png', imgoutpil)
imgout = MyInterpolation(img, 2, 3, 4)
cv2.imwrite('out-my-lanzcos3.png', imgout)
psnr = compare_psnr(gt, imgout)
print(psnrpil, psnr)
if __name__ == '__main__':
np.set_printoptions(precision=6, suppress=True)
TestImage()结果图如下，
和PIL模块里的bicubic插值对比，容差设置的为1，可以看到差异很小，当然不知道哪里原因，还是有些差异，基本上也符合预期。
还进行了PSNR对比，结果如下，可以看到，基本上和PIL模块还是相似的。
展望
这几种插值方法在图像处理中确实应用较多，根据不同的应用场景，使用不同的方法，窗口越大，虽然效果可以增加，但运算量和成本也就越高。
当然，这几种方法，哪怕是最好的插值算法lanzcos也是不能完全满足实际需要的，有的场景效果要求更高。
可以考虑在此基础上进行改进，比如，边缘方向插值，在强边缘处，考虑沿着边缘方向插值，这样可以更好的保护边缘。
这篇directional cubic convolution interpolation[4]可以作为参考，不过我也没有找到一种通用的边缘插值方法，就是把这些插值方式直接和边缘插值结合起来，例如，45°直线，就降插值核旋转45°之类的，目前没有找到相应的方法。
现在深度学习的方法在SR方法可以取的很不错的效果，比这种要好很多，但这些插值方法在实际中还是有很多需要。
参考资料
[1] Cubic Convolution Interpolation for Digital Image Processing
[2] https://en.wikipedia.org/wiki/Lanczos_resampling
[3] https://blog.csdn.net/q923714892/article/details/118343736
[4] image zooming using directional cubic convolution interpolation
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）