三次插值ClampedB样条曲线Matlab代码实现
时间:2022-09-03 07:00:01
三次插值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