B样条插值算法

K阶B样条插值应用非常广，其中函数性质也是对称的，通过矩阵求逆，很容易得到系数矩阵。从而得到任意点的值，以及n阶导数。

K阶B样条函数Ω(k,x)的递推性质：Ω(k,x)=1/k*(x+(k+1)/2)*Ω(k-1,x+0.5)-1/k*(x-(k+1)/2)*Ω(k-1,x-0.5)

K阶B样条函数的导数性质：Ω‘(k,x)=Ω(k-1,x+0.5)-Ω(k-1,x-0.5)

在求解系数矩阵的时候需要额外的K-1个方程条件，一般情况下设置首尾一阶导数值或是二阶导数值。具体的条件可以自己确定，从而得到完整的方阵。

代码很简单，如下：

/**

* K阶B样条函数 Ω(k,x)
* 
* Ω(k,x)的导数性质：Ω‘(k,x)=Ω(k-1,x+0.5)-Ω(k-1,x-0.5)
* 
* Ω(k,x)的递推性质：Ω(k,x)=1/k*(x+(k+1)/2)*Ω(k-1,x+0.5)-
1/k*(x-(k+1)/2)*Ω(k-1,x-0.5)
* 根据递推性质即可求得函数
* @param x 
* @param k B样条阶数
* @param n n阶导数,n=0时得到y值
* @return
*/
private static double getOkx(double x,double k,int n){
if(n==0){
if(k==0){
double abs=Math.abs(x);
if(abs==0.5)
return  0.5f;
else if(abs>0.5)
return 0;
else
return 1;
}else{

return 1/k*(x+(k+1)/2)*getOkx( x+0.5, k-1,0)-
1/k*(x-(k+1)/2)*getOkx( x-0.5, k-1,0)
;
}
}else{
return getOkx( x+0.5f, k-1,n-1)-getOkx( x-0.5, k-1,n-1);
}

}

/**
* 根据 S(x)=∑(C*Ω(k,x))可得：
* ∑(C*Ω(k,0))=y[0]
* ∑(C*Ω(k,1))=y[1]
* ∑(C*Ω(k,2))=y[2]
* ......
* ∑(C*Ω(k,n))=y[n]
* 
* 
* 以及加入条件方程组，得到N+K矩阵
* 
* A[i]=∑Ω(k,i),V=y
* 
* A*C=V --> C=A逆*V
* 
* @return C
*/
public double[] getCi(){
double[][]A=new double[N+K][N+K];
double[][]V=new double[N+K][1];
for(int i=0;i<=N;i++){
for(int j=0;j<N+K;j++){
A[i][j]=getOkx(i-j+(K-1)/2,K,0);
}
V[i][0]=y[i];
}

//还有K-1个条件(n阶导数为0)
//自定义K-1个条件方程
//这里是x0的导数=y0d,xn的导数=ynd;

int m=0;
for(int i=N+1;i<N+K;i++){
if(m==0){
for(int j=0;j<N+K;j++){
A[i][j]=getOkx(0-j+(K-1)/2,K,1);//x0的一阶导数系数
}
V[i][0]=y0d;
}else if(m==1){
for(int j=0;j<N+K;j++){
A[i][j]=getOkx(N-j+(K-1)/2,K,1);//xn的一阶导数导数系数
}
V[i][0]=ynd;
}else{
for(int j=0;j<N+K;j++){
A[i][j]=getOkx(0-j+(K-1)/2,K,m);
}
V[i][0]=0;
}
m++;
}

double[][]A1=inverseMatrix(A);//矩阵求逆
double [][]t=times(A1, V);//矩阵相乘
double []C=new double[N+K];
for(int i=0;i<t.length;i++){
C[i]=t[i][0];
}
return C;//返回系数

}

/**
* 获取曲线上的点S(x)=∑(C*Ω(k,x))
* @param x （x,y）
* @param C B样条函数系数数组
* @param n n阶导数,n=0时得到S(x)值
* @return S(x), S'(x),S''(x)....
*/
public static double getY(double x,double[]C,int n){
double rst=0;
if(x<0){
return 0;
}
int size=0;
for(int i=(int) x;i<C.length;i++){
double temp=getOkx(x-i+(K-1)/2,K,n);
rst+=C[i]*temp;
size++;
if(size==K+1){
break;
}
}
return rst;
}

完整的例子见：http://download.csdn.net/detail/zhq1314zhq/9746788