1、大連理工大學微分幾何作業學院:姓名:學號:授課教師:電子郵箱:完成時間:機械工程學院王曉明2012年6月25日1/18非均勻有理B樣條（NURBS曲線和曲面1,非均勻有理B樣條曲線（3階）NURB珈線提供了對標準解析幾何和自由曲線、曲面的統一數學描述方法,它可通過調整控制頂點和權因子,方便地改變曲面形狀，同時也可方便地轉換成對應Bezier曲面。給定n+1個控制點及權因子,則k階（k-1）次NURB釉線的表達式為:r(x)nNi,k(x)WRi0nNi,k(x)Wii0其中：Ni,k（u）為非均勻有理B樣條基函數，按照deBoor-Cox遞推公式，基函數滿足一下關系:M,k(x)01,當Uiu

2、Ui1其它M,k(x)(xUi)Ni,k1(x)(Uikx)Ni01(x)1.1.曲線的建立Uik1UiUikUi1本文選取四個控制頂點P1,P2,P3,P4,已經四個權因子W,W2,%,WA,通過Matlab計算其基函數，并作出曲線。P控制頂點：0.51；權因子：w20.513；1.52/18非均勻節點：u12456.57.589.3由Matlab繪制曲線如下:-控制多邊形A-B-樣條曲線在每個點上標注其Frenet標架如圖:控制多邊形*B-樣條曲線alphabetagama3/18在每個點上標注其曲率、撓率如圖:控制多邊形控制多邊形B-樣條曲線-撓率2.非均勻有理B樣條曲面(3階)假設名&

3、amp;定R3空間中，齊次坐標下的m1n1個點R(X。,Zj,Wj)(WjXij,WijYj，WjZij,Wij),i0,1,m;j=0,1,n,w。0,則參數曲面：mnR(u,v)(X(u,v),Y(u,v),Z(u,v),w(u,v)2Bm(u)B：(v),u0,v01 0j0稱為矩形域0,1x0,1上的mxn次有理B樣條曲面。其中巳稱為控制頂點，w。稱為權因子。NURBS®面的表達式為:nR(u,v)Wi,jR,jNi,k(u)Nj(v)j0nWi,jNi,k(u)N"(v)i01.1. 曲面計算的基本公式曲面第一基本形式:Frurv,Grvi;4/18曲面第二基本形

4、式:Lnruu,Mnruv,Nnrwo高斯曲率計算公式:2LNMEGF2平均曲率計算公式:NE2MFLG2EGF2曲面主曲率計算公式：kiH4HK,k2H4HK1.2. 曲面的建立選才?16個頂點構造B樣條曲面，如表1所示：表116個頂點坐標數據廳Pxyz1-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.967

5、2.967.9112-6.1581.5711.9513-16.3453.371.0414-12.8163.545.5715-7.6173.4110.0416-1.4581.7213.675/1820.513w權因子：0.310.8211.510.80.5010.51計算結果如圖:圖2B樣條曲面平均曲率H圖3E值圖4F值圖5G值圖6L值6/18圖9k1值圖10k2值3.參考文獻1王省富.樣條函數及其應用M.西北工業大學出版社，1989,9,西安.2張化光，劉鑫蕊，孫秋野.MATLAB/SIMULINK實用教程M.人民郵電出版社，2009,3,北京.3JohnH.Mathews,KurtisD.F

6、ink.數值方法（MATLAB版）.電子工業出版社，2002,6,北京.7/18附錄(Matlab源程序)：非均勻有理B樣條曲線(3階)functionNURBS3clear;clc;symsrP=010;000.5;101;111.5;w=2,0.5,1,3;%插入點權重%w=2,2,2,2;%權重相等時u=12456.57.589.3;%非均勻節點(構造樣條基函數)k=3;%三階非均勻有理B樣條曲線n=size(P,1);%控制點數目nu=size(u,2);%非均勻節點數ifnu=n+k+1M=100;r1=zeros(M,3);r2=zeros(M,3);r=zeros(M,3);X=

7、linspace(u(k+1),u(n+1),M);forl=1:M%構造非均勻有理B樣條基函數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+18/18N(i,j)=(x-u(i)/(u(i+卜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);

8、%分子分母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的維數是否滿足n+k+1');endaxisequal獷階導xt1=diff(r(:,1)',1);yt1=diff(r(:,2)',1

9、);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);9/18fori=1:size(ft2,1)n1=ft1(i,:);n

10、2=ft2(i,:);ka=sqrt(sum(cross(n1,n2).*cross(n1,n2);kb=sqrt(sum(n1.*n1);k(i)=ka/kbA3;end%fc點上表示曲率大小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-樣條曲

11、線','曲率)；%率計算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/kaA2;end%在點上表示撓率大小figure(3)plot3(P(:,1),P(:,2),P(:,3),'-ko',r(:,1),r(:,2),r(:,3),'-r.');holdon;gridon;axisequalforj=2:2:M-310/18quiver3(r(j,1),r(j,2),r(j,3)

12、,0,0,t(j),b');endlegend('控制多邊形','B-樣條曲線',撓率)；%fc點上表示其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');

13、beta=ft2(j,:)/(k(j)*norm(ft2(j,:)A2);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','

14、beta','gama');非均勻有理B樣條曲面(3階)P1,P2,P3,F4,已經四個權因子w,w，wi,W4,通過Matlab計算其基函數，并作出曲線。Clear；clcfid=fopen('Data.dat','r');k=1;s=fscanf(fid,'%s',1);%讀字符串1維num=0;%NUMB行數Dotnum=zeros(1);%保存每個NUM勺點的個數Dot=cell(1);%保存每個NUM勺點的坐標whilek=111/18ifstrcmp(s,'NUM')a=fscanf(fid,&

15、#39;%d',2);%讀整型數據2維num=num+1;k=1;elsek=0;break;endb=fscanf(fid,'%f,3*a(2);%讀浮點型數據3*a(2)維Dotnum(num)=a(2);b=reshape(b,3,a(2);%將b轉化為3行19列Dotnum=b'%clearaclearbend%hold;view(3);s=fscanf(fid,'%s',1);%讀字符串1維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=(

16、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*

17、u;bbb1=(6-6*v);bbb2=(-12+18*v);bbb3=(6-18*v);bbb4=6*v;12/18forh=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

18、=1:num-3forj=1:Dotnum(num)-3xc=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),

19、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);13/18Doti+2(j,3),Doti+2(j+1,3),Doti+2(j+2,3),Do

20、ti+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),shadinginterp;colormap(jet)%彩色顯示endendfori=1:num-3forj=1:Dotnum(num)-3xc=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)

21、,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

22、+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);forh=1:length(u)fork=1:length(v)14/18A=a(h,:)*xc*b(k,:)

23、'%曲面上點的坐標賦值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,

24、:)*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)=do

25、t(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)A2-K(h,k);k2(h,k)=H(h,k)-sqrt(H(h,k)A2-K(h,k);

26、%»%bl=1000;if(H(h,k)=NaN)break;15/18elseAD(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,

27、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;16/18fori=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;colormap(jet)%彩色顯示holdon;gridon;fori=1:size(rA,1)forj