好例子网_非均匀B作业-反算.doc

非均匀B样条控制顶点反算

已知曲线数据如下，利用非均匀B样条方法反算控制顶点，再取△t=0.02时，计算该曲线上的点，一阶和二阶矢量，并绘制出点与一阶单位切矢量、点与二阶单位切矢量的图形，打印出计算程序清单。

34.450764.897992.7423

53.213362.910378.7443

72.617458.538964.6679

90.492253.872651.7908

107.626250.252740.0300

124.311548.538229.7099

142.917447.891420.1672

165.059346.123311.5884

187.763341.41425.4394

210.385132.38901.3837

232.500919.6438-0.7049

253.29474.2544-1.1379

程序：

clear

pt=load(data.txt);

m=length(pt);

n=m+1;

k=3;

U=zeros(1,n+k+2);%节点矢量

%********节点矢量************%

x=pt(:,1);

y=pt(:,2);

z=pt(:,3);

temp=zeros(1,m-1);

fori=1:m-1%m个数据点，利用积累弦长法

temp(i)=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2+(z(i+1)-z(i))^2);

end

fori=1:k+1%前k+1个节点为0

U(i)=0;

end

fori=n+k-1:n+k+2%后k+1个节点为1

U(i)=1;

end

fori=k+1:m+k-1%m个数据点，m=n-1,内节点为m-2个，U(k+1)作为初始值

U(i+1)=U(i)+temp(i-k);%m-1段弦长，U（k+1）到U（m+k）共m个节点

end

fori=k+1:m+k

U(i)=U(i)/U(m+k);

end

%***************反算n+1个控制点**************

dU=zeros(1,n+k+2);%△U

dpt1=[-0.64890.70790.2789];%给定首数据点切矢

dptm=[0.12270.7740-0.6212];%给定末数据点切矢

A=zeros(n-1);

E=zeros(n-1,3);

fori=k+1:n+k+1

dU(i)=U(i+1)-U(i);

end

A(1,1)=1;%切矢条件b1=1,c1=a1=0,e1

E(1,:)=pt(1,:)+(dU(4)/3)*dpt1;

A(n-1,n-1)=1;%切矢条件cn-1=an-1=0,bn-1=0

E(n-1,:)=pt(m,:)-(dU(n+1)/3)*dptm;

fori=2:n-2%a,b,c,e的值

A(i,i-1)=dU(i+3).^2/(dU(i+1)+dU(i+2)+dU(i+3));

A(i,i)=dU(i+3)*(dU(i+1)+dU(i+2))/(dU(i+1)+dU(i+2)+dU(i+3))+...

dU(i+2)*(dU(i+3)+dU(i+4))/(dU(i+2)+dU(i+3)+dU(i+4));

A(i,i+1)=dU(i+2).^2/(dU(i+2)+dU(i+3)+dU(i+4));

E(i,:)=(dU(i+2)+dU(i+3))*pt(i,:);

end

D=A\E;

D=[pt(1,:);D;pt(m,:)];%加上首末端点，控制顶点比数据点多两个

[s,t]=size(D);

%矩阵行数由1开头，故除k和l以外的参数都加1,n本身已加1

dt=0.02;

k=3;

P=[];

dP=[];

ddP=[];

symsx;

%定义节点间的递增变量为变量x

u_0(1:k+1)=0;

fort=k+2:s

u_0(t)=u_0(t-1)+1/(s-k);

end

u_0(s+1:s+k+1)=1;

%求出初始的节点值

%r=0,数据点：

fori=k+1:s

u=u_0;

d=sym(D);

%每次迭代后将d恢复初值

forl=1:k

forj=i-k:i-l

alpha(j)=(x-u(j+l))/(u(j+k