曲面加工几何学

1、大连理工大学微分几何作业学 院： 机械工程学院 姓 名： 学 号： 授课教师： 王晓明 电子邮箱： 完成时间：2012年6月25日非均匀有理B样条（NURBS）曲线和曲面1. 非均匀有理B样条曲线（3阶）NURBS曲线提供了对标准解析几何和自由曲线、曲面的统一数学描述方法，它可通过调整控制顶点和权因子，方便地改变曲面形状，同时也可方便地转换成对应Bezier曲面。给定n+1个控制点及权因子，则k阶（k-1）次NURBS曲线的表达式为：其中：为非均匀有理B样条基函数，按照deBoor-Cox递推公式，基函数满足一下关系：其它1.1. 曲线的建立本文选取四个控制顶点P1,，P2,，P3,，P4,，

2、已经四个权因子w1，w2，w4，w4，通过Matlab计算其基函数，并作出曲线。控制顶点：每个坐标有三个值 立体的；权因子：；非均匀节点：；由Matlab绘制曲线如下：在每个点上标注其Frenet标架如图：在每个点上标注其曲率、挠率如图： 2. 非均匀有理B样条曲面（3阶）假设给定R3空间中，齐次坐标下的个点，则参数曲面：称为矩形域0,1x0,1上的mxn次有理B样条曲面。其中称为控制顶点，称为权因子。NURBS曲面的表达式为：2.1. 曲面计算的基本公式曲面第一基本形式：曲面第二基本形式：高斯曲率计算公式：平均曲率计算公式：曲面主曲率计算公式： ， 2.2. 曲面的建立选择16个顶点构造B样

3、条曲面，如表1所示：表1 16个顶点坐标数据序号xyz1-24.8849.76-4.352-23.660.560.153-19.4771.915.174-11.9981.099.785-22.8950.75-3.016-21.0361.451.447-16.6572.446.348-9.4781.3510.729-20.2551.91-1.2910-17.6662.433.1811-12.9672.967.9112-6.1581.5711.9513-16.3453.371.0414-12.8163.545.5715-7.6173.4110.0416-1.4581.7213.67权因子：；计算结

4、果如图： 图1 B样条曲面高斯曲率K 图2 B样条曲面平均曲率H 图3 E值 图4 F值 图5 G值 图6 L值 图7 L值 图8 N值 图9 k1值 图10 k2值3. 参考文献1 王省富.样条函数及其应用M. 西北工业大学出版社，1989,9，西安.2 张化光，刘鑫蕊，孙秋野.MATLAB/SIMULINK实用教程M.人民邮电出版社，2009，3，北京.3 John H. Mathews, Kurtis D. Fink.数值方法（MATLAB版）.电子工业出版社，2002,6，北京.附录（Matlab源程序）：l 非均匀有理B样条曲线（3阶）function NURBS3clear;clc

5、;syms rP= 0 1 0; 0 0 0.5; 1 0 1; 1 1 1.5;w=2,0.5,1,3; %插入点权重% w=2,2,2,2; %权重相等时u=1 2 4 5 6.5 7.5 8 9.3;%非均匀节点（构造样条基函数）k=3; %三阶非均匀有理B样条曲线n=size(P,1); %控制点数目nu=size(u,2); %非均匀节点数if nu=n+k+1 M=100; r1=zeros(M,3); r2=zeros(M,3); r=zeros(M,3); X=linspace(u(k+1),u(n+1),M); for l=1:M %构造非均匀有理B样条基函数 x=X(l);

