1、基本思路:namespacesp/樣條曲線曲面域名空間static doubleeps = 1.0e-8;/精 度設定精度void SetEps(double Eps)if(Eps b|a -b;/ k次規范B樣條基函數double N(int i,int k,double su,double u)if(k=0)if(u-su0 =sui&u 0)if(MBiger(sui+k-sui,eps)&MBiger(sui+k+1-sui+1,eps)return (u-sui)*N(i,k-1,su,u)/(sui+k-sui)+(sui+k+1-u)*N(i+1,k-1,su,u)/(sui+k

2、+1-sui+1);else if(!MBiger(sui+k-sui,eps)&MBiger(sui+k+1-sui+1,eps)return (sui+k+1-u)*N(i+1,k-1,su,u)/(sui+k+1-sui+1);else if(MBiger(sui+k-sui,eps)&!MBiger(sui+k+1-sui+1,eps)return (u-sui)*N(i,k-1,su,u)/(sui+k-sui);elsereturn0.0;elsereturn0.0;/B樣條曲線double BSL1(int k,double su,double d,int n,double u)

3、int i;double s=0.0;for(i=0;i n;i+)s=s+di*N(i,k,su,u);return s;德布爾算法計算B樣條曲線inline double alpha(int l,int j,int k,double su,double u)if(suj+k+1-l-suj -eps) return 0.0;return (u-suj)/(suj+k+1-l-suj);double d1(int l,int j,int k,double su,double d,double u)if(l =0) return dj;double a=alpha(l,j,k,su,u);re

4、turn (1-a)*d1(l-1,j-1,k,su,d,u)+a*d1(l-1,j,k,su,d,u);double BSL(int k,double su,double d,int n,double u)int i;for(i=1;i n+k;i+)if(u k&u=sui) break;i-； return d1(k,i,k,su,d,u);德布爾算法計算B樣條曲線的導矢double d2(int l,int j,int k,double su,double d,double u)if(l =0) return dj;double z=suj+k+1-l-suj;if(z -eps) r

5、eturn 0.0;return (k+1-l)*(d2(l-1,j,k,su,d,u)-d2(l-1,j-1,k,su,d,u)/z;double GetDerValue(int k,double su,double d,int n,int r,double u) int i,j;double s=0.0;for(i=1;i n+k;i+)if(u k&u=sui) break;i-;for(j=i-k+r;j =i;j+)s+=d2(r,j,k,su,d,u)*N(j,k-r,su,u); return s;/B樣條曲面double BSS(int k,double su,int l,do

6、uble sv,double *d,int m,double u,int n,double v) int i,j;double s=0.0;for(i=0;i m;i+)for(j=0;j n;j+)s=s+dij*N(i,k,su,u)*N(j,l,sv,v);return s;/NURBS 曲線double NURBSL(int k,double su,double d,double w,int n,double u) int i,t=k;double s1=0.0,s2=0.0,wN;for(i=0;i n+1;i+)if(u =sui)t=i-1;break;if(t k)t=k;fo

7、r(i=t-k;i eps|s2 -eps)return s1/s2;else return 0.0;/NURBS 曲面double NURBSS(int k,double su,int l,double sv,double *d,double *w,int m,double u,int n,double v) int i,j,tu=k,tv=l;double s1=0.0,s2=0.0,wN;for(i=0;i m+1;i+)if(u =sui)tu=i-1;break;for(j=0;j n+1;j+) if(v =svj)tv=j-1;break;if(tu k)tu=k;if(tv l

8、)tv=l;for (i=tu-k;i =tu;i+)for(j=tv-l;j eps|s2 -eps)return s1/s2;else return 0.0; /均勻B樣條曲線double UBSL(int k,double d,int n,double u)double s,*su=new doublen+k+1;int i;for (i=0;i n+k+1;i+)sui=double(i)/(n+k);s=BSL(k,su,d,n,u);delete su;return s; /準均勻B樣條曲線doubleQUBSL(intk,double d,intn,double u)double

9、s,*su=newdoublen+k+1;int i;for(i=0;i k;i+)sui=0.0;for(i=k;i =n;i+) sui=double(i-k)/(n-k);for(i=n+1;i n+k+1;i+)sui=1.0;s=BSL(k,su,d,n,u);delete su;return s;/均勻B樣條曲面double UBSS(int k,int l,double *d,int m,double u,int n,double v)double s,*su=new doublem+k+1,*sv=new doublen+l+1;int i;for( i=0;i m+k+1;i

10、+) sui=double(i)/(m+k);for( i=0;i n+l+1;i+)svi=double(i)/(n+l);s=BSS(k,su,l,sv,d,m,u,n,v);delete su;delete sv;return s; /準均勻B樣條曲面double QUBSS(int k,int l,double *d,int m,double u,int n,double v)double s,*su=new doublem+k+1,*sv=new doublen+l+1;int i;for(i=0;i k;i+)sui=0.0;for(i=k;i =m;i+) sui=double(

11、i-k)/(m-k);for(i=m+1;i m+k+1;i+)sui=1.0;for(i=0;i l;i+)svi=0.0;for(i=l;i =n;i+) svi=double(i-l)/(n-l);for(i=n+1;i n+l+1;i+)svi=1.0;s=BSS(k,su,l,sv,d,m,u,n,v);delete su;delete sv; return s; 非均勻三次B樣條Riesenfeld方法void RiBSL(double *dx,double *dy,int n,double u,double *x,double *y) int i;double *su=new doublen+3;double *l=new doublen;double L=0.0,l1;l0=0.0;for(i=1;i n;i+)li=sqrt(pow2(dx