数值计算 --- B样条函数（B-spline）

松下J27

已于 2024-07-29 15:03:55 修改

阅读量1.5w

点赞数 91

分类专栏： Numerical Methods Digital Image Processing 文章标签： B样条 b-spline B spline 样条函数

于 2021-09-27 19:38:35 首次发布

本文链接：https://blog.csdn.net/daduzimama/article/details/120515540

版权

Digital Image Processing 同时被 2 个专栏收录

12 篇文章 133 订阅

订阅专栏

Numerical Methods

2 篇文章 13 订阅

订阅专栏

B样条函数（B-spline）

B-spline就是Basis Spline的简称：

术语B-spline一词，是由罗马尼亚裔美国数学家Isaac Jacob Schoenberg创造的（就是下图的这个老爷爷，我在wiki百科上找到了他的照片）。一个n阶的B-spline函数，是一个变量为x的n-1次分段多项式函数f(x)。分段函数之间，段与段之间的断点被称为knots。

B-spline的定义：（虽然我也很不想一上来就给个定义）

B-spline函数(上图中的f(t))，可以看成是一系列控制点 $P_{i}$ (Control points)和权重函数 $N_{i,D}$ (Weight function)的线性组合。简而言之，B-spline曲线就是每个控制点乘以与之对应的权重函数的加权和。控制点是我们后期指定的，而权重函数是我们预先定义好的，只跟阶数D有关，不会因为控制点的变化而改变。

权重函数 $\large N_{i,D}(t)$

（注意：以下只讨论权重函数N(t)，暂时还不涉及控制点P）

权重函数 $N_{i,D}(t)$ 的定义：

当阶数D=1时：

函数 $N_{i,D}(t)$ 的定义，可以参考Cox–de Boor的递归公式。对于给定的knots，当D=1时，一阶权重函数 $N_{i,1}(t)$ 可以定义为：

（公式1）

当阶数D>1时：

有了一阶权重函数 $N_{i,1}(t)$ 以后，其他所有阶数(D)大于1的高阶权重函数都可以按照如下递归表达式来定义的：

（公式2）

公式1和公式2合称为Cox-de Boor recursion formula。

且有如下性质：

1，阶数为D的权重函数 $N_{i,D}(t)$ 的，是关于t的D-1次多项式。

2，Positivity: 当t在此区间 $t_{i}< t< t_{i+D}$ 内时，函数值一定不为0，即 $N_{i,D}(t)>0$ 。

3，Local support: 前面的Positivity定义了高阶权重函数在哪里有值，而这里的local support则是定义了高阶权重函数在哪里没有值。和一阶权重函数一样，对于区间[ $t_{i}$ ~ $t_{i+D}$ ]之外的t， $N_{i,D}(t)=0$ 。
4，Partition of Unity: 在区间[ $t_{i}$ ~ $t_{i+1}$ ]内，所有非0阶权重函数的和为1。(这个性质我自己也没看懂)

5，Recursion: 递归性，这一性质由公式2决定。

例如，当D=2时：

$\LARGE \mathbf N_{i,2}(t)=\frac{(t-t_{i})N_{i,1}(t)}{t_{i+1}-t_{i}}+\frac{(t_{i+2}-t)N_{i+1,1}(t)}{t_{i+2}-t_{i+1}},t_{i}\leqslant t < t_{i+2}$

其中 $\LARGE N_{i,1}(t)$ 和 $\LARGE N_{i+1,1}(t)$ 都是已知的一阶权重函数。

6，一个D阶B-spline函数是由一些D-1阶多项式曲线连起来的(参考性质1)。每段函数之间的断点/连接点(breakpoints)，被称为knot，用非递减序列 $t_{0}\leqslant t_{1}\leqslant t_{2}...\leqslant t_{m}$ 表示，这些knot的合集被称为knot vector list，共m+1个knot。

$T=(t_{0},t_{1},...,t_{m})$

knot与knot之间的半开半闭区间[ $t_{i}$ , $t_{i+1}$ )被称为knot span或i-th knot span，意思就是两个knot之间的跨度。如果knot vector list中，所有的knot都是等间隔的，即 $t_{1}-t_{0}=t_{2}-t_{1}=t_{3}-t_{2}=....=t_{m}-t_{m-1}$ ，则我们说T是uniform的，否则，称之为non-uniform。(本文中讨论的都是uniform B-spline)

7，一但knot vector list确定了以后，对于阶数为D的B-spline函数，每段函数与函数之间在断点knot处保证 $C^{D-2}$ 阶连续。（注：函数的连续性和knot的选择方式有关）