6、 N=zeros(n+k,k+1); for i=1:n+k if x>=u(i)&&x<u(i+1) N(i,1)=1; end end for j=2:k+1 for i=1:n+k-j+1 N(i,j)=(x-u(i)/(u(i+j-1)-u(i)*N(i,j-1)+(u(i+j)-x)/(u(i+j)-u(i+1)*N(i+1,j-1); end end for m1=1:n for m2=1:3 r1(l,m2)=r1(l,m2)+ w(m1)*N(m1,k+1)*P(m1,m2); %分子 end r2(l,:)=r2(l,:)+w(m1)*N(m1,k

7、+1);%分母 end end r=r1./r2; plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');hold on;grid on;axis equal legend('控制多边形', 'B-样条曲线');else disp('输入错误，请核查u的维数是否满足n+k+1');endaxis equal%一阶导xt1=diff(r(:,1)',1);yt1=diff(r(:,2)',1);zt1=diff(r(:,3)'

8、,1);ft1=xt1;yt1;zt1'%2阶导xt2=diff(r(:,1)',2);yt2=diff(r(:,2)',2);zt2=diff(r(:,3)',2);ft2=xt2;yt2;zt2'%3阶导xt3=diff(r(:,1)',3);yt3=diff(r(:,2)',3);zt3=diff(r(:,3)',3);ft3=xt3;yt3;zt3'%曲率计算k=zeros(size(ft2,1),1);for i=1:size(ft2,1) n1=ft1(i,:); n2=ft2(i,:); ka=sqrt(su

