【图像处理】最近邻插值、双线性插值和双三次插值
但是大家只是在用这个算法，很少关注这个算法的实现细节：插值算法是如何工作的。上图效果是最近邻法的计算过程示意图，由上图可见，最近邻法不需要计算只需要寻找原图中对应的点，所以最近邻法速度最快，但是会破坏原图像中像素的渐变关系，原图像中的像素点的值是渐变的，但是在新图像中局部破坏了这种渐变关系。最近邻法实际上是不需要计算新图像矩阵中点的数值的，直接找到原图像中对应的点，将数值赋值给新图像矩阵中的点，根
目录
1.什么是插值
Interpolation is a method of constructing new data points within the range of a discrete set of known data points. Image interpolation refers to the“guess”of intensity values at missing locations.
图片放大是图像处理中的一个特别基础的操作。在几乎每一个图片相关的项目中，从传统图像处理到深度学习，都有应用。生活里，和朋友通过微信传张图片，从图片发出，到朋友收到图片，查看图片，都会数次的的改变图像的尺寸，从而用到这个算法。但是大家只是在用这个算法，很少关注这个算法的实现细节：插值算法是如何工作的。
简单来说，插值指利用已知的点来“猜”未知的点，图像领域插值常用在修改图像尺寸的过程，由旧的图像矩阵中的点计算新图像矩阵中的点并插入，不同的计算过程就是不同的插值算法。
下图是自己实现双线性插值的效果
2.常用的插值算法
插值算法有很多种，这里列出关联比较密切的三种：
最近邻法(Nearest Interpolation)：计算速度最快，但是效果最差。
双线性插值(Bilinear Interpolation)：双线性插值是用原图像中4(2*2)个点计算新图像中1个点，效果略逊于双三次插值，速度比双三次插值快，属于一种平衡美，在很多框架中属于默认算法。
双三次插值(Bicubic interpolation)：双三次插值是用原图像中16(4*4)个点计算新图像中1个点，效果比较好，但是计算代价过大。
3.最近邻法(Nearest Interpolation)
双线性插值法由原图中4个点计算新图中的1个点，在介绍计算过程前，需要先了解如何找到这4个点，双线性插值法找寻4个点的方式和最近邻法相似，这里顺带的了解下最近邻法的计算流程。
最近邻法实际上是不需要计算新图像矩阵中点的数值的，直接找到原图像中对应的点，将数值赋值给新图像矩阵中的点，根据对应关系找到原图像中的对应的坐标，这个坐标可能不是整数，这时候找最近的点进行插值。对应关系如下：
变量含义如下：
import os
import numpy as np
import cv2
def nearest_factor(image, factor_w, factor_h):
height, width, _ = image.shape
height_new = int(height * factor_h)
width_new = int(width * factor_w)
image_new = np.zeros((height_new, width_new, 3), dtype="uint8")
for i in range(height_new):
for j in range(width_new):
ori_i = round(i / factor_h)
ori_j = round(j / factor_w)
if ori_i >= height:
ori_i = height - 1
if ori_j >= width:
ori_j = width - 1
image_new[i, j, :] = image[ori_i, ori_j, :]
return image_new
def nearest_size(image, height_new, width_new):
height, width, _ = image.shape
image_new = np.zeros((height_new, width_new, 3), dtype="uint8")
factor_h = height_new / height
factor_w = width_new / width
for i in range(height_new):
for j in range(width_new):
ori_i = round(i / factor_h)
ori_j = round(j / factor_w)
if ori_i >= height:
ori_i = height - 1
if ori_j >= width:
ori_j = width - 1
image_new[i, j, :] = image[ori_i, ori_j, :]
return image_new
if __name__ == "__main__":
image_path = r"1.bmp"
image = cv2.imread(image_path, -1)
new_image = nearest_factor(image, 0.33333, 0.333333)
cv2.imwrite(os.path.join("nearest_factor.bmp"), new_image)3.1总结
上图效果是最近邻法的计算过程示意图，由上图可见，最近邻法不需要计算只需要寻找原图中对应的点，所以最近邻法速度最快，但是会破坏原图像中像素的渐变关系，原图像中的像素点的值是渐变的，但是在新图像中局部破坏了这种渐变关系。
3.2双线性插值对应关系
双线性插值的对应公式和前面的最近邻法一样，不一样的是根据对应关系不再是找最近的1个点，而是找最近的4个点，如下图所示。
这里可能还会有点小疑问，如果根据对应关系找到原图中的点不是在不同的点之间，而是跟原图像中的点重合，那该如何找4个点？关于这个问题要看一下双线性插值的计算公式，双线性插值实际上是从2个方向一共进行了3次单线性插值，咱们先了解单线性插值的计算方式。
4.单线性插值
已知中P1点和P2点，坐标分别为(x1, y1)、(x2, y2)，要计算 [x1, x2] 区间内某一位置 x 在直线上的y值
根据初中的知识，2点求一条直线公式(这是双线性插值所需要的唯一的基础公式)
经过简单整理成下面的格式：
这里没有写成经典的AX+B的形式，因为这种形式从权重的角度更好理解。
首先看分子，分子可以看成x与x1和x2的距离作为权重，这也是很好理解的，P点与P1、P2点符合线性变化关系，所以P离P1近就更接近P1，反之则更接近P2。
现在再把公式中的分式看成一个整体，原式可以理解成y1与y2是加权系数，如何理解这个加权，要返回来思考一下，咱们先要明确一下根本的目的：咱们现在不是在求一个公式，而是在图像中根据2个点的像素值求未知点的像素值。这样一个公式是不满足咱们写代码的要求的。
现在根据实际的目的理解，就很好理解这个加权了，y1与y2分别代表原图像中的像素值，上面的公式可以写成如下形式：
5.双线性插值
已知Q11(x1,y1)、Q12(x1,y2)、Q21(x2,y1)、Q22(x2,y2)，求其中点P(x,y)的值。
前面介绍过双线性插值是分别在两个方向计算了共3次单线性插值，如图所示，先在x方向求2次单线性插值，获得R1(x, y1)、R2(x, y2)两个临时点，再在y方向计算1次单线性插值得出P(x, y)（实际上调换2次轴的方向先y后x也是一样的结果）。
1.x方向单线性插值 直接带入前一步单线性插值最后的公式
2.y方向单线性插值
将第一步结果带入第二步
回顾一下上面双线性插值对应关系的图，不难发现，在计算中有这样的关系：
那么上面的公式中的分母全都为0，如下：
在有些资料中，会写成权重的形式，上面的展开式是下面的权重表达式的正确求法
这种权重的表达式也不难理解，观察一下可以发现每个点的权重都和待求点和对角点的距离有关，比如：f(Q11)的权重与f(Q22)的坐标有关。
5.1遗留问题
上面遗留了个小问题，就是当对应关系公式带入后跟原图像中的点发生重合后怎么选取剩下3个点，根据上面的公式可以发现，无论我们怎么选取，其实其余3点的权重都至少1项为0，所以不论我们怎么取剩下的3个点，对最终的结果都不会产生影响。
5.2总结
到此双线性插值的计算已经完成了，但是我们现在只能说是勉强的完成双线性插值算法，为什么这么说，现在的计算方式有比较大的问题，看下面的内容。
6.双线性插值的优化
原始公式：
双线性插值的对应关系看似比较清晰，但还是有2个问题，我画图片进行说明。
根据坐标系的不同，产生的结果不同 这张图是左上角为坐标系原点的情况，我们可以发现最左边x=0的点都会有概率直接复制到目标图像中（至少原点肯定是这样），而且就算不不和原图像中的点重合，也相当于进行了1次单线性插值（仔细想想为什么，带入到权重公式中会发现结果）
现在这张图是右上角为坐标系原点的情况，我们可以发现最右面的点都会有概率直接复制到目标图像中（至少原点肯定是这样），而且就算不不和原图像中的点重合，也相当于进行了1次单线性插值。这样如果我们采用不用的坐标系产生的结果是不一样的，而且无论我们采用什么坐标系，最左侧和最右侧（最上侧和最下侧）的点是不“公平的”，这是第一个问题。
整体的图像相对位置会发生变化看下面这张图，左侧是原图像(33)，右侧是目标图像(55)，原图像的几何中心点是(1, 1)，目标图像的几何中心点是(2, 2)，根据对应关系，目标图像的几何中心点对应的原图像的位置是(1.2, 1.2)，如图所示，那么问题来了，目标图像的原点(0, 0)点和原始图像的原点是重合的，但是目标图像的几何中心点相对于原始图像的几何中心点偏右下，那么整体图像的位置会发生偏移，为什么这样说，其实图像是由1个个的像素点组成，单纯说1个像素点是没有太大的意义的，1个像素点跟相邻像素点的值的渐变或者突变形成图像颜色的渐变或者边界，所以参与计算的点相对都往右下偏移会产生相对的位置信息损失。这是第二个问题。
6.1计算机视觉中的蝴蝶效应
其实对于咱们人眼，上面的2个问题都不会产生太大的结果，但是对于现在的基于学习的学习类算法，计算机通过卷积神经网络提取图像中的深层信息，在这个过程中，我们人眼难以发现的变化也许会发生想象之外的变化，所以不要小看这两个问题
6.2解决方案
几何中心点重合对应公式：
再带入到上面的情况，可以发现问题解决了，到此，才算完成完整的双线性插值法，当然如果这样计算发现结果跟OpenCV的结果不一样，是因为OpenCV还进行了很多速度上的优化，比如用整形计算代替浮点数计算。
import math
import numpy as np
import cv2
import os
def bilinear_factor(image, factor_w, factor_h):
height, width, _ = image.shape
height_new = int(height * factor_h)
width_new = int(width * factor_w)
image_new = np.zeros((height_new, width_new, 3), dtype="uint8")
for i in range(height_new):
for j in range(width_new):
ori_i = (i + 0.5) / factor_h - 0.5
ori_j = (j + 0.5) / factor_w - 0.5
if ori_i < 0:
ori_i = 0
if ori_j < 0:
ori_j = 0
if ori_i >= height - 1:
ori_i = height - 2
if ori_j >= width - 1:
ori_j = width - 2
ori_i_int = math.floor(ori_i)
ori_j_int = math.floor(ori_j)
ori_i_float = ori_i - ori_i_int
ori_j_float = ori_j - ori_j_int
if ori_i_int +1 == height or ori_j_int +1 == width:
image_new[i, j, :] = image[ori_i_int, ori_j_int, :]
else:
image_new[i, j, :] = (1 - ori_i_float) * (1 - ori_j_float) * image[ori_i_int, ori_j_int, :] + \
ori_i_float * (1 - ori_j_float) * image[ori_i_int + 1, ori_j_int, :] + \
(1 - ori_i_float) * ori_j_float * image[ori_i_int, ori_j_int + 1, :] + \
ori_i_float * ori_j_float * image[ori_i_int + 1, ori_j_int + 1, :]
return image_new
if __name__ == "__main__":
image_path = r"1.bmp"
image = cv2.imread(image_path, -1)
new_image = bilinear_factor(image, 0.33333, 0.33333)
cv2.imwrite(os.path.join("bilinear_factor.bmp"), new_image)
print("end")7. 双三次插值
7.1理论
双三次插值也被称为三线性插值，立方卷积插值，三次卷积插值等。
该算法利用待采样点周围16个点的灰度值作三次插值，不仅考虑到4 个直接相邻点的灰度影响，而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果，但也导致了运算量的急剧增加。
假设源图像A大小为mn，缩放K倍后的目标图像B的大小为MN，即K=M/m。A的每一个像素点是已知的，B是未知的，我们想要求出目标图像B中每一像素点(X,Y)的值，必须先找出像素(X,Y)在源图像A中对应的像素(x,y)，再根据源图像A距离像素(x,y)最近的16个像素点作为计算目标图像B(X,Y)处像素值的参数，利用BiCubic基函数求出16个像素点的权重，图B像素(X,Y)的值就等于16个像素点的加权叠加。
根据比例关系x/X=m/M=1/K，我们可以得到B(X,Y)在A上的对应坐标为A(x,y)=A(X/K,Y/K)。如图所示P点就是目标图像B在(X,Y)处对应于源图像A中的位置，P的坐标位置会出现小数部分，所以我们假设 P的坐标为P(x+u,y+v)，其中x,y分别表示整数部分，u,v分别表示小数部分（蓝点到a11方格中红点的距离）。那么我们就可以得到如图所示的最近16个像素的位置，在这里用a(i,j)(i,j=0,1,2,3)来表示，如上图。
BiCubic函数：
BiCubic基函数是一维的，而像素是二维的，所以我们将像素点的行与列分开计算。
BiCubic函数中的参数x表示该像素点到P点的距离，例如a00距离P(x+u,y+v)的距离为(1+u,1+v)，因此a00的横坐标权重i_0=W(1+u)，纵坐标权重j_0=W(1+v)。
重点：距离是权重
a00对B(X,Y)的贡献值为：（a00像素值）* i_0* j_0。
因此，P点横坐标权重分别为W(1+u)，W(u)，W(1-u)，W(2-u)；纵坐标权重分别为W(1+v)，W(v)，W(1-v)，W(2-v)。
写成矩阵形式，a也可以不取-0.5，但是原作者的建议是-0.5。
如将a设置为-1，S（x）是三次插值核函数，权重函数，可由如下公式近似：
7.2python实现
import numpy as np
import cv2
# 三次卷积的核函数
def S(x):
if abs(x) <= 1:
y = 1 - 2 * np.power(x, 2) + abs(np.power(x, 3))
elif abs(x) > 1 and abs(x) < 2:
y = 4 - 8 * abs(x) + 5 * np.power(x, 2) - abs(np.power(x, 3))
else:
y = 0
return y
# 三次卷积插值
def bicubic_interpolation(src, dst_shape):
# 获取原图维度
src_height, src_width = src.shape[0], src.shape[1]
# 计算新图维度 注意channel数要想同
dst_height, dst_width, channels = dst_shape[0], dst_shape[1], dst_shape[2]
dst = np.zeros(shape=(dst_height, dst_width, channels), dtype=np.uint8)
for dst_x in range(dst_height):
for dst_y in range(dst_width):
# 寻找源图像对应坐标
src_x = (dst_x + 0.5) * (src_width / dst_width) - 0.5
src_y = (dst_y + 0.5) * (src_width / dst_width) - 0.5
if src_x < 0:
src_x = 0
if src_y < 0:
src_y = 0
i, j = int(src_x), int(src_y)
u, v = src_x - i, src_y - j
# 边界条件
x1 = min(max(0, i - 1), src_height - 4)
x2 = x1 + 4
y1 = min(max(0, j - 1), src_width - 4)
y2 = y1 + 4
# 计算双三次插值
A = np.array([S(u + 1), S(u), S(u - 1), S(u - 2)])
C = np.array([S(v + 1), S(v), S(v - 1), S(v - 2)])
B = src[x1:x2, y1:y2]
f0 = [A @ B[..., i] @ C.T for i in range(channels)]
f1 = np.stack(f0)
f = np.clip(f1, 0, 255) # 处理一下越界的数据
# 插值
dst[dst_x, dst_y, :] = f.astype(np.uint8)
return dst
# 三次卷积插值，改进边界
def bicubic_interpolation_boundary(src, dst_shape):
# 获取原图维度
src_height, src_width = src.shape[0], src.shape[1]
# 计算新图维度 注意channel数要想同
dst_height, dst_width, channels = dst_shape[0], dst_shape[1], dst_shape[2]
src_pad = np.pad(src, ((2, 2), (2, 2), (0, 0)), 'symmetric')
dst = np.zeros(shape=(dst_height, dst_width, channels), dtype=np.uint8)
for dst_x in range(dst_height):
for dst_y in range(dst_width):
# 寻找源图像对应坐标
src_x = (dst_x + 0.5) * (src_width / dst_width) - 0.5
src_y = (dst_y + 0.5) * (src_width / dst_width) - 0.5
i, j = int(src_x), int(src_y)
u, v = src_x - i, src_y - j
x1 = i - 1 + 2
x2 = x1 + 4
y1 = j - 1 + 2
y2 = y1 + 4
# 计算双三次插值
A = np.array([S(u + 1), S(u), S(u - 1), S(u - 2)])
C = np.array([S(v + 1), S(v), S(v - 1), S(v - 2)])
B = src_pad[x1:x2, y1:y2]
f0 = [A @ B[..., i] @ C.T for i in range(channels)]
f1 = np.stack(f0)
f = np.clip(f1, 0, 255) # 处理一下越界的数据
# 插值
dst[dst_x, dst_y, :] = f.astype(np.uint8)
return dst
if __name__ == "__main__":
root_path = r"1.bmp"
img = cv2.imread(root_path, -1)
interpolated_img = bicubic_interpolation_boundary(img, (719, 1279, 3))
cv2.imwrite(r"cubic3.bmp", interpolated_img)
参考
开放原子开发者工作坊旨在鼓励更多人参与开源活动，与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动，如meetup、训练营等，主打技术交流，干货满满，真诚地邀请各位开发者共同参与！
更多推荐
0
0- 0
已为社区贡献27条内容
所有评论(0)