MATLAB---构造一个插值三次样条曲线

function InterpCubicSplineCurv()
%本程序的功能是构造一个插值三次样条曲线

N = 12;
SamplPs = CollectSPFergusonCi2(N);%采集型值点
x=SamplPs(1,:);
y=SamplPs(2,:);
sT = [0,1/3];
eT = [0,-1/3];

s = AccuArcL(SamplPs');%采用累加弦长参数化方法对数据点进行参数化
%也就是为每个点指定一个参数

mxs = Intrepolation3Spline(s,x,sT(1),eT(1));%计算切矢量在x方向的分量
mys = Intrepolation3Spline(s,y,sT(2),eT(2));%计算切矢量在y方向的分量

mxss=[sT(1) mxs' eT(1)];
myss=[sT(2) mys' eT(2)];

n=length(x)-1;
k=1;
for i=1:n
    smx=mxss(i);
    emx=mxss(i+1);
    smy=myss(i);
    emy=myss(i+1);
    h=s(i+1)-s(i);
     %在每段中采集21个点
    for j=1:21
        u=(j-1)/20;
        xf(k)=[1 u u*u u*u*u]*[1 0 0 0;0 0 1 0;-3 3 -2 -1;2 -2 1 1]*[x(i) x(i+1) h*smx h*emx]';
        yf(k)=[1 u u*u u*u*u]*[1 0 0 0;0 0 1 0;-3 3 -2 -1;2 -2 1 1]*[y(i) y(i+1) h*smy h*emy]';
        %采集的点放入整体数组xf和yf
        k=k+1;
    end
    %在每段中采集21个点
end
%对计算出的点进行密化采样

figure%开一个新的图形窗口
hold on
plot(xf,yf,'lineWidth',2)
plot(x,y,'k*')
hold off
box on%坐标轴显示为包围盒的形式
axis equal

function ms = Intrepolation3Spline(X,Y,sm,em)
%采用插值参数三次样条为每个型值点计算切矢量

n=length(X)-2;%注意到首末端点要添加边界条件的
CoeffMat=zeros(n,n);%预定义线性方程组的系数矩阵
Eqb=zeros(1,n);%预定义线性方程组的常数向量
n=n+2;
for i=2:n-1
    hi=X(i)-X(i-1);
    hiq=X(i+1)-X(i);
    lmdi=hiq/(hi+hiq);
    miui=1-lmdi;
    ci=3*(lmdi*(Y(i)-Y(i-1))/hi+miui*(Y(i+1)-Y(i))/hiq);
    j=i-1;
    if(j==1)
        CoeffMat(j,1)=2;
        CoeffMat(j,2)=miui;
        Eqb(j)=ci-lmdi*sm;
    elseif(j==n-2)
        CoeffMat(j,n-3)=lmdi;
        CoeffMat(j,n-2)=2;
        Eqb(j)=ci-miui*em;
    else
        CoeffMat(j,j-1)=lmdi;
        CoeffMat(j,j)=2;
        CoeffMat(j,j+1)=miui;
        Eqb(j)=ci;
    end
end
%为矩阵CoeffMat和常数向量赋值Eqb
ms=CoeffMat\Eqb';%求解线性方程组得到切矢量

function SamplPs = CollectSPFergusonCi2(N)
%N表示样点个数

xita = [0,pi/27,pi/20,pi/18,pi/5,pi/3,pi/2,3*pi/4,3.5*pi/4,pi];
%根据参数方程计算圆周上的点
SamplPs(1,:) = cos(xita);%x的坐标放在矩阵的第一行
SamplPs(2,:) = sin(xita);%y的坐标放在矩阵的第二行
%根据参数方程计算圆周上的点
% plot(SamplPs(1,:),SamplPs(2,:),'k*')%测试程序正确性用的
% axis equal;