举个例子，简要的说明一下上述性质。对于一阶权重函数 $N_{i,1}(t)$ 而言，阶数D=1，每段函数都是一个D-1=0阶多项式。当i=0时，有一阶权重函数 $N_{0,1}(t)$ ，对应knot span[0，1)。在[0，1)之间的函数值不为0，在这之外函数值都等于0。

权重函数 $N_{i,D}(t)$ 定义的简化：

为了方便计算，上面定义的权重函数可以做进一步的简化。一般情况下，会把knot vector list也就是函数之间的断点 $t_{0},t_{1},t_{2},...,t_{m}$ ，简化成一组非递减的正整数。简而言之，就是把原来的 $t_{0},t_{1},t_{2},...,t_{m}$ 变成非递减的正整数0,1,2,...m。

（简化后的公式1）

（简化后的公式2）

为了更好的说明当D>1时的Cox–de Boor递归公式，我这里借鉴了一张别人网站里的插图，为了后续描述方便，我这里暂且称它为Cox–de Boor递归流程图（见参考文献1）。

首先，让我们把目光放在图中三角形的最左上角的一个小三角形 $N_{0,1}$ , $N_{1,1}$ 和 $N_{0,2}$ ，他所要表达的是，要想求出二阶函数 $N_{0,2}$ ，我们需要先知道一阶函数 $N_{0,1}$ 和 $N_{1,1}$ 。这就是说，要想知道上图中左起第二列的所有二阶函数 $N_{0,2}$ ， $N_{1,2}$ ，。。。 $N_{4,2}$ ，必须先知道他们所对应的左起第一列的一阶函数 $N_{0,1}$ ， $N_{1,1}$ 。。。 $N_{5,1}$ 。依此类推，当我们要计算六阶函数 $N_{0,6}$ 时，首先需要先知道五阶函数 $N_{0,5}$ 和 $N_{1,5}$ ，其次我们需要相应的所有四阶函数，三阶函数。。。直到一阶函数。

从公式的角度看也一样，根据简化后的公式，我们发现要想得到阶数为D的第i段函数 $N_{i,D}(t)$ ，需要先知道阶数为D-1的第i段函数 $N_{i,D-1}(t)$ 和阶数为D-1的第i+1段函数 $N_{i+1,D-1}(t)$ 。以此类推，要想知道D-1阶的第i段函数 $N_{i,D-1}(t)$ 和第i+1段函数 $N_{i+1,D-1}(t)$ ，我们又分别需要 $N_{i,D-2}(t),N_{i+1,D-2}(t)$ 和 $N_{i+1,D-2}(t),N_{i+2,D-2}(t)$ ，以此类推直到对应的一阶权重函数。

下面我们由低到高，详细探讨一下几个不同阶数的权重函数 $N_{i,D}(t)$ 。我这里再次强调，到目前为止都没有讨论过控制点 $P_{i}$ ，而只讨论权重函数，目的是为了避免混淆对权重函数 $N_{i,D}(t)$ 和控制点 $P_{i}$ ，这两个B-spline中最重要的概念的理解。因为，所谓的B-spline就是这两个概念的线性组合，只有分别理解他们的各自的作用，才能更好的理解B-spline。

一阶权重函数 $N_{i,1}(t)$ :

一阶权重函数，D=1，t=[i,i+1)，有效区间为i+1。

首先，我们知道权重函数都有一个有效区间，这里我暂且把他称之为t的定义域，毕竟函数在定义域之外是有值的，且值为0。

这里我们统一以7个knots为例，knot vector list T={0,1,2,3,4,5,6}，每段函数的knot span为1。对于第0段函数，也就是当i=0时，t的定义域为[0,1)，函数 $N_{0,1}(u)$ 的断点knot为0，1。对于最后一段函数，也就是第5段函数，此时，i=5，u=[5,6)，函数两端的端点knot为5，6。(注意，此时还没有用到递归函数。)

一阶权重函数 $N_{i,1}(t)$ 的图像:

注意：相对于第i段权重函数的图像而言，第i+1段权重函数只不过是对第i段函数向右移动一位后的结果，以此类推。

Matlab code:

%% B-spline demo for CSDN
% Created: early fall, 2022. (2022/09/02)
% Author: Z.Zhu, J27
% Copy Rights Reserved.

close all
clear all


