曲面加工几何学.docx

曲面加工几何学.docx

《曲面加工几何学.docx》由会员分享，可在线阅读，更多相关《曲面加工几何学.docx（21页珍藏版）》请在冰豆网上搜索。

曲面加工几何学

大连理工大学

微分几何作业

学院：

机械工程学院

姓名：

学号：

授课教师：

王晓明

电子邮箱：

完成时间：

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,1]x[0,1]上的mxn次有理B样条曲面。

其中

称为控制顶点，

称为权因子。

NURBS曲面的表达式为：

2.1.曲面计算的基本公式

曲面第一基本形式：

曲面第二基本形式：

高斯曲率计算公式：

平均曲率计算公式：

曲面主曲率计算公式：

，

2.2.曲面的建立

选择16个顶点构造B样条曲面，如表1所示：

表116个顶点坐标数据

序号

x

y

z

1

-24.88

49.76

-4.35

2

-23.6

60.56

0.15

3

-19.47

71.91

5.17

4

-11.99

81.09

9.78

5

-22.89

50.75

-3.01

6

-21.03

61.45

1.44

7

-16.65

72.44

6.34

8

-9.47

81.35

10.72

9

-20.25

51.91

-1.29

10

-17.66

62.43

3.18

11

-12.96

72.96

7.91

12

-6.15

81.57

11.95

13

-16.34

53.37

1.04

14

-12.81

63.54

5.57

15

-7.61

73.41

10.04

16

-1.45

81.72

13.67

权因子：

；

计算结果如图：

图1B样条曲面高斯曲率K图2B样条曲面平均曲率H

图3E值图4F值

图5G值图6L值

图7L值图8N值

图9k1值图10k2值

3.参考文献

[1]王省富.样条函数和其应用[M].西北工业大学出版社，1989,9，西安.

[2]张化光，刘鑫蕊，孙秋野.MATLAB/SIMULINK实用教程[M].人民邮电出版社，2009，3，北京.

[3]JohnH.Mathews,KurtisD.Fink.数值方法（MATLAB版）.电子工业出版社，2002,6，北京.

附录（Matlab源程序）：

●非均匀有理B样条曲线（3阶）

functionNURBS3

clear;clc;symsr

P=[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+1

M=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%构造非均匀有理B样条基函数

x=X（l）;

N=zeros（n+k,k+1）;

fori=1:

n+k

ifx>=u（i）&&x

N（i,1）=1;

end

end

forj=2:

k+1

fori=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

form1=1:

n

form2=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.'）;holdon;gridon;axisequal

legend（'控制多边形','B-样条曲线'）;

elsedisp（'输入错误，请核查u的维数是否满足n+k+1'）;

end

axisequal

%一阶导

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;axisequal

forj=1:

2:

M-3

quiver3（r（j,1）,r（j,2）,r（j,3）,0,0,k（j）,'b'）;

end

legend（'控制多边形','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;axisequal

forj=2:

2:

M-3

quiver3（r（j,1）,r（j,2）,r（j,3）,0,0,t（j）,'b'）;

end

legend（'控制多边形','B-样条曲线','挠率'）;

%%在点上表示其Frenetic标架

figure（4）

plot3（P（:

1）,P（:

2）,P（:

3）,'-ko',r（:

1）,r（:

2）,r（:

3）,'-r.'）;holdon;gridon;axisequal

forj=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'）;

end

legend（'控制多边形','B-样条曲线','alpha','beta','gama'）;

●非均匀有理B样条曲面（3阶）

P1,，P2,，P3,，P4,，已经四个权因子w1，w2，w4，w4，通过Matlab计算其基函数，并作出曲线。

Clear；clc

fid=fopen（'Data.dat','r'）;

k=1;

s=fscanf（fid,'%s',1）;%读字符串1维

num=0;%NUM的行数

Dotnum=zeros

（1）;%保存每个NUM的点的个数

Dot=cell

（1）;%保存每个NUM的点的坐标

whilek==1

ifstrcmp（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列

Dot{num}=b';%

cleara

clearb

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;

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）];

end

fork=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-3

forj=1:

Dotnum（num）-3

xc=[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）%彩色显示

end

end

fori=1:

num-3

forj=1:

Dotnum（num）-3

xc=[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;

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（[AAD（h,k）],[BBD（h,k）],[CCD（h,k）],'Color','b'）;

end

end

end

end

end

figure

（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'）;

end

end

figure（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'）;

end

end

figure（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'）;

end

end

figure（5）;

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,L（i,j）,'b'）;

end

end

figure（6）;

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,M（i,j）,'b'）;

end

end

figure（7）;

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,N（i,j）,'b'）;

end

end

figure（8）;

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,k1（i,j）,'b'）;

end

end

figure（9）;

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,k2（i,j）,'m'）;

end

end

figure（10）;

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,K（i,j）,'m'）;

end

end