图像处理中的三次卷积插值(Cubic Convolution)

2021年6月22日 99次阅读来源: struggler_1990

　　图像插值常用于二维图像中的缩放，以及图像畸变校正。由于图像数据量一般都比较大，本文主要参考Cubic Convolution Interpolation for Digital Image Processing一文，实现并分析文中提出的三次卷积插值算法。三次卷积插值算法的精度介于现行插值和三次样条算法之间，但他才用卷积的方式，减少了计算量，而且在两个方向的插值是可分的，是高档图像软件中图像缩放常用的算法。

1. 三次卷积算法

　　由于该算法在图像两个维度上是可分的，这里只介绍一维的三次卷积插值方法，在另一个维度上用同样的方法在做一次即可。下面介绍一维函数的三次卷积插值。
　　假设一维空间中的函数为 f f ，我们的插值函数为 g g ，那么对于每个插值点 xk x k ，有 g(xk)=f(xk) g ( x k ) = f ( x k ) 。对于等间距采样的的数据，许多插值函数可以表示为如下形式：
　　

g(x)=∑kcku(x−xkh)(1) (1) g ( x ) = ∑ k c k u ( x − x k h )

　　在式

(1) ( 1 ) 中，

h h 表示采样间距，

xk x k 表示插值点，

u u 表示插值卷积用的

kernel k e r n e l ，

g g 为插值函数，

ck c k 为依赖于采样点数据的参数，他们必须满足

g(xk)=f(xk) g ( x k ) = f ( x k ) 的条件。

1.1 三次卷积 kernel k e r n e l 的基本公式

　　三次卷积 kernel k e r n e l 满足一下几个条件：
　　1. 对于三次卷积插值，卷积插值内核是定义在子区间 (−2,−1) ( − 2 , − 1 ) ， (−1,0) ( − 1 , 0 ) ， (0,1) ( 0 , 1 ) ， (1,2) ( 1 , 2 ) ，在区间 (−2,2) ( − 2 , 2 ) 外，插值kernel为0。因此，式 (1) ( 1 ) 中用于计算插值函数的数据点减为4个。
　　2. 卷积核必须是对称函数，结合上一个条件，卷积核 u u 可以写成如下形式：
　　

