BSpline曲线逼近

在逼近问题中，曲线的误差界限 E E E和要拟合的数据一起被输入。一般预先并不知道需要多少个控制点才能达到预期的精度 E E E，因此逼近一般都需要通过迭代来实现。

给定一个确定的控制点数目，比如 n n n，构造曲线逼近给定的数据。有很多方法可以完成这一任务。例如，可建立一个非线性优化问题，以控制点、参数值 u ‾ k \overline{u}_k uk、节点甚至权值作为未知量，而要最小化的目标函数通常是某种形式的误差的度量，例如误差的平方和或最大偏差。在下面将介绍的曲线逼近方法中，只有固定数目的控制点是未知的，并且它们通过线性最小二乘技术来求解。

最小二乘曲线逼近

为了避免非线性问题，我们设定权值为1，并预先计算好数据点的参数值和节点矢量，然后，建立线性最小二乘问题来求解未知控制点。假定 p ⩾ 1 , n ⩾ p p\geqslant 1,n\geqslant p p⩾1,n⩾p，并且给定数据点 Q 0 , ⋯ , Q m ( m ⩾ n ) Q_0,\cdots,Q_m \left(m\geqslant n\right) Q0,⋯,Qm(m⩾n)。我们试图寻求一条 p p p次BSpline曲线
C ( u ) = ∑ i = 0 n N i , p ( u ) P i , u ∈ [ 0 , 1 ] C\left(u\right)=\sum^n_{i=0}N_{i,p}\left(u\right)P_i,\quad u\in\left[0,1\right] C(u)=i=0∑nNi,p(u)Pi,u∈[0,1]
满足条件

∙ Q 0 = C ( 0 ) = P 0 , Q m = C ( 1 ) = P n \bullet\quad Q_0=C(0)=P_0,\ Q_m=C(1)=P_n ∙Q0=C(0)=P0, Qm=C(1)=Pn；

∙ \bullet\quad ∙其余的数据点 Q k Q_k Qk在最小二乘的意义下被逼近，即
∑ k = 1 m − 1 ∣ Q k − C ( u ‾ k ) ∣ 2 \sum^{m-1}_{k=1}\left|Q_k-C(\overline{u}_k)\right|^2 k=1∑m−1∣Qk−C(uk)∣2
关于 n + 1 n+1 n+1个变量 P i P_i Pi达到最小； { u ‾ k } \{\overline{u}_k\} {uk}是预先计算好的参数值。

