三次插值ClampedB样条曲线Matlab代码实现

在代码实现过程中，与原博客中的矩阵计算存在误差，主要体现在边缘值对控制点的系数求解上。系数求解为0、1/6、4/6、1/6和0…,代码和自己的手算解为0，1/4，7/12，1/6…

计算复杂，如有误望指正。

结果如下：

代码如下：

clc;clear;close; %输入值 l 1个 points=[4，4]；[5，6]；[7，7].[9，9]； plot(points(:,1),points(，2)or','MarkerSize',10); hold on 样条曲线%三次 k=m-n-1=3 B样条曲线的时间节点比型值点多6个 m 1=l 1 6 tArr=getTimeArray(points);   控制节点按照准均匀节点和插值节点计算 n 1=m-3=l 1 2个 controlPoint=getControlPoint(tArr,points);    插值样条曲线根据控制点和时间数组绘制 for times=0:0.002:1     b=getB(times,tArr);     p=(b.')*controlPoint;     plot(p(1,1),p(1,2),'.b');     hold on  end  %获得准均匀时间节点数组 function tArr=getTimeArray(points)     num=size(points,1);     tArr=zeros(1,num 6);     for i=1:1:num 6         if i<=4             tArr(1,i)=0;         end         if i>=num 3             tArr(1,i)=1;         end         if i>4 && i<num 3             tArr(1,i)=1/(num-1)*(i-4);         end     end  end  %%%根据时间节点和型值点获得控制点 l 1个 %根据非纽结边界条件 第一点与第二点的三次求导相同 最后一点与倒数第二点三次求导相同 function controlPoint=getControlPoint(tArr,points)     %型节点=系数矩阵*控制点     %型节点的数量     len=size(points,1);     %控制点     controlPoint=zeros(len 2,2);     %系数矩阵     bMatrix=zeros(len 2,len 2);     for i=1:1:len 2         %t=0三阶导数表达式         if i==1             t1=tArr(1,4);             t2=tArr(1,5);             vec1=getB3(t1,tArr).';             vec2=getB3(t2,tArr).';             bMatrix(i,:)=vec1-vec2;         end         if i>1&&i<len 2             t=tArr(1,i 2);             vec=getB(t,tArr).';             bMatrix(i,:)=vec;         end         if i==len 2             t1=tArr(1,len 1);             t2=tArr(1,len 2);             vec1=getB3(t1,tArr).';             vec2=getB3(t2,tArr).';             bMatrix(i,:)=vec1-vec2;         end     end     z=zeros(len 2,2);     for i=1:1:len 2         if i==1 ||i==len 2             z(i,1)=0;             z(i,2)=0;         end         if i>1&&i<len 2         z(i,:)=points(i-1,:);         end     end     bm=bMatrix^(-1);     controlPoint=bm*z; end  %%获取t时刻控制点系数 B0,n 1(t)-Bn 1,n 1(t) function bt=getB(t,tArr)     m=size(tArr,2)-1;     bt=zeros(m-3,1);     %t=1     if t==1         bt(m-3,1)=1;         return     end     if t==0         bt(1,1)=1;         return     end     初始位置的计算%     loc=0;     for i=m:-1:1         if t-tArr(1,i)>=0             loc=i;             break         end     end     %直接代入公式 进行系数求解      if loc>=1 && loc<=m-3         bt(loc,1)=((t-tArr(1,loc))^3)/((tArr(1,loc 3)-tArr(1,loc))*(tArr(1,loc 2)-tArr(1,loc))*(tArr(1,loc 1)-tArr(1,loc)));     end     if loc-1>=1 && loc-1<=m-3         bt1=(t-tArr(1,loc-1))*(t-tArr(1,loc-1))*(tArr(1,loc 1)-t);         bt11=((tArr(1,loc 2)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc)));         bt2=((t-tArr(1,loc-1))*(tArr(1,loc 2)-t)*(t-tArr(1,loc)));         bt21=((tArr(1,loc 2)-tArr(1,loc-1))*(tArr(1,loc 2)-tArr(1,loc))*(tArr(1,loc 1)-tArr(1,loc)));         bt3=(tArr(1,loc 3)-t)*(t-tArr(1,loc))*(t-tArr(1,loc));         bt31=((tArr(1,loc 3)-tArr(1,loc))*(tArr(1,loc 2)-tArr(1,loc))*(tArr(1,loc 1)-tArr(1,loc)));         bt(loc-1,1)=bt1/bt11 bt2/bt21 bt3/bt31;     end     if loc-2>=1 && loc-2<=m-3         bt21=(t-tArr(1,loc-2))*((tArr(1,loc 1)-t)^2);         bt211=(tArr(1,loc 1)-tArr(1,loc-2))*(tArr(1,loc 1)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc));         bt22=(tArr(1,loc 2)-t)*(t-tArr(1,loc-1))*(tArr(1,loc 1)-t);         bt221=(tArr(1,loc 2)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc));         bt23=(tArr(1,loc 2)-t)*(tArr(1,loc 2)-t)*(t-tArr(1,loc));         bt231=(tArr(1,loc 2)-tArr(1,loc-1))*(tArr(1,loc 2)-tArr(1,loc))*(tArr(1,loc 1)-tArr(1,loc));         bt(loc-2,1)=bt21/bt211 bt22/bt221 bt23/bt231;     end     if loc-3>=1 && loc-3<=m-3         bt(loc-3,1)=((tArr(1,loc 1)-t)^3)/((tArr(1,loc 1)-tArr(1,loc-2))*(tArr(1,loc 1)-tArr(1,loc-1))*(tArr(1,loc 1)-tArr(1,loc)));     end end  function bt3=getB3(t,tArr)     m=size(tArr,2)-1;     bt3=zeros(m-3,1);     初始位置的计算%     loc=0;     for i=m:-1:1         if t-tArr(1,i)>=0             loc=i;             break         end     end     %直接代入公式 求解系数      if loc>=1 && loc<=m-3         bt3(loc,1)=6/((tArr(1,loc 3)-tArr(1,loc))*(tArr(1,loc 2)-tArr(1,loc))*(tArr(1,loc 1)-tArr(1,loc)));     end     if loc-1>=1 && loc-1<=m-3         bt1=-6;         bt1=((tArr(1,loc+2)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc)));
        bt2=-6;
        bt21=((tArr(1,loc+2)-tArr(1,loc-1))*(tArr(1,loc+2)-tArr(1,loc))*(tArr(1,loc+1)-tArr(1,loc)));
        bt33=-6;
        bt331=((tArr(1,loc+3)-tArr(1,loc))*(tArr(1,loc+2)-tArr(1,loc))*(tArr(1,loc+1)-tArr(1,loc)));
        bt3(loc-1,1)=bt1/bt11+bt2/bt21+bt33/bt331;
    end
    if loc-2>=1 && loc-2<=m-3
        bt21=6;
        bt211=(tArr(1,loc+1)-tArr(1,loc-2))*(tArr(1,loc+1)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc));
        bt22=6;
        bt221=(tArr(1,loc+2)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc));
        bt23=6;
        bt231=(tArr(1,loc+2)-tArr(1,loc-1))*(tArr(1,loc+2)-tArr(1,loc))*(tArr(1,loc+1)-tArr(1,loc));
        bt3(loc-2,1)=bt21/bt211+bt22/bt221+bt23/bt231;
    end
    if loc-3>=1 && loc-3<=m-3
        bt3(loc-3,1)=-6/((tArr(1,loc+1)-tArr(1,loc-2))*(tArr(1,loc+1)-tArr(1,loc-1))*(tArr(1,loc+1)-tArr(1,loc)));
    end

end