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