在此我们要强调，生成的曲线一般不精确地通过数据点 Q k Q_k Qk，并且，一般 C ( u ‾ k ) C(\overline{u}_k) C(uk)不是曲线上与 Q k Q_k Qk最接近的点。令
R k = Q k − N 0 , p ( u ‾ k ) Q 0 − N n , p ( u ‾ k ) Q m , k = 1 , ⋯ , m − 1 R_k=Q_k-N_{0,p}\left(\overline{u}_k\right)Q_0-N_{n,p}\left(\overline{u}_k\right)Q_m,\quad k=1,\cdots,m-1 Rk=Qk−N0,p(uk)Q0−Nn,p(uk)Qm,k=1,⋯,m−1
然后令
f = ∑ k = 1 m − 1 ∣ Q k − C ( u ‾ k ) ∣ 2 = ∑ k = 1 m − 1 ∣ R k − ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ∣ 2 = ∑ k = 1 m − 1 ( R k − ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) ( R k − ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) = ∑ k = 1 m − 1 [ R k ⋅ R k − 2 ∑ i = 1 n − 1 N i , p ( u ‾ k ) ( R k ⋅ P i ) + ( ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) ⋅ ( ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) ] \begin{aligned} f=& \sum^{m-1}_{k=1}|Q_k-C(\overline{u}_k)|^2 =\sum^{m-1}_{k=1}\left|R_k-\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right|^2 \\ =& \sum^{m-1}_{k=1}\left(R_k-\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right) \left(R_k-\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right) \\ =& \sum^{m-1}_{k=1}\left[R_k\cdot R_k -2\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)(R_k\cdot P_i) +\left(\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right) \cdot\left(\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right)\right] \end{aligned} f===k=1∑m−1∣Qk−C(uk)∣2=k=1∑m−1∣∣∣∣∣Rk−i=1∑n−1Ni,p(uk)Pi∣∣∣∣∣2k=1∑m−1(Rk−i=1∑n−1Ni,p(uk)Pi)(Rk−i=1∑n−1Ni,p(uk)Pi)k=1∑m−1[Rk⋅Rk−2i=1∑n−1Ni,p(uk)(Rk⋅Pi)+(i=1∑n−1Ni,p(uk)Pi)⋅(i=1∑n−1Ni,p(uk)Pi)]
f f f是关于 n − 1 n-1 n−1个变量 P 1 , ⋯ , P n − 1 P_1,\cdots,P_{n-1} P1,⋯,Pn−1的标量值函数。现我们应用标准的线性最小二乘拟合技术，为使目标函数 f f f最小，我们令 f f f关于 n − 1 n-1 n−1个未知控制点 P l P_l Pl的偏导都等于零。它的第 l l l 个偏导为
∂ f ∂ P l = ∑ k = 1 m − 1 ( − 2 N l , p ( u ‾ k ) R k + 2 N l , p ( u ‾ k ) ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) \dfrac{\partial f}{\partial P_l} =\sum^{m-1}_{k=1}\left(-2N_{l,p}(\overline{u}_k)R_k +2N_{l,p}(\overline{u}_k)\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i\right) ∂Pl∂f=k=1∑m−1(−2Nl,p(uk)Rk+2Nl,p(uk)i=1∑n−1Ni,p(uk)Pi)
这意味着
− ∑ k = 1 m − 1 N l , p ( u ‾ k ) R k + ∑ k = 1 m − 1 ∑ i = 1 n − 1 N l , p ( u ‾ k ) N i , p ( u ‾ k ) P i = 0 -\sum^{m-1}_{k=1}N_{l,p}(\overline{u}_k)R_k +\sum^{m-1}_{k=1}\sum^{n-1}_{i=1}N_{l,p}(\overline{u}_k)N_{i,p}(\overline{u}_k)P_i =0 −k=1∑m−1Nl,p(uk)Rk+k=1∑m−1i=1∑n−1Nl,p(uk)Ni,p(uk)Pi=0
于是
∑ i = 1 n − 1 ( ∑ k = 1 m − 1 N l , p ( u ‾ k ) N i , p ( u ‾ k ) ) P i = ∑ k = 1 m − 1 N l , p ( u ‾ k ) R k \sum^{n-1}_{i=1}\left(\sum^{m-1}_{k=1}N_{l,p}(\overline{u}_k)N_{i,p}(\overline{u}_k)\right)P_i =\sum^{m-1}_{k=1}N_{l,p}(\overline{u}_k)R_k i=1∑n−1(k=1∑m−1Nl,p(uk)Ni,p(uk))Pi=k=1∑m−1Nl,p(uk)Rk
这给出了一个以控制顶点 P 1 , ⋯ , P n − 1 P_1,\cdots,P_{n-1} P1,⋯,Pn−1为未知量的线性方程。让 l = 1 , 2 , ⋯ , n − 1 l=1,2,\cdots,n-1 l=1,2,⋯,n−1，则得到含 n − 1 n-1 n−1个未知量和 n − 1 n-1 n−1个方程的线性方程组
( N T N ) P = R \left(N^TN\right)P=R (NTN)P=R
这里， N N N是由标量组成的 ( m − 1 ) × ( n − 1 ) \left(m-1\right)\times\left(n-1\right) (m−1)×(n−1)的矩阵
N = [ N 1 , p ( u ‾ 1 ) ⋯ N n − 1 , p ( u ‾ 1 ) ⋮ ⋮ N 1 , p ( u ‾ m − 1 ) ⋯ N n − 1 , p ( u ‾ m − 1 ) ] \begin{aligned} N= \begin{bmatrix} N_{1,p}(\overline{u}_1) & \cdots & N_{n-1,p}(\overline{u}_1) \\ \vdots & & \vdots \\ N_{1,p}(\overline{u}_{m-1}) & \cdots & N_{n-1,p}(\overline{u}_{m-1}) \end{bmatrix} \end{aligned} N=⎣⎢⎡N1,p(u1)⋮N1,p(um−1)⋯⋯Nn−1,p(u1)⋮Nn−1,p(um−1)⎦⎥⎤
R R R是由 n − 1 n-1 n−1个点组成的列向量：
R = [ N 1 , p ( u ‾ 1 ) R 1 + ⋯ + N 1 , p ( u ‾ m − 1 ) R m − 1 ⋮ N n − 1 , p ( u ‾ 1 ) R 1 + ⋯ + N n − 1 , p ( u ‾ m − 1 ) R m − 1 ] \begin{aligned} R= \begin{bmatrix} N_{1,p}(\overline{u}_1)R_1 +\cdots +N_{1,p}(\overline{u}_{m-1})R_{m-1} \\ \vdots \\ N_{n-1,p}(\overline{u}_1)R_1 +\cdots +N_{n-1,p}(\overline{u}_{m-1})R_{m-1} \end{bmatrix} \end{aligned} R=⎣⎢⎡N1,p(u1)R1+⋯+N1,p(um−1)Rm−1⋮Nn−1,p(u1)R1+⋯+Nn−1,p(um−1)Rm−1⎦⎥⎤
并且
P = [ P 1 ⋮ P n − 1 ] \begin{aligned} P= \begin{bmatrix} P_1 \\ \vdots \\ P_{n-1} \end{bmatrix} \end{aligned} P=⎣⎢⎡P1⋮Pn−1⎦⎥⎤
为了求解线性方程组计算控制点，需要知道各数据点的参数值 { u ‾ k } \{\overline{u}_k\} {uk}和节点矢量 U = { u 0 , u 1 , ⋯ , u n + p + 1 } U=\{u_0,u_1,\cdots,u_{n+p+1}\} U={u0,u1,⋯,un+p+1}。

