B 样条曲线拟合算法

一、B 样条曲线拟合的概念

1、B 样条曲线定义

一条 $p$ 次 B 样条曲线的定义为
$$C(u)=\sum _{i=0}^n{N}_{i,p}(u)P_i,a\le u\le b\text{\hspace{1em}(1)}$$

其中 $P_i$ 是曲线的控制顶点，$N_{i,p}(u)$ 是定义在非周期节点矢量 $U$ 上的 $p$ 次 B 样条基函数，其值按照下式求取。
$$\begin{array}{rl}{N}_{i,0}(u)& =\{\begin{array}{ll}1,& u_i\le u\le u_{+1},\\ 0,& else\end{array}\\ N_{i,p}(u)& =\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u)\end{array}\text{\hspace{1em}(2)}$$

节点矢量 $U$ 中 $a$ 和 $b$ 的重复度都为 $p+1$，除非特别声明，通常取 $a=0,b=1$
$$U=(\underset{p+1}{\underbrace{a,\cdots ,a}},u_{p+1},\cdots ,u_{m-p-1},\underset{p+1}{\underbrace{b,\cdots ,b}})\text{\hspace{1em}(3)}$$

2、B 样条曲线拟合

如果数据点 Q = { q j : j = 1 , 2 , ⋯   , r } Q=\{q_j:j=1,2,\cdots,r\} Q={qj​:j=1,2,⋯,r} 落在 B 样条曲线上，其则应该满足式（1），即
$$q_j(u_j)=\sum _{j=0}^n{N}_{j,p}(u_j)P_j,j=1,2,\cdots ,r\text{\hspace{1em}(4)}$$

式（4）的矩阵形式为 $Q=NP$ 其中 $Q$ 是包含 $r$ 个数据点的 $r\times 3$ 矩阵，$N$ 是 $r\times n$ 的 B 样条基函数系数矩阵，$P$ 是包含 $n$ 个未知控制顶点的 $n\times 3$ 矩阵。若对数据点 $Q$ 进行参数化，然后确定曲线的节点矢量，则可以求出基函数系数矩阵 $N$，进而求出控制顶点 $P$。可见，通过落在曲线上的点，可反求出曲线的有关参数，并可定义出一条 B 样条曲线，该过程就是所谓的插值。

作为曲线拟合的另一种情况，曲线逼近不要求曲线通过所有的给定数据点，而是所得曲线与给定数据点的偏差落在指定范围内。

二、初始拟合曲线的获取

先对二维数据进行预处理，包括去除重复点、过滤曲率过大的尖锐噪声点和重采样点；然后对采样点的曲率近似求解并分析，确定轮廓关键点；在对关键点进行参数化，并计算节点矢量和控制顶点，最后完成对关键点进行插值。

2.1、有序点列的预处理

首先去除数据点中重复点和孤立噪声点，这方面的成熟方法很多，在此不再赘述。然后对原始点列进行等距重采样处理，其原理如图所示，原始数据点为 { A , B , ⋯   , G } \{A,B,\cdots,G\} {A,B,⋯,G}，重采样后的点列为 { A ′ , B ′ , ⋯   , G ′ } \{A^\prime,B^\prime,\cdots,G^\prime\} {A′,B′,⋯,G′}，设依次连接 $n$ 个原始点所得的线段总长为 $d_n$。现重新布局各点，并使各相关点之间满足如下关系：
$$\mid AB\mid +\mid B{B}^{\prime }\mid =\mid B^{\prime }{C}^{\prime }\mid =\cdots =\mid F^{\prime }F\mid +\mid FG\mid =d_n/n$$

若将原始点的连线（即图中的虚线）近似看成是弧线，而相邻两个重采样点连线（即图中的实现）可近似看作是等弧长所对应的弦长，因此新采样点是近似等距离分布的，这将有助于得到光滑的拟合曲线。

2.2、轮廓关键点的确定

待拟合点列中的某些特定点会对逼近曲线形状产生重要影响，此处称之为轮廓关键点。而曲率的分布可以反映采样点列的总体和局部特征，本文采取基于曲率分布的轮廓关键点遴选方法。

2.2.1、曲率求解与分析