u(s)=⎧⎩⎨⎪⎪⎪⎪⎪⎪1A1|s|3+B1|s|2+C1|s|+D1A2|s|3+B1|s|2+C2|s|+D20s=00<|s|<11<|s|<21=|s|　or　2≤|s|.　　(2) (2) u ( s ) = { 1 s = 0 A 1 | s | 3 + B 1 | s | 2 + C 1 | s | + D 1 0 < | s | < 1 A 2 | s | 3 + B 1 | s | 2 + C 2 | s | + D 2 1 < | s | < 2 0 1 = | s | 　 o r 　 2 ≤ | s | . 　　

　　由于

h h 为采样间隔，可知

xj−xk=(j−k)h x j − x k = ( j − k ) h ,代入

(1) ( 1 ) 式得，

g(x)=∑kcku(j−k)(3) (3) g ( x ) = ∑ k c k u ( j − k )

　　3. 由

g(xk)=f(xk) g ( x k ) = f ( x k ) ，及条件2可知

cj=f(xj) c j = f ( x j ) 。这样不需要像三次样条算法那样，通过求解方程组来获得参数

ck c k ，大大较少计算量。

　　4. 插值kernel必须是连续的，而且一阶导数连续。通过这个条件，将

|s|=0,1,2 | s | = 0 , 1 , 2 代入

u(s),u′(s) u ( s ) , u ′ ( s ) 的解析式，可以获得

7 7 个方程。

　　通过求解方程组，

(2) ( 2 ) 式可写为：

u(s)=⎧⎩⎨⎪⎪⎪⎪⎪⎪1(a+2)|s|3−(a+3)|s|2+1a|s|3−5a|s|2+8a|s|−4a0s=00<|s|<11<|s|<21=|s|　or　2≤|s|.　　(4) (4) u ( s ) = { 1 s = 0 ( a + 2 ) | s | 3 − ( a + 3 ) | s | 2 + 1 0 < | s | < 1 a | s | 3 − 5 a | s | 2 + 8 a | s | − 4 a 1 < | s | < 2 0 1 = | s | 　 o r 　 2 ≤ | s | . 　　　　

　　上式只是一个简单的代数方程组求解，有兴趣的可以参考原文。

1.2 三次卷积 kernel k e r n e l 中 a a 的计算

　　对于某个需要获得值的插值点 x x ，它肯定位于两个样本点 xj x j 和 xj+1 x j + 1 之间。令 s=(x−xj)/h s = ( x − x j ) / h ，有 (x−xk)/h=(x−xj+xj−xk)/h=s+j−k ( x − x k ) / h = ( x − x j + x j − x k ) / h = s + j − k 。因此，式 (1) ( 1 ) 可写为：
　　

g(x)=∑kcku(s+j−k).(5) (5) g ( x ) = ∑ k c k u ( s + j − k ) .

且由于

u u 在区间

(−2,2) ( − 2 , 2 ) ，

0<s<1 0 < s < 1 ,式

(5) ( 5 ) 简化为：

g(x)=cj−1u(s+1)+cju(s)+cj+1u(s−1)+cj+2u(s−2)(6) (6) g ( x ) = c j − 1 u ( s + 1 ) + c j u ( s ) + c j + 1 u ( s − 1 ) + c j + 2 u ( s − 2 )

通过泰勒公式以及关系式

cj=f(xj) c j = f ( x j ) ，可以将

(6) ( 6 ) 式写为：

g(x)=−(2a+1)[2hf′(xj)+h2f′′(xj)]s3+[(6a+3)hf′(xj+(4a+3)h2f′′(xj/2]s2−2ahf′(xj)s+f(xj)+o(h3)(7) (7) g ( x ) = − ( 2 a + 1 ) [ 2 h f ′ ( x j ) + h 2 f ″ ( x j ) ] s 3 + [ ( 6 a + 3 ) h f ′ ( x j + ( 4 a + 3 ) h 2 f ″ ( x j / 2 ] s 2 − 2 a h f ′ ( x j ) s + f ( x j ) + o ( h 3 )

函数

f(x) f ( x ) 的泰勒展开式为：

f(x)=f(xj)+shf′(xj)+s2h2f′′(xj)/2+o(h3)(8) (8) f ( x ) = f ( x j ) + s h f ′ ( x j ) + s 2 h 2 f ″ ( x j ) / 2 + o ( h 3 )

式

(7) ( 7 ) 减

(8) ( 8 ) 可得：

f(x)−g(x)=(2a+1)[2hf′(xj+h2f′′(xj)]s3−(2a+1)[3hf′(xj+h2f′′(xj)]s2+(2a+1)shf′(xj)+o(h3)(9) (9) f ( x ) − g ( x ) = ( 2 a + 1 ) [ 2 h f ′ ( x j + h 2 f ″ ( x j ) ] s 3 − ( 2 a + 1 ) [ 3 h f ′ ( x j + h 2 f ″ ( x j ) ] s 2 + ( 2 a + 1 ) s h f ′ ( x j ) + o ( h 3 )

令

a=−12 a = − 1 2 ,则

(9) ( 9 ) 式为：

f(x)−g(x)=o(h3)(10) (10) f ( x ) − g ( x ) = o ( h 3 )

公式

(10) ( 10 ) 表明，当

a=−12 a = − 1 2 时，

g(x) g ( x ) 与

f(x) f ( x ) 之间的误差，以正比于

h3 h 3 速率趋近于0，即

g(x) g ( x ) 的误差为

o(h3) o ( h 3 ) 。将

a=−12 a = − 1 2 代入

(4) ( 4 ) 有：

u(s)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪132|s|3−52|s|2+1−12|s|3+52|s|2+4|s|+20s=00<|s|<11<|s|<21=|s|　or　2≤|s|.　　(11) (11) u ( s ) = { 1 s = 0 3 2 | s | 3 − 5 2 | s | 2 + 1 0 < | s | < 1 − 1 2 | s | 3 + 5 2 | s | 2 + 4 | s | + 2 1 < | s | < 2 0 1 = | s | 　 o r 　 2 ≤ | s | . 　　　　

1.3 边界条件

　　当 j=0或n−1(假设有n个数据点) j = 0 或 n − 1 ( 假设有 n 个数据点 ) 时，式 (6) ( 6 ) 中的 c−1或cn c − 1 或 c n 无法通过采样的数据点获得，因此需要估计 c−1 c − 1 和 cn c n 的值。
　　为了保证 g(x) g ( x ) 在边界上也能以 o(h3) o ( h 3 ) 接近 f(x) f ( x ) ，通过类似1.2节中的方法可以得出：
　　

　　c−1　　cn==f(x2)−3f(x1)+3f(x0)f(xn−3)−3f(xn−2)+3f(xn−1)　　(12)(13) (12) 　　 c − 1 = f ( x 2 ) − 3 f ( x 1 ) + 3 f ( x 0 ) (13) 　　 c n = f ( x n − 3 ) − 3 f ( x n − 2 ) + 3 f ( x n − 1 ) 　　

　　具体推导过程可以参考原文。

　　具体的图像缩放三次卷积插值过程可以总结为如下过程：

　　1. 设计和计算卷积

kernel k e r n e l （本文为式

(11) ( 11 ) ）。

　　2. 补全插值的边界采样点（本文方法为公式

(12)/(13) ( 12 ) / ( 13 ) ）。

　　3. 计算

x x 方向待插值的点

xi x i ，如果为图像缩放。可以用

xi=x′i/τ x i = x i ′ / τ 的方法计算待插值点，其中

x′i x i ′ 代表缩放后的像素坐标，

xi x i 为缩放前的坐标。

　　4. 通过每个

xi x i 计算

xj x j 和

s s ，然后通过

(11) ( 11 ) 式和

(6) ( 6 ) 式计算获得

g(xi) g ( x i ) 的值。

　　5. 在

y y 方向同样做步骤1-4的过程，即可完成图像的三次卷积插值。

2. 双三次卷积插值的 matlab m a t l a b 实现

　　双三次卷积插值的Matlab代码实现，其cubic_conv.m函数文件为：

function [ scaled_image ] = cubic_conv( input_image, scaled_sz)
% CUBIC_CONV Implement cubic convolution along x direction,
% then along y direction.
% input_image Input image.
% scaled_sz Size of output image.
% scaled_image Output image.
    scaled_image_x = cubic_conv_x(input_image, scaled_sz(2));
    scaled_image = cubic_conv_x(scaled_image_x', scaled_sz(1));
    scaled_image = scaled_image';
end

function [ scaled_image_x ] = cubic_conv_x( input_image, scaled_sz_x )
    %1. Compute the coordinates to be interpolated
    org_sz = [size(input_image,1), size(input_image,2)];
    x_org = (0:scaled_sz_x-1)/(scaled_sz_x-1)*(org_sz(2)-1);
    x = zeros([1, scaled_sz_x, 4]);
    x_org_dec = x_org - floor(x_org);
    x(1,:,1) = cubic_kernel(x_org_dec + 1);
    x(1,:,2) = cubic_kernel(x_org_dec);
    x(1,:,3) = cubic_kernel(1 - x_org_dec);
    x(1,:,4) = cubic_kernel(2 - x_org_dec);
    X = repmat(x, org_sz(1), 1, 1);

    %2. Compute the boudary data
    col_neg_1 = input_image(:,3) - 3*input_image(:,2) + 3*input_image(:,1);
    col_N_plus_1 = input_image(:,end-2) - 3*input_image(:,end-1) + 3*input_image(:,end);
    extent_image = [col_neg_1, input_image, col_N_plus_1];

    x_int = floor(x_org);
    x_int(x_int >= org_sz(2)-1) = org_sz(2)-2;
    % x_int_to_ext_idx = x_int + 1;
    %3. Convolution for the interpolated data
    conv_x = cat(3,extent_image(:,x_int+1), extent_image(:,x_int+2),...
                 extent_image(:,x_int+3), extent_image(:,x_int+4));
    scaled_image_x = dot(X, conv_x, 3);
    scaled_image_x(:,1) = input_image(:,1);
    scaled_image_x(:,end) = input_image(:,end);
end

cubic_conv_x调用的cubic_kernel.m函数文件为：

function [ kernel_res ] = cubic_kernel( s )
    kernel_res = zeros(size(s));
    kernel_res(abs(s)<eps) = 1;
    kernel_res(abs(s-1)<eps) = 0;
    kernel_res(abs(s-2)<eps) = 0;
    kernel_res((s>0 & s < 1)) = kernel_0_1(s(s>0 & s < 1));
    kernel_res((s>1 & s < 2)) = kernel_1_2(s(s>1 & s < 2));
end

function [ kernel_res ] = kernel_0_1( s )
    s2 = s.*s;
    s3 = s2.*s;
    kernel_res = 1.5*s3 - 2.5*s2 + 1;
end

function [ kernel_res ] = kernel_1_2( s )
    s2 = s.*s;
    s3 = s2.*s;
    kernel_res = -0.5*s3 + 2.5*s2 - 4*s + 2;
end

3. 和别的插值算法比较

　　本文主要比较 cubic convolution c u b i c c o n v o l u t i o n 和最近邻插值、双线性插值算法的精度、频谱响应以及计算效率。
　　精度
　　 cubic convolution c u b i c c o n v o l u t i o n 能够完美重建二次多项式，而线性插值只能无误差重建一次多项式，最近邻插值只能重建常数图像。而增加插值的采样数据点，还能提高 cubic convolution c u b i c c o n v o l u t i o n 的精度到重建三次多项式（具体方法可以参考原文）。
　　频响
　　由于 cubic convolution c u b i c c o n v o l u t i o n 本质上是卷积过程，因此可以计算其频响曲线，从而分析其高频丢失以及走样情况。 cubic convolution c u b i c c o n v o l u t i o n 和最近邻插值、双线性插值算法的 kernel k e r n e l 如图1，频响曲线如图2。虚线为理想的频响曲线， nyquist n y q u i s t 频率以上的频响为0， nyquist n y q u i s t 频率以下的频响都为1，既能保持图像的锐度，又不发生走样。从图2可以发现，最近邻插值会在高频出现走样，线性插值在中低频部分相对有较大的频响损失。而 cubic convolution c u b i c c o n v o l u t i o n 的频响最接近理想响应曲线，其对图像的缩放最能保持图像原样，而不发生走样。
　　效率
　　最近邻插值、线性插值以及 cubic convolution c u b i c c o n v o l u t i o n 算法，其计算复杂度是越来越高的。但是 cubic convolution c u b i c c o n v o l u t i o n 在效果和计算效率上能做到比较好的平衡，而且其实现本质是卷积，利于硬件实现，所以能很好的进行加速，所以高档软件以及印刷常用 cubic convolution c u b i c c o n v o l u t i o n 算法。
　　　

《图像处理中的三次卷积插值(Cubic Convolution)》

　　图1 不同算法的kernel

　
《图像处理中的三次卷积插值(Cubic Convolution)》

　　图2 不同算法kernel的频响曲线

4. 插值案例

　　由于 cubic convolution c u b i c c o n v o l u t i o n 本质上是卷积过程其中，虚线为理想的频响曲线， nyquist n y q u i s t 频率以上的频响为0， nyquist n y q u i s t 频率以下的频响都为1，即保持图像的锐度，又不发生走样。从图2可以发现，最近邻插值会在高频出现走样，线性插值在中低频部分相对有较大的频响损失。而 cubic convolution c u b i c c o n v o l u t i o n 的频响最接近理想响应曲线，其对图像的缩放最能保持图像原样，而不发生走样。
　　效率
　　最近邻插值、线性插值以及 cubic convolution c u b i c c o n v o l u t i o n 算法，其计算复杂度是越来越高的。但是 cubic convolution c u b i c c o n v o l u t i o n 在效果和计算效率上能做到比较好的平衡，而且其实现本质是卷积，利于硬件实现，所以能很好的进行加速，所以高档软件以及印刷常用 cubic convolution c u b i c c o n v o l u t i o n 算法。
　　对二维连续函数 f(x,y)=sin(0.5r2) f ( x , y ) = s i n ( 0.5 r 2 ) ，在区间 [0,xmax]∗[y,ymax] [ 0 , x m a x ] ∗ [ y , y m a x ] 上进行采样获得一张分辨率为 64∗64 64 ∗ 64 的小图 I0 I 0 ，其采样间隔为 (h0x,h0y) ( h 0 x , h 0 y ) 。然后在同样的函数区间 [0,xmax]∗[y,ymax] [ 0 , x m a x ] ∗ [ y , y m a x ] 上，用采样间隔为 (h1x,h1y) ( h 1 x , h 1 y ) 的方式，获得分辨率为 350∗336 350 ∗ 336 图像 Iorg I o r g 。
　　然后用 I0 I 0 ，分别通过最近邻、双线性和 cubic convolution c u b i c c o n v o l u t i o n 算法插值获得 Inearest I n e a r e s t 、 Ibilinear I b i l i n e a r 和 Icubic I c u b i c 。
　　《图像处理中的三次卷积插值(Cubic Convolution)》
　　