>>Drew的主页--->个人兴趣--------->NURBS算法
应一些网友的要求，drew在这里提供以前在学校写的非均匀有理B样条（NURBS）曲线曲面造型的主要算法函数源码，供感兴趣的朋友参考。
由于这是多年前在学校所做的程序，为方便编程，坐标数据主要采用三维数组形式。当时使用计算机内存很小（好像是4M内存）,所以所取控制点数有一定的限制(15个），其它未列的函数是一些数据结构，内存分配和空间转换的函数，和算法无关，可以不必考虑。由于这些源码函数是从drew的论文程序中摘录出来的，所以需要做一些修改才能使用，朋友们可以参考和修改这些算法函数，加入自己的程序。
关于NURBS，曲线曲面造型算法最好的中文资料，drew推荐北航施法中教授编写的《计算机辅助几何设计与非均匀有理B样条(CAGD&NURBS)》.
//哈特利(Hartley)-贾德(Judd,1978)的弦长参数化算法求节点矢量U,用于由控制点正算曲线曲面
//coeff:控制顶点。
//v1,v2:U向，W向的节点序列。 count,countj:控制点行数和列数。void Hartley_Judd(coeff,v1,v2,counti,countj)
float coeff[15][15][3];
float v1[20],v2[20];
int counti,countj;
{
int i,s,j,k=DEGREE_J,n1=counti-1,n2=countj-1;
float l1[20],l2[20];
float ll01,ll11,ll02,ll12;
for(i=1;i<=n1;i++)
l1[i]=(float)sqrt((coeff[i][0][0]-coeff[i-1][0][0])*
(coeff[i][0][0]-coeff[i-1][0][0])+
(coeff[i][0][1]-coeff[i-1][0][1])*
(coeff[i][0][1]-coeff[i-1][0][1])+
(coeff[i][0][2]-coeff[i-1][0][2])*
(coeff[i][0][2]-coeff[i-1][0][2]));
for(i=1;i<=n2;i++)
l2[i]=(float)sqrt((coeff[0][i][0]-coeff[0][i-1][0])*
(coeff[0][i][0]-coeff[0][i-1][0])+
(coeff[0][i][1]-coeff[0][i-1][1])*
(coeff[0][i][1]-coeff[0][i-1][1])+
(coeff[0][i][2]-coeff[0][i-1][2])*
(coeff[0][i][2]-coeff[0][i-1][2]));
ll01=0;
ll11=0;
for(i=0;i<=k;i++) v1[i]=0;
for(i=k+1;i<=n1+1;i++)
{
for(j=i-k;j<=i-1;j++)
ll01=ll01+l1[j];
for(s=k+1;s<=n1+1;s++)
for(j=s-k;j<=s-1;j++)
ll11=ll11+l1[j];
v1[i]=ll01/ll11+v1[i-1];
}
for(i=n1+1;i<=n1+k+1;i++) v1[i]=1;
if(v1[2]==1) v1[2]=0;
ll02=0;
ll12=0;
for(i=0;i<=k;i++) v2[i]=0;
for(i=k+1;i<=n2+1;i++)
{
for(j=i-k;j<=i-1;j++)
ll02=ll02+l2[j];
for(s=k+1;s<=n2+1;s++)
for(j=s-k;j<=s-1;j++)
ll12=ll12+l2[j];
v2[i]=ll02/ll12+v2[i-1];
}
for(i=n2+1;i<=n2+k+1;i++) v2[i]=1;
if(v2[2]==1) v2[2]=0;
}
// 伯姆(Boehm)节点插入算法：
// coeff:原控制顶点.v:原节点序列.uu:新插入节点.counti:原控制点数.d:新控制顶点
void boehm(coeff,v,uu,counti,d)
float coeff[15][3],v[20],d[15][3],uu;
int counti;
{ int r=0,g=0,i,j;
float a[15],vd[20];
for(i=0;i<=counti+DEGREE_J;i++)
{ if(uu==v[i]) r++;
else if(uu>=v[i]) g++;
}
g--;
for(i=0;i<=g-DEGREE_J;i++)
for(j=0;j<=2;j++)
d[i][j]=coeff[i][j];
for(i=g-DEGREE_J+1;i<=g-r;i++)
for(j=0;j<=2;j++)
{ a[i]=(float)((uu-v[i])/(v[i+DEGREE_J]-v[i]));
d[i][j]=(1-a[i])*coeff[i-1][j]+a[i]*coeff[i][j];
}
for(i=g-r+1;i<=counti;i++)
for(j=0;j<=2;j++)
d[i][j]=coeff[i-1][j];
for(i=0;i<=g;i++)
vd[i]=v[i];
vd[g+1]=(float)(uu);
for(i=g+1;i<=counti+DEGREE_J+1;i++)
vd[i+1]=v[i];
for(i=0;i<=counti+DEGREE_J+1;i++) v[i]=vd[i];
}//德布尔-考克斯递推公式 for 重节点
void deBoor_2D(float uw,float n[15][4],int count,float v[20])
{ int i,k=0;
for(i=0;i<=count+1;i++)
{ if(uw>=v[i]&&uw<=v[i+1]) n[i][0]=1;
else n[i][0]=0;
}
for(k=1;k<=DEGREE;k++)
for(i=0;i<=count+1;i++)
{ if((v[i+k]-v[i])==0)
{ if((v[i+k+1]-v[i+1])==0) n[i][k]=0;
else
n[i][k]=((v[i+k+1]-uw)*n[i+1][k-1])/(v[i+k+1]-v[i+1]);
}
else if((v[i+k+1]-v[i+1])==0)
n[i][k]=((uw-v[i])*n[i][k-1])/(v[i+k]-v[i]);
else
n[i][k]=(((uw-v[i])*n[i][k-1])/(v[i+k]-v[i]))
+(((v[i+k+1]-uw)*n[i+1][k-1])/(v[i+k+1]-v[i+1]));
}
}
//由给定曲线，曲面型值点反算控制点及节点矢量算法和处理重节点情况算法。
//cff:给定型值点数组. ww:权因子数组。 coeff:反算出的控制顶点。
//v1,v2:U向，W向的节点序列。 count,countj:控制点行数和列数。
void Counter_count(cff,ww,coeff,v1,v2,count,countj)
float cff[15][15][3],coeff[15][15][3],v1[20],v2[20],ww[15][15];
int count,countj;
{ float vu[20],bb[35],bb1[35],n[15][4],nn[15][4];
float bx[15],by[15],bz[15];
float u=0,p0_=HEAD_QIE,pn_=END_QIE,cc=0,ccc;
int i,j,k,k1,kkk;
DEGREE=3;
DEGREE_I=3;
DEGREE_J=3;
CURVE=1;
vu[0]=0;
cc=0;
for(i=1;i<=count-1;i++)
cc=cc+(float)sqrt((cff[i][0][0]-cff[i-1][0][0])
*(cff[i][0][0]-cff[i-1][0][0])
+(cff[i][0][1]-cff[i-1][0][1])
*(cff[i][0][1]-cff[i-1][0][1])
+(cff[i][0][2]-cff[i-1][0][2])
*(cff[i][0][2]-cff[i-1][0][2]));
for(i=1;i<=count-1;i++)
vu[i]=vu[i-1]+(float)sqrt((cff[i][0][0]-cff[i-1][0][0])
*(cff[i][0][0]-cff[i-1][0][0])
+(cff[i][0][1]-cff[i-1][0][1])
*(cff[i][0][1]-cff[i-1][0][1])
+(cff[i][0][2]-cff[i-1][0][2])
*(cff[i][0][2]-cff[i-1][0][2]))/cc;
for(i=0;i<=3;i++) v1[i]=0;
for(i=4;i<=count-1;i++) v1[i]=vu[i-3];
for(i=count;i<=count+3;i++) v1[i]=1;
if(count==2)
{ for(i=0;i<=3;i++) v1[i]=0;
for(i=4;i<=7;i++) v1[i]=1.1F;
}
for(k=0;k<=count-1;k++)
{ vu[0]=0;
cc=0;
for(k1=1;k1<=countj-1;k1++)
cc=cc+(float)sqrt((cff[k][k1][0]-cff[k][k1-1][0])
*(cff[k][k1][0]-cff[k][k1-1][0])
+(cff[k][k1][1]-cff[k][k1-1][1])
*(cff[k][k1][1]-cff[k][k1-1][1])
+(cff[k][k1][2]-cff[k][k1-1][2])
*(cff[k][k1][2]-cff[k][k1-1][2]));
for(k1=1;k1<=countj-1;k1++)
vu[k1]=vu[k1-1]+(float)sqrt((cff[k][k1][0]-cff[k][k1-1][0])
*(cff[k][k1][0]-cff[k][k1-1][0])
+(cff[k][k1][1]-cff[k][k1-1][1])
*(cff[k][k1][1]-cff[k][k1-1][1])
+(cff[k][k1][2]-cff[k][k1-1][2])
*(cff[k][k1][2]-cff[k][k1-1][2]))/cc;
for(k1=0;k1<=3;k1++) v2[k1]=0;
for(k1=4;k1<=countj+2;k1++) v2[k1]=vu[k1-3];
for(k1=countj+2;k1<=countj+5;k1++) v2[k1]=1;
for(i=0;i<=2;i++)
{ coeff[k][0][i]=cff[k][0][i];
coeff[k][countj+1][i]=cff[k][countj-1][i];
}
for(i=0;i<=2;i++)
{ coeff[k][1][i]=p0_*ww[k][0]/(2*ww[k][1])+coeff[k][0][i];
coeff[k][countj][i]=coeff[k][countj+1][i]
-pn_*ww[k][countj+1]/(2*ww[k][countj]);
}
deBoor_2D(v2[4],n,countj+2,v2);
cc=Nurbs_bzer1(ww,u,v2[4],v1,v2,1,countj+2);
bb[0]=ww[k][2]*n[2][3]/cc;
bb[1]=ww[k][3]*n[3][3]/cc;
for(i=0;i<=2;i++)
cff[k][1][i]=cff[k][1][i]-coeff[k][1][i]*ww[k][1]*n[1][3]/cc;
deBoor_2D(v2[countj+1],n,countj+2,v2);
cc=Nurbs_bzer1(ww,u,v2[countj+1],v1,v2,1,countj+2);
bb[(countj-4)*3+2]=ww[k][countj-2]*n[countj-2][3]/cc;
bb[(countj-4)*3+3]=ww[k][countj-1]*n[countj-1][3]/cc;
for(i=0;i<=2;i++)
cff[k][countj-2][i]=cff[k][countj-2][i]-coeff[k][countj][i]
*ww[k][countj]*n[countj][3]/cc;
i=0;
j=2;
while((j+3)<=countj)
{ deBoor_2D(v2[j+3],n,countj+2,v2);
cc=Nurbs_bzer1(ww,u,v2[j+3],v1,v2,1,countj+2);
bb[i+2]=ww[k][j]*n[j][3]/cc;
bb[i+3]=ww[k][j+1]*n[j+1][3]/cc;
bb[i+4]=ww[k][j+2]*n[j+2][3]/cc;
i=i+3;
j++;
}
j=1;
i=0;
loop: while(i<=countj-2)
{ if(cff[k][i][0]==cff[k][i+1][0]&&cff[k][i][1]==cff[k][i+1][1]
&&cff[k][i][2]==cff[k][i+1][2])
{
j++;
if(cff[k][i+1][0]==cff[k][i+2][0]
&&cff[k][i+1][1]==cff[k][i+2][1]
&&cff[k][i+1][2]==cff[k][i+2][2])
j++;
}
if(j>1) break;
i++;
}
kkk=i;
if(j>=2)
{ for(k1=0;k1<=2;k1++) cff[k][kkk][k1]=1;
deBoor_2D(v2[kkk+3],n,countj+2,v2);
cc=Nurbs_bzer1(ww,u,v2[kkk+3],v1,v2,1,countj+2);
nn[kkk][3]=3*(n[kkk][2]/(v2[kkk+3]-v2[kkk])
-n[kkk+1][2]/(v2[kkk+4]-v2[kkk+1]));
nn[kkk+1][3]=3*(n[kkk+1][2]/(v2[kkk+4]-v2[kkk+1])
-n[kkk+2][2]/(v2[kkk+5]-v2[kkk+2]));
nn[kkk+2][3]=3*(n[kkk+2][2]/(v2[kkk+5]-v2[kkk+2])
-n[kkk+3][2]/(v2[kkk+6]-v2[kkk+3]));
ccc=0;
for(k1=kkk;k1<=kkk+2;k1++) ccc=ccc+ww[k][k1]*nn[k1][3];
bb[(kkk-2)*3+2]=ww[k][kkk]*(nn[kkk][3]*cc-n[kkk][3]*ccc)/(cc*cc);
bb[(kkk-2)*3+3]=ww[k][kkk+1]*(nn[kkk+1][3]*cc-n[kkk+1][3]*ccc)/(cc*cc);
bb[(kkk-2)*3+4]=ww[k][kkk+2]*(nn[kkk+2][3]*cc-n[kkk+1][3]*ccc)/(cc*cc);
}
i=i+j;
if(i<=countj-2) { j=1; goto loop; }
for(i=0;i<=countj-3;i++)
{ bx[i]=cff[k][i+1][0];
by[i]=cff[k][i+1][1];
bz[i]=cff[k][i+1][2];
}
for(i=0;i<=(countj-4)*3+3;i++) bb1[i]=bb[i];
run_run(bb1,countj-2,4+3*(countj-4),bx);
for(i=0;i<=(countj-4)*3+3;i++) bb1[i]=bb[i];
run_run(bb1,countj-2,4+3*(countj-4),by);
run_run(bb,countj-2,4+3*(countj-4),bz);
for(i=2;i<=countj-1;i++)
{ coeff[k][i][0]=bx[i-2];
coeff[k][i][1]=by[i-2];
coeff[k][i][2]=bz[i-2];
}
}
CURVE=0;
}
//NURBS算法，由控制顶点，权因子，节点序列等计算插值点集
void cal_surf(SURF *ps,HDC hdc)
{ GRID *pg;
float coeffx[15][15],coeffy[15][15],coeffz[15][15];
int m,n,i,k,j,mark=0;
float u,w,nn;
o=0;p=0;
DEGREE_I=ps->degree_u;
DEGREE_J=ps->degree_w;
for(i=0;i<=ps->nums_u-1;i++)
for(j=0;j<=ps->nums_w-1;j++)
{ coeffx[i][j]=ps->pro[i][j].x;
coeffy[i][j]=ps->pro[i][j].y;
coeffz[i][j]=ps->pro[i][j].z;
}
m=ps->sgridn_u;n=ps->sgridn_w;
if(ROTA_FLG==0)
{ ALLOC(ps->pgrid,GRID);
pg=ps->pgrid;
}
if(CURVE==1)
{ if(NC_STYLE==1) { NC_STYLE=0; mark=1; } m=1; }
for(i=0;i<m+1;i++)
{ u=(float)(i*1./m);
p=0;
for(k=0;k<n+1;k++)
{ w=(float)(k*1./n);
if(NC_STYLE!=1)
{ nn=Nurbs_bzer1(ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,ps->nums_w);
x_[o][p]=L_S_bzer(coeffx,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
y_[o][p]=L_S_bzer(coeffy,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
z_[o][p]=L_S_bzer(coeffz,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
}
else
{ nn=Nurbs_bzer1(ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,ps->nums_w);
x_[o][p]=L_S_bzer(coeffx,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
y_[o][p]=L_S_bzer(coeffy,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
z_[o][p]=L_S_bzer(coeffz,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
ps->nums_w)/nn;
}
if(ROTA_FLG==0)
{ ALLOC(pg->next,GRID);
pg=pg->next;
// pg->value_u=u; pg->value_w=w;
pg->value_g.x=x_[o][p];
pg->value_g.y=y_[o][p];
pg->value_g.z=z_[o][p];
}
if(MAIN_W==1)
{ if(p==0) MoveTo(hdc,(int)(x_[o][p]),(int)(y_[o][p]));
else LineTo(hdc,(int)(x_[o][p]),(int)(y_[o][p]));
}
if(LEFT_W==1)
{ if(p==0) MoveTo(hdc,(int)(z_[o][p]+320),(int)(y_[o][p]));
else LineTo(hdc,(int)(z_[o][p]+320),(int)(y_[o][p]));
}
if(DOWN_W==1)
{ if(p==0) MoveTo(hdc,(int)(x_[o][p]),(int)(z_[o][p]+220));
else LineTo(hdc,(int)(x_[o][p]),(int)(z_[o][p]+220));
}
p++;
}
o++;
}
if(ROTA_FLG==0) pg->next=NULL;
if(CURVE!=1)
{ if(MAIN_W==1)
for(i=0;i<=p-1;i++)
for(j=0;j<=o-1;j++)
{
if(j==0) MoveTo(hdc,(int)(x_[j][i]),(int)(y_[j][i]));
LineTo(hdc,(int)(x_[j][i]),(int)(y_[j][i]));
}
if(LEFT_W==1)
for(i=0;i<=p-1;i++)
for(j=0;j<=o-1;j++)
{ if(j==0) MoveTo(hdc,(int)(z_[j][i]+320),(int)(y_[j][i]));
LineTo(hdc,(int)(z_[j][i]+320),(int)(y_[j][i]));
}
if(DOWN_W==1)
for(i=0;i<=p-1;i++)
for(j=0;j<=o-1;j++)
{ if(j==0) MoveTo(hdc,(int)(x_[j][i]),(int)(z_[j][i]+220));
LineTo(hdc,(int)(x_[j][i]),(int)(z_[j][i]+220));
}
}
if(WATCH_3==1)
{ for(i=0;i<4;i++)
for(j=0;j<4;j++)
{ ff[i][j]=0;
if(i==j) ff[i][j]=0.5f;
}
ff[3][0]=(float)((1-0.25)*320);
ff[3][1]=(float)((1-0.25)*220);
ff[3][2]=0.25f; ff[3][3]=1;
affine3(x_,y_,z_);
for(i=0;i<4;i++)
for(j=0;j<4;j++)
{ ff[i][j]=0;
if(i==j) ff[i][j]=1;
}
ff[3][0]=-240; ff[3][1]=-140; ff[3][2]=0; ff[3][3]=1;
affine3(x_,y_,z_); //XOY
carton_draw(x_,y_,hdc);
ff[3][0]=320;ff[3][1]=0;ff[3][2]=480;
affine3(x_,y_,z_);
carton_draw(z_,y_,hdc);
ff[3][0]=-320;ff[3][1]=0;ff[3][2]=-150;
affine3(x_,y_,z_);
carton_draw(x_,z_,hdc);
ff[3][0]=-190;ff[3][1]=620;ff[3][2]=0;
affine3(x_,y_,z_);
rotate_any(320.,220.,0.,200.,0.,0.,45*PI/180);
affine3(x_,y_,z_);
rotate_any(320.,220.,0.,0.,290.,0.,45*PI/180);
affine3(x_,y_,z_);
carton_draw(x_,y_,hdc);
}
if(mark==1) NC_STYLE=1;
}
//追赶法求解方程组
//b:系数数组。 n:方程个数。 m:系数个数。 d:未知数
void run_run(b,n,m,d)
int n,m;
float b[],d[];
{ int k,j;
float s;
if(m!=(3*n-2)) return;
for(k=0;k<=n-2;k++)
{ j=3*k;
s=b[j];
if(fabs(s)+1.0==1.0)
{ MessageBox(NULL,"计算错误","系统消息",MB_OK|MB_ICONEXCLAMATION);
return;
}
b[j+1]=b[j+1]/s;
d[k]=d[k]/s;
b[j+3]=b[j+3]-b[j+2]*b[j+1];
d[k+1]=d[k+1]-b[j+2]*d[k];
}
s=b[3*n-3];
if(fabs(s)+1.0==1.0) return;
d[n-1]=d[n-1]/s;
for(k=n-2;k>=0;k--)
d[k]=d[k]-b[3*k+1]*d[k+1];
}
float Nurbs_bzer1(ww,u,w,v1,v2,counti,countj)
float ww[15][15],v1[20],v2[20];
float u,w;
int counti,countj;
{ float n2[15][4],n1[15][4];
float c;
int i,j;
DEGREE=DEGREE_I;
deBoor_2D(u,n2,counti,v1);
DEGREE=DEGREE_J;
deBoor_2D(w,n1,countj,v2);
c=0;
if(CURVE==1) n2[0][DEGREE_I]=1;
for(i=0;i<=counti-1;i++)
for(j=0;j<=countj-1;j++)
c=c+ww[i][j]*n2[i][DEGREE_I]*n1[j][DEGREE_J];
return(c);
}