Nurbs中文资料.docx

应一些网友的要求，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 coeff15153;float v120,v220;int counti,countj; int i,s,j,k=DEGREE_J,n1=counti-1,n2=countj-1; float l120,l220; float ll01,ll11,ll02,ll12;for(i=1;i=n1;i+) l1i=(float)sqrt(coeffi00-coeffi-100)* (coeffi00-coeffi-100)+ (coeffi01-coeffi-101)* (coeffi01-coeffi-101)+ (coeffi02-coeffi-102)* (coeffi02-coeffi-102);for(i=1;i=n2;i+) l2i=(float)sqrt(coeff0i0-coeff0i-10)* (coeff0i0-coeff0i-10)+ (coeff0i1-coeff0i-11)* (coeff0i1-coeff0i-11)+ (coeff0i2-coeff0i-12)* (coeff0i2-coeff0i-12);ll01=0;ll11=0;for(i=0;i=k;i+) v1i=0;for(i=k+1;i=n1+1;i+) for(j=i-k;j=i-1;j+) ll01=ll01+l1j; for(s=k+1;s=n1+1;s+) for(j=s-k;j=s-1;j+) ll11=ll11+l1j; v1i=ll01/ll11+v1i-1;for(i=n1+1;i=n1+k+1;i+) v1i=1;if(v12=1) v12=0;ll02=0;ll12=0;for(i=0;i=k;i+) v2i=0;for(i=k+1;i=n2+1;i+) for(j=i-k;j=i-1;j+) ll02=ll02+l2j; for(s=k+1;s=n2+1;s+) for(j=s-k;j=s-1;j+) ll12=ll12+l2j; v2i=ll02/ll12+v2i-1;for(i=n2+1;i=n2+k+1;i+) v2i=1; if(v22=1) v22=0;/ 伯姆(Boehm)节点插入算法：/ coeff:原控制顶点.v:原节点序列.uu:新插入节点.counti:原控制点数.d:新控制顶点void boehm(coeff,v,uu,counti,d)float coeff153,v20,d153,uu;int counti; int r=0,g=0,i,j; float a15,vd20; for(i=0;i=vi) g+; g-;for(i=0;i=g-DEGREE_J;i+) for(j=0;j=2;j+) dij=coeffij;for(i=g-DEGREE_J+1;i=g-r;i+) for(j=0;j=2;j+) ai=(float)(uu-vi)/(vi+DEGREE_J-vi); dij=(1-ai)*coeffi-1j+ai*coeffij; for(i=g-r+1;i=counti;i+) for(j=0;j=2;j+) dij=coeffi-1j;for(i=0;i=g;i+) vdi=vi;vdg+1=(float)(uu);for(i=g+1;i=counti+DEGREE_J+1;i+) vdi+1=vi;for(i=0;i=counti+DEGREE_J+1;i+) vi=vdi;/德布尔-考克斯递推公式 for 重节点void deBoor_2D(float uw,float n154,int count,float v20) int i,k=0; for(i=0;i=vi&uw=vi+1) ni0=1; else ni0=0;for(k=1;k=DEGREE;k+) for(i=0;i=count+1;i+) if(vi+k-vi)=0) if(vi+k+1-vi+1)=0) nik=0; else nik=(vi+k+1-uw)*ni+1k-1)/(vi+k+1-vi+1); else if(vi+k+1-vi+1)=0) nik=(uw-vi)*nik-1)/(vi+k-vi); else nik=(uw-vi)*nik-1)/(vi+k-vi) +(vi+k+1-uw)*ni+1k-1)/(vi+k+1-vi+1); /由给定曲线，曲面型值点反算控制点及节点矢量算法和处理重节点情况算法。/cff:给定型值点数组. ww:权因子数组。 coeff:反算出的控制顶点。/v1,v2:U向，W向的节点序列。 count,countj:控制点行数和列数。void Counter_count(cff,ww,coeff,v1,v2,count,countj)float cff15153,coeff15153,v120,v220,ww1515;int count,countj; float vu20,bb35,bb135,n154,nn154; float bx15,by15,bz15; 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; vu0=0; cc=0; for(i=1;i=count-1;i+) cc=cc+(float)sqrt(cffi00-cffi-100) *(cffi00-cffi-100) +(cffi01-cffi-101) *(cffi01-cffi-101) +(cffi02-cffi-102) *(cffi02-cffi-102); for(i=1;i=count-1;i+) vui=vui-1+(float)sqrt(cffi00-cffi-100) *(cffi00-cffi-100) +(cffi01-cffi-101) *(cffi01-cffi-101) +(cffi02-cffi-102) *(cffi02-cffi-102)/cc; for(i=0;i=3;i+) v1i=0; for(i=4;i=count-1;i+) v1i=vui-3; for(i=count;i=count+3;i+) v1i=1; if(count=2) for(i=0;i=3;i+) v1i=0; for(i=4;i=7;i+) v1i=1.1F; for(k=0;k=count-1;k+) vu0=0; cc=0; for(k1=1;k1=countj-1;k1+) cc=cc+(float)sqrt(cffkk10-cffkk1-10) *(cffkk10-cffkk1-10) +(cffkk11-cffkk1-11) *(cffkk11-cffkk1-11) +(cffkk12-cffkk1-12) *(cffkk12-cffkk1-12); for(k1=1;k1=countj-1;k1+) vuk1=vuk1-1+(float)sqrt(cffkk10-cffkk1-10) *(cffkk10-cffkk1-10) +(cffkk11-cffkk1-11) *(cffkk11-cffkk1-11) +(cffkk12-cffkk1-12) *(cffkk12-cffkk1-12)/cc; for(k1=0;k1=3;k1+) v2k1=0; for(k1=4;k1=countj+2;k1+) v2k1=vuk1-3; for(k1=countj+2;k1=countj+5;k1+) v2k1=1; for(i=0;i=2;i+) coeffk0i=cffk0i; coeffkcountj+1i=cffkcountj-1i; for(i=0;i=2;i+) coeffk1i=p0_*wwk0/(2*wwk1)+coeffk0i; coeffkcountji=coeffkcountj+1i -pn_*wwkcountj+1/(2*wwkcountj); deBoor_2D(v24,n,countj+2,v2); cc=Nurbs_bzer1(ww,u,v24,v1,v2,1,countj+2); bb0=wwk2*n23/cc; bb1=wwk3*n33/cc; for(i=0;i=2;i+) cffk1i=cffk1i-coeffk1i*wwk1*n13/cc; deBoor_2D(v2countj+1,n,countj+2,v2); cc=Nurbs_bzer1(ww,u,v2countj+1,v1,v2,1,countj+2); bb(countj-4)*3+2=wwkcountj-2*ncountj-23/cc; bb(countj-4)*3+3=wwkcountj-1*ncountj-13/cc; for(i=0;i=2;i+) cffkcountj-2i=cffkcountj-2i-coeffkcountji *wwkcountj*ncountj3/cc; i=0; j=2; while(j+3)=countj) deBoor_2D(v2j+3,n,countj+2,v2); cc=Nurbs_bzer1(ww,u,v2j+3,v1,v2,1,countj+2); bbi+2=wwkj*nj3/cc; bbi+3=wwkj+1*nj+13/cc; bbi+4=wwkj+2*nj+23/cc; i=i+3; j+; j=1; i=0;loop: while(i1) break; i+; kkk=i; if(j=2) for(k1=0;k1=2;k1+) cffkkkkk1=1; deBoor_2D(v2kkk+3,n,countj+2,v2); cc=Nurbs_bzer1(ww,u,v2kkk+3,v1,v2,1,countj+2); nnkkk3=3*(nkkk2/(v2kkk+3-v2kkk) -nkkk+12/(v2kkk+4-v2kkk+1); nnkkk+13=3*(nkkk+12/(v2kkk+4-v2kkk+1) -nkkk+22/(v2kkk+5-v2kkk+2); nnkkk+23=3*(nkkk+22/(v2kkk+5-v2kkk+2) -nkkk+32/(v2kkk+6-v2kkk+3); ccc=0; for(k1=kkk;k1=kkk+2;k1+) ccc=ccc+wwkk1*nnk13; bb(kkk-2)*3+2=wwkkkk*(nnkkk3*cc-nkkk3*ccc)/(cc*cc); bb(kkk-2)*3+3=wwkkkk+1*(nnkkk+13*cc-nkkk+13*ccc)/(cc*cc); bb(kkk-2)*3+4=wwkkkk+2*(nnkkk+23*cc-nkkk+13*ccc)/(cc*cc); i=i+j; if(i=countj-2) j=1; goto loop; for(i=0;i=countj-3;i+) bxi=cffki+10; byi=cffki+11; bzi=cffki+12; for(i=0;i=(countj-4)*3+3;i+) bb1i=bbi; run_run(bb1,countj-2,4+3*(countj-4),bx); for(i=0;i=(countj-4)*3+3;i+) bb1i=bbi; run_run(bb1,countj-2,4+3*(countj-4),by); run_run(bb,countj-2,4+3*(countj-4),bz); for(i=2;idegree_u; DEGREE_J=ps-degree_w; for(i=0;inums_u-1;i+) for(j=0;jnums_w-1;j+) coeffxij=ps-proij.x; coeffyij=ps-proij.y; coeffzij=ps-proij.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;im+1;i+) u=(float)(i*1./m); p=0; for(k=0;kwww,u,w,ps-v_u,ps-v_w,ps-nums_u,ps-nums_w); x_op=L_S_bzer(coeffx,ps-www,u,w,ps-v_u,ps-v_w,ps-nums_u, ps-nums_w)/nn; y_op=L_S_bzer(coeffy,ps-www,u,w,ps-v_u,ps-v_w,ps-nums_u, ps-nums_w)/nn; z_op=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_op=L_S_bzer(coeffx,ps-www,w,u,ps-v_u,ps-v_w,ps-nums_u, ps-nums_w)/nn; y_op=L_S_bzer(coeffy,ps-www,w,u,ps-v_u,ps-v_w,ps-nums_u, ps-nums_w)/nn; z_op=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_op; pg-value_g.y=y_op; pg-value_g.z=z_op; if(MAIN_W=1) if(p=0) MoveTo(hdc,(int)(x_op),(int)(y_op); else LineTo(hdc,(int)(x_op),(int)(y_op); if(LEFT_W=1) if(p=0) MoveTo(hdc,(int)(z_op+320),(int)(y_op); else LineTo(hdc,(int)(z_op+320),(int)(y_op); if(DOWN_W=1) if(p=0) MoveTo(hdc,(int)(x_op),(int)(z_op+220); else LineTo(hdc,(int)(x_op),(int)(z_op+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_ji),(int)(y_ji); LineTo(hdc,(int)(x_ji),(int)(y_ji); 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_ji+320),(int)(y_ji); LineTo(hdc,(int)(z_ji+320),(int)(y_ji); 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_ji),(int)(z_ji+220); LineTo(hdc,(int)(x_ji),(int)(z_ji+220); if(WATCH_3=1) for(i=0;i4;i+) for(j=0;j4;j+) ffij=0; if(i=j) ffij=0.5f; ff30=(float)(1-0.25)*320); ff31=(float)(1-0.25)*220); ff32=0.25f; ff33=1; affine3(x_,y_,z_); for(i=0;i4;i+) for(j=0;j4;j+) ffij=0; if(i=j) ffij=1; ff30=-240; ff31=-140; ff32=0; ff33=1; affine3(x_,y_,z_); /XOY carton_draw(x_,y_,hdc); ff30=320;ff31=0;ff32=480; affine3(x_,y