python非均匀三次样条函数拟合 python 三次样条插值

文章目录

• 三次样条插值

• 基本过程
• 代码验证

• INTER_AREA

• 基本逻辑
• FAST_AREA

• 验证代码：

• AREA

• 验证代码

三次样条插值

接上一篇，在opencv中，对于放大图像，还有一个参数是：INTER_CUBIC。这个枚举值代表了另外一种插值方法：三次样条插值。

这个方法念上去挺拗口，其实逻辑上和前一篇讲到的双线性插值的总体过程是差不多的，只是在具体计算目标图像的像素值时，使用的计算方法不一样。

参考：

双线性插值

基本过程

• 对于每个目标图像中的点，通过
  $srcX=dstX(src width /dst width )+0.5(src width /dst width −1)$
  $srcY=dstY(src height /dst height )+0.5(src width /dst width −1)$
  两个公式来计算在原图中的坐标，这个和双线性插值是一样的。
• 接下来就是通过什么样的逻辑来计算目标像素的问题，在双线性插值中，实际上是利用了映射坐标点 (srcX, srcY)所在的周边的四个原图像素点，利用这四个点的像素值在X和Y上做了三次插值计算，就可以得到目标图像像素点的像素值。
• 在线性插值算法中，把点到点之间的像素值关系考虑成线性的关系。而在三次样条插值算法中，目标像素的像素值取决于映射点(srcX, srcY)周边的4 * 4的区域。如图：

目标图像中的像素值，需要根据点P周边的16个像素点的像素值来进行加权计算。根据上一步，位置的映射关系已经确定了，那么重点就在于确定每个点的权重了。

• 权重计算是根据一个三次函数来计算的，估计这就是被成为三次样条的原因吧，函数为：