下面先求出采样点的曲率。离散点曲率半径的计算比较复杂，这里采用近似法求解，即第 $i$ 个数据点处的曲率半径 ${\rho }_i$，可近似为经过 3 个相邻数据点 $d_{i-1},d_i$ 和 $d_{i+1}$ 的圆弧半径，其数值用式（5）求出：
$${\rho }_i=\frac{\mid \overrightarrow{d_i{d}_{i+1}}\mid \mid \overrightarrow{d_i{d}_{i-1}}\mid \mid \overrightarrow{d_i{d}_{i+1}}-\overrightarrow{d_i{d}_{i-1}}\mid }{2\mid \overrightarrow{d_i{d}_{i+1}}\times \overrightarrow{d_i{d}_{i-1}}\mid }\text{\hspace{1em}(5)}$$

因此可得第 $i$ 个数据点处得曲率 $k_i=1/rh{o}_i$。

2.2.2、拟定轮廓关键点

不同于常见的曲率极值法，本文通过曲率临界值来确定轮廓关键点，即在分析曲率分布图的基础上，设定曲率临界值 $[k]$，将曲率超过 $[k]$ 的部分数据点初定为轮廓关键点。如果不是封闭曲线，则需同时将两个端点直接确定为关键点。

对于 $[k]$ 值的选取，力求既能反映曲线的大致形状，又尽可能减少关键点数量。较小的 $[k]$ 值将导致出现较多的控制顶点，而过大的 $[k]$ 值则容易失去对曲线形状的控制。基于上述考虑，此处建议取曲率平均值 $k_{\mathrm{avg}}$ 作为临界值 $[k]$。

需要说明的是，一方面，由于数据点列轮廓通常会有凹凸变化，而外凸点和内凹点的曲率值较大，故此时以 $[k]=k_{\mathrm{avg}}$ 为条件确定的关键点数目通常小于数据点总数的一半，因而可有效压缩控制顶点数量。另一方面，由于 $[k]$ 值的选取而造成的逼近误差，可以通过后续的局部优化得到补偿，因而曲线拟合的预设精度可以得到保证。

2.3、初始拟合 B 样条曲线

2.3.1、关键点的参数化

计算量最小的是均匀参数化，但该方法易造成曲线打圈相交的现象。另一种是弦长参数化法，但当数据点急转弯变化时该方法得到的曲线光滑性较差。为使曲线充分反映型值点的分布情况，此处采用向心参数化方法。对于 $n+1$ 个关键点 $Q_k(k=0,1,\cdots ,n)$，记
$$d=\sum _{k=1}^n\sqrt{\mid Q_k-Q_{k-1}\mid }$$

则 ${\overline{u}}_0=0,{\overline{u}}_n=1$，而
$${\overline{u}}_k=\overline{(}u{)}_{k-1}+\sqrt{\mid Q_k-Q_{k-1}\mid }/d,k=1,2,\cdots ,n-1\text{\hspace{1em}(6)}$$

2.3.3、节点矢量的计算

从理论上讲，节点矢量可以等距分布，即 $u_0=u_1=\cdots =u_p=0,u_{m-p}=u_{m-(p-1)}=\cdots =u_m=1,u_{j+p}=1/(n-p+1),j=1,2,\cdots ,n-p$。但是如果讲将其与式（6）联合使用，方程组（4）可能变成奇异方程组。为避免这种情况出现，此处采用取平均值的方法，中间的节点矢量用式（7）求出：
$$u_{j+p}=\frac{1}{p}\sum _{i=j}^{j+p-1}{\overline{u}}_i,j=1,2,\cdots ,n-p\text{\hspace{1em}(7)}$$

一方面，由于式（6）中考虑了轮廓关键点对应的参数值 ${\overline{u}}_i$，可以反映关键点的分布；另一方面，联立式（6）和（7），可以确保方程组（4）的系数矩阵为正定的带状矩阵，当 $\mid i-k\mid \ge p$ 时，$N_{i,p}({\overline{u}}_k)=0$。故可以采用不选主元的高斯消去法求解，减少了算法的计算量。

2.3.3、控制顶点的反算

根据式（2）和（3），通过参数 ${\overline{u}}_k$ 和节点矢量 $U$ 求出基函数 $N_{i,p}({\overline{u}}_k)$，控制顶点 $P_i$ 是式（4）中唯一未知量。