求教关于曲面拟合程序，计算公式

5爱football

步骤和公式见附件，谢谢高手帮忙。

5爱football

5爱football

此程序是曲线的，好像有点错误，请大牛们指教，多谢！最好给个曲面重构的
%X:原始资料,d:控制顶点
%n:数据条数,k:B样条的次数
%
%
X=load('data.txt');
[n,numy]=size(X); %得数据维数;
k=3;
%弦长参数化
u(k+1)=0;
for i=1:(n-1)
    u(k+i+1)=u(k+i+1)=u(k+i)+abs((X(i+1,1)-X(i,1))^2+(X(i+1,2)-X(i,2))^2+(X(i+1,3)-X(i,3))^2)^(1/2);
end
%====规范参数化=========================
temp=u(n+k);
for i=(k):(n+k) %(k+1):(n-1+k)足够了
    u(i)=u(i)/temp;
end
%首末节点重复k+1个
for i=1:(k+1)
    u(i)=0;
    u(n-1+k+i)=1;
end
%% ===============================
%A:方程系数------采用递推算法-----------------------
A=zeros(n+2);
%控制多形与曲线相切,切曲线图端点为第1个控制点重合
%A(1,1)=1;A(1,2)=-2;A(1,3)=1;A(n+2,n)=1;A(n+2,n+1)=-2;A(n+2,n+2)=1;
%控制多形与曲线相切,切曲线图端点和第1,2个控制点重合
A(1,1)=1;A(1,2)=-1;A(n+2,n+1)=-1;A(n+2,n+2)=1;
A(2,2)=1;A(n+1,n+1)=1;
for i=3:n
    for j=0:2
        A(i,i+j-1)=Bbase(i+j-1,3,u,u(2+i));
    end
end
%e:方程右边.
e=0;
for i=1:numy
    e(n+2,i)=0;
end
for i=1:n
    e(i+1,:)=X(i,:);
end
%得到控制点,
d=inv(A)*e;
%clear A;
%clear e;
%===画出图形================================
hold on
%原始数据,红色,点
plot(X(:,1),X(:,2),'r.');
%控制多边形,蓝色,线
plot(d(:,1),d(:,2),'b');
%===插值B样条曲线============================
x=0;y=0;down=0;
for j=1:(n-1)
    uu=(u(j+3)):0.005:u(j+4);
    for kk=1:length(uu)
        temp=j+4;
        down=down+1;
        x(down)=d(j,1)*Bbase(temp-4,3,u,uu(kk))+d(j+1,1)*Bbase(temp-3,3,u,uu(kk))+d(j+2,1)*Bbase(temp-2,3,u,uu(kk))+d(j+3,1)*Bbase(temp-1,3,u,uu(kk));
        y(down)=d(j,2)*Bbase(temp-4,3,u,uu(kk))+d(j+1,2)*Bbase(temp-3,3,u,uu(kk))+d(j+2,2)*Bbase(temp-2,3,u,uu(kk))+d(j+3,2)*Bbase(temp-1,3,u,uu(kk));
    end
end
plot(x,y,'g');
%clear uu x y;
%===输入增加的节点并得到该点数据==============
fprintf('请输入一个%d-%d之间的',min(X(:,1)),max(X(:,1)) );
Knot=input('一个补充节点:');
%Knot=25;
if Knot<min(X(:,1)) || Knot>max(X(:,1))
    fprintf('你的输入有误,超过了范围.\n');
    return;
end
%转换为节点
knot=(Knot-min(X(:,1)))/(max(X(:,1))-min(X(:,1)));
kn=find(u>knot,1,'first');
newu=u;
newu(kn)=knot;
newu(kn+1:length(u)+1)=u(kn:length(u));
%===得到插入节点的数据======================
inknot=DeBoor(d,u,knot,kn-1,3);
plot(inknot(1),inknot(2),'m.');
%==生成新的控制序列dd,Deboor算法=============
dd=d;
dd(n+3,:)=0;
dd(kn:n+3,:)=dd(kn-1:n+2,:);
for i=kn-k:(kn-1)
    dd(i,:)=DeBoor(d,u,knot,i,1);
end
%新控制多边形,蓝色,线
plot(dd(:,1),dd(:,2),'b--');
%===画出修正后的图=========================
x=0;y=0;down=0;
for j=1:n
    uu=(newu(j+3)):0.005:newu(j+4);
    for kk=1:length(uu)
        temp=j+4;
        down=down+1;
        x(down)=dd(j,1)*Bbase(temp-4,3,newu,uu(kk))+dd(j+1,1)*Bbase(temp-3,3,newu,uu(kk))+dd(j+2,1)*Bbase(temp-2,3,newu,uu(kk))+dd(j+3,1)*Bbase(temp-1,3,newu,uu(kk));
        y(down)=dd(j,2)*Bbase(temp-4,3,newu,uu(kk))+dd(j+1,2)*Bbase(temp-3,3,newu,uu(kk))+dd(j+2,2)*Bbase(temp-2,3,newu,uu(kk))+dd(j+3,2)*Bbase(temp-1,3,newu,uu(kk));
    end
end
plot(x,y,'r-.');
hold off

%****************第i段k次B样条基,Deboor递推递归算法******************************
function result = Bbase(i,k,u,t)
%第i段k次B样条基,Deboor递推递归算法
%t为变量,u(i)<=t<u(i+1),k=0时result=1;
if k==0
    if (u(i)<=t && t<u(i+1)) %注意t=u(i+1)) 时的情况,要用t<=u(i+1);
        result=1;
        return;
    else
        result=0;
        return;
    end
end
    if u(i+k)-u(i)==0
        alpha=0;
    else
        alpha=(t-u(i))/(u(i+k)-u(i));
    end
    if u(i+k+1)-u(i+1)==0
         beta=0;
    else
        beta=(u(i+k+1)-t)/(u(i+k+1)-u(i+1));
    end
result=alpha*Bbase(i,k-1,u,t)+beta*Bbase(i+1,k-1,u,t);
end

%************递归,德布尔递推算法***************************************
function y=DeBoor(d,u,kn,i,L)
%控制点序列,u节点序列,kn插入的节点,u(i)<kn<u(i+1),L递推层次
% k B样条次数
k=3;
if L==0
     y=d(i,:);
     return;
end
temp=(kn-u(i))/(u(i+k+1-L)-u(i));
y=(1-temp)*DeBoor(d,u,kn,i-1,L-1)+temp*DeBoor(d,u,kn,i,L-1);