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

图像插值常用于二维图像中的缩放,以及图像畸变校正。由于图像数据量一般都比较大,本文主要参考Cubic Convolution Interpolation for Digital Image Processing一文,实现并分析文中提出的三次卷积插值算法。三次卷积插值算法的精度介于现行插值和三次样条算法之间,但他才用卷积的方式,减少了计算量,而且在两个方向的插值是可分的,是高档图像软件中图像缩放常用的算法。
#1. 三次卷积算法
  由于该算法在图像两个维度上是可分的,这里只介绍一维的三次卷积插值方法,在另一个维度上用同样的方法在做一次即可。下面介绍一维函数的三次卷积插值。
  假设一维空间中的函数为 f f f,我们的插值函数为 g g g,那么对于每个插值点 x k x_{k} xk,有 g ( x k ) = f ( x k ) g(x_{k})=f(x_{k}) g(xk)=f(xk)。对于等间距采样的的数据,许多插值函数可以表示为如下形式:
   g ( x ) = ∑ k c k u ( x − x k h ) (1) g(x)=\sum_{k}c_{k}u(\frac{x-x_{k}}{h}) \tag{1} g(x)=kcku(hxxk)(1)
  在式 ( 1 ) (1) (1)中, h h h表示采样间距, x k x_{k} xk表示插值点, u u u表示插值卷积用的 k e r n e l kernel kernel g g g为插值函数, c k c_{k} ck为依赖于采样点数据的参数,他们必须满足 g ( x k ) = f ( x k ) g(x_{k})=f(x_{k}) g(xk)=f(xk)的条件。
##1.1 三次卷积 k e r n e l kernel kernel的基本公式
  三次卷积 k e r n e l kernel kernel满足一下几个条件:
  1. 对于三次卷积插值,卷积插值内核是定义在子区间 ( − 2 , − 1 ) (-2,-1) (2,1) ( − 1 , 0 ) (-1,0) (1,0) ( 0 , 1 ) (0,1) (0,1) ( 1 , 2 ) (1,2) (1,2),在区间 ( − 2 , 2 ) (-2,2) (2,2)外,插值kernel为0。因此,式 ( 1 ) (1) (1)中用于计算插值函数的数据点减为4个。
  2. 卷积核必须是对称函数,结合上一个条件,卷积核 u u u可以写成如下形式:
      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 ∣ .       (2)   u(s)=\begin{cases} 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| or 2\leq|s|. \tag{2}   \end{cases}      u(s)= 1A1s3+B1s2+C1s+D1A2s3+B1s2+C2s+D20s=00<s<11<s<21=s or 2s∣.    (2)
  由于 h h h为采样间隔,可知 x j − x k = ( j − k ) h x_{j}-x_{k}=(j-k)h xjxk=(jk)h,代入 ( 1 ) (1) (1)式得,
      g ( x ) = ∑ k c k u ( j − k )    (3)   g(x)=\sum_{k}c_{k}u(j-k) \tag{3}      g(x)=kcku(jk)  (3)
  3. 由 g ( x k ) = f ( x k ) g(x_{k})=f(x_{k}) g(xk)=f(xk),及条件2可知 c j = f ( x j ) c_{j}=f(x_{j}) cj=f(xj)。这样不需要像三次样条算法那样,通过求解方程组来获得参数 c k c_{k} ck,大大较少计算量。
  4. 插值kernel必须是连续的,而且一阶导数连续。通过这个条件,将 ∣ s ∣ = 0 , 1 , 2 |s|=0,1,2 s=0,1,2代入 u ( s ) , u ′ ( s ) u(s),u'(s) u(s),u(s)的解析式,可以获得 7 7 7个方程。
  通过求解方程组, ( 2 ) (2) (2)式可写为:
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 ∣ .    (4) u(s)=\begin{cases} 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| or 2\leq|s|. \tag{4}   \end{cases} u(s)= 1(a+2)s3(a+3)s2+1as35as2+8as4a0s=00<s<11<s<21=s or 2s∣.  (4)  
  上式只是一个简单的代数方程组求解,有兴趣的可以参考原文。
##1.2 三次卷积 k e r n e l kernel kernel a a a的计算
  对于某个需要获得值的插值点 x x x,它肯定位于两个样本点 x j x_{j} xj x j + 1 x_{j+1} xj+1之间。令 s = ( x − x j ) / h s=(x-x_{j})/h s=(xxj)/h,有 ( x − x k ) / h = ( x − x j + x j − x k ) / h = s + j − k (x-x_{k})/h=(x-x_{j}+x_{j}-x_{k})/h=s+j-k (xxk)/h=(xxj+xjxk)/h=s+jk。因此,式 ( 1 ) (1) (1)可写为:
      g ( x ) = ∑ k c k u ( s + j − k ) .    (5)   g(x)=\sum_{k}c_{k}u(s+j-k). \tag{5}      g(x)=kcku(s+jk).  (5)
