曲面加工幾何學

大連理工大學微分幾何作業(yè)學院：機械工程學院姓名：學號：授課教師：王曉明電子郵箱：完成時間：2012年6月25日非均勻有理B樣條〔NURBS〕曲線和曲面非均勻有理B樣條曲線〔3階〕NURBS曲線提供了對標準解析幾何和自由曲線、曲面的統(tǒng)一數(shù)學描述方法，它可通過調(diào)整控制頂點和權(quán)因子，方便地改變曲面形狀，同時也可方便地轉(zhuǎn)換成對應Bezier曲面。給定n+1個控制點及權(quán)因子，那么k階〔k-1〕次NURBS曲線的表達式為：其中：為非均勻有理B樣條基函數(shù)，按照deBoor-Cox遞推公式，基函數(shù)滿足一下關(guān)系：其它曲線的建立本文選取四個控制頂點P1,，P2,，P3,，P4,，已經(jīng)四個權(quán)因子w1，w2，w4，w4，通過Matlab計算其基函數(shù)，并作出曲線?？刂祈旤c：每個坐標有三個值立體的；權(quán)因子：；每個坐標有三個值立體的非均勻節(jié)點：；由Matlab繪制曲線如下：在每個點上標注其Frenet標架如圖：在每個點上標注其曲率、撓率如圖：：非均勻有理B樣條曲面〔3階〕假設(shè)給定R3空間中，齊次坐標下的個點，那么參數(shù)曲面：稱為矩形域[0,1]x[0,1]上的mxn次有理B樣條曲面。其中稱為控制頂點，稱為權(quán)因子。NURBS曲面的表達式為：曲面計算的根本公式曲面第一根本形式：曲面第二根本形式：高斯曲率計算公式：平均曲率計算公式：曲面主曲率計算公式：，曲面的建立選擇16個頂點構(gòu)造B樣條曲面，如表1所示：表116個頂點坐標數(shù)據(jù)序號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權(quán)因子：；計算結(jié)果如圖：圖1B樣條曲面高斯曲率K圖2B樣條曲面平均曲率H圖3E值圖4F值圖5G值圖6L值圖7L值圖8N值圖9k1值圖10k2值參考文獻[1]王省富.樣條函數(shù)及其應用[M].西北工業(yè)大學出版社，1989,9，西安.[2]張化光，劉鑫蕊，孫秋野.MATLAB/SIMULINK實用教程[M].人民郵電出版社，2009，3，北京.[3]JohnH.Mathews,KurtisD.Fink.數(shù)值方法〔MATLAB版〕.電子工業(yè)出版社，2002,6，北京.附錄〔Matlab源程序〕：非均勻有理B樣條曲線〔3階〕functionNURBS3clear;clc;symsrP=[010;000.5;101;111.5];w=[2,0.5,1,3];%插入點權(quán)重%w=[2,2,2,2];%權(quán)重相等時u=[12456.57.589.3];%非均勻節(jié)點〔構(gòu)造樣條基函數(shù)〕k=3;%三階非均勻有理B樣條曲線n=size(P,1);%控制點數(shù)目nu=size(u,2);%非均勻節(jié)點數(shù)ifnu==n+k+1M=100;r1=zeros(M,3);r2=zeros(M,3);r=zeros(M,3);X=linspace(u(k+1),u(n+1),M);forl=1:M%構(gòu)造非均勻有理B樣條基函數(shù)x=X(l);N=zeros(n+k,k+1);fori=1:n+kifx>=u(i)&&x<u(i+1)N(i,1)=1;endendforj=2:k+1fori=1:n+k-j+1N(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);endendform1=1:nform2=1:3r1(l,m2)=r1(l,m2)+w(m1)*N(m1,k+1)*P(m1,m2);%分子endr2(l,:)=r2(l,:)+w(m1)*N(m1,k+1);%分母endendr=r1./r2;plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');holdon;gridon;axisequallegend('控制多邊形','B-樣條曲線');elsedisp('輸入錯誤，請核查u的維數(shù)是否滿足n+k+1');endaxisequal%一階導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);fori=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/kb^3;end%在點上表示曲率大小figure(2)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');holdon;gridon;axisequalforj=1:2:M-3quiver3(r(j,1),r(j,2),r(j,3),0,0,k(j),'b');endlegend('控制多邊形','B-樣條曲線','曲率');%%%撓率計算t=zeros(size(ft3,1),1);fori=1:size(ft3,1)nn1=ft1(i,:);nn2=ft2(i,:);nn3=ft3(i,:);ta=det([nn1;nn2;nn3]);t(i)=ta/ka^2;end%在點上表示撓率大小figure(3)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');holdon;gridon;axisequalforj=2:2:M-3quiver3(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.');holdon;gridon;axisequalforj=1:M-2alpha=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');非均勻有理B樣條曲面〔3階〕P1,，P2,，P3,，P4,，已經(jīng)四個權(quán)因子w1，w2，w4，w4，通過Matlab計算其基函數(shù)，并作出曲線。Clear；clcfid=fopen('Data.dat','r');k=1;s=fscanf(fid,'%s',1);%讀字符串1維num=0;%NUM的行數(shù)Dotnum=zeros(1);%保存每個NUM的點的個數(shù)Dot=cell(1);%保存每個NUM的點的坐標whilek==1ifstrcmp(s,'NUM')a=fscanf(fid,'%d',2);%讀整型數(shù)據(jù)2維num=num+1;k=1;elsek=0;break;endb=fscanf(fid,'%f',3*a(2));%讀浮點型數(shù)據(jù)3*a(2)維Dotnum(num)=a(2);b=reshape(b,3,a(2));%將b轉(zhuǎn)化為3行19列Dot{num}=b';%clearaclearbs=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;forh=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)];endfork=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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fori=1:num-3forj=1:Dotnum(num)-3xc=[Dot{i}(j,1),Dot{i}(j+1,1),Dot{i}(j+2,1),Dot{i}(j+3,1);Dot{i+1}(j,1),Dot{i+1}(j+1,1),Dot{i+1}(j+2,1),Dot{i+1}(j+3,1);Dot{i+2}(j,1),Dot{i+2}(j+1,1),Dot{i+2}(j+2,1),Dot{i+2}(j+3,1);Dot{i+3}(j,1),Dot{i+3}(j+1,1),Dot{i+3}(j+2,1),Dot{i+3}(j+3,1);];yc=[Dot{i}(j,2),Dot{i}(j+1,2),Dot{i}(j+2,2),Dot{i}(j+3,2);Dot{i+1}(j,2),Dot{i+1}(j+1,2),Dot{i+1}(j+2,2),Dot{i+1}(j+3,2);Dot{i+2}(j,2),Dot{i+2}(j+1,2),Dot{i+2}(j+2,2),Dot{i+2}(j+3,2);Dot{i+3}(j,2),Dot{i+3}(j+1,2),Dot{i+3}(j+2,2),Dot{i+3}(j+3,2);];zc=[Dot{i}(j,3),Dot{i}(j+1,3),Dot{i}(j+2,3),Dot{i}(j+3,3);Dot{i+1}(j,3),Dot{i+1}(j+1,3),Dot{i+1}(j+2,3),Dot{i+1}(j+3,3);Dot{i+2}(j,3),Dot{i+2}(j+1,3),Dot{i+2}(j+2,3),Dot{i+2}(j+3,3);Dot{i+3}(j,3),Dot{i+3}(j+1,3),Dot{i+3}(j+2,3),Dot{i+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),shadinginterp;colormap(jet)%彩色顯示endendfori=1:num-3forj=1:Dotnum(num)-3xc=[Dot{i}(j,1),Dot{i}(j+1,1),Dot{i}(j+2,1),Dot{i}(j+3,1);Dot{i+1}(j,1),Dot{i+1}(j+1,1),Dot{i+1}(j+2,1),Dot{i+1}(j+3,1);Dot{i+2}(j,1),Dot{i+2}(j+1,1),Dot{i+2}(j+2,1),Dot{i+2}(j+3,1);Dot{i+3}(j,1),Dot{i+3}(j+1,1),Dot{i+3}(j+2,1),Dot{i+3}(j+3,1);];yc=[Dot{i}(j,2),Dot{i}(j+1,2),Dot{i}(j+2,2),Dot{i}(j+3,2);Dot{i+1}(j,2),Dot{i+1}(j+1,2),Dot{i+1}(j+2,2),Dot{i+1}(j+3,2);Dot{i+2}(j,2),Dot{i+2}(j+1,2),Dot{i+2}(j+2,2),Dot{i+2}(j+3,2);Dot{i+3}(j,2),Dot{i+3}(j+1,2),Dot{i+3}(j+2,2),Dot{i+3}(j+3,2);];zc=[Dot{i}(j,3),Dot{i}(j+1,3),Dot{i}(j+2,3),Dot{i}(j+3,3);Dot{i+1}(j,3),Dot{i+1}(j+1,3),Dot{i+1}(j+2,3),Dot{i+1}(j+3,3);Dot{i+2}(j,3),Dot{i+2}(j+1,3),Dot{i+2}(j+2,3),Dot{i+2}(j+3,3);Dot{i+3}(j,3),Dot{i+3}(j+1,3),Dot{i+3}(j+2,3),Dot{i+3}(j+3,3);];forh=1:length(u)fork=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;elseAD(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([AAD(h,k)],[BBD(h,k)],[CCD(h,k)],'Color','b');endendendendendfigure(2);surf(rA,rB,rC),shadinginterp;colormap(jet)%彩色顯示holdon;gridon;fori=1:size(rA,1)forj=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),shadinginterp;colormap(jet)%彩色顯示holdon;gridon;fori=1:size(rA,1)forj=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),shadinginterp;colormap(jet)%彩色顯示holdon;gridon;fori=1:size(rA,1)forj=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),shadinginterp;color