图像双立方插值——C实现

原理:详细参见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,fyfyfxyfxy已知,那么,单位正方形内部任意点的插值可以表示为:
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=03j=03aijxiyj
因此插值问题变成求解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=03j=03aij.
根据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=13j=03aiji,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=03j=13aijj.
四个点上的二阶混合偏导数(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=13j=13aijij.
上面的四个混合偏导数还可以表示为:
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=13j=03aijixi1yj,py(x,y)=i=03j=13aijxijyj1,pxy(x,y)=i=13j=13aijixi1jyj1.
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]} A1=1032000030962064003200000096006400000000309620640000000000960064012100000363024200110000003300220000000003630242000000000033002200001032206410320000003200640032000000001032103200000000003200320000012102420121000000110022001100000000012101210000000000110011
另外,上面的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)=1100011101020103a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a331000111101000123,
等式右边两边左乘 [ 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=1032003201210011f(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)1000001033212211,

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]a00a10a20a30a01a11a21a31a02a12a22a32a03a13a23a331yy2y3.

差分计算可以采用下面公式,参考: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(xh,y)2kf(x,y+k)f(x,yk)h2f(x+h,y)2f(x,y)+f(xh,y)k2f(x,y+k)2f(x,y)+f(x,yk)4hkf(x+h,y+k)f(x+h,yk)f(xh,y+k)+f(xh,yk)
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 &lt; ∣ x ∣ &lt; 2 , 0 otherwise , {\displaystyle W(x)={\begin{cases}(a+2)|x|^{3}-(a+3)|x|^{2}+1&amp;{\text{for }}|x|\leq 1,\\a|x|^{3}-5a|x|^{2}+8a|x|-4a&amp;{\text{for }}1&lt;|x|&lt;2,\\0&amp;{\text{otherwise}},\end{cases}}} W(x)=(a+2)x3(a+3)x2+1ax35ax2+8ax4a0for x1,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&amp;t&amp;t^{2}&amp;t^{3}\\\end{bmatrix}}{\begin{bmatrix}0&amp;2&amp;0&amp;0\\-1&amp;0&amp;1&amp;0\\2&amp;-5&amp;4&amp;-1\\-1&amp;3&amp;-3&amp;1\\\end{bmatrix}}{\begin{bmatrix}f_{-1}\\f_{0}\\f_{1}\\f_{2}\\\end{bmatrix}}} p(t)=21[1tt2t3]0121205301430011f1f0f1f2
对图像进行插值时,可以先进行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 b1=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,b1,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可以参看本人上一篇关于二次线性插值的博客。

效果图:
原图:
在这里插入图片描述
放大图:
在这里插入图片描述

立方插值(Bilinear Interpolation)是一种在二维空间中对图像或栅格数据进行平滑近似的方法。在Python中,你可以使用NumPy库来方便地实现这个算法。以下是一个简单的Python实现: ```python import numpy as np def bilinear_interpolation(x, y, grid, values): """ 立方插值函数 :param x: 横坐标 :param y: 纵坐标 :param grid: 细分网格的行和列索引 :param values: 格子中的值数组 :return: 插值后的值 """ # 计算坐标对应的四个邻居的索引 x0, y0 = np.floor((grid[:, 0] - x)).astype(int), np.floor((grid[0, :] - y)).astype(int) x1, y1 = x0 + 1, y0 + 1 # 如果坐标在边界之外,处理边界情况 x0[x0 < 0] = 0 y0[y0 < 0] = 0 x1[x1 >= grid.shape] = grid.shape - 1 y1[y1 >= grid.shape] = grid.shape - 1 # 计算权重(每个邻居的权重是其距离目标位置的平方和的倒数) weights_x0 = (1 - (x - x0)) ** 2 weights_x1 = (1 - (x1 - x)) ** 2 weights_y0 = (1 - (y - y0)) ** 2 weights_y1 = (1 - (y1 - y)) ** 2 # 计算插值值 interpolated_value = np.sum(weights_x0 * weights_y0 * values[y0, x0]) + \ np.sum(weights_x0 * weights_y1 * values[y1, x0]) + \ np.sum(weights_x1 * weights_y0 * values[y0, x1]) + \ np.sum(weights_x1 * weights_y1 * values[y1, x1]) return interpolated_value # 示例用法 grid = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) values = np.random.rand(3, 3) # 假设这里有一张3x3的随机数值网格 x, y = 1.5, 2.5 # 需要插值的坐标 result = bilinear_interpolation(x, y, grid, values) ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值