9、m(cross(n1,n2).*cross(n1,n2); kb=sqrt(sum(n1.*n1); k(i)=ka/kb3;end%在点上表示曲率大小figure(2)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');hold on;grid on;axis equalfor j=1:2:M-3 quiver3(r(j,1),r(j,2),r(j,3),0,0,k(j),'b');endlegend('控制多边形', 'B-样条曲线',

10、9;曲率');%挠率计算t=zeros(size(ft3,1),1);for i=1:size(ft3,1) nn1=ft1(i,:); nn2=ft2(i,:); nn3=ft3(i,:); ta=det(nn1;nn2;nn3); t(i)=ta/ka2;end%在点上表示挠率大小figure(3)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');hold on;grid on;axis equalfor j=2:2:M-3 quiver3(r(j,1),r(j,2),r(j,3

11、),0,0,t(j),'b');endlegend('控制多边形', 'B-样条曲线','挠率');% %在点上表示其Frenetic标架figure(4)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');hold on;grid on;axis equalfor j=1:M-2 alpha=ft1(j,:)/norm(ft1(j,:); quiver3(r(j,1),r(j,2),r(j,3),alpha(1),alpha(2

12、),alpha(3),'b'); beta=ft2(j,:)/(k(j)*norm(ft2(j,:)2); beta=beta/norm(beta); quiver3(r(j,1),r(j,2),r(j,3),beta(1),beta(2),beta(3),'g'); gama=cross(alpha,beta)/norm(cross(alpha,beta); quiver3(r(j,1),r(j,2),r(j,3),gama(1),gama(2),gama(3),'m');endlegend('控制多边形', 'B-样条

13、曲线','alpha','beta','gama');l 非均匀有理B样条曲面（3阶）P1,，P2,，P3,，P4,，已经四个权因子w1，w2，w4，w4，通过Matlab计算其基函数，并作出曲线。Clear；clcfid=fopen('Data.dat','r');k=1;s=fscanf(fid,'%s',1);%读字符串 1维num=0; %NUM的行数Dotnum=zeros(1); %保存每个NUM的点的个数Dot=cell(1); %保存每个NUM的点的坐标while k=1 i

14、f strcmp(s,'NUM') a=fscanf(fid,'%d',2);%读整型数据 2维 num=num+1; k=1; else k=0; break; end b=fscanf(fid,'%f',3*a(2);%读浮点型数据 3*a(2)维 Dotnum(num)=a(2); b=reshape(b,3,a(2);%将b转化为3行 19列 Dotnum=b' % clear a clear b s = fscanf(fid,'%s',1);%读字符串 1维end%hold;view(3);u=0:0.08:1;

15、v=0:0.04:1;a1=(1-3*u+3*u.*u-u.*u.*u);a2=(3*u-6*u.*u+3*u.*u.*u);a3=(3*u.*u-3*u.*u.*u);a4=u.*u.*u;b1=(1-3*v+3*v.*v-v.*v.*v);b2=(3*v-6*v.*v+3*v.*v.*v);b3=(3*v.*v-3*v.*v.*v);b4=v.*v.*v;aa1=(-3+6*u-3*u.*u);aa2=(3-12*u+9*u.*u);aa3=(6*u-9*u.*u);aa4=3*u.*u;bb1=(-3+6*v-3*v.*v);bb2=(3-12*v+9*v.*v);bb3=(6*v-9*

16、v.*v);bb4=3*v.*v;aaa1=(6-6*u);aaa2=(-12+18*u);aaa3=(6-18*u);aaa4=6*u;bbb1=(6-6*v);bbb2=(-12+18*v);bbb3=(6-18*v);bbb4=6*v;for h=1:length(u) a(h,:)=a1(h) a2(h) a3(h) a4(h); aa(h,:)=aa1(h) aa2(h) aa3(h) aa4(h); aaa(h,:)=aaa1(h) aaa2(h) aaa3(h) aaa4(h);endfor k=1:length(v) b(k,:)=b1(k) b2(k) b3(k) b4(k)

17、; bb(k,:)=bb1(k) bb2(k) bb3(k) bb4(k); bbb(k,:)=bbb1(k) bbb2(k) bbb3(k) bbb4(k);end%for i=1:num-3for j=1:Dotnum(num)-3 xc=Doti(j,1),Doti(j+1,1),Doti(j+2,1),Doti(j+3,1);Doti+1(j,1),Doti+1(j+1,1),Doti+1(j+2,1),Doti+1(j+3,1); Doti+2(j,1),Doti+2(j+1,1),Doti+2(j+2,1),Doti+2(j+3,1);Doti+3(j,1),Doti+3(j+1,

18、1),Doti+3(j+2,1),Doti+3(j+3,1); yc=Doti(j,2),Doti(j+1,2),Doti(j+2,2),Doti(j+3,2);Doti+1(j,2),Doti+1(j+1,2),Doti+1(j+2,2),Doti+1(j+3,2); Doti+2(j,2),Doti+2(j+1,2),Doti+2(j+2,2),Doti+2(j+3,2);Doti+3(j,2),Doti+3(j+1,2),Doti+3(j+2,2),Doti+3(j+3,2); zc=Doti(j,3),Doti(j+1,3),Doti(j+2,3),Doti(j+3,3);Doti+1

19、(j,3),Doti+1(j+1,3),Doti+1(j+2,3),Doti+1(j+3,3); Doti+2(j,3),Doti+2(j+1,3),Doti+2(j+2,3),Doti+2(j+3,3);Doti+3(j,3),Doti+3(j+1,3),Doti+3(j+2,3),Doti+3(j+3,3); A=a*xc*b'%曲面上点的坐标赋值 B=a*yc*b' C=a*zc*b' rA=A;rB=B;rC=C; surf(A,B,C),shading interp;colormap(jet) %彩色显示 endendfor i=1:num-3 for j=1

20、:Dotnum(num)-3 xc=Doti(j,1),Doti(j+1,1),Doti(j+2,1),Doti(j+3,1);Doti+1(j,1),Doti+1(j+1,1),Doti+1(j+2,1),Doti+1(j+3,1); Doti+2(j,1),Doti+2(j+1,1),Doti+2(j+2,1),Doti+2(j+3,1);Doti+3(j,1),Doti+3(j+1,1),Doti+3(j+2,1),Doti+3(j+3,1); yc=Doti(j,2),Doti(j+1,2),Doti(j+2,2),Doti(j+3,2);Doti+1(j,2),Doti+1(j+1,

21、2),Doti+1(j+2,2),Doti+1(j+3,2); Doti+2(j,2),Doti+2(j+1,2),Doti+2(j+2,2),Doti+2(j+3,2);Doti+3(j,2),Doti+3(j+1,2),Doti+3(j+2,2),Doti+3(j+3,2); zc=Doti(j,3),Doti(j+1,3),Doti(j+2,3),Doti(j+3,3);Doti+1(j,3),Doti+1(j+1,3),Doti+1(j+2,3),Doti+1(j+3,3); Doti+2(j,3),Doti+2(j+1,3),Doti+2(j+2,3),Doti+2(j+3,3);D

22、oti+3(j,3),Doti+3(j+1,3),Doti+3(j+2,3),Doti+3(j+3,3); for h=1:length(u) for k=1:length(v) A=a(h,:)*xc*b(k,:)'%曲面上点的坐标赋值 B=a(h,:)*yc*b(k,:)'C=a(h,:)*zc*b(k,:)' Au=aa(h,:)*xc*b(k,:)' Av=a(h,:)*xc*bb(k,:)' Bu=aa(h,:)*yc*b(k,:)' Bv=a(h,:)*yc*bb(k,:)' Cu=aa(h,:)*zc*b(k,:)'

23、 Cv=a(h,:)*zc*bb(k,:)' Auu=aaa(h,:)*xc*b(k,:)' Avv=a(h,:)*xc*bbb(k,:)' Buu=aaa(h,:)*yc*b(k,:)' Bvv=a(h,:)*yc*bbb(k,:)' Cuu=aaa(h,:)*zc*b(k,:)' Cvv=a(h,:)*zc*bbb(k,:)' Auv=aa(h,:)*xc*bb(k,:)' Buv=aa(h,:)*yc*bb(k,:)' Cuv=aa(h,:)*zc*bb(k,:)' Sz=Au,Bu,Cu; Sf=Av,Bv

24、,Cv; Sff=Avv,Bvv,Cvv; Sfz=Auv,Buv,Cuv; Szz=Auu,Buu,Cuu; N=cross(Sz,Sf); n=N/norm(N,2); E(h,k)=dot(Sz,Sz);F(h,k)=dot(Sz,Sf);G(h,k)=dot(Sf,Sf);L(h,k)=dot(Szz,n);M(h,k)=dot(Sfz,n);N(h,k)=dot(Sff,n); f2(h,k)=L(h,k)*N(h,k)-M(h,k)*M(h,k); f1(h,k)=E(h,k)*G(h,k)-F(h,k)*F(h,k); K(h,k)=f2(h,k)/f1(h,k); %高斯曲率

25、H(h,k)=L(h,k)*G(h,k)-2*M(h,k)*F(h,k)+N(h,k)*E(h,k)/(f1(h,k);%平均曲率 k1(h,k)=H(h,k)+sqrt(H(h,k)2-K(h,k); k2(h,k)=H(h,k)-sqrt(H(h,k)2-K(h,k); %高斯曲率 bl=1000; if (H(h,k)=NaN) break; else AD(h,k)=A+bl*1.0/H(h,k)*n(:,1); BD(h,k)=B+bl*1.0/H(h,k)*n(:,2); CD(h,k)=C+bl*1.0/H(h,k)*n(:,3); line(A AD(h,k),B BD(h,k

26、),C CD(h,k),'Color','b'); end end end endendfigure(2);surf(rA,rB,rC),shading interp;colormap(jet) %彩色显示hold on;grid on;for i=1:size(rA,1) for j=1:size(rA,2) quiver3(rA(i,j),rB(i,j),rC(i,j),0,0,E(i,j),'b'); endendfigure(3);surf(rA,rB,rC),shading interp;colormap(jet) %彩色显示hold

27、on;grid on;for i=1:size(rA,1) for j=1:size(rA,2) quiver3(rA(i,j),rB(i,j),rC(i,j),0,0,F(i,j)/20,'b'); endendfigure(4);surf(rA,rB,rC),shading interp;colormap(jet) %彩色显示hold on;grid on;for i=1:size(rA,1) for j=1:size(rA,2) quiver3(rA(i,j),rB(i,j),rC(i,j),0,0,G(i,j)/20,'b'); endendfigure(5);surf(rA,rB,rC),shading interp;colormap(jet) %彩色显示hold on;grid on;for i=1:size(rA,1) for j=1:size(rA,2) quiver3(rA(i,j),rB(i,j),rC(i,j)