大連理工大學(xué)微分幾何作業(yè)學(xué) 院： 機(jī)械工程學(xué)院 姓 名： 學(xué) 號： 授課教師： 王曉明 電子郵箱： 完成時間：2012年6月25日非均勻有理B樣條（NURBS）曲線和曲面1. 非均勻有理B樣條曲線（3階）NURBS曲線提供了對標(biāo)準(zhǔn)解析幾何和自由曲線、曲面的統(tǒng)一數(shù)學(xué)描述方法，它可通過調(diào)整控制頂點和權(quán)因子，方便地改變曲面形狀，同時也可方便地轉(zhuǎn)換成對應(yīng)Bezier曲面。給定n+1個控制點及權(quán)因子，則k階（k-1）次NURBS曲線的表達(dá)式為：其中：為非均勻有理B樣條基函數(shù)，按照deBoor-Cox遞推公式，基函數(shù)滿足一下關(guān)系：其它1.1. 曲線的建立本文選取四個控制頂點P1,，P2,，P3,，P4,，已經(jīng)四個權(quán)因子w1，w2，w4，w4，通過Matlab計算其基函數(shù)，并作出曲線?？刂祈旤c：每個坐標(biāo)有三個值 立體的；權(quán)因子：；非均勻節(jié)點：；由Matlab繪制曲線如下：在每個點上標(biāo)注其Frenet標(biāo)架如圖：在每個點上標(biāo)注其曲率、撓率如圖： 2. 非均勻有理B樣條曲面（3階）假設(shè)給定R3空間中，齊次坐標(biāo)下的個點，則參數(shù)曲面：稱為矩形域0,1x0,1上的mxn次有理B樣條曲面。其中稱為控制頂點，稱為權(quán)因子。NURBS曲面的表達(dá)式為：2.1. 曲面計算的基本公式曲面第一基本形式：曲面第二基本形式：高斯曲率計算公式：平均曲率計算公式：曲面主曲率計算公式： ， 2.2. 曲面的建立選擇16個頂點構(gòu)造B樣條曲面，如表1所示：表1 16個頂點坐標(biāo)數(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é)果如圖： 圖1 B樣條曲面高斯曲率K 圖2 B樣條曲面平均曲率H 圖3 E值 圖4 F值 圖5 G值 圖6 L值 圖7 L值 圖8 N值 圖9 k1值 圖10 k2值3. 參考文獻(xiàn)1 王省富.樣條函數(shù)及其應(yīng)用M. 西北工業(yè)大學(xué)出版社，1989,9，西安.2 張化光，劉鑫蕊，孫秋野.MATLAB/SIMULINK實用教程M.人民郵電出版社，2009，3，北京.3 John H. Mathews, Kurtis D. Fink.數(shù)值方法（MATLAB版）.電子工業(yè)出版社，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; %插入點權(quán)重% w=2,2,2,2; %權(quán)重相等時u=1 2 4 5 6.5 7.5 8 9.3;%非均勻節(jié)點（構(gòu)造樣條基函數(shù)）k=3; %三階非均勻有理B樣條曲線n=size(P,1); %控制點數(shù)目nu=size(u,2); %非均勻節(jié)點數(shù)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 %構(gòu)造非均勻有理B樣條基函數(shù) 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的維數(shù)是否滿足n+k+1);endaxis equal%一階導(dǎo)xt1=diff(r(:,1),1);yt1=diff(r(:,2),1);zt1=diff(r(:,3),1);ft1=xt1;yt1;zt1;%2階導(dǎo)xt2=diff(r(:,1),2);yt2=diff(r(:,2),2);zt2=diff(r(:,3),2);ft2=xt2;yt2;zt2;%3階導(dǎo)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標(biāo)架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,，已經(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的點的坐標(biāo)while k=1 if strcmp(s,NUM) a=fscanf(fid,%d,2);%讀整型數(shù)據(jù) 2維 num=num+1; k=1; else k=0; break; end b=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列 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;%曲面上點的坐標(biāo)賦值 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,:);%曲面上點的坐標(biāo)賦值 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),