Data=[6,6,6,6,1,1];
%% 一阶B-spline中的权重函数Ni,1(u),ti为正整数的定义域
%n=6, i=0~6,共7个点, each knot span=order=1,总共可以生成7-1=6段曲线。
%N0,1(u),u=[t0~t1)=[0~1)
freq=0.01;
ti_D1=0:freq:6;
step=1/freq;
%第一段函数
u0=ti_D1(1:step);
N01=ones(1,step);
figure;
subplot(3,1,1);
plot(u0,N01);
xlabel('0<=u<1')
legend('N01(u)')
axis([ 0 6 0 2])
%第二段函数
u1=ti_D1(step+1:2*step);
N11=ones(1,step);
subplot(3,1,2);
plot(u1,N11)
xlabel('1<=u<2')
legend('N11(u)')
axis([ 0 6 0 2])
%第六段函数
u5=ti_D1(5*step+1:6*step);
N51=ones(1,step);
subplot(3,1,3);
plot(u5,N51)
xlabel('5<=u<6')
legend('N51(u)')
axis([ 0 6 0 2])

二阶权重函数 $N_{i,2}(t)$ :

二阶权重函数，D=2，t=[i,i+2)，第i段函数的有效区间为i+2。

同样是7个knots，knot vector list T={0,1,2,3,4,5,6}，每段函数的knot span=D=2。

根据Cox–de Boor的递归公式，一阶权重函数的值是分区间给定的，而对于二阶权重函数的计算，需要根据递归函数来调用一阶权重函数（下图中简化后的公式2）。

首先，我们把目光移到公式中的阶数D，当我们要计算 $N_{i,D}(t)$ ，我们预先要知道 $N_{i,D-1}(u)$ 和 $N_{i+1,D-1}(u)$ 。也就是说，要想得到第i段2阶权重函数，需要预先知道第i段1阶权重函数和阶数同为1的第i+1段权重函数。

同样，现在我们分别把i=0，1，5和D=2代入简化后的公式2，分别考察这些函数。

当i=0时：有第0段函数 $N_{0,2}(t)$ ，定义域t=[0,2)。而用于合成他的两个一阶函数 $N_{0,1}(t)$ 和 $N_{1,1}(t)$ 的定义域分别是[0,1)和[1,2)。

为了更好的理解该函数，可以借助之前的Cox–de Boor递归流程图：

这幅图指出，当t=[0,1)时，仅仅只有函数 $N_{0,1}(t)$ 对函数 $N_{0,2}(t)$ 起作用，也就是式中的前半部分 $tN_{0,1}(t)$ 。

此时， $N_{0,1}(t)=1$ 推出 $N_{0,2}(t)=t$ 。

同理，当t=[1,2)时，仅仅只有函数 $N_{1,1}(t)$ 对函数 $N_{0,2}(t)$ 起作用，也就是式中的后半部分 $(2-t)N_{1,1}(t)$ 。

此时， $N_{1,1}(t)=1$ 推出 $N_{0,2}(t)=2-t$ 。

当i=1时：有第1段函数 $N_{1,2}(t)$ ，定义域t=[1,3)：

当i=4时：有第4段函数 $N_{4,2}(t)$ ，定义域t=[4,6)：

（注意：当D=2时，由于只有0~6个knots，i的最大取值范围不会超过4，i的最大取值范围等于knot list中的最大值6-D=4）

二阶权重函数 $N_{i,2}(t)$ 的图像:

注意：和一阶权重函数的图像类似，二阶权重函数的图像也可以通过把前一段函数整体向右移动一位得到。这一点也可以由于，函数自身的变化看出来。 $N_{0,2}(t)$ 中前半段曲线t向右移动一位，得到t-1，正好等于 $N_{1,2}(t)$ 的前半段。同理， $N_{0,2}(t)$ 的后半段曲线2-t向右移动一位，得到2-(t-1)=3-t，也正好等于 $N_{1,2}(t)$ 的后半段。

Matlab code:

