图像插值常用于二维图像中的缩放,以及图像畸变校正。由于图像数据量一般都比较大,本文主要参考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)=k∑cku(hx−xk)(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)=⎩
⎨
⎧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)
由于
h
h
h为采样间隔,可知
x
j
−
x
k
=
(
j
−
k
)
h
x_{j}-x_{k}=(j-k)h
xj−xk=(j−k)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)=k∑cku(j−k) (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)∣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)
上式只是一个简单的代数方程组求解,有兴趣的可以参考原文。
##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=(x−xj)/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
(x−xk)/h=(x−xj+xj−xk)/h=s+j−k。因此,式
(
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)=k∑cku(s+j−k). (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)=cj−1u(s+1)+cju(s)+cj+1u(s−1)+cj+2u(s−2)(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]s2−2ahf′(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)=⎩
⎨
⎧123∣s∣3−25∣s∣2+1−21∣s∣3+25∣s∣2−4∣s∣+20s=00<∣s∣<11<∣s∣<21=∣s∣ or 2≤∣s∣. (11)
##1.3 边界条件
当
j
=
0
或
n
−
1
(
假设有
n
个数据点
)
j=0或n-1(假设有n个数据点)
j=0或n−1(假设有n个数据点)时,式
(
6
)
(6)
(6)中的
c
−
1
或
c
n
c_{-1}或c_{n}
c−1或cn无法通过采样的数据点获得,因此需要估计
c
−
1
c_{-1}
c−1和
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图](https://i-blog.csdnimg.cn/blog_migrate/dc3c0ae8c19b736b918fa4aea89e0482.jpeg)
![图2 不同算法kernel的频响曲线](https://i-blog.csdnimg.cn/blog_migrate/644d9d2d0d3a4c497519dbc9369f049f.jpeg)
#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
64∗64的小图
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
350∗336图像
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。
![三次卷积插值图像](https://i-blog.csdnimg.cn/blog_migrate/b425ab6588539a26665fd7ee429850f1.jpeg)
![原函数采样图像](https://i-blog.csdnimg.cn/blog_migrate/a2e8faf76e8e99f5bf086857840cb987.jpeg)
可以看出 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](https://i-blog.csdnimg.cn/blog_migrate/fa427fda858a37eb78d0a47ecbd941b8.jpeg)
![双线性error](https://i-blog.csdnimg.cn/blog_migrate/8ca0b74a04c4a333113e9efbfbbae838.jpeg)
![三次卷积error](https://i-blog.csdnimg.cn/blog_migrate/db7772d1f0662dacf55e572563889a19.jpeg)
通过插值图像与原始函数采样的图像的差error,可以看出三次卷积插值误差最小。
#5. 代码链接
1. CSDN资源: cubic convolution/bilinear/nearest image resize matlab
2. Github: cubic_convolution