nurbs曲线

[w47553112]

运行该程序时会报错：this statement is incomplete
不知道什么地方有问题 麻烦哪位高手帮忙看看 如果有绘制NURBS曲线的程序能否借鉴参考下（最好用泰勒级数展开的方法进行插值的）

而且程序运行不会出现三维NURBS曲线

clear;
U=[0 0 0 0 0.354 0.5 0.646 1 1 1 1]
X=[0 -150 -150 0 150 150 0]
Y=[0 -150 150 0 -150 150 0]
Z=[0 0 0 0 0 0 0]
W=[1 25 25 1 25 25 1]
%U=input('input U:') %输入节点向量
%X=input('please input X:') %输入控制顶点X坐标
%Y=input('please input Y:') %输入控制顶点Y坐标
%Z=input('please input Z:') %输入控制顶点Z坐标
%W=input('please input W:') %输入权因子
%Uz=U;%读取NURBS参数:节点矢量，控制顶点，权因子
%Xz=X;
%Yz=Y;
%Zz=Z;
%Wz=W;
%T=0.001;hm=0.001;am=4900;F=18*1000/60;A=2400;tag=1;
%J=48000;fs=0;T1=A/J;T3=T1;f1=0.5*J*T1*T1;T2=(F-fs-2*f1)/A;f2=f1+A*T2;all=T1/T+T2/T+T3/T+1;e=0;

syms b %b就为参数u
n=length(W); %n个控制顶点



u(1)=0;
dfx=(b-X(1))^3./(X(4)-X(1))./(X(3)-X(1))./(X(4)-X(1));
dfy=(b-Y(1))^3./(Y(4)-Y(1))./(Y(3)-Y(1))./(Y(4)-Y(1));
dfz=(b-Z(1))^3./(Z(4)-Z(1))./(Z(3)-Z(1))./(Z(4)-Z(1));
ddx=diff(dfx,'b');
ddy=diff(dfy,'b');
ddz=diff(dfz,'b');
b=u(1);
cx=eval(ddx);
cy=eval(ddy);
cz=eval(ddz);
u(2)=u(1)+3./sqrt((cx^2)+(cy^2)+(cz^2));
b=u(2);
j=3;

%求下一个插补点



for i=1:n
WX(i)=W(i).*X(i);
WY(i)=W(i).*Y(i);
WZ(i)=W(i).*Z(i);
end

while (b<1)
%求B样条基函数Ni,k(b)
for i=1:n

if (b>=U(i))&&(b<U(i+1))
N3(i)=(b-U(i))^3/(U(i+3)-U(i))/(U(i+2)-U(i))/(U(i+1)-U(i)); %N3(i)代表N i,3
N2(i)=(b-U(i))^2/(U(i+2)-U(i))/(U(i+1)-U(i));
N1(i)=(b-U(i))/(U(i+1)-U(i));

else if (b>=U(i+1))&&(b<U(i+2))
N3(i)=(b-U(i))^2*(U(i+2)-b)/(U(i+3)-U(i))/(U(i+2)-U(i))/(U(i+2)-U(i+1))+(b-U(i))*(U(i+3)-b)*(b-U(i+1))/(U(i+3)-U(i))/(U(i+3)-U(i+1))/(U(i+2)-U(i+1))+(b-U(i+1))^2*(U(i+4)-b)/(U(i+4)-U(i+1))/(U(i+3)-U(i+1))/(U(i+2)-U(i+1));
N2(i)=(b-U(i))*(U(i+2)-b)/(U(i+2)-U(i))/(U(i+2)-U(i+1))+(b-U(i+1))*(U(i+3)-b)/(U(i+3)-U(i+1))/(U(i+2)-U(i+1));
N1(i)=(U(i+2)-b)/(U(i+2)-U(i+1));

else if (b>=U(i+2))&&(b<U(i+3))
N3(i)=(U(i+3)-b)^2*(b-U(i))/(U(i+3)-U(i))/(U(i+3)-U(i+1))/(U(i+3)-U(i+2))+(b-U(i+1))*(U(i+3)-b)*(U(i+4)-b)/(U(i+3)-U(i+1))/(U(i+3)-U(i+2))/(U(i+4)-U(i+1))+(U(i+4)-b)^2*(b-U(i+2))/(U(i+4)-U(i+1))/(U(i+4)-U(i+2))/(U(i+3)-U(i+2));
N2(i)=(U(i+3)-b)^2/(U(i+3)-U(i+1))/(U(i+3)-U(i+2));
N1(i)=0;

else if (b>=U(i+3))&&(b<U(i+4))
N3(i)=(U(i+4)-b)^3/(U(i+4)-U(i+1))/(U(i+4)-U(i+2))/(U(i+4)-U(i+3));
N2(i)=0;
N1(i)=0;

else
N3(i)=0;
N2(i)=0;
N1(i)=0;

end
end

x=(WX*N3')./(W*N3');
y=(WY*N3')./(W*N3');
z=(WZ*N3')./(W*N3');

for i=1:n
M3(i)=3*(N2(i)./(U(i+3)-U(i))-N2(i+1)./(U(i+4)-U(i+1)));%M3(i)是N3(i)对b的导数
end

dX=((W*N3')*(M3*WX')-(W*M3')*(N3*WX'))/(W*N3')^2; %dX是x对b的一阶导数
dY=((W*N3')*(M3*WY')-(W*M3')*(N3*WY'))/(W*N3')^2;
dZ=((W*N3')*(M3*WZ')-(W*M3')*(N3*WZ'))/(W*N3')^2;


u(j)=b+3./sqrt(dX^2+dY^2+dZ^2);%3处应该是L
b=u(j);
j=j+1;

plot3(x,y,z)
hold on
clear x y z
end