%% 二阶B-spline中的权重函数Ni,2(u)
%n=6, i=0~6,共7个点, each knot span=order=2,总共可以生成7-2=5段曲线。
freq=0.01;
ti_D2=0:freq:6;
step=1/freq;
%第一段函数
%N0,2(u),u=[t0~t2)=[0~2)
u0=ti_D2(1:2*step);
N02=zeros(1,2*step);
%前半部分[0,1)
N02(1:step)=u0(1:step);
%后半部分[1,2)
N02(step+1:2*step)=2-u0(step+1:2*step);
figure;
subplot(3,1,1);
plot(u0,N02);
xlabel('0<=u<2')
legend('N02(u)')
axis([ 0 6 0 2])
%第二段函数
%N1,2(u),u=[t1~t3)=[1~3)
u1=ti_D2(step+1:3*step);
N12=zeros(1,2*step);
%前半部分[1,2)
N12(1:step)=u1(1:step)-1;
%后半部分[2,3)
N12(step+1:2*step)=3-u1(step+1:2*step);
subplot(3,1,2);
plot(u1,N12);
xlabel('1<=u<3')
legend('N12(u)')
axis([ 0 6 0 2])
%第五段函数
%N4,2(u),u=[t4~t6)=[4~6)
u4=ti_D2(4*step+1:6*step);
N42=zeros(1,2*step);
%前半部分[1,2)
N42(1:step)=u4(1:step)-4;
%后半部分[2,3)
N42(step+1:2*step)=6-u4(step+1:2*step);
subplot(3,1,3);
plot(u4,N42);
xlabel('4<=u<6')
legend('N42(u)')
axis([ 0 6 0 2])

三阶权重函数 $N_{i,3}(t)$ ：

三阶权重函数，D=3，t=[i,i+3)，第i段函数的有效区间为i+3。

总共7个knots，knot vector list T={0,1,2,3,4,5,6}，每段函数的knot span=D=3。

按照Cox–de Boor递归公式，要求三阶的权重函数 $N_{i,3}(t)$ ，首先我们要先知道二阶权重函数 $N_{i,2}(t)$ ，而这些函数都已经在前面求出来了。和前面一样，我们把i=0，1，3和D=3逐一代入简化后的公式2，看一下三阶权重函数。

1，当i=0时，t=[0,3)，第0段函数 $N_{0,3}(t)$ ：

我们可以看到，和前面一样，为了求出 $N_{0,3}(t)$ ，我们需要先知道 $N_{0,2}(t)$ 和 $N_{1,2}(t)$ （前面已经算好了）。

此时，我们要特别注意函数 $N_{0,3}(t)$ 的定义域t=[0,3)，用于合成他的子函数之一 $N_{0,2}(t)$ 的定义域是[0,2)，继而我们还会发现，对于 $N_{0,2}(t)$ 而言，用于合成他的子函数之一 $N_{0,1}(t)$ 的定义域是[0,1)，在这以外都是0。也就是说，对于这个递归分段函数，我们要特别注意他的有效区间，或者说是每个子函数的有效定义域。（这在Cox–de Boor递归流程图中会更加直观，我们分三个区间讨论）

1.1，当i=0时，对t=[0,1)而言：

这时只有 $N_{0,2}(t)$ 在此处有定义，其他都为0。因此，只有 $N_{0,2}(t)$ 对函数 $N_{0,3}(t)$ 起作用，此时， $N_{0,2}(t)$ 在[0,1)上等于t， $N_{1,2}(t)$ 在[0,1)上等于0，得到 $N_{0,3}(t)=t^{2}/2$ 。而这一点，在我们之前的推导结果中也显而易见。

（上图中画红线处）

1.2，当i=0时，对t=[1,2)而言：

此时， $N_{0,2}(t)$ 和 $N_{1,2}(t)$ 都有定义，他们同时对函数 $N_{0,3}(t)$ 起作用，此时， $N_{0,2}(t)$ 在[1,2)上等于2-t， $N_{1,2}(t)$ 在[1,2)上等于t-1，得到 $N_{0,3}(t)=(-2t^{2}+6t-3)/2$ 。

（上图中画红线处）

1.3，当i=0时，对t=[2,3)而言：

t=[2,3)时的情形和t=[0,1)时的情形很类似，此时，只有 $N_{1,2}(t)$ 在此处有定义，其他都为0。因此，这时只有函数 $N_{1,2}(t)$ 对函数 $N_{0,3}(t)$ 起作用。此时， $N_{1,2}(t)$ 等于3-t， $N_{0,2}(t)$ 等于0，得到 $N_{0,3}(t)=(3-t)^{2}/2$ 。

（上图中画红线处）

2，当i=1时，t=[1,4)，第1段函数 $N_{1,3}(t)$ ：

3，当i=3时，t=[3,6)，第3段函数 $N_{3,3}(t)$ ：

（注意：当D=3时，由于只有0~6个knots，i的最大取值范围不会超过3，i的最大取值范围等于knot list中的最大值6-D=3）

三阶权重函数 $N_{i,3}(t)$ 的图像:

Matlab code:

