三次插值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 && i1&&i1&&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