原理：详细参见https://en.wikipedia.org/wiki/Bicubic_interpolation
在数学上，双立方插值是2个变量函数在方形网格上立方插值的一个扩展，在图像处理中，双立方插值考虑插值像素点周围16个像素，当不在乎算法执行速度时，一般选择双立方插值，相对最近邻域插值和二次线性插值，双立方插值具有较少的插值畸变（artifacts ，不知道翻译成什么好）。
双立方插值与最近邻域插值和二次线性插值的比较如下图：
假设在(0,0),(1,0),(0,1),(1,1)单位正方形四个顶点函数值f 及其偏导数
f
x
,
f
y
f
y
和
f
x
y
f
x
y
f_{x}, {\displaystyle f_{y}} f_y 和 {\displaystyle f_{xy}} f_{xy}
fx,fyfy和fxyfxy已知，那么，单位正方形内部任意点的插值可以表示为：
p
(
x
,
y
)
=
∑
i
=
0
3
∑
j
=
0
3
a
i
j
x
i
y
j
p(x,y) = \sum\limits_{i=0}^3 \sum_{j=0}^3 a_{ij} x^i y^j
p(x,y)=i=0∑3j=0∑3aijxiyj
因此插值问题变成求解16个参数的问题，四个点
(
0
,
0
)
,
(
1
,
0
)
,
(
0
,
1
)
,
和
(
1
,
1
)
(0,0), {\displaystyle (1,0)} , {\displaystyle (0,1)}, 和{\displaystyle } (1,1)
(0,0),(1,0),(0,1),和(1,1)代入上式可得：
f
(
0
,
0
)
=
p
(
0
,
0
)
=
a
00
,
f
(
1
,
0
)
=
p
(
1
,
0
)
=
a
00
+
a
10
+
a
20
+
a
30
,
f
(
0
,
1
)
=
p
(
0
,
1
)
=
a
00
+
a
01
+
a
02
+
a
03
,
f
(
1
,
1
)
=
p
(
1
,
1
)
=
∑
i
=
0
3
∑
j
=
0
3
a
i
j
.
{\displaystyle f(0,0)=p(0,0)=a_{00},}\newline {\displaystyle f(1,0)=p(1,0)=a_{00}+a_{10}+a_{20}+a_{30},} \newline {\displaystyle f(0,1)=p(0,1)=a_{00}+a_{01}+a_{02}+a_{03},}\newline {\displaystyle f(1,1)=p(1,1)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=0}^{3}a_{ij}.}
f(0,0)=p(0,0)=a00,f(1,0)=p(1,0)=a00+a10+a20+a30,f(0,1)=p(0,1)=a00+a01+a02+a03,f(1,1)=p(1,1)=i=0∑3j=0∑3aij.
根据4个点x和y方向上的偏导数，可得8个方程：
f
x
(
0
,
0
)
=
p
x
(
0
,
0
)
=
a
10
,
f
x
(
1
,
0
)
=
p
x
(
1
,
0
)
=
a
10
+
2
a
20
+
3
a
30
,
f
x
(
0
,
1
)
=
p
x
(
0
,
1
)
=
a
10
+
a
11
+
a
12
+
a
13
,
f
x
(
1
,
1
)
=
p
x
(
1
,
1
)
=
∑
i
=
1
3
∑
j
=
0
3
a
i
j
i
,
f
y
(
0
,
0
)
=
p
y
(
0
,
0
)
=
a
01
,
f
y
(
1
,
0
)
=
p
y
(
1
,
0
)
=
a
01
+
a
11
+
a
21
+
a
31
,
f
y
(
0
,
1
)
=
p
y
(
0
,
1
)
=
a
01
+
2
a
02
+
3
a
03
,
f
y
(
1
,
1
)
=
p
y
(
1
,
1
)
=
∑
i
=
0
3
∑
j
=
1
3
a
i
j
j
.
{\displaystyle f_{x}(0,0)=p_{x}(0,0)=a_{10},}\newline {\displaystyle f_{x}(1,0)=p_{x}(1,0)=a_{10}+2a_{20}+3a_{30},}\newline {\displaystyle f_{x}(0,1)=p_{x}(0,1)=a_{10}+a_{11}+a_{12}+a_{13},}\newline {\displaystyle f_{x}(1,1)=p_{x}(1,1)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}i,}\newline {\displaystyle f_{y}(0,0)=p_{y}(0,0)=a_{01},}\newline {\displaystyle f_{y}(1,0)=p_{y}(1,0)=a_{01}+a_{11}+a_{21}+a_{31},} \newline{\displaystyle f_{y}(0,1)=p_{y}(0,1)=a_{01}+2a_{02}+3a_{03},}\newline {\displaystyle f_{y}(1,1)=p_{y}(1,1)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}a_{ij}j.}\newline
fx(0,0)=px(0,0)=a10,fx(1,0)=px(1,0)=a10+2a20+3a30,fx(0,1)=px(0,1)=a10+a11+a12+a13,fx(1,1)=px(1,1)=i=1∑3j=0∑3aiji,fy(0,0)=py(0,0)=a01,fy(1,0)=py(1,0)=a01+a11+a21+a31,fy(0,1)=py(0,1)=a01+2a02+3a03,fy(1,1)=py(1,1)=i=0∑3j=1∑3aijj.
四个点上的二阶混合偏导数（mixed partial derivative）可得四个方程：
f
x
y
(
0
,
0
)
=
p
x
y
(
0
,
0
)
=
a
11
,
f
x
y
(
1
,
0
)
=
p
x
y
(
1
,
0
)
=
a
11
+
2
a
21
+
3
a
31
,
f
x
y
(
0
,
1
)
=
p
x
y
(
0
,
1
)
=
a
11
+
2
a
12
+
3
a
13
,
f
x
y
(
1
,
1
)
=
p
x
y
(
1
,
1
)
=
∑
i
=
1
3
∑
j
=
1
3
a
i
j
i
j
.
{\displaystyle f_{xy}(0,0)=p_{xy}(0,0)=a_{11},}\newline {\displaystyle f_{xy}(1,0)=p_{xy}(1,0)=a_{11}+2a_{21}+3a_{31},} \newline {\displaystyle f_{xy}(0,1)=p_{xy}(0,1)=a_{11}+2a_{12}+3a_{13},}\newline {\displaystyle f_{xy}(1,1)=p_{xy}(1,1)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=1}^{3}a_{ij}ij.}
fxy(0,0)=pxy(0,0)=a11,fxy(1,0)=pxy(1,0)=a11+2a21+3a31,fxy(0,1)=pxy(0,1)=a11+2a12+3a13,fxy(1,1)=pxy(1,1)=i=1∑3j=1∑3aijij.
上面的四个混合偏导数还可以表示为：
p
x
(
x
,
y
)
=
∑
i
=
1
3
∑
j
=
0
3
a
i
j
i
x
i
−
1
y
j
,
p
y
(
x
,
y
)
=
∑
i
=
0
3
∑
j
=
1
3
a
i
j
x
i
j
y
j
−
1
,
p
x
y
(
x
,
y
)
=
∑
i
=
1
3
∑
j
=
1
3
a
i
j
i
x
i
−
1
j
y
j
−
1
.
{\displaystyle p_{x}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}ix^{i-1}y^{j},}\newline {\displaystyle p_{y}(x,y)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}a_{ij}x^{i}jy^{j-1},}\newline {\displaystyle p_{xy}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=1}^{3}a_{ij}ix^{i-1}jy^{j-1}.}
px(x,y)=i=1∑3j=0∑3aijixi−1yj,py(x,y)=i=0∑3j=1∑3aijxijyj−1,pxy(x,y)=i=1∑3j=1∑3aijixi−1jyj−1.
将
a
i
j
a_{ij}
aij组合成一个向量：
α
=
[
a
00
a
10
a
20
a
30
a
01
a
11
a
21
a
31
a
02
a
12
a
22
a
32
a
03
a
13
a
23
a
33
]
T
\alpha=\left[\begin{matrix}a_{00}&a_{10}&a_{20}&a_{30}&a_{01}&a_{11}&a_{21}&a_{31}&a_{02}&a_{12}&a_{22}&a_{32}&a_{03}&a_{13}&a_{23}&a_{33}\end{matrix}\right]^T
α=[a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a33]T令：
x
=
[
f
(
0
,
0
)
f
(
1
,
0
)
f
(
0
,
1
)
f
(
1
,
1
)
f
x
(
0
,
0
)
f
x
(
1
,
0
)
f
x
(
0
,
1
)
f
x
(
1
,
1
)
f
y
(
0
,
0
)
f
y
(
1
,
0
)
f
y
(
0
,
1
)
f
y
(
1
,
1
)
f
x
y
(
0
,
0
)
f
x
y
(
1
,
0
)
f
x
y
(
0
,
1
)
f
x
y
(
1
,
1
)
]
T
,
{\displaystyle x=\left[{\begin{matrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&f_{x}(0,0)&f_{x}(1,0)&f_{x}(0,1)&f_{x}(1,1)&f_{y}(0,0)&f_{y}(1,0)&f_{y}(0,1)&f_{y}(1,1)&f_{xy}(0,0)&f_{xy}(1,0)&f_{xy}(0,1)&f_{xy}(1,1)\end{matrix}}\right]^{T},}
x=[f(0,0)f(1,0)f(0,1)f(1,1)fx(0,0)fx(1,0)fx(0,1)fx(1,1)fy(0,0)fy(1,0)fy(0,1)fy(1,1)fxy(0,0)fxy(1,0)fxy(0,1)fxy(1,1)]T,
上面的16个式子可以写成
A
α
=
x
A\alpha=x
Aα=x,求解这个线性方程组可得
a
i
j
a_{ij}
aij。
A
A
A的逆矩阵如下：
A
−
1
=
[
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
−
3
3
0
0
−
2
−
1
0
0
0
0
0
0
0
0
0
0
2
−
2
0
0
1
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
−
3
3
0
0
−
2
−
1
0
0
0
0
0
0
0
0
0
0
2
−
2
0
0
1
1
0
0
−
3
0
3
0
0
0
0
0
−
2
0
−
1
0
0
0
0
0
0
0
0
0
−
3
0
3
0
0
0
0
0
−
2
0
−
1
0
9
−
9
−
9
9
6
3
−
6
−
3
6
−
6
3
−
3
4
2
2
1
−
6
6
6
−
6
−
3
−
3
3
3
−
4
4
−
2
2
−
2
−
2
−
1
−
1
2
0
−
2
0
0
0
0
0
1
0
1
0
0
0
0
0
0
0
0
0
2
0
−
2
0
0
0
0
0
1
0
1
0
−
6
6
6
−
6
−
4
−
2
4
2
−
3
3
−
3
3
−
2
−
1
−
2
−
1
4
−
4
−
4
4
2
2
−
2
−
2
2
−
2
2
−
2
1
1
1
1
]
{\displaystyle A^{-1}=\left[{\begin{matrix}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\-3&3&0&0&-2&-1&0&0&0&0&0&0&0&0&0&0\\2&-2&0&0&1&1&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&1&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0&0&0&1&0&0&0\\0&0&0&0&0&0&0&0&-3&3&0&0&-2&-1&0&0\\0&0&0&0&0&0&0&0&2&-2&0&0&1&1&0&0\\-3&0&3&0&0&0&0&0&-2&0&-1&0&0&0&0&0\\0&0&0&0&-3&0&3&0&0&0&0&0&-2&0&-1&0\\9&-9&-9&9&6&3&-6&-3&6&-6&3&-3&4&2&2&1\\-6&6&6&-6&-3&-3&3&3&-4&4&-2&2&-2&-2&-1&-1\\2&0&-2&0&0&0&0&0&1&0&1&0&0&0&0&0\\0&0&0&0&2&0&-2&0&0&0&0&0&1&0&1&0\\-6&6&6&-6&-4&-2&4&2&-3&3&-3&3&-2&-1&-2&-1\\4&-4&-4&4&2&2&-2&-2&2&-2&2&-2&1&1&1&1\end{matrix}}\right]}
A−1=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡10−320000−309−620−64003−2000000−96006−40000000030−96−206−400000000009−600−6401−2100000−36−302−4200−110000003−300−220000000003−630−24−20000000000−33002−2000010−32−206−410−320000003−200−64003−200000000−103−210−320000000000−32003−2000001−210−24−201−21000000−11002−200−11000000000−12−101−2100000000001−100−11⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
另外，上面的16个方程的简明表示：
[
f
(
0
,
0
)
f
(
0
,
1
)
f
y
(
0
,
0
)
f
y
(
0
,
1
)
f
(
1
,
0
)
f
(
1
,
1
)
f
y
(
1
,
0
)
f
y
(
1
,
1
)
f
x
(
0
,
0
)
f
x
(
0
,
1
)
f
x
y
(
0
,
0
)
f
x
y
(
0
,
1
)
f
x
(
1
,
0
)
f
x
(
1
,
1
)
f
x
y
(
1
,
0
)
f
x
y
(
1
,
1
)
]
=
[
1
0
0
0
1
1
1
1
0
1
0
0
0
1
2
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
1
0
0
0
1
1
1
0
1
0
2
0
1
0
3
]
,
{\displaystyle {\begin{bmatrix}f(0,0)&f(0,1)&f_{y}(0,0)&f_{y}(0,1)\\f(1,0)&f(1,1)&f_{y}(1,0)&f_{y}(1,1)\\f_{x}(0,0)&f_{x}(0,1)&f_{xy}(0,0)&f_{xy}(0,1)\\f_{x}(1,0)&f_{x}(1,1)&f_{xy}(1,0)&f_{xy}(1,1)\end{bmatrix}}={\begin{bmatrix}1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&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&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3\end{bmatrix}},}
⎣⎢⎢⎡f(0,0)f(1,0)fx(0,0)fx(1,0)f(0,1)f(1,1)fx(0,1)fx(1,1)fy(0,0)fy(1,0)fxy(0,0)fxy(1,0)fy(0,1)fy(1,1)fxy(0,1)fxy(1,1)⎦⎥⎥⎤=⎣⎢⎢⎡1100011101020103⎦⎥⎥⎤⎣⎢⎢⎡a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a33⎦⎥⎥⎤⎣⎢⎢⎡1000111101000123⎦⎥⎥⎤,
等式右边两边左乘
[
1
0
0
0
1
1
1
1
0
1
0
0
0
1
2
3
]
{\begin{bmatrix}1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&3\end{bmatrix}}
⎣⎢⎢⎡1100011101020103⎦⎥⎥⎤的逆矩阵，然后右乘
[
1
1
0
0
0
1
1
1
0
1
0
2
0
1
0
3
]
{\begin{bmatrix}1&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3\end{bmatrix}}
⎣⎢⎢⎡1000111101000123⎦⎥⎥⎤的逆矩阵，可得：
[
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
0
0
0
0
0
1
0
−
3
3
−
2
−
1
2
−
2
1
1
]
[
f
(
0
,
0
)
f
(
0
,
1
)
f
y
(
0
,
0
)
f
y
(
0
,
1
)
f
(
1
,
0
)
f
(
1
,
1
)
f
y
(
1
,
0
)
f
y
(
1
,
1
)
f
x
(
0
,
0
)
f
x
(
0
,
1
)
f
x
y
(
0
,
0
)
f
x
y
(
0
,
1
)
f
x
(
1
,
0
)
f
x
(
1
,
1
)
f
x
y
(
1
,
0
)
f
x
y
(
1
,
1
)
]
[
1
0
−
3
2
0
0
3
−
2
0
1
−
2
1
0
0
−
1
1
]
,
{\displaystyle {\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&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end{bmatrix}}{\begin{bmatrix}f(0,0)&f(0,1)&f_{y}(0,0)&f_{y}(0,1)\\f(1,0)&f(1,1)&f_{y}(1,0)&f_{y}(1,1)\\f_{x}(0,0)&f_{x}(0,1)&f_{xy}(0,0)&f_{xy}(0,1)\\f_{x}(1,0)&f_{x}(1,1)&f_{xy}(1,0)&f_{xy}(1,1)\end{bmatrix}}{\begin{bmatrix}1&0&-3&2\\0&0&3&-2\\0&1&-2&1\\0&0&-1&1\end{bmatrix}},}
⎣⎢⎢⎡a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a33⎦⎥⎥⎤=⎣⎢⎢⎡10−32003−201−2100−11⎦⎥⎥⎤⎣⎢⎢⎡f(0,0)f(1,0)fx(0,0)fx(1,0)f(0,1)f(1,1)fx(0,1)fx(1,1)fy(0,0)fy(1,0)fxy(0,0)fxy(1,0)fy(0,1)fy(1,1)fxy(0,1)fxy(1,1)⎦⎥⎥⎤⎣⎢⎢⎡10000010−33−2−12−211⎦⎥⎥⎤,
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⎦⎥⎥⎤.
差分计算可以采用下面公式，参考：https://en.wikipedia.org/wiki/Finite_difference
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可以参看本人上一篇关于二次线性插值的博客。
效果图：
原图：
放大图：
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