%% 三阶B-spline中的权重函数Ni,2(u)
%n=6, i=0~6,共7个点, each knot span=order=3,总共可以生成7-3=4段曲线。
freq=0.01;
ti_D3=0:freq:6;
step=1/freq;
%第一段函数
%N0,3(u),u=[t0~t3)=[0~3)
u0=ti_D3(1:3*step);
N03=zeros(1,6*step+1);
N03(1:step)=u0(1:step).^2./2;%0~1
N03(step+1:2*step)=(-2*u0(step+1:2*step).^2+6*u0(step+1:2*step)-3)/2;%1~2
N03(2*step+1:3*step)=(3-u0(2*step+1:3*step)).^2/2;%2~3
figure;
subplot(4,1,1);
plot(ti_D3,N03);
xlabel('0<=u<3')
legend('N03(u)')
axis([ 0 6 0 2])
%第二段函数
%N1,3(u),u=[t1~t4)=[1~4)
u1=ti_D3(step+1:4*step);
N13=zeros(1,6*step+1);
N13(step+1:2*step)=(u1(1:step)-1).^2./2;%1~2
N13(2*step+1:3*step)=(-2*u1(step+1:2*step).^2+10*u1(step+1:2*step)-11)/2;%2~3
N13(3*step+1:4*step)=(4-u1(2*step+1:3*step)).^2/2;%3~4
subplot(4,1,2);
plot(ti_D3,N13);
xlabel('1<=u<4')
legend('N13(u)')
axis([ 0 6 0 2])
%第三段函数
%N2,3(u),u=[t2~t5)=[2~5)
u2=ti_D3(2*step+1:5*step);
N23=zeros(1,6*step+1);
N23(2*step+1:3*step)=(u2(1:step)-2).^2./2;%2~3
N23(3*step+1:4*step)=(-2*u2(step+1:2*step).^2+14*u2(step+1:2*step)-23)/2;%3~4
N23(4*step+1:5*step)=(5-u2(2*step+1:3*step)).^2/2;%4~5
subplot(4,1,3);
plot(ti_D3,N23);
xlabel('2<=u<5')
legend('N23(u)')
axis([ 0 6 0 2])
%第四段函数
%N3,3(u),u=[t3~t6)=[3~6)
u3=ti_D3(3*step+1:6*step);
N33=zeros(1,6*step+1);
N33(3*step+1:4*step)=(u3(1:step)-3).^2./2;%1~2
N33(4*step+1:5*step)=(-2*u3(step+1:2*step).^2+18*u3(step+1:2*step)-39)/2;%2~3
N33(5*step+1:6*step)=(6-u3(2*step+1:3*step)).^2/2;%3~4
subplot(4,1,4);
plot(ti_D3,N33);
xlabel('3<=u<6')
legend('N33(u)')
axis([ 0 6 0 2])

关于权重函数的两点总结：

1，在文章的开头处，我提到了权重函数 $N_{i,D}$ 的几点特征，其中的Positivity非常重要，他指出了权重函数的有效区间，也就是定义域t=[i,i+D)。函数在[i,i+D)以内的值一定不为0，而在这以外都是0。这一点我们依然可以借助Cox–de Boor递归流程图看出来。

比如说，我们要看一下函数 $N_{1,4}(t)$ 的定义域，我们可以沿着 $N_{1,4}(t)$ 的两条放射线朝着图像左边一直追溯到底。他所夹着的区域，就是他的定义域，或者说是有效区间，也叫非零区间。同时，我们还能在这个夹角区间内，找到 $N_{1,4}(t)$ 的子函数，以及子函数的子函数所对应的定义域。

2，对于给定的knot span[i,i+1)，如何才能知道，生成哪些权重函数的时候会用到它？或者说，他这个区间所对应的最低阶的一阶函数 $N_{i,1}(t)$ ，可以生成哪些高阶函数？要回答这个问题，我们需要在Cox–de Boor递归流程图中，沿着 $N_{i,1}(t)$ 的两条放射线向右追溯。比如说 $N_{1,4}(t)$ ，他所能生成的高阶函数，如下图所示：

带有控制点 $\large P_{i}$ (Control Point)的B-Spline Curve

（注意：现在开始讨论控制点CP的应用）

现在让我们再回到之前对B-Spline的定义，之前我们已经详细介绍了权重函数N，现在我们着重讨论控制点P，以及加入控制点P后最终生成的B-Spline曲线。

