三次非均匀B样条曲面的绘制

%%%%%%%%%%三次非均匀B样条曲面的绘制

clc

clf

clear

A=input('请输入控制点：')%16个控制顶点要求以列向量表示点的坐标[2 3 4 5 2 3 4 5 2 3 4 5 2 3 4 5;2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5;1 2 3 1 3 4 6 2 7 5 6 2 9 1 3 4];

n=input('请输入所需绘制B样条曲线的次数：');

Vx=[A(1,1:4);A(1,5:8);A(1,9:12);A(1,13:16)];%将控制顶点的x坐标放到维数为4的方阵中

V y=[A(2,1:4);A(2,5:8);A(2,9:12);A(2,13:16)];%将控制顶点的y坐标放到维数为4的方阵中Vz=[A(3,1:4);A(3,5:8);A(3,9:12);A(3,13:16)];%将控制顶点的z坐标放到维数为4的方阵中plot3(A(1,:),A(2,:),A(3,:),'o');%输出控制顶点

hold on

syms u

syms w

P_uwx=[0];

P_uwy=[0];

P_uwz=[0];

for i=1:1:4;

for j=1:1:4

u0=0;u1=0;u2=0;u3=0;u4=1;u5=1;u6=1;u7=1;%给定4个控制顶点所取的节点矢量

B1=(u4-u3-u)^3/((u4-u3)*(u4-u2)*(u4-u1));

B2=((u3-u1+u)*(u4-u3-u)^2)/((u4-u3)*(u4-u2)*(u4-u1))+((u5-u3-u)*(u4-u3-u)*(u3-u2+u))/((u5-u 2)*(u4-u3)*(u4-u2))+((u5-u3-u)^2*u)/((u5-u3)*(u5-u2)*(u4-u2));

B3=((u4-u3-u)*(u3-u2+u)^2)/((u5-u2)*(u4-u3)*(u4-u2))+((u5-u3-u)*u*(u3-u2+u))/((u5-u3)*(u5-u2)*(u4-u3))+((u6-u3-u)*u^2)/((u6-u3)*(u5-u3)*(u4-u3));

B4=u^3/((u6-u3)*(u5-u3)*(u4-u3));

B=[B1;B2;B3;B4];%三次非均匀B样条基函数矩阵（自由曲线曲面造型技术139页）

w0=0;w1=0;w2=0;w3=0;w4=1;w5=1;w6=1;w7=1;%给定4个控制顶点所取的节点矢量

C1=(w4-w3-w)^3/((w4-w3)*(w4-w2)*(w4-w1));

C2=((w3-w1+w)*(w4-w3-w)^2)/((w4-w3)*(w4-w2)*(w4-w1))+((w5-w3-w)*(w4-w3-w)*(w3-w2 +w))/((w5-w2)*(w4-w3)*(w4-w2))+((w5-w3-w)^2*w)/((w5-w3)*(w5-w2)*(w4-w2));

C3=((w4-w3-w)*(w3-w2+w)^2)/((w5-w2)*(w4-w3)*(w4-w2))+((w5-w3-w)*w*(w3-w2+w))/((w 5-w3)*(w5-w2)*(w4-w3))+((w6-w3-w)*w^2)/((w6-w3)*(w5-w3)*(w4-w3));

C4=w^3/((w6-w3)*(w5-w3)*(w4-w3));

C=[C1;C2;C3;C4];%三次非均匀B样条基函数矩阵

Puwx=Vx(i,j)* B(i)*C(j);%控制顶点的x坐标与两个方向的基函数做乘积

Puwy=V y(i,j)* B(i)*C(j);

Puwz=Vz(i,j)* B(i)*C(j);

P_uwx=P_uwx+Puwx;%得到曲面的x坐表关于u,w的表达式（自由曲线曲面造型技术146页）

P_uwy=P_uwy+Puwy;%得到曲面的y坐表关于u,w的表达式

P_uwz=P_uwz+Puwz;%得到曲面的z坐表关于u,w的表达式

end

end

%%%%%%%%%%%%%%%%%%%%%%%%曲面的绘制

X1=zeros(11);Y1=zeros(11);Z1=zeros(11);%形成11维0矩阵

i=1;j=1;

for m=0:1/10:1;

for n=0:1/10:1;

X=subs(P_uwx,{u,w},{m,n});%P=subs(P,{u,w},{m,n})就是把P表达式中所有u,w都用具体的m,n值代替

Y=subs(P_uwy,{u,w},{m,n});

Z=subs(P_uwz,{u,w},{m,n});

X1(i,j)=X;%将X的值放入11维的方阵

Y1(i,j)=Y;

Z1(i,j)=Z;

j=j+1;

end

j=1;

i=i+1;

end

grid;

surf(X1,Y1,Z1);

hold on

输出结果：

图1.三次非均匀曲面图