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

图像插值常用于二维图像中的缩放，以及图像畸变校正。由于图像数据量一般都比较大，本文主要参考Cubic Convolution Interpolation for Digital Image Processing一文，实现并分析文中提出的三次卷积插值算法。三次卷积插值算法的精度介于现行插值和三次样条算法之间，但他才用卷积的方式，减少了计算量，而且在两个方向的插值是可分的，是高档图像软件中图像缩放常用的算法。
#1. 三次卷积算法
　　由于该算法在图像两个维度上是可分的，这里只介绍一维的三次卷积插值方法，在另一个维度上用同样的方法在做一次即可。下面介绍一维函数的三次卷积插值。
　　假设一维空间中的函数为$f$，我们的插值函数为$g$，那么对于每个插值点$x_k$，有$g(x_k)=f(x_k)$。对于等间距采样的的数据，许多插值函数可以表示为如下形式：
　　$$g(x)=\sum _k{c}_{k}u(\frac{x-x_k}{h})\text{\hspace{1em}(1)}$$
　　在式$(1)$中，$h$表示采样间距，$x_k$表示插值点，$u$表示插值卷积用的$kernel$，$g$为插值函数，$c_k$为依赖于采样点数据的参数，他们必须满足$g(x_k)=f(x_k)$的条件。
##1.1 三次卷积$kernel$的基本公式
　　三次卷积$kernel$满足一下几个条件：
　　1. 对于三次卷积插值，卷积插值内核是定义在子区间$(-2,-1)$，$(-1,0)$，$(0,1)$，$(1,2)$，在区间$(-2,2)$外，插值kernel为0。因此，式$(1)$中用于计算插值函数的数据点减为4个。
　　2. 卷积核必须是对称函数，结合上一个条件，卷积核$u$可以写成如下形式：
　　$$u(s)=\{\begin{array}{llc}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|or2\le |s|.\end{array}\text{\hspace{1em}(2)}$$
　　由于$h$为采样间隔，可知$x_j-x_k=(j-k)h$,代入$(1)$式得，
　　$$g(x)=\sum _k{c}_{k}u(j-k)\text{\hspace{1em}(3)}$$
　　3. 由$g(x_k)=f(x_k)$，及条件2可知$c_j=f(x_j)$。这样不需要像三次样条算法那样，通过求解方程组来获得参数$c_k$，大大较少计算量。
　　4. 插值kernel必须是连续的，而且一阶导数连续。通过这个条件，将$|s|=0,1,2$代入$u(s),u^{\prime }(s)$的解析式，可以获得$7$个方程。
　　通过求解方程组，$(2)$式可写为：
$$u(s)=\{\begin{array}{llc}1& & s=0\\ (a+2)|s{|}^3-(a+3)|s{|}^2+1& & 0<|s|<1\\ a|s{|}^3-5a|s{|}^2+8a|s|-4a& & 1<|s|<2\\ 0& & 1=|s|or2\le |s|.\end{array}\text{\hspace{1em}(4)}$$　　
　　上式只是一个简单的代数方程组求解，有兴趣的可以参考原文。
##1.2 三次卷积$kernel$中$a$的计算
　　对于某个需要获得值的插值点$x$，它肯定位于两个样本点$x_j$和$x_{j+1}$之间。令$s=(x-x_j)/h$，有$(x-x_k)/h=(x-x_j+x_j-x_k)/h=s+j-k$。因此，式$(1)$可写为：
　　$$g(x)=\sum _k{c}_{k}u(s+j-k).\text{\hspace{1em}(5)}$$
且由于$u$在区间$(-2,2)$，$0<s<1$,式$(5)$简化为：
$$g(x)=c_{j-1}u(s+1)+c_{j}u(s)+c_{j+1}u(s-1)+c_{j+2}u(s-2)\text{\hspace{1em}(6)}$$
通过泰勒公式以及关系式$c_j=f(x_j)$，可以将$(6)$式写为：
$$g(x)=-(2a+1)[2h{f}^{\prime }(x_j)+h^2{f}^{\prime \prime }(x_j)]s^3+[(6a+3)h{f}^{\prime }(x_j+(4a+3)h^2{f}^{\prime \prime }(x_j/2]s^2-2ah{f}^{\prime }(x_j)s+f(x_j)+o(h^3)\text{\hspace{1em}(7)}$$
函数$f(x)$的泰勒展开式为：
$$f(x)=f(x_j)+sh{f}^{\prime }(x_j)+s^2{h}^2{f}^{\prime \prime }(x_j)/2+o(h^3)\text{\hspace{1em}(8)}$$
式$(7)$减$(8)$可得：
$$f(x)-g(x)=(2a+1)[2h{f}^{\prime }(x_j+h2f^{\prime \prime }(x_j)]s^3-(2a+1)[3h{f}^{\prime }(x_j+h^2{f}^{\prime \prime }(x_j)]s^2+(2a+1)sh{f}^{\prime }(x_j)+o(h^3)\text{\hspace{1em}(9)}$$

令$a=-\frac{1}{2}$,则$(9)$式为：
$$f(x)-g(x)=o(h^3)\text{\hspace{1em}(10)}$$
公式$(10)$表明，当$a=-\frac{1}{2}$时，$g(x)$与$f(x)$之间的误差，以正比于$h^3$速率趋近于0，即$g(x)$的误差为$o(h^3)$。将$a=-\frac{1}{2}$代入$(4)$有：
$$u(s)=\{\begin{array}{llc}1& & s=0\\ \frac{3}{2}|s{|}^3-\frac{5}{2}|s{|}^2+1& & 0<|s|<1\\ -\frac{1}{2}|s{|}^3+\frac{5}{2}|s{|}^2-4|s|+2& & 1<|s|<2\\ 0& & 1=|s|or2\le |s|.\end{array}\text{\hspace{1em}(11)}$$　　
##1.3 边界条件
　　当$j=0或n-1(假设有n个数据点)$时，式$(6)$中的$c_{-1}或c_n$无法通过采样的数据点获得，因此需要估计$c_{-1}$和$c_n$的值。
　　为了保证$g(x)$在边界上也能以$o(h^3)$接近$f(x)$，通过类似1.2节中的方法可以得出：
　　KaTeX parse error: No such environment: eqnarray at position 10: 　　\begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲ 　　c_{-1}&=&f(x…
　　具体推导过程可以参考原文。
　　具体的图像缩放三次卷积插值过程可以总结为如下过程：
　　1. 设计和计算卷积$kernel$（本文为式$(11)$）。
　　2. 补全插值的边界采样点（本文方法为公式$(12)/(13)$）。
　　3. 计算$x$方向待插值的点$x_i$，如果为图像缩放。可以用$x_i=x_i^{\prime }/\tau$的方法计算待插值点，其中$x_i^{\prime }$代表缩放后的像素坐标，$x_i$为缩放前的坐标。
　　4. 通过每个$x_i$计算$x_j$和$s$，然后通过$(11)$式和$(6)$式计算获得$g(x_i)$的值。
　　5. 在$y$方向同样做步骤1-4的过程，即可完成图像的三次卷积插值。
#2. 双三次卷积插值的$matlab$实现
　　双三次卷积插值的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. 和别的插值算法比较
　　本文主要比较$cubicconvolution$和最近邻插值、双线性插值算法的精度、频谱响应以及计算效率。
　　精度
　　$cubicconvolution$能够完美重建二次多项式，而线性插值只能无误差重建一次多项式，最近邻插值只能重建常数图像。而增加插值的采样数据点，还能提高$cubicconvolution$的精度到重建三次多项式（具体方法可以参考原文）。
　　频响
　　由于$cubicconvolution$本质上是卷积过程，因此可以计算其频响曲线，从而分析其高频丢失以及走样情况。$cubicconvolution$和最近邻插值、双线性插值算法的$kernel$如图1，频响曲线如图2。虚线为理想的频响曲线，$nyquist$频率以上的频响为0，$nyquist$频率以下的频响都为1，既能保持图像的锐度，又不发生走样。从图2可以发现，最近邻插值会在高频出现走样，线性插值在中低频部分相对有较大的频响损失。而$cubicconvolution$的频响最接近理想响应曲线，其对图像的缩放最能保持图像原样，而不发生走样。
　　效率
　　最近邻插值、线性插值以及$cubicconvolution$算法，其计算复杂度是越来越高的。但是$cubicconvolution$在效果和计算效率上能做到比较好的平衡，而且其实现本质是卷积，利于硬件实现，所以能很好的进行加速，所以高档软件以及印刷常用$cubicconvolution$算法。
　　　

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

#4. 插值案例
　　由于$cubicconvolution$本质上是卷积过程其中，虚线为理想的频响曲线，$nyquist$频率以上的频响为0，$nyquist$频率以下的频响都为1，即保持图像的锐度，又不发生走样。从图2可以发现，最近邻插值会在高频出现走样，线性插值在中低频部分相对有较大的频响损失。而$cubicconvolution$的频响最接近理想响应曲线，其对图像的缩放最能保持图像原样，而不发生走样。
　　效率
　　最近邻插值、线性插值以及$cubicconvolution$算法，其计算复杂度是越来越高的。但是$cubicconvolution$在效果和计算效率上能做到比较好的平衡，而且其实现本质是卷积，利于硬件实现，所以能很好的进行加速，所以高档软件以及印刷常用$cubicconvolution$算法。
　　对二维连续函数$f(x,y)=sin(0.5r^2)$，在区间$[0,x_{max}]\ast [y,y_{max}]$上进行采样获得一张分辨率为$64\ast 64$的小图$I_0$，其采样间隔为$(h_{0x},h_{0y})$。然后在同样的函数区间$[0,x_{max}]\ast [y,y_{max}]$上，用采样间隔为$(h_{1x},h_{1y})$的方式，获得分辨率为$350\ast 336$图像$I_{org}$。
　　然后用$I_0$，分别通过最近邻、双线性和$cubicconvolution$算法插值获得$I_{nearest}$、$I_{bilinear}$和$I_{cubic}$。
　　
　　

图3 最近邻插值结果（左）和线性插值结果（右）
　　
　　 图4 三次卷积插值结果（左）和函数重采样结果（右）
　　可以看出 $cubicconvlution$最接近原图，下面给出了三个算法和 $I_{org}$的差。
　　
　　 图5 最近邻error（左）和双线性error（右）
　
　　 图6 三次卷积插值error
　　通过插值图像与原始函数采样的图像的差error，可以看出三次卷积插值误差最小。
#5. 代码链接
　　1. CSDN资源： cubic convolution/bilinear/nearest image resize matlab
　　2. Github: cubic_convolution