双三次插值（卷积插值）bicubic interplotation

给几个传送门：http://www.tinaja.com/glib/pixintpl.pdf

http://en.wikipedia.org/wiki/Bicubic_interpolation

http://luoweifu.iteye.com/blog/1726068

http://www.codeproject.com/Articles/236394/Bi-Cubic-and-Bi-Linear-Interpolation-with-GLSL

[        16*16的系数矩阵T        ]*[a00,a10,a20,a30,a01,a11,a21,a31,a02,a12,a22,a32,a03,a13,a23,a33]^T=[f00,f10,f01,f11,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*16的系数矩阵F ]*[p00,p01,p02,p03,p10,p11,p12,p13,p20,p21,p22,p23,p30,p31,p32,p33]^T

p00到p33是附近的16个点的像素值（原图 Origin picture）

T*A=F*P

求A（a00,a10,a20……）只要求T^-1*F *P，求T^-1*F 即可以求出A和P的关系 （a00到a33与16个像素值的关系）

再回代入

可以求出v(x,y) x为[0,1]的浮点小数,y同理

T矩阵＝

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 2 3 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0
0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 2 0 0 0 3 0 0 0
0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 2 3 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 2 0 0 0 3 0 0
0 0 0 0 0 1 2 3 0 2 4 6 0 3 6 9

a00
a10
a20
a30
a01
a11
a21
a31
这样的顺序

T^-1为

F矩阵为

 

计算T^-1*F的代码如下

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;
double a[100][100],b[100][100],c[100][100];
int da[16],db[16],dc[16],dd[16],bel[16][16];
int main()
{
freopen("1.in","r",stdin);
freopen("11.out","w",stdout);
for (int i=0;i<16;i++)
for (int j=0;j<16;j++)
scanf("%lf",&a[i][j]);
for (int i=0;i<16;i++)
for (int j=0;j<16;j++)
scanf("%lf",&b[i][j]);
for (int i=0;i<16;i++)
for (int j=0;j<16;j++)
for (int k=0;k<16;k++)
{
double tmp=a[i][k]*b[k][j];
c[i][j]=c[i][j]+tmp;
}
printf("%.2lf %.2lf\n",a[0][0],b[0][5]);

printf("\n");
for (int i=0;i<16;i++)
{
for (int j=0;j<16;j++)
printf("%.2lf ",b[i][j]);
printf("\n");
}
printf("above was b:\n");
for (int i=0;i<16;i++)
{
for (int j=0;j<16;j++)
printf("%.2lf ",c[i][j]);
printf("\n");
}
int cnt=0;
for (int i=0;i<4;i++)
{
for (int j=0;j<4;j++)
{
da[cnt]=j;
db[cnt]=i;
bel[j][i]=cnt;
cnt++;
}
}
cnt=0;
for (int i=0;i<4;i++)
{
for (int j=0;j<4;j++)
{
dc[cnt]=i;
dd[cnt]=j;
cnt++;
}
}
for (int i=0;i<16;i++) printf("%d %d\n",da[i],db[i]);
for (int i=0;i<16;i++)
{
printf("a %d%d = ",dc[i],dd[i]);
for (int j=0;j<16;j++)
{
if (c[bel[dc[i]][dd[i]]][j])
{
printf("%.2lfp[%d][%d]+",c[bel[dc[i]][dd[i]]][j],dc[j],dd[j]);
}
}
printf("\n");
}
printf("%.2lf %.2lf %.2lf\n",c[5][0],a[5][12],b[12][0]);
return 0;
}

系数为

a 00 = 1.00*p[1][1]
a 01 = -0.50*p[1][0]+0.50*p[1][2]
a 02 = 1.00*p[1][0]-2.50*p[1][1]+2.00*p[1][2]-0.50*p[1][3]
a 03 = -0.50*p[1][0]+1.50*p[1][1]-1.50*p[1][2]+0.50*p[1][3]
a 10 = -0.50*p[0][1]+0.50*p[2][1]
a 11 = 0.25*p[0][0]-0.25*p[0][2]-0.25*p[2][0]+0.25*p[2][2]
a 12 = -0.50*p[0][0]+1.25*p[0][1]-1.00*p[0][2]+0.25*p[0][3]+0.50*p[2][0]-1.25*p[2][1]+1.00*p[2][2]-0.25*p[2][3]
a 13 = 0.25*p[0][0]-0.75*p[0][1]+0.75*p[0][2]-0.25*p[0][3]-0.25*p[2][0]+0.75*p[2][1]-0.75*p[2][2]+0.25*p[2][3]
a 20 = 1.00*p[0][1]-2.50*p[1][1]+2.00*p[2][1]-0.50*p[3][1]
a 21 = -0.50*p[0][0]+0.50*p[0][2]+1.25*p[1][0]-1.25*p[1][2]-1.00*p[2][0]+1.00*p[2][2]+0.25*p[3][0]-0.25*p[3][2]
a 22 = 1.00*p[0][0]-2.50*p[0][1]+2.00*p[0][2]-0.50*p[0][3]-2.50*p[1][0]+6.25*p[1][1]-5.00*p[1][2]+1.25*p[1][3]+2.00*p[2][0]-5.00*p[2][1]+4.00*p[2][2]-1.00*p[2][3]-0.50*p[3][0]+1.25*p[3][1]-1.00*p[3][2]+0.25*p[3][3]
a 23 = -0.50*p[0][0]+1.50*p[0][1]-1.50*p[0][2]+0.50*p[0][3]+1.25*p[1][0]-3.75*p[1][1]+3.75*p[1][2]-1.25*p[1][3]-1.00*p[2][0]+3.00*p[2][1]-3.00*p[2][2]+1.00*p[2][3]+0.25*p[3][0]-0.75*p[3][1]+0.75*p[3][2]-0.25*p[3][3]
a 30 = -0.50*p[0][1]+1.50*p[1][1]-1.50*p[2][1]+0.50*p[3][1]
a 31 = 0.25*p[0][0]-0.25*p[0][2]-0.75*p[1][0]+0.75*p[1][2]+0.75*p[2][0]-0.75*p[2][2]-0.25*p[3][0]+0.25*p[3][2]
a 32 = -0.50*p[0][0]+1.25*p[0][1]-1.00*p[0][2]+0.25*p[0][3]+1.50*p[1][0]-3.75*p[1][1]+3.00*p[1][2]-0.75*p[1][3]-1.50*p[2][0]+3.75*p[2][1]-3.00*p[2][2]+0.75*p[2][3]+0.50*p[3][0]-1.25*p[3][1]+1.00*p[3][2]-0.25*p[3][3]
a 33 = 0.25*p[0][0]-0.75*p[0][1]+0.75*p[0][2]-0.25*p[0][3]-0.75*p[1][0]+2.25*p[1][1]-2.25*p[1][2]+0.75*p[1][3]+0.75*p[2][0]-2.25*p[2][1]+2.25*p[2][2]-0.75*p[2][3]-0.25*p[3][0]+0.75*p[3][1]-0.75*p[3][2]+0.25*p[3][3]