原理：详细参见https://en.wikipedia.org/wiki/Bicubic_interpolation
在数学上，双立方插值是2个变量函数在方形网格上立方插值的一个扩展，在图像处理中，双立方插值考虑插值像素点周围16个像素，当不在乎算法执行速度时，一般选择双立方插值，相对最近邻域插值和二次线性插值，双立方插值具有较少的插值畸变（artifacts ，不知道翻译成什么好）。
双立方插值与最近邻域插值和二次线性插值的比较如下图：
p
(
x
,
y
)
p(x,y)
p(x,y)插值可由下式给出：
p
(
x
,
y
)
=
[
1
x
x
2
x
3
]
[
a
00
a
01
a
02
a
03
a
10
a
11
a
12
a
13
a
20
a
21
a
22
a
23
a
30
a
31
a
32
a
33
]
[
1
y
y
2
y
3
]
.
{\displaystyle p(x,y)={\begin{bmatrix}1&x&x^{2}&x^{3}\end{bmatrix}}{\begin{bmatrix}a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}\end{bmatrix}}{\begin{bmatrix}1\\y\\y^{2}\\y^{3}\end{bmatrix}}.}
p(x,y)=[1xx2x3]⎣⎢⎢⎡a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a33⎦⎥⎥⎤⎣⎢⎢⎡1yy2y3⎦⎥⎥⎤.
差分计算可以采用下面公式，参考：
f
x
(
x
,
y
)
≈
f
(
x
+
h
,
y
)
−
f
(
x
−
h
,
y
)
2
h
f
y
(
x
,
y
)
≈
f
(
x
,
y
+
k
)
−
f
(
x
,
y
−
k
)
2
k
f
x
x
(
x
,
y
)
≈
f
(
x
+
h
,
y
)
−
2
f
(
x
,
y
)
+
f
(
x
−
h
,
y
)
h
2
f
y
y
(
x
,
y
)
≈
f
(
x
,
y
+
k
)
−
2
f
(
x
,
y
)
+
f
(
x
,
y
−
k
)
k
2
f
x
y
(
x
,
y
)
≈
f
(
x
+
h
,
y
+
k
)
−
f
(
x
+
h
,
y
−
k
)
−
f
(
x
−
h
,
y
+
k
)
+
f
(
x
−
h
,
y
−
k
)
4
h
k
{\displaystyle {\begin{aligned}f_{x}(x,y)&\approx {\frac {f(x+h,y)-f(x-h,y)}{2h}}\\f_{y}(x,y)&\approx {\frac {f(x,y+k)-f(x,y-k)}{2k}}\\f_{xx}(x,y)&\approx {\frac {f(x+h,y)-2f(x,y)+f(x-h,y)}{h^{2}}}\\f_{yy}(x,y)&\approx {\frac {f(x,y+k)-2f(x,y)+f(x,y-k)}{k^{2}}}\\f_{xy}(x,y)&\approx {\frac {f(x+h,y+k)-f(x+h,y-k)-f(x-h,y+k)+f(x-h,y-k)}{4hk}}\end{aligned}}}
fx(x,y)fy(x,y)fxx(x,y)fyy(x,y)fxy(x,y)≈2hf(x+h,y)−f(x−h,y)≈2kf(x,y+k)−f(x,y−k)≈h2f(x+h,y)−2f(x,y)+f(x−h,y)≈k2f(x,y+k)−2f(x,y)+f(x,y−k)≈4hkf(x+h,y+k)−f(x+h,y−k)−f(x−h,y+k)+f(x−h,y−k)
h
=
1
h=1
h=1
另外一种实现如下：
上面的实现对每个
p
(
x
,
y
)
p(x,y)
p(x,y)都要求解一个线性方程，可以采用卷积方式达到同样的目的，卷积核如下：
W
(
x
)
=
{
(
a
+
2
)
∣
x
∣
3
−
(
a
+
3
)
∣
x
∣
2
+
1
for
∣
x
∣
≤
1
,
a
∣
x
∣
3
−
5
a
∣
x
∣
2
+
8
a
∣
x
∣
−
4
a
for
1
<
∣
x
∣
<
2
,
0
otherwise
,
{\displaystyle W(x)={\begin{cases}(a+2)|x|^{3}-(a+3)|x|^{2}+1&{\text{for }}|x|\leq 1,\\a|x|^{3}-5a|x|^{2}+8a|x|-4a&{\text{for }}1<|x|<2,\\0&{\text{otherwise}},\end{cases}}}
W(x)=⎩⎪⎨⎪⎧(a+2)∣x∣3−(a+3)∣x∣2+1a∣x∣3−5a∣x∣2+8a∣x∣−4a0for ∣x∣≤1,for 1<∣x∣<2,otherwise,
a
=
−
0.5
a=-0.5
a=−0.5。
上式可以表示成矩阵形式：
p
(
t
)
=
1
2
[
1
t
t
2
t
3
]
[
0
2
0
0
−
1
0
1
0
2
−
5
4
−
1
−
1
3
−
3
1
]
[
f
−
1
f
0
f
1
f
2
]
{\displaystyle p(t)={\tfrac {1}{2}}{\begin{bmatrix}1&t&t^{2}&t^{3}\\\end{bmatrix}}{\begin{bmatrix}0&2&0&0\\-1&0&1&0\\2&-5&4&-1\\-1&3&-3&1\\\end{bmatrix}}{\begin{bmatrix}f_{-1}\\f_{0}\\f_{1}\\f_{2}\\\end{bmatrix}}}
p(t)=21[1tt2t3]⎣⎢⎢⎡0−12−120−53014−300−11⎦⎥⎥⎤⎣⎢⎢⎡f−1f0f1f2⎦⎥⎥⎤
对图像进行插值时，可以先进行x方向插值，然后再进行y方向插值。
b
−
1
=
p
(
t
x
,
f
(
−
1
,
−
1
)
,
f
(
0
,
−
1
)
,
f
(
1
,
−
1
)
,
f
(
2
,
−
1
)
)
,
b
0
=
p
(
t
x
,
f
(
−
1
,
0
)
,
f
(
0
,
0
)
,
f
(
1
,
0
)
,
f
(
2
,
0
)
)
,
b
1
=
p
(
t
x
,
f
(
−
1
,
1
)
,
f
(
0
,
1
)
,
f
(
1
,
1
)
,
f
(
2
,
1
)
)
,
b
2
=
p
(
t
x
,
f
(
−
1
,
2
)
,
f
(
0
,
2
)
,
f
(
1
,
2
)
,
f
(
2
,
2
)
)
,
p
(
x
,
y
)
=
p
(
t
y
,
b
−
1
,
b
0
,
b
1
,
b
2
)
.
{\displaystyle b_{-1}=p(t_{x},f_{(-1,-1)},f_{(0,-1)},f_{(1,-1)},f_{(2,-1)}),}\newline {\displaystyle b_{0}=p(t_{x},f_{(-1,0)},f_{(0,0)},f_{(1,0)},f_{(2,0)}),} \newline {\displaystyle b_{1}=p(t_{x},f_{(-1,1)},f_{(0,1)},f_{(1,1)},f_{(2,1)}),} \newline {\displaystyle b_{2}=p(t_{x},f_{(-1,2)},f_{(0,2)},f_{(1,2)},f_{(2,2)}),}\newline {\displaystyle p(x,y)=p(t_{y},b_{-1},b_{0},b_{1},b_{2}).}\newline
b−1=p(tx,f(−1,−1),f(0,−1),f(1,−1),f(2,−1)),b0=p(tx,f(−1,0),f(0,0),f(1,0),f(2,0)),b1=p(tx,f(−1,1),f(0,1),f(1,1),f(2,1)),b2=p(tx,f(−1,2),f(0,2),f(1,2),f(2,2)),p(x,y)=p(ty,b−1,b0,b1,b2).
算法实现
方法一：
static inline void matrix_mul(float c[4][4], float a[4][4], float b[4][4])
{
int n, k;
for (n = 0; n < 4; n++)
{
for (k = 0; k < 4; k++)
{
c[n][k] = a[n][0] * b[0][k] + a[n][1] * b[1][k] + a[n][2] * b[2][k] + a[n][3] * b[3][k];
}
}
}
void img_resize_using_bicubic(uint8_2d* src, uint8_2d* dst)
{
int src_rows, src_cols;
int dst_rows, dst_cols;
int i, j;
if (!src || !dst)return;
uint8_t** src_arr;
uint8_t** dst_arr;
float xratio;
float yratio;
float x, y;
int xi, yi;
int xminus1, yminus1;
int xadd1, yadd1;
int xadd2, yadd2;
float u, v;
float u2, u3;
float v2, v3;
int val;
src_arr = src->arr;
dst_arr = dst->arr;
src_rows = src->rows;
src_cols = src->cols;
dst_rows = dst->rows;
dst_cols = dst->cols;
float m1[4][4] = { { 1, 0, 0, 0 },
{ 0, 0, 1, 0},
{-3, 3, -2,-1},
{ 2,-2, 1, 1} };
float m0[4][4] = { { 1, 0, -3, 2 },
{ 0, 0, 3, -2 },
{ 0, 1, -2, 1 },
{ 0, 0, -1, 1 } };
float d[4][4];
float a[4][4];
xratio = (float)src_cols / (float)dst_cols;
yratio = (float)src_rows / (float)dst_rows;
int v00, v01, v02, v03;
int v10, v11, v12, v13;
int v20, v21, v22, v23;
int v30, v31, v32, v33;
float tmp[4];
for (i = 0; i < dst_rows; i++)
{
for (j = 0; j < dst_cols; j++)
{
x = xratio*j;
y = yratio*i;
xi = (int)x;
yi = (int)y;
u = x - xi;
v = y - yi;
xminus1 = xi - 1;
yminus1 = yi - 1;
xadd1 = xi + 1;
yadd1 = yi + 1;
xadd2 = xi + 2;
yadd2 = yi + 2;
if (xi >= src_cols)xi = src_cols - 1;
if (yi >= src_rows)yi = src_rows - 1;
if (xminus1 < 0)xminus1 = 0;
if (yminus1 < 0)yminus1 = 0;
if (xadd1 >= src_cols)xadd1 = src_cols - 1;
if (yadd1 >= src_rows)yadd1 = src_rows - 1;
if (xadd2 >= src_cols)xadd2 = src_cols - 1;
if (yadd2 >= src_rows)yadd2 = src_rows - 1;
v00 = src_arr[yminus1][xminus1];
v01 = src_arr[yminus1][xi];
v02 = src_arr[yminus1][xadd1];
v03 = src_arr[yminus1][xadd2];
v10 = src_arr[yi][xminus1];
v11 = src_arr[yi][xi];
v12 = src_arr[yi][xadd1];
v13 = src_arr[yi][xadd2];
v20 = src_arr[yadd1][xminus1];
v21 = src_arr[yadd1][xi];
v22 = src_arr[yadd1][xadd1];
v23 = src_arr[yadd1][xadd2];
v30 = src_arr[yadd2][xminus1];
v31 = src_arr[yadd2][xi];
v32 = src_arr[yadd2][xadd1];
v33 = src_arr[yadd2][xadd2];
d[0][0] = v11;//f00
d[0][1] = v21;//f01
d[1][0] = v12;//f10
d[1][1] = v22;//f11
//d[0][2] = (v21-v01)/2.0;//fy00
//d[0][3] = (v31-v11)/2.0;//fy01
//d[1][2] = (v22-v02)/2.0;//fy10;
//d[1][3] = (v32-v12)/2.0;//fy11;
//d[2][0] = (v12-v10)/2.0;//fx00;
//d[2][1] = (v22-v20)/2.0;//fx01;
//d[3][0] = (v13-v11)/2.0;//fx10;
//d[3][1] = (v23-v21)/2.0;//fx11;
//这里采用sobel算子求一阶偏导数，为了避免引入一些噪声
d[0][2] = (v20 + 2*v21+v22 - v00 - 2*v01- v02) / 8.0;//fy00
d[0][3] = (v30 + 2*v31+v32 - v10 - 2*v11- v12) / 8.0;//fy01
d[1][2] = (v21 + 2*v22+v23 - v01 - 2*v02- v03) / 8.0;//fy10;
d[1][3] = (v31 + 2*v32+v33 - v11 - 2*v12- v13) / 8.0;//fy11;
d[2][0] = (v02 + 2*v12+v22 - v00 - 2*v10 - v20) / 8.0;//fx00;
d[2][1] = (v12 + 2*v22+v32 - v10 - 2*v20 - v30) / 8.0;//fx01;
d[3][0] = (v03 + 2*v13+v23 - v01 - 2*v11 - v21) / 8.0;//fx10;
d[3][1] = (v13 + 2*v23+v33 - v11 - 2*v21 - v31) / 8.0;//fx11;
d[2][2] = (v22+v00-v02-v20) / 4.0; //fxy00
d[2][3] = (v32+v10-v30-v12) / 4.0; //fxy01;
d[3][2] = (v23+v01-v03-v21) / 4.0; //fxy10;
d[3][3] = (v33+v11-v13-v31) / 4.0; //fxy11
matrix_mul(a, d, m0);
matrix_mul(d, m1, a);
u2 = u*u;
u3 = u2*u;
v2 = v*v;
v3 = v2*v;
tmp[0] = d[0][0] + d[0][1] * v + d[0][2] * v2 + d[0][3] * v3;
tmp[1] = d[1][0] + d[1][1] * v + d[1][2] * v2 + d[1][3] * v3;
tmp[2] = d[2][0] + d[2][1] * v + d[2][2] * v2 + d[2][3] * v3;
tmp[3] = d[3][0] + d[3][1] * v + d[3][2] * v2 + d[3][3] * v3;
val = tmp[0] + tmp[1] * u + tmp[2] * u2 + tmp[3] * u3;
if (val < 0)val = 0;
if (val > 255)val = 255;
dst_arr[i][j] = val;
}
}
}
方法二：
static inline void matrix_mul1(float c[4], float a[4][4], float b[4])
{
int n;
for (n = 0; n < 4; n++)
{
c[n] = a[n][0] * b[0] + a[n][1] * b[1] + a[n][2] * b[2] + a[n][3] * b[3];
}
}
void img_resize_using_bicubic1(uint8_2d* src, uint8_2d* dst)
{
int src_rows, src_cols;
int dst_rows, dst_cols;
int i, j;
if (!src || !dst)return;
uint8_t** src_arr;
uint8_t** dst_arr;
float xratio;
float yratio;
float x, y;
int xi, yi;
int xminus1, yminus1;
int xadd1, yadd1;
int xadd2, yadd2;
float u, v;
float u2, u3;
float v2, v3;
int val;
src_arr = src->arr;
dst_arr = dst->arr;
src_rows = src->rows;
src_cols = src->cols;
dst_rows = dst->rows;
dst_cols = dst->cols;
float m[4][4] = { { 0, 2, 0, 0},
{ -1, 0, 1, 0},
{ 2, -5, 4, -1},
{ -1, 3, -3, 1}};
xratio = (float)src_cols / (float)dst_cols;
yratio = (float)src_rows / (float)dst_rows;
float tmp[4];
float pix[4];
float b[4];
for (i = 0; i < dst_rows; i++)
{
for (j = 0; j < dst_cols; j++)
{
x = xratio*j;
y = yratio*i;
xi = (int)x;
yi = (int)y;
u = x - xi;
v = y - yi;
xminus1 = xi - 1;
yminus1 = yi - 1;
xadd1 = xi + 1;
yadd1 = yi + 1;
xadd2 = xi + 2;
yadd2 = yi + 2;
if (xi >= src_cols)xi = src_cols - 1;
if (yi >= src_rows)yi = src_rows - 1;
if (xminus1 < 0)xminus1 = 0;
if (yminus1 < 0)yminus1 = 0;
if (xadd1 >= src_cols)xadd1 = src_cols - 1;
if (yadd1 >= src_rows)yadd1 = src_rows - 1;
if (xadd2 >= src_cols)xadd2 = src_cols - 1;
if (yadd2 >= src_rows)yadd2 = src_rows - 1;
u2 = u*u;
u3 = u2*u;
v2 = v*v;
v3 = v2*v;
pix[0] = src_arr[yminus1][xminus1];
pix[1] = src_arr[yminus1][xi];
pix[2] = src_arr[yminus1][xadd1];
pix[3] = src_arr[yminus1][xadd2];
matrix_mul1(tmp, m, pix);
b[0] = 0.5*(tmp[0] + tmp[1] * u + tmp[2] * u2 + tmp[3] * u3);
pix[0] = src_arr[yi][xminus1];
pix[1] = src_arr[yi][xi];
pix[2] = src_arr[yi][xadd1];
pix[3] = src_arr[yi][xadd2];
matrix_mul1(tmp, m, pix);
b[1] = 0.5*(tmp[0] + tmp[1] * u + tmp[2] * u2 + tmp[3] * u3);
pix[0] = src_arr[yadd1][xminus1];
pix[1] = src_arr[yadd1][xi];
pix[2] = src_arr[yadd1][xadd1];
pix[3] = src_arr[yadd1][xadd2];
matrix_mul1(tmp, m, pix);
b[2] = 0.5*(tmp[0] + tmp[1] * u + tmp[2] * u2 + tmp[3] * u3);
pix[0] = src_arr[yadd2][xminus1];
pix[1] = src_arr[yadd2][xi];
pix[2] = src_arr[yadd2][xadd1];
pix[3] = src_arr[yadd2][xadd2];
matrix_mul1(tmp, m, pix);
b[3] = 0.5*(tmp[0] + tmp[1] * u + tmp[2] * u2 + tmp[3] * u3);
matrix_mul1(tmp, m, b);
val = 0.5*(tmp[0] + tmp[1] * v + tmp[2] * v2 + tmp[3] * v3);
if (val < 0)val = 0;
if (val > 255)val = 255;
dst_arr[i][j] = val;
}
}
}
测试代码：
void test_bicubic_resize(dai_image* img, float factor)
{
uint8_2d* r = NULL;
uint8_2d* g = NULL;
uint8_2d* b = NULL;
uint8_2d* r1 = NULL;
uint8_2d* g1 = NULL;
uint8_2d* b1 = NULL;
dai_image* img1 = NULL;
if (!img)return;
split_img_data(img, &r, &g, &b);
int w, h;
int w1, h1;
w = img->width;
h = img->height;
w1 = factor*w;
h1 = factor*h;
r1 = create_uint8_2d(h1, w1);
g1 = create_uint8_2d(h1, w1);
b1 = create_uint8_2d(h1, w1);
if (factor < 1.0)
{
real_2d* kernel = create_gaussian_kernel(1.4);
real_2d* r2 = uint8_to_real(r->data, r->cols, r->rows);
gaussian_blur(r2, kernel, 0);
real_2d* g2 = uint8_to_real(g->data, g->cols, g->rows);
gaussian_blur(g2, kernel, 0);
real_2d* b2 = uint8_to_real(b->data, b->cols, b->rows);
gaussian_blur(b2, kernel, 0);
destroy_uint8_2d(&r);
destroy_uint8_2d(&g);
destroy_uint8_2d(&b);
r = real_to_uint8(r2);
g = real_to_uint8(g2);
b = real_to_uint8(b2);
}
img_resize_using_bicubic(r, r1);
img_resize_using_bicubic(g, g1);
img_resize_using_bicubic(b, b1);
merge_img_data(r1, g1, b1, &img1);
if (img1)
{
img1->type = EJPEG;
dai_save_image("resize_bicubic.jpg",img1);
dai_destroy_image(&img1);
}
destroy_uint8_2d(&r);
destroy_uint8_2d(&g);
destroy_uint8_2d(&b);
destroy_uint8_2d(&r1);
destroy_uint8_2d(&g1);
destroy_uint8_2d(&b1);
}
对于scale小于1.0的情况，需要先进行低通滤波然后再进行亚采样插值，split_img_data和merge_img_data可以参看本人上一篇关于二次线性插值的博客。
效果图：
原图：
因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- 7swz.com 版权所有 赣ICP备2024042798号-8
违法及侵权请联系：TEL:199 18 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务