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

本文代码根据插值样条曲线原理实现。

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

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

结果如下：

代码如下：

clc;clear;close;
%输入型值点 l+1个
points=[[4,4];[5,6];[7,7];[8,8.5];[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:1b=getB(times,tArr);p=(b.')*controlPoint;plot(p(1,1),p(1,2),'.b');hold onend%获取准均匀时间节点数组
function tArr=getTimeArray(points)num=size(points,1);tArr=zeros(1,num+6);for i=1:1:num+6if i<=4tArr(1,i)=0;endif i>=num+3tArr(1,i)=1;endif i>4 && i<num+3tArr(1,i)=1/(num-1)*(i-4);endendend%%根据时间节点与型值点获取控制点 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==1t1=tArr(1,4);t2=tArr(1,5);vec1=getB3(t1,tArr).';vec2=getB3(t2,tArr).';bMatrix(i,:)=vec1-vec2;endif i>1&&i<len+2t=tArr(1,i+2);vec=getB(t,tArr).';bMatrix(i,:)=vec;endif i==len+2t1=tArr(1,len+1);t2=tArr(1,len+2);vec1=getB3(t1,tArr).';vec2=getB3(t2,tArr).';bMatrix(i,:)=vec1-vec2;endendz=zeros(len+2,2);for i=1:1:len+2if i==1 ||i==len+2z(i,1)=0;z(i,2)=0;endif i>1&&i<len+2z(i,:)=points(i-1,:);endendbm=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=1if t==1bt(m-3,1)=1;returnendif t==0bt(1,1)=1;returnend%计算初始位置loc=0;for i=m:-1:1if t-tArr(1,i)>=0loc=i;breakendend%直接代入公式 进行系数求解if loc>=1 && loc<=m-3bt(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)));endif loc-1>=1 && loc-1<=m-3bt1=(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;endif loc-2>=1 && loc-2<=m-3bt21=(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;endif loc-3>=1 && loc-3<=m-3bt(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
endfunction bt3=getB3(t,tArr)m=size(tArr,2)-1;bt3=zeros(m-3,1);%计算初始位置loc=0;for i=m:-1:1if t-tArr(1,i)>=0loc=i;breakendend%直接代入公式 进行系数求解if loc>=1 && loc<=m-3bt3(loc,1)=6/((tArr(1,loc+3)-tArr(1,loc))*(tArr(1,loc+2)-tArr(1,loc))*(tArr(1,loc+1)-tArr(1,loc)));endif loc-1>=1 && loc-1<=m-3bt1=-6;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=-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;endif loc-2>=1 && loc-2<=m-3bt21=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;endif loc-3>=1 && loc-3<=m-3bt3(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)));endend

三次插值ClampedB样条曲线Matlab代码实现相关推荐

非均匀三次B样条曲线插值实现及MATLAB代码
这篇博客跟我上一篇博客<均匀三次B样条曲线插值实现及MATLAB代码>的内容有点像,只是在基函数的计算上不同,造成均匀/非均匀的区别. 参考资料: [1](这个PPT讲得很通俗,但对于多插 ...
均匀三次B样条曲线插值实现及MATLAB代码
参考资料: [1](这个PPT讲得很通俗,但对于多插值点分段曲线的内容漏讲了一个知识点)三次周期B样条曲线的算法 - 百度文库 (baidu.com) [2](这个介绍只有两个插值点的三次B样条曲线, ...
C++编写任意次clampedB样条曲线（曲线分别与第一个控制点和最后一个控制点的第一边和最后一边相切）
我们知道,B样条曲线根据其起始点与控制点的关系可分为:openB样条曲线,clampedB样条曲线,closedB样条曲线.OpenB样条曲线不会与第一个控制点和最后一个控制点的第一条边和最后一条边相 ...
材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）
材料力学求解器-刚架与桁架杆系的计算机求解(附matlab代码) 1 刚架的计算机求解 1.1位移法与刚度矩阵 1.2 matlab程序 2 桁架的计算机求解材料力学是一门非常成熟的学科,里面有大量 ...
一种用于模拟电晕放电的高效半拉格朗日算法（Matlab代码实现）
目录摘要 1 概述 2 数学模型与方法 3 讲解 3.1 测试1 3.2 测试2 3.3 测试3 3.4 测试4 4 Matlab代码实现摘要提出了一种无需通量校正的高效电晕放电模拟算法.称为 ...
基于神经网络多项式插值的图像超分辨重构研究-附Matlab代码
⭕⭕ 目录 ⭕⭕ ✳️ 一.引言 ✳️ 二.基于单帧图像的超分辨率重构技术 ✳️ 2.1 最近邻域插值法 ✳️ 2.2 双线性插值法 ✳️ 2.3 双三次插值法(Keys'插值) ✳️ 三.神经网络 ...
气象数据Grib格式解析的Python代码和Matlab代码
以.grb/.grb1/.grb2为扩展名的都是气象数据,气象数据中可以存储多个内容,如云量.雪深.气压.风速等内容,或者具有时间序列的云量等.这些文件不可以直接打开成图片,若想直观地查看grib数据 ...
龙格-库塔法(runge-kutta)matlab代码及含义,龙格-库塔法(Runge-Kutta)matlab代码及含义...
龙格-库塔法(Runge-Kutta)matlab代码及含义龙格-库塔法(Runge-Kutta) 数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭 ...
arima模型matlab代码_PSTR面板平滑转换模型简介（附Matlab代码分享）
写论文的时候用到的~相关的资料太少了,做一些简单内容和资料的分享.(PSTR模型的Matlab代码分享在最后)本文主要为简单理论和粗暴实操~ 有用的话可以点个赞哟(知乎小白卑微求赞) 嘻嘻下面进入正题 ...

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

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

三次插值ClampedB样条曲线Matlab代码实现相关推荐

最新文章

热门文章