数据点的参数化

弦长参数化法计算 { u ‾ k } \{\overline{u}_k\} {uk}：令 d d d为总弦长
d = ∑ k = 1 m ∣ Q k − Q k − 1 ∣ d=\sum^m_{k=1}|Q_k-Q_{k-1}| d=k=1∑m∣Qk−Qk−1∣
则 u ‾ 0 = 0 , u ‾ m = 1 \overline{u}_0=0,\ \overline{u}_m=1 u0=0, um=1
u ‾ k = u ‾ k − 1 + ∣ Q k − Q k − 1 ∣ d , k = 1 , 2 , ⋯ , m − 1 \overline{u}_k =\overline{u}_{k-1}+\dfrac{|Q_k-Q_{k-1}|}{d},\quad k=1,2,\cdots,m-1 uk=uk−1+d∣Qk−Qk−1∣,k=1,2,⋯,m−1

节点配置

根据端点插值与曲线定义域要求，节点矢量 U = [ u 0 , u 1 , ⋯ , u n + p + 1 ] U=\left[u_0,u_1,\cdots,u_{n+p+1}\right] U=[u0,u1,⋯,un+p+1] 常采用定义域两端节点为 p + 1 p+1 p+1重的重节点端点条件，也即固支条件。于是有 u 0 = u 1 = ⋯ = u p = 0 u_0=u_1=\cdots=u_p=0 u0=u1=⋯=up=0和 u n + 1 = u n + 2 = ⋯ = u n + p + 1 = 1 u_{n+1}=u_{n+2}=\cdots=u_{n+p+1}=1 un+1=un+2=⋯=un+p+1=1。现在的问题是怎样确定所含内节点数和内节点值，也即怎样进行节点配置以确定节点矢量。

节点的选择对生成的逼近曲线的形状有很大的影响，不合理的节点矢量可能导致不可接受的形状。因此，BSpline曲线逼近成功的关键在于尽可能给出合适的节点配置。皮格尔和蒂勒推荐采用如下的方法：
u p + j = { 1 p ∑ i = j j + p − 1 u ‾ i , m = n ( 1 − α ) u ‾ i − 1 + α u ‾ i , m > n ( j = 1 , 2 , ⋯ , n − p ) 其中： i = i n t ( j c ) , c = m + 1 n − p + 1 , α = j c − i \begin{aligned} u_{p+j}= \begin{cases} \dfrac{1}{p}\sum\limits^{j+p-1}_{i=j}\overline{u}_i,\quad m=n \\ (1-\alpha)\overline{u}_{i-1}+\alpha\overline{u}_i,\quad m>n\ (j=1,2,\cdots,n-p) \\ 其中：i=int(jc),\ c=\dfrac{m+1}{n-p+1},\ \alpha=jc-i \end{cases} \end{aligned} up+j=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧p1i=j∑j+p−1ui,m=n(1−α)ui−1+αui,m>n (j=1,2,⋯,n−p)其中：i=int(jc), c=n−p+1m+1, α=jc−i
此即用于 m = n m=n m=n时的平均技术即AVG技术和用于 m > n m>n m>n时的节点配置技术 ( ( (the knot placement technique，KTP技术 ) ) )。可见这一套方法包括插值 ( m = n ) (m=n) (m=n)和逼近 ( m > n ) (m>n) (m>n)两种情况，这里把插值 ( m = n ) (m=n) (m=n)作为特殊情况进行特殊处理，内节点值取顺序 p p p个参数值的算术平均。而在逼近 ( m > n ) (m>n) (m>n)情况下，采用另一种不同的平均方法即 ( m + 1 ) / ( n − p + 1 ) (m+1)/(n-p+1) (m+1)/(n−p+1)和据此决定的 α \alpha α对顺序两参数值的加权平均，据此决定的内节点值保证了定义域每个节点区间包含至少一个 u ‾ \overline{u} u。德布尔表明在这样确定内节点的条件下，矩阵 ( N T N ) \left(N^TN\right) (NTN)对角占优，有一个小于 p + 1 p+1 p+1的带宽。 ( N T N ) \left(N^TN\right) (NTN)是正定的和情况良好的，这给出了一个稳定方程组，可由高斯消元法求解。