假定控制点的数量为n+1，他们分别是控制点 $P_{0}$ ， $P_{1}$ ，。。。 $P_{n}$ 。权重函数的阶数为D，那么在生成B-spline时所需的knot断点数m，必须满足：

m=n+D+1

也就是说，变量t所在的knot vector list为T={0，1，。。。m-1}。

我们以三阶B-spline函数为例，D=3，给定6个控制点[2,6,3,8,1,1]，即n=5，共需m=5+3+1=9个knots。首先，按照我们上面的推导，也就是knots={0~6}共7个knots的函数推导，我们已经得到了N03，N13,N23和N33。

根据带有控制点的计算公式：

首先，我们分别计算函数f(0)，f(1)，f(2)，f(3)。其中，f(0)=P0xN03，f(1)=P1xN13，f(2)=P2xN23，f(3)=P3xN33，f(4)=P4xN43，f(5)=P5xN53。共六个控制点，每个控制点乘以对应的权重函数，就能得到当前点所生成的曲线。他们的函数图像如下所示：

这里我们要注意两点:

1，其实这里除了控制点的值不同，每个权重函数都相同，只是他们所对应的定义域不同，即，曲线相同而位置不同。

2，根据我们之前的推导，我们现在只有N03~N33的函数。因此，我们暂时无法使用最后两个控制点P4和P5，因为我们暂时还不知道权重函数N43和N53。

3，上面我在给出控制点时，只给了具体的值，而没有给出每个控制点所对应的位置。在实际应用中需要特别注意。

接下来，只需把他们都加到一起，就能得到最终的B-Spline Curve：

f(t)=f(0)+f(1)+f(2)+f(3)=P0xN03+P1xN13+P2xN23+P3xN33

下图是最终生成的B-Spline Curve，注意他在每个knot断点处都是连续的。根据权重函数的定义，他的连续性为 $C^{D-2}$ ，也就是说最终生成的Curve在每个断点处都连续 $C^{0}$ ，且，在断点处的导数也连续 $C^{1}$ 。

对于三阶B-Spline而言，正如我们前面说的，每个控制点CP所乘的权重函数其实是一模一样的，只不过是在不同的定义域定义的，这也直接决定了输出函数的位置。也就是说，CP值的改变会极大地影响整个曲线的形状。现在我们把控制点的值由[2,6,3,8,1,1]改成[6,6,6,6,1,1]，我们再来看看最终生成的Curve，并和之前的曲线做一个比较（他们使用的是相同的权重函数，只是使用了不同的控制点）。注意：最后两个控制点P4=1,P5=1在这里用不到。

Matlab code:

%% add Control Points 
%P(u)=PN
Px0=Data(1)*N03;%P(0)=CP0*N03,domain=0~3
Px1=Data(2)*N13;%P(1)=CP1*N13,domain=1~4
Px2=Data(3)*N23;%P(2)=CP2*N23,domain=2~5
Px3=Data(4)*N33;%P(3)=CP3*N33,domain=3~6

figure;
subplot(4,1,1);
plot(ti_D3,Px0);
xlabel('0<=u<3')
legend('P(0)=CP0*N03')
axis([ 0 6 0 10])
subplot(4,1,2);
plot(ti_D3,Px1);
xlabel('1<=u<4')
legend('P(1)=CP1*N13')
axis([ 0 6 0 10])
subplot(4,1,3);
plot(ti_D3,Px2);
xlabel('2<=u<5')
legend('P(2)=CP2*N23')
axis([ 0 6 0 10])
subplot(4,1,4);
plot(ti_D3,Px3);
xlabel('3<=u<6')
legend('P(3)=CP3*N33')
axis([ 0 6 0 10])

%Bspline = P(0) + P(1) + P(2) + ...
BsplineY=Px0+Px1+Px2+Px3;

figure;plot(ti_D3,BsplineY)
axis([ 0 6 0 10])
hold on
stem(0:5,Data)
legend('B-spline','Control Points')
title("B-Spline generated with CP0~CP3.Note:the last two CP4 and CP5 are not used!");

（全文完）

2022年8月4日晚，上海

作者 --- 松下J27

1，Matlab code added, 2022/09/02下班后

2，对文章的排版做了部分修改，2022/10/28

3，对文中B-spline函数的定义进行了修改（修改了原有的插图），2023/07/31

4，对文章进行了大改，2023/08/03

参考文献(鸣谢)：

1，B-spline Basis Functions: Definition

2，1.4.1 B-splines

3，https://en.wikipedia.org/wiki/Isaac_Jacob_Schoenberg

诗词赏析：