kqx.m

采用德布尔算法进行B样条的生成

%%%%  n为控制点数减1，m为数据点数减1，u为反算的节点矢量，q为数据点
function d=kqx(n,u,q,m)
[I,J]=size(q);
if I==2
for i=2:n-2
    d(i:i+3)=u(i+2:i+5)-u(i+1:i+4);
    a(i)=d(i+2)^2/(d(i)+d(i+1)+d(i+2));
    b(i)=d(i+2)*(d(i)+d(i+1))/(d(i)+d(i+1)+d(i+2))+d(i+1)*(d(i+2)+d(i+3))/(d(i+1)+d(i+2)+d(i+3));
    c(i)=d(i+1)^2/(d(i+1)+d(i+2)+d(i+3));
    e1(i)=(d(i+1)+d(i+2))*q(1,i);
    e2(i)=(d(i+1)+d(i+2))*q(2,i);
end
%a(1)=1;b(1)=0;c(1)=0;
%e1(1)=d1(1);e2(1)=d2(1);
%a(n-1)=0;b(n-1)=0;c(n-1)=1;
%e1(n-1)=d1(n-1);e2(n-1)=d2(n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%自由端点条件
a(1)=2-d(3)*d(4)/(d(3)+d(4))^2;
b(1)=d(3)/(d(3)+d(4))*(d(4)/(d(3)+d(4))-d(3)/(d(3)+d(4)+d(5)));
c(1)=d(3)*d(3)/((d(3)+d(4))*(d(3)+d(4)+d(5)));
e1(1)=q(1,1)+q(1,2);
e2(1)=q(2,1)+q(2,2);
a(n-1)=-d(n)^2/(d(n-1)+d(n))/(d(n-2)+d(n-1)+d(n));
b(n-1)=d(n)/(d(n-1)+d(n))*(d(n)/(d(n-2)+d(n-1)+d(n))-d(n-1)/(d(n-1)+d(n)));
c(n-1)=d(n-1)*d(n)/(d(n-1)+d(n))^2-2;
e1(n-1)=-q(1,m+1)-q(1,m);
e2(n-1)=-q(2,m+1)-q(2,m);

A(1,1)=a(1);A(1,2)=b(1);A(1,3)=c(1);
j=1;
for i=1:n-3
    j=j+1;
    A(j,i)=a(j);
    A(j,i+1)=b(j);
    A(j,i+2)=c(j);
    
end
A(n-1,n-3)=a(n-1);
A(n-1,n-2)=b(n-1);
A(n-1,n-1)=c(n-1);

ep=0.001;
[d1,k]=Jacobi(A,e1,ep);
d1=[q(1,1);d1;q(1,m+1)];
[d2,k]=Jacobi(A,e2,ep);
d2=[q(2,1);d2;q(2,m+1)];

d=[d1';d2'];
%for k=1:n+1
%    plot(d(1,k),d(2,k),'kp');
%end
%plot(d(1,:),d(2,:),'r');
else 
  for i=1:n-2
    d(i:i+3)=u(i+2:i+5)-u(i+1:i+4);
    a(i)=d(i+2)^2/(d(i)+d(i+1)+d(i+2));
    b(i)=d(i+2)*(d(i)+d(i+1))/(d(i)+d(i+1)+d(i+2))+d(i+1)*(d(i+2)+d(i+3))/(d(i+1)+d(i+2)+d(i+3));
    c(i)=d(i+1)^2/(d(i+1)+d(i+2)+d(i+3));
    e1(i)=(d(i+1)+d(i+2))*q(1,i);
    e2(i)=(d(i+1)+d(i+2))*q(2,i);
    e3(i)=(d(i+1)+d(i+2))*q(3,i);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%自由端点条件
a(1)=2-d(3)*d(4)/(d(3)+d(4))^2;
b(1)=d(3)/(d(3)+d(4))*(d(4)/(d(3)+d(4))-d(3)/(d(3)+d(4)+d(5)));
c(1)=d(3)*d(3)/((d(3)+d(4))*(d(3)+d(4)+d(5)));
e1(1)=q(1,1)+q(1,2);
e2(1)=q(2,1)+q(2,2);
e3(1)=q(3,1)+q(3,2);
a(n-1)=-d(n)^2/(d(n-1)+d(n))/(d(n-2)+d(n-1)+d(n));
b(n-1)=d(n)/(d(n-1)+d(n))*(d(n)/(d(n-2)+d(n-1)+d(n))-d(n-1)/(d(n-1)+d(n)));
c(n-1)=d(n-1)*d(n)/(d(n-1)+d(n))^2-2;
e1(n-1)=-q(1,m+1)-q(1,m);
e2(n-1)=-q(2,m+1)-q(2,m);
e3(n-1)=-q(3,m+1)-q(3,m);

A(1,1)=a(1);A(1,2)=b(1);A(1,3)=c(1);
j=1;
for i=1:n-3
    j=j+1;
    A(j,i)=a(j);
    A(j,i+1)=b(j);
    A(j,i+2)=c(j);
    
end
A(n-1,n-3)=a(n-1);
A(n-1,n-2)=b(n-1);
A(n-1,n-1)=c(n-1);

ep=0.001;
[d1,k]=Jacobi(A,e1,ep);
d1=[q(1,1);d1;q(1,m+1)];
[d2,k]=Jacobi(A,e2,ep);
d2=[q(2,1);d2;q(2,m+1)];
[d3,k]=Jacobi(A,e3,ep);
d3=[q(3,1);d3;q(3,m+1)];
d=[d1';d2';d3'];
end