• 公式中的a在opencv 4.5中，取值为-0.75。
• 至于为什么要用这么个函数去计算，有兴趣的同学可以去了解一下，什么一阶二阶导数都存在且连续什么的，我是不太搞得明白。
• 函数中的x代表的是上图中的16个原图像像素点到映射点P的距离。
• 距离是在X和Y的方向上分开计算。比如目标图像中的像素点$(5, 5)$，经过上面的坐标映射后，在原图中P点的坐标为$(1.33 1.33)$，那么可以确定16个点的X坐标和Y坐标的列表都是：$(0, 1, 2, 3)$，总共组成16个点，
• 那么(1, 1)这个原图的像素点到P点的距离就是
  $(dx = abs(1 - 1.33), dy = abs(1 - 1.33)$
  再把这个值 0.33 带入到上面的公式中去进行计算，可以得到这个点的W(dx)和W(dy)，那么在整个计算公式中，像素点(1, 1)的影响值就是：
  $f(1, 1) * W(d_x) * W(d_y)$
• 把16个像素点都计算一遍，再全部相加，就是目标图像像素点的像素值(如果计算出来为负数，则赋值为0)

代码验证

我的验证逻辑是，使用opencv的resize函数，使用调试模式观察输出的img数据，然后再根据自己的逻辑写一段代码，比较输出是否相同。

img_large = cv2.resize(img, dst_size, 3.0, 3.0, interpolation=cv2.INTER_CUBIC)

调试模式得到的结果：

1. 首先是W(x)的这个权重函数，从opencv的源码中抠出来的，做了一点修改：

def interpolateCubic(x):

    A = -0.75

    if x < 1:
        result = ((A + 2)*x - (A + 3))*x*x + 1
    elif 1 < x < 2:
        # result = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A
        result = ((A * x - 5 * A) * x + 8 * A) * x - 4 * A
    else:
        result = 0

    return result

2. 坐标的映射，因为是一个demo代码，所以是写死一个像素点，我用的是坐标 (5, 5)，从上图中可以看到，像素值为24。

test_x = 5
test_y = 5

# 目标图像坐标往原图映射
sx = (test_x + 0.5) * 1/3 - 0.5
sy = (test_y + 0.5) * 1 / 3 - 0.5

2. 计算距离，为了更方便阅读，代码写的比较啰嗦

src_x = [0,0,0,0]

src_x[0] = math.floor(sx)
src_x[1] = math.floor(sx) - 1
src_x[2] = math.ceil(sx)
src_x[3] = math.ceil(sx) + 1

src_y = [0,0,0,0]
src_y[0] = math.floor(sy) -1
src_y[1] = math.floor(sy)
src_y[2] = math.ceil(sy)
src_y[3] = math.ceil(sy) + 1

2. 根据距离，计算权重，X和Y的是分开计算的。

src_x_weight = [0,0,0,0]
src_x_weight[0] = interpolateCubic(math.fabs(src_x[0] - sx))
src_x_weight[1] = interpolateCubic(math.fabs(src_x[1] - sx))
src_x_weight[2] = interpolateCubic(math.fabs(src_x[2] - sx))
src_x_weight[3] = interpolateCubic(math.fabs(src_x[3] - sx))

src_y_weight = [0,0,0,0]
src_y_weight[0] = interpolateCubic(math.fabs(src_y[0] - sy))
src_y_weight[1] = interpolateCubic(math.fabs(src_y[1] - sy))
src_y_weight[2] = interpolateCubic(math.fabs(src_y[2] - sy))
src_y_weight[3] = interpolateCubic(math.fabs(src_y[3] - sy))

2. 加权计算得到目标像素值

dst_pix = img[src_y[0], src_x[0]] * src_x_weight[0] * src_y_weight[0] \
          + img[src_y[0], src_x[1]] * src_x_weight[1] * src_y_weight[0] \
          + img[src_y[0], src_x[2]] * src_x_weight[2] * src_y_weight[0] \
          + img[src_y[0], src_x[3]] * src_x_weight[3] * src_y_weight[0] \
          + img[src_y[1], src_x[0]] * src_x_weight[1] * src_y_weight[1] \
          + img[src_y[1], src_x[1]] * src_x_weight[1] * src_y_weight[1] \
          + img[src_y[1], src_x[2]] * src_x_weight[2] * src_y_weight[1] \
          + img[src_y[1], src_x[3]] * src_x_weight[3] * src_y_weight[1] \
          + img[src_y[2], src_x[0]] * src_x_weight[1] * src_y_weight[2] \
          + img[src_y[2], src_x[1]] * src_x_weight[1] * src_y_weight[2] \
          + img[src_y[2], src_x[2]] * src_x_weight[2] * src_y_weight[2] \
          + img[src_y[2], src_x[3]] * src_x_weight[3] * src_y_weight[2] \
          + img[src_y[3], src_x[0]] * src_x_weight[1] * src_y_weight[3] \
          + img[src_y[3], src_x[1]] * src_x_weight[1] * src_y_weight[3] \
          + img[src_y[3], src_x[2]] * src_x_weight[2] * src_y_weight[3] \
          + img[src_y[3], src_x[3]] * src_x_weight[3] * src_y_weight[3]

print(round(dst_pix))

最终得到的结果是 24。和opencv里面的逻辑保持一直。

2. 还可以计算(4, 4)，(5, 6)，(6, 6)的像素值，都是一致的。当然，源码中还有一些边界处理的过程，我这里就不阐述了，有兴趣的可以去翻一翻opencv的源码。

INTER_AREA

这个插值算法主要用于图像缩小用，如果是放大，会使用类似线性插值的方式去实现，在opencv代码里有描述：

// true "area" interpolation is only implemented for the case (scale_x >= 1 && scale_y >= 1).
        // In other cases it is emulated using some variant of bilinear interpolation

至于具体是什么逻辑，不是这篇文章所关注的，也就是说如果是放大图像的话，就不要使用INTER_AREA这个枚举值。

基本逻辑

和这个插值算法名称一样，算法的基本逻辑就是根据缩小的倍数来确定将整个源图像划分成多少个区域。

• 比如10 * 10的区域，在X和Y的方向上都缩小2倍的话，总体上就是缩小4倍。那么算法就会以scaleX * scaleY，也就是2 * 2为一个区域，将整个原始图像划分成25个5 * 5区域。每个区域的像素平均值就作为新图像的像素值。

用excel简单画了个图，大概就是这么个概念，那么根据缩小的倍数是否能将源图分成完整的区域，也就是width/scalex和height/scaleY是否是整数，在opencv中是分成了两种处理方式：

• FAST_AREA
• AREA

FAST_AREA

如果width/scalex和height/scaleY都是整数的话，比如上面提到的10 * 10的区域在X和Y的方向上都缩小到原来的一半，那么就是分成25个区域，每个新区域都包含了原来的4个像素点，那么新像素点就是这4个像素点的平均值。

验证代码：

def demoShrink():
    img = np.zeros((10, 10, 3), dtype="uint8")
    cv2.rectangle(img, (2, 2), (7, 7), (255, 255, 255), 1)

    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    dst_size = (round(img.shape[0] * 0.5), round(img.shape[1] * 0.5))
    img_shrink = cv2.resize(img, dst_size, 0.5, 0.5, interpolation=cv2.INTER_AREA)

    print(img)
    print(img_shrink)

输出结果：

img:

[[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 255 255 255 255 255 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 255 255 255 255 255 0 0]
[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0]]

img_shrink：

[[ 0 0 0 0 0]
[ 0 191 128 191 0]
[ 0 128 0 128 0]
[ 0 191 128 191 0]
[ 0 0 0 0 0]]

在这个例子里，img_shrink的(1， 1)对应的AREA就是由(2, 2), (2, 3), (3, 2), (3, 3)组成的区域，算下来的平均值就是191，其他的部分可以照这个方法进行计算。

AREA

如果width/scalex和height/scaleY不是整数的话，起始逻辑上也差不多。比如我们将上面的10 * 10的区域分成 9份，也就是scaleX = 3，scaleY = 3。就无法被整除了。

为什么说逻辑是一样的呢，看看刚才那张图，我们看一下X方向上的划分，每个区域应该就是10/3个像素：

比如说第一个区域，源图在X方向上的第四个像素src(4, 0)，可以这里理解这个像素的位置，这个像素值有1/3在目标图的第一个区域内，那么就把他的像素值的1/3计算进去，也就是说，如果不是整除的话，计算逻辑和可以整除是一样的，都是区域内所有像素的平均值，只不过边上像素只是像素值的一部分，也就是说，目标图的第一个像素dst(0 ,0)的计算公式是：
$dst(0,0) = \frac{src(0, 0) + src(0, 1) + src(0, 2) + src(1, 0) + src(1, 1) + src(1, 2) + src(2, 0) + src(2, 1) + src(2, 2) + src(0, 3)/scaleY + src(1, 3)/scaleY + src(2, 3)/scaleY + src(3, 0)/scaleY + src(3, 1)/scaleY + src(3, 2)/scaleY + src(3, 3)/scaleY/scaleX}{scaleX * scaleY}$
将上面的数据代入之后可以得到第一个数字是38.25，截断成像素就是38。

验证代码

将scale换成0.3

def demoShrink():
    img = np.zeros((10, 10, 3), dtype="uint8")
    cv2.rectangle(img, (2, 2), (7, 7), (255, 255, 255), 1)

    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    dst_size = (round(img.shape[0] * 0.3), round(img.shape[1] * 0.3))
    img_shrink = cv2.resize(img, dst_size, 0.3, 0.3, interpolation=cv2.INTER_AREA)

    print(img)
    print(img_shrink)

结果和我们的计算一致：

img：

[[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 255 255 255 255 255 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 0 0 0 0 255 0 0]
[ 0 0 255 255 255 255 255 255 0 0]
[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0]]

img_shrink:

[[38 76 38]
[76 0 76]
[38 76 38]]