下面的这个算法自己看的不是很懂，是图像配准算法中进行仿射变换之后的插值，不知道这个插值算法到底是怎么实现的，知道的人能帮忙解释解释吗？谢谢啦
* \说明:
* 该函数根据邻近位置的数值进行插值。
*此图象的大小为rectNewImg
* unsigned char CImageMatch::CalSpline(unsigned char *pUnchCorr, double x, double y)
{ double ret=0, Cx, Cy;
double Temp; for(int i=0;i <4;i++)
for(int j=0;j <4;j++)
{
Temp = pUnchCorr[i*4 + j];
if(fabs(y-i) <1)
Cy = 1-2*fabs(y-i)*fabs(y-i)+fabs(y-i)*fabs(y-i)*fabs(y-i);
if(fabs(y-i)>=1)
Cy = 4-8*fabs(y-i)+5*fabs(y-i)*fabs(y-i)-fabs(y-i)*fabs(y-i)*fabs(y-i);
if(fabs(x-j) <1)
Cx = 1-2*fabs(x-j)*fabs(x-j)+fabs(x-j)*fabs(x-j)*fabs(x-j);
if(fabs(x-j)>=1)
Cx = 4-8*fabs(x-j)+5*fabs(x-j)*fabs(x-j)-fabs(x-j)*fabs(x-j)*fabs(x-j);
ret += Temp*Cy*Cx;
}
if(ret <0)
ret=0;
if(ret>255)
ret=255; return (unsigned char)ret; }
* \说明:
* 该函数根据邻近位置的数值进行插值。
*此图象的大小为rectNewImg
* unsigned char CImageMatch::CalSpline(unsigned char *pUnchCorr, double x, double y)
{ double ret=0, Cx, Cy;
double Temp; for(int i=0;i <4;i++)
for(int j=0;j <4;j++)
{
Temp = pUnchCorr[i*4 + j];
if(fabs(y-i) <1)
Cy = 1-2*fabs(y-i)*fabs(y-i)+fabs(y-i)*fabs(y-i)*fabs(y-i);
if(fabs(y-i)>=1)
Cy = 4-8*fabs(y-i)+5*fabs(y-i)*fabs(y-i)-fabs(y-i)*fabs(y-i)*fabs(y-i);
if(fabs(x-j) <1)
Cx = 1-2*fabs(x-j)*fabs(x-j)+fabs(x-j)*fabs(x-j)*fabs(x-j);
if(fabs(x-j)>=1)
Cx = 4-8*fabs(x-j)+5*fabs(x-j)*fabs(x-j)-fabs(x-j)*fabs(x-j)*fabs(x-j);
ret += Temp*Cy*Cx;
}
if(ret <0)
ret=0;
if(ret>255)
ret=255; return (unsigned char)ret; }
pUnchCorr是x,y周围的16个点的颜色值矩阵。
这个函数就是要综合这16个点，插值求出x,y点的颜色值，这个颜色值实际上就是个加权平均。
循环里的Temp就是i,j点的颜色值，cx, cy就是这个点的权重。
fabs(y-i)和fabs(x-j)就是(x,y)到每个(i,j)的距离，权重就是由这个距离决定的。
权重的公式，反正前辈已经算好了。返回值ret就是所有16个点的加权平均，也就是我们要的x,y点的颜色
确定变换后的坐标是由下面的式子求出来的，可总感觉如果这样计算的话，tx,ty的坐标都在（1，1）附近，这又是为什么？// 确定变换的坐标
ty = ty - (int)ty + 1;
tx = tx - (int)tx + 1;
附部分程序：
for(int i=0;i<rectNewImg.bottom-rectNewImg.top;i++)
for(int j=0;j<rectNewImg.right-rectNewImg.left;j++)
{
tx = pDbAffPara[0*3 + 0]*(j+rectNewImg.left) +
pDbAffPara[0*3 + 1]*(i+rectNewImg.top) + pDbAffPara[0*3 + 2];
ty = pDbAffPara[1*3 + 0]*(j+rectNewImg.left) +
pDbAffPara[1*3 + 1]*(i+rectNewImg.top) + pDbAffPara[1*3 + 2]; for(int m=(int)ty-1;m<=(int)ty+2;m++)
for(int n=(int)tx-1;n<=(int)tx+2;n++)
{
if(m<0||m>=nSamplImgHeight||n<0||n>=nSamplImgWidth)
pUnchSect[(m-(int)ty+1)*4 + (n-(int)tx+1)] = 0;
else
pUnchSect[(m-(int)ty+1)*4 + (n-(int)tx+1)] =
m_pDibSamp->m_lpImage[(nSamplImgHeight-m-1)*nSampSaveWidth + n];
} // 确定变换的坐标
ty = ty - (int)ty + 1;
tx = tx - (int)tx + 1; // 确定变换后此坐标的数值
pUnchAftAffSamp[i*nNewImgWidth + j] = CalSpline(pUnchSect,tx,ty);//pUnchSect插值的点
tx = tx - (int)tx + 1;
这个变换是把任意一个点按整数平移到，(1,1)和(2,2)之间。这样如果找最近的16个点，肯定是以这个格子为中心的一个3*3的九宫格，坐标从(0,0)到(3,3)。这样正好和16个点的矩阵坐标相匹配。这样在CalSpline函数里坐标数据就不用再平移了，节省运算量吧。
typedef enum _condition
{//边界条件
NATURAL,//自然边界
CLAMPED,//夹持边界
NOTaKNOT//非扭结边界
}condition;typedef struct _coefficient
{//对应每个区间段的基函数对应的系数
double a3;
double b2;
double c1;
double d0;
}coefficient;typedef struct _splinePoint
{
coefficient *coe;//存储两点间多项式系数
double *xCoordinate;//被插值值对应(x,y)坐标
double *yCoordinate;
int num;//对应插值坐标点数 double f0;//对应区间的边界条件
double fn; condition con;//何种边界条件
}splinePoint;
上面对应变量声明！下面为插值函数
已经试用过，效果很好！
int Spline(splinePoint *point)
{//三次样条插值
double *x = (*point).xCoordinate;
double *y = (*point).yCoordinate;
int n = (*point).num - 1;
coefficient *coe = (*point).coe;
condition con = (*point).con;
double *h, *d;
double *a, *b, *c, *f, *m;
double temp;
int i;
h = (double *)malloc( n * sizeof(double) ); /* 0,1--(n-1),n */
d = (double *)malloc( n * sizeof(double) ); /* 0,1--(n-1),n */
a = (double *)malloc( (n + 1) * sizeof(double) ); /* 特别使用,1--(n-1), n */
b = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1), n */
c = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1),特别使用 */
f = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1), n */
m = b;
if(f == NULL)
{
free(h);
free(d);
free(a);
free(b);
free(c);
return 1;
}
/* 计算 h[] d[] */
for (i = 0; i < n; i++)
{
h[i] = x[i + 1] - x[i]; //存储每段之间的h值
d[i] = (y[i + 1] - y[i]) / h[i];//d[i]表示点的导数，除去最后一个点
} /**********************
** 初始化系数增广矩阵
**********************/
a[0] = (*point).f0;//置零
c[n] = (*point).fn;
/* 计算 a[] b[] d[] f[] */
switch(con)
{
case NATURAL://自然边界
f[0] = a[0];//边界为零
f[n] = c[n]; a[0] = 0;
c[n] = 0; c[0] = 0;
a[n] = 0; b[0] = 1;
b[n] = 1;
break;
case CLAMPED://夹持边界
f[0] = 6 * (d[0] - a[0]);
f[n] = 6 * (c[n] - d[n - 1]); a[0] = 0;
c[n] = 0; c[0] = h[0];
a[n] = h[n - 1]; b[0] = 2 * h[0];
b[n] = 2 * h[n - 1];
break;
case NOTaKNOT://非扭结边界
f[0] = 0;
f[n] = 0; a[0] = h[0];
c[n] = h[n - 1]; c[0] = -(h[0] + h[1]);
a[n] = -(h[n - 2] + h[n - 1]); b[0] = h[1];
b[n] = h[n - 2];
break;
} for (i = 1; i < n; i++)
{
a[i] = h[i - 1];//存放x段的距离
b[i] = 2 * (h[i - 1] + h[i]);
c[i] = h[i];
f[i] = 6 * (d[i] - d[i - 1]);
}
/* for (i = 0; i < n+1; i++) printf("%f\n", c[i]); */
/*************
** 高斯消元
*************/
/* 第0行到第(n-3)行的消元 */
for(i = 0; i <= n - 3; i++)
{
/* 选主元 */
if( fabs(a[i + 1]) > fabs(b[i]) )
{
temp = a[i + 1]; a[i + 1] = b[i]; b[i] = temp;//a[i + 1]与b[i]对换
temp = b[i + 1]; b[i + 1] = c[i]; c[i] = temp;//b[i + 1]与c[i]对换
temp = c[i + 1]; c[i + 1] = a[i]; a[i] = temp;//c[i + 1]与a[i]对换
temp = f[i + 1]; f[i + 1] = f[i]; f[i] = temp;//f[i + 1]与f[i]对换
}
temp = a[i + 1] / b[i];
a[i + 1] = 0;
b[i + 1] = b[i + 1] - temp * c[i];
c[i + 1] = c[i + 1] - temp * a[i];
f[i + 1] = f[i + 1] - temp * f[i];
}
/* 最后3行的消元 */
/* 选主元 */
if( fabs(a[n - 1]) > fabs(b[n - 2]) )
{
temp = a[n - 1]; a[n - 1] = b[n - 2]; b[n - 2] = temp;
temp = b[n - 1]; b[n - 1] = c[n - 2]; c[n - 2] = temp;
temp = c[n - 1]; c[n - 1] = a[n - 2]; a[n - 2] = temp;
temp = f[n - 1]; f[n - 1] = f[n - 2]; f[n - 2] = temp;
}
/* 选主元 */
if( fabs(c[n]) > fabs(b[n - 2]) )
{
temp = c[n]; c[n] = b[n - 2]; b[n - 2] = temp;
temp = a[n]; a[n] = c[n - 2]; c[n - 2] = temp;
temp = b[n]; b[n] = a[n - 2]; a[n - 2] = temp;
temp = f[n]; f[n] = f[n - 2]; f[n - 2] = temp;
}
/* 第(n-1)行消元 */
temp = a[n - 1] / b[n - 2];
a[n - 1] = 0;
b[n - 1] = b[n - 1] - temp * c[n - 2];
c[n - 1] = c[n - 1] - temp * a[n - 2];
f[n - 1] = f[n - 1] - temp * f[n - 2];
/* 第n行消元 */
temp = c[n] / b[n - 2];
c[n] = 0;
a[n] = a[n] - temp * c[n - 2];
b[n] = b[n] - temp * a[n - 2];
f[n] = f[n] - temp * f[n - 2];
/* 选主元 */
if( fabs(a[n]) > fabs(b[n - 1]) )
{
temp = a[n]; a[n] = b[n - 1]; b[n - 1] = temp;
temp = b[n]; b[n] = c[n - 1]; c[n - 1] = temp;
temp = f[n]; f[n] = f[n - 1]; f[n - 1] = temp;
}
/* 最后一次消元 */
temp = a[n] / b[n-1];
a[n] = 0;
b[n] = b[n] - temp * c[n - 1];
f[n] = f[n] - temp * f[n - 1];
/* 回代求解 m[] */
m[n] = f[n] / b[n];
m[n - 1] = (f[n - 1] - c[n - 1] * m[n]) / b[n-1];
for(i = n - 2; i >= 0; i--)
{
m[i] = ( f[i] - (m[i + 2] * a[i] + m[i + 1] * c[i]) ) / b[i];
printf("%f",m[i]);
// printf("\n");
}
/* for (i = 0; i < n+1; i++) printf("%f\n", m[i]); */
free(a);
free(c);
free(f);
/************
** 放置参数
************/
// Y(j+1) - Yj ( 2*Mj + M(j+1) ) Mj M(j+1) - Mj
//Sj(x) = Yj + [ ------------- - -----------------*Hj ]*( x - Xj) + ----(x - Xj)^2 + ------------(x - Xj)^3
// Hj 6 2 6Hj
//
//Hj = X(j+1) - Xj x属于[ Xj , X(j+1) ],j = ( 0 , 1 , 2 ,...,n-1 )
//
//
//a3,b2,c1,d0分别对应三次基函数，二次基函数，一次基函数，零次基函数的系数，直接代入即可。注意各段区间内使用的系数
//
//
for(i = 0; i < n; i++)
{
coe[i].a3 = (m[i + 1] - m[i]) / (6 * h[i]);
coe[i].b2 = m[i] / 2;
coe[i].c1 = d[i] - (h[i] / 6) * (2 * m[i] + m[i + 1]);//d[i] = (y[i + 1] - y[i]) / h[i]
coe[i].d0 = y[i];
}
free(h);
free(d);
free(b);
return n + 1;
}