曲面加工几何学.doc

大连理工大学微分几何作业学 院： 机械工程学院 姓 名： 学 号： 授课教师： 王晓明 电子邮箱： 完成时间：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,，已经四个权因子w1，w2，w4，w4，通过Matlab计算其基函数，并作出曲线。控制顶点：每个坐标有三个值 立体的；权因子：；非均匀节点：；由Matlab绘制曲线如下：在每个点上标注其Frenet标架如图：在每个点上标注其曲率、挠率如图： 2. 非均匀有理B样条曲面（3阶）假设给定R3空间中，齐次坐标下的个点，则参数曲面：称为矩形域0,1x0,1上的mxn次有理B样条曲面。其中称为控制顶点，称为权因子。NURBS曲面的表达式为：2.1. 曲面计算的基本公式曲面第一基本形式：曲面第二基本形式：高斯曲率计算公式：平均曲率计算公式：曲面主曲率计算公式： ， 2.2. 曲面的建立选择16个顶点构造B样条曲面，如表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权因子：；计算结果如图： 图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;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); N=zeros(n+k,k+1); for i=1:n+k if x=u(i)&xu(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+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),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(sum(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-样条曲线,曲率);%挠率计算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),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),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-样条曲线,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 if 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;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*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); 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,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(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: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,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);Doti+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,:); 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,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); %高斯曲率 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),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 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),