一、B 样条曲线拟合的概念
1、B 样条曲线定义
一条
p
p
p 次 B 样条曲线的定义为
C
(
u
)
=
∑
i
=
0
n
N
i
,
p
(
u
)
P
i
,
a
≤
u
≤
b
(1)
C(u)=\sum_{i=0}^nN_{i,p}(u)P_i,\quad a\leq u\leq b\tag{1}
C(u)=i=0∑nNi,p(u)Pi,a≤u≤b(1)
其中
P
i
P_i
Pi 是曲线的控制顶点,
N
i
,
p
(
u
)
N_{i,p}(u)
Ni,p(u) 是定义在非周期节点矢量
U
\pmb U
U 上的
p
p
p 次 B 样条基函数,其值按照下式求取。
N
i
,
0
(
u
)
=
{
1
,
u
i
≤
u
≤
u
+
1
,
0
,
e
l
s
e
N
i
,
p
(
u
)
=
u
−
u
i
u
i
+
p
−
u
i
N
i
,
p
−
1
(
u
)
+
u
i
+
p
+
1
−
u
u
i
+
p
+
1
−
u
i
+
1
N
i
+
1
,
p
−
1
(
u
)
(2)
\begin{aligned} N_{i,0}(u)&=\begin{cases}1,\quad &u_i\leq u\leq u_{+1},\\ 0,&else \end{cases}\\[3ex] 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)\tag{2} \end{aligned}
Ni,0(u)Ni,p(u)={1,0,ui≤u≤u+1,else=ui+p−uiu−uiNi,p−1(u)+ui+p+1−ui+1ui+p+1−uNi+1,p−1(u)(2)
节点矢量
U
\pmb U
U 中
a
a
a 和
b
b
b 的重复度都为
p
+
1
p+1
p+1,除非特别声明,通常取
a
=
0
,
b
=
1
a=0,b=1
a=0,b=1
U
=
(
a
,
⋯
,
a
⏟
p
+
1
,
u
p
+
1
,
⋯
,
u
m
−
p
−
1
,
b
,
⋯
,
b
⏟
p
+
1
)
(3)
U=\bigg(\underbrace{a,\cdots,a}_{p+1},u_{p+1},\cdots,u_{m-p-1},\underbrace{b,\cdots,b}_{p+1}\bigg)\tag{3}
U=(p+1
a,⋯,a,up+1,⋯,um−p−1,p+1
b,⋯,b)(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
)
=
∑
j
=
0
n
N
j
,
p
(
u
j
)
P
j
,
j
=
1
,
2
,
⋯
,
r
(4)
q_j(u_j)=\sum_{j=0}^nN_{j,p}(u_j)P_j,\quad j=1,2,\cdots,r\tag{4}
qj(uj)=j=0∑nNj,p(uj)Pj,j=1,2,⋯,r(4)
式(4)的矩阵形式为 Q = N P \pmb Q=\pmb N\pmb P Q=NP 其中 Q \pmb Q Q 是包含 r r r 个数据点的 r × 3 r\times 3 r×3 矩阵, N \pmb N N 是 r × n r\times n r×n 的 B 样条基函数系数矩阵, P \pmb P P 是包含 n n n 个未知控制顶点的 n × 3 n\times 3 n×3 矩阵。若对数据点 Q \pmb Q Q 进行参数化,然后确定曲线的节点矢量,则可以求出基函数系数矩阵 N \pmb N N,进而求出控制顶点 P \pmb P 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
n
n 个原始点所得的线段总长为
d
n
d_n
dn。现重新布局各点,并使各相关点之间满足如下关系:
∣
A
B
∣
+
∣
B
B
′
∣
=
∣
B
′
C
′
∣
=
⋯
=
∣
F
′
F
∣
+
∣
F
G
∣
=
d
n
/
n
\mid AB\mid+\mid BB^\prime\mid\ =\ \mid B^\prime C^\prime\mid\ =\cdots=\ \mid F^\prime F\mid+\mid FG\mid\ =\ d_n/n
∣AB∣+∣BB′∣ = ∣B′C′∣ =⋯= ∣F′F∣+∣FG∣ = dn/n
若将原始点的连线(即图中的虚线)近似看成是弧线,而相邻两个重采样点连线(即图中的实现)可近似看作是等弧长所对应的弦长,因此新采样点是近似等距离分布的,这将有助于得到光滑的拟合曲线。
2.2、轮廓关键点的确定
待拟合点列中的某些特定点会对逼近曲线形状产生重要影响,此处称之为轮廓关键点。而曲率的分布可以反映采样点列的总体和局部特征,本文采取基于曲率分布的轮廓关键点遴选方法。
2.2.1、曲率求解与分析
下面先求出采样点的曲率。离散点曲率半径的计算比较复杂,这里采用近似法求解,即第
i
i
i 个数据点处的曲率半径
ρ
i
\rho_i
ρi,可近似为经过 3 个相邻数据点
d
i
−
1
,
d
i
d_{i-1},d_i
di−1,di 和
d
i
+
1
d_{i+1}
di+1 的圆弧半径,其数值用式(5)求出:
ρ
i
=
∣
d
i
d
i
+
1
→
∣
∣
d
i
d
i
−
1
→
∣
∣
d
i
d
i
+
1
→
−
d
i
d
i
−
1
→
∣
2
∣
d
i
d
i
+
1
→
×
d
i
d
i
−
1
→
∣
(5)
\rho_i=\frac{\mid\overrightarrow{d_id_{i+1}}\mid\mid\overrightarrow{d_id_{i-1}}\mid\mid\overrightarrow{d_id_{i+1}}-\overrightarrow{d_id_{i-1}}\mid}{2\mid \overrightarrow{d_id_{i+1}}\times\overrightarrow{d_id_{i-1}}\mid}\tag{5}
ρi=2∣didi+1×didi−1∣∣didi+1∣∣didi−1∣∣didi+1−didi−1∣(5)
因此可得第 i i i 个数据点处得曲率 k i = 1 / r h o i k_i=1/rho_i ki=1/rhoi。
2.2.2、拟定轮廓关键点
不同于常见的曲率极值法,本文通过曲率临界值来确定轮廓关键点,即在分析曲率分布图的基础上,设定曲率临界值 [ k ] [k] [k],将曲率超过 [ k ] [k] [k] 的部分数据点初定为轮廓关键点。如果不是封闭曲线,则需同时将两个端点直接确定为关键点。
对于 [ k ] [k] [k] 值的选取,力求既能反映曲线的大致形状,又尽可能减少关键点数量。较小的 [ k ] [k] [k] 值将导致出现较多的控制顶点,而过大的 [ k ] [k] [k] 值则容易失去对曲线形状的控制。基于上述考虑,此处建议取曲率平均值 k a v g k_{\rm {avg}} kavg 作为临界值 [ k ] [k] [k]。
需要说明的是,一方面,由于数据点列轮廓通常会有凹凸变化,而外凸点和内凹点的曲率值较大,故此时以 [ k ] = k a v g [k]=k_{\rm{avg}} [k]=kavg 为条件确定的关键点数目通常小于数据点总数的一半,因而可有效压缩控制顶点数量。另一方面,由于 [ k ] [k] [k] 值的选取而造成的逼近误差,可以通过后续的局部优化得到补偿,因而曲线拟合的预设精度可以得到保证。
2.3、初始拟合 B 样条曲线
2.3.1、关键点的参数化
计算量最小的是均匀参数化,但该方法易造成曲线打圈相交的现象。另一种是弦长参数化法,但当数据点急转弯变化时该方法得到的曲线光滑性较差。为使曲线充分反映型值点的分布情况,此处采用向心参数化方法。对于
n
+
1
n+1
n+1 个关键点
Q
k
(
k
=
0
,
1
,
⋯
,
n
)
Q_k(k=0,1,\cdots,n)
Qk(k=0,1,⋯,n),记
d
=
∑
k
=
1
n
∣
Q
k
−
Q
k
−
1
∣
d=\sum_{k=1}^n\sqrt{\mid Q_k-Q_{k-1}\mid}
d=k=1∑n∣Qk−Qk−1∣
则
u
‾
0
=
0
,
u
‾
n
=
1
\overline{u}_0=0,\overline{u}_n=1
u0=0,un=1,而
u
‾
k
=
(
‾
u
)
k
−
1
+
∣
Q
k
−
Q
k
−
1
∣
/
d
,
k
=
1
,
2
,
⋯
,
n
−
1
(6)
\overline{u}_k=\overline(u)_{k-1}+\sqrt{\mid Q_k-Q_{k-1}\mid}/d\ ,\quad k=1,2,\cdots,n-1\tag{6}
uk=(u)k−1+∣Qk−Qk−1∣/d ,k=1,2,⋯,n−1(6)
2.3.3、节点矢量的计算
从理论上讲,节点矢量可以等距分布,即
u
0
=
u
1
=
⋯
=
u
p
=
0
,
u
m
−
p
=
u
m
−
(
p
−
1
)
=
⋯
=
u
m
=
1
,
u
j
+
p
=
1
/
(
n
−
p
+
1
)
,
j
=
1
,
2
,
⋯
,
n
−
p
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
u0=u1=⋯=up=0,um−p=um−(p−1)=⋯=um=1,uj+p=1/(n−p+1), j=1,2,⋯,n−p。但是如果讲将其与式(6)联合使用,方程组(4)可能变成奇异方程组。为避免这种情况出现,此处采用取平均值的方法,中间的节点矢量用式(7)求出:
u
j
+
p
=
1
p
∑
i
=
j
j
+
p
−
1
u
‾
i
,
j
=
1
,
2
,
⋯
,
n
−
p
(7)
u_{j+p}=\frac{1}{p}\sum_{i=j}^{j+p-1}\overline{u}_i\ ,\quad j=1,2,\cdots,n-p\tag{7}
uj+p=p1i=j∑j+p−1ui ,j=1,2,⋯,n−p(7)
一方面,由于式(6)中考虑了轮廓关键点对应的参数值 u ‾ i \overline{u}_i ui,可以反映关键点的分布;另一方面,联立式(6)和(7),可以确保方程组(4)的系数矩阵为正定的带状矩阵,当 ∣ i − k ∣ ≥ p \mid i-k\mid\geq p ∣i−k∣≥p 时, N i , p ( u ‾ k ) = 0 N_{i,p}(\overline{u}_k)=0 Ni,p(uk)=0。故可以采用不选主元的高斯消去法求解,减少了算法的计算量。
2.3.3、控制顶点的反算
根据式(2)和(3),通过参数 u ‾ k \overline{u}_k uk 和节点矢量 U \pmb U U 求出基函数 N i , p ( u ‾ k ) N_{i,p}(\overline{u}_k) Ni,p(uk),控制顶点 P i P_i Pi 是式(4)中唯一未知量。