且由于 u u u在区间 ( − 2 , 2 ) (-2,2) (2,2) 0 < s < 1 0<s<1 0<s<1,式 ( 5 ) (5) (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 ) (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) \tag{6} g(x)=cj1u(s+1)+cju(s)+cj+1u(s1)+cj+2u(s2)(6)
通过泰勒公式以及关系式 c j = f ( x j ) c_{j}=f(x_{j}) cj=f(xj),可以将 ( 6 ) (6) (6)式写为:
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 ) (7) g(x)=-(2a+1)[2hf'(x_{j})+h^2f''(x_{j})]s^3 +[(6a+3)hf'(x_{j}+(4a+3)h^2f''(x_{j}/2]s^2 -2ahf'(x_{j})s+f(x_{j})+o(h^3) \tag{7} g(x)=(2a+1)[2hf(xj)+h2f′′(xj)]s3+[(6a+3)hf(xj+(4a+3)h2f′′(xj/2]s22ahf(xj)s+f(xj)+o(h3)(7)
函数 f ( x ) f(x) f(x)的泰勒展开式为:
f ( x ) = f ( x j ) + s h f ′ ( x j ) + s 2 h 2 f ′ ′ ( x j ) / 2 + o ( h 3 ) (8) f(x)=f(x_{j})+shf'(x_{j})+s^2h^2f''(x_{j})/2+o(h^3) \tag{8} f(x)=f(xj)+shf(xj)+s2h2f′′(xj)/2+o(h3)(8)
( 7 ) (7) (7) ( 8 ) (8) (8)可得:
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 ) (9) f(x)-g(x)=(2a+1)[2hf'(x_{j}+h2f''(x_{j})]s^3 - (2a+1)[3hf'(x_{j}+h^2f''(x_{j})]s^2 +(2a+1)shf'(x_{j})+o(h^3) \tag{9} 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)

a = − 1 2 a=-\frac{1}{2} a=21,则 ( 9 ) (9) (9)式为:
f ( x ) − g ( x ) = o ( h 3 ) (10) f(x)-g(x)=o(h^3) \tag{10} f(x)g(x)=o(h3)(10)
公式 ( 10 ) (10) (10)表明,当 a = − 1 2 a=-\frac{1}{2} a=21时, g ( x ) g(x) g(x) f ( x ) f(x) f(x)之间的误差,以正比于 h 3 h^3 h3速率趋近于0,即 g ( x ) g(x) g(x)的误差为 o ( h 3 ) o(h^3) o(h3)。将 a = − 1 2 a=-\frac{1}{2} a=21代入 ( 4 ) (4) (4)有:
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 ∣ .    (11) u(s)=\begin{cases} 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| or 2\leq|s|. \tag{11}   \end{cases} u(s)= 123s325s2+121s3+25s24∣s+20s=00<s<11<s<21=s or 2s∣.  (11)  
##1.3 边界条件
  当 j = 0 或 n − 1 ( 假设有 n 个数据点 ) j=0或n-1(假设有n个数据点) j=0n1(假设有n个数据点)时,式 ( 6 ) (6) (6)中的 c − 1 或 c n c_{-1}或c_{n} c1cn无法通过采样的数据点获得,因此需要估计 c − 1 c_{-1} c1 c n c_{n} cn的值。
  为了保证 g ( x ) g(x) g(x)在边界上也能以 o ( h 3 ) o(h^3) o(h3)接近 f ( x ) f(x) 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. 设计和计算卷积 k e r n e l kernel kernel(本文为式 ( 11 ) (11) (11))。
  2. 补全插值的边界采样点(本文方法为公式 ( 12 ) / ( 13 ) (12)/(13) (12)/(13))。
  3. 计算 x x x方向待插值的点 x i x_{i} xi,如果为图像缩放。可以用 x i = x i ′ / τ x_{i}=x'_{i}/\tau xi=xi/τ的方法计算待插值点,其中 x i ′ x'_{i} xi代表缩放后的像素坐标, x i x_{i} xi为缩放前的坐标。
  4. 通过每个 x i x_{i} xi计算 x j x_{j} xj s s s,然后通过 ( 11 ) (11) (11)式和 ( 6 ) (6) (6)式计算获得 g ( x i ) g(x_{i}) g(xi)的值。
  5. 在 y y y方向同样做步骤1-4的过程,即可完成图像的三次卷积插值。
#2. 双三次卷积插值的 m a t l a b matlab 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. 和别的插值算法比较
  本文主要比较 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution和最近邻插值、双线性插值算法的精度、频谱响应以及计算效率。
  精度
   c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution能够完美重建二次多项式,而线性插值只能无误差重建一次多项式,最近邻插值只能重建常数图像。而增加插值的采样数据点,还能提高 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution的精度到重建三次多项式(具体方法可以参考原文)。
  频响
  由于 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution本质上是卷积过程,因此可以计算其频响曲线,从而分析其高频丢失以及走样情况。 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution和最近邻插值、双线性插值算法的 k e r n e l kernel kernel如图1,频响曲线如图2。虚线为理想的频响曲线, n y q u i s t nyquist nyquist频率以上的频响为0, n y q u i s t nyquist nyquist频率以下的频响都为1,既能保持图像的锐度,又不发生走样。从图2可以发现,最近邻插值会在高频出现走样,线性插值在中低频部分相对有较大的频响损失。而 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution的频响最接近理想响应曲线,其对图像的缩放最能保持图像原样,而不发生走样。
  效率
  最近邻插值、线性插值以及 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution算法,其计算复杂度是越来越高的。但是 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution在效果和计算效率上能做到比较好的平衡,而且其实现本质是卷积,利于硬件实现,所以能很好的进行加速,所以高档软件以及印刷常用 c u b i c   c o n v o l u t i o n cubic\ convolution cubic convolution算法。
   

图1 不同算法的kernel图

  
图1 不同算法的kernel

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

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

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

图3 最近邻插值结果(左)和线性插值结果(右)

   三次卷积插值图像 原函数采样图像
  
图4 三次卷积插值结果(左)和函数重采样结果(右)

  可以看出 c u b i c   c o n v l u t i o n cubic\ convlution cubic convlution最接近原图,下面给出了三个算法和 I o r g I_{org} Iorg的差。
   最近邻error 双线性error
  
图5 最近邻error(左)和双线性error(右)

 
三次卷积error

  
图6 三次卷积插值error

  通过插值图像与原始函数采样的图像的差error,可以看出三次卷积插值误差最小。
#5. 代码链接
  1. CSDN资源: cubic convolution/bilinear/nearest image resize matlab
  2. Github: cubic_convolution

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值