B-Spline(11):样条曲线拟合-光顺逼近

曲线拟合包含两个方面,插值(interpolation)和逼近(approximation)。用于曲线拟合的离散点通常不具有非常高的精度,直接插值得到的曲线可能不满足“光顺(fair)”要求,本节的目标是介绍光顺的定义,以及给出一种满足“光顺”要求的最小二乘逼近方法。

1)光顺的定义

曲线光顺其实是NURBS曲线研究的一支非常重要的分支,作为一个非数学专业的学生,我们不能完全给出非常严谨的数学定义,仅能从应用方面介绍光顺的意义及应用。首先,我已经介绍了曲线的插值是令样条曲线严格的穿过所有的型值点(Fit Point),这样求出来的曲线,因为给定点的精度问题,通常具有曲率变化幅度特别大的问题。比如下图,虚线表示秋曲线的“曲率伴随曲线”,可以看到,曲线的曲率变化是非常剧烈的(但是,曲线的曲率仍然是连续的,不要混淆),这种现象,称之为曲线不光顺。

Poliakoff J F, Wong Y K, Thomas P D. An automated curve fairing algorithm for cubic B -spline curves[J]. Journal of Computational & Applied Mathematics, 1999, 102(1):73-85.

为什么会产生这样的问题呢?可以用下图解释,假设理想的曲线应该是圆(灰色),但是中间的黑点存在偏差的情况下,曲线对应的部分会有一个不必要的“弯折”,表现在曲率上,会出现这部分曲率剧烈变大的情况。

假设样条曲线是一条钢丝(就像最早的造船业使用的样条一样),被上面的四个黑色的“duck”约束,中间的黑点部分就产生了多余的弹性势能,因此,广泛的曲线“光顺”的定义就是曲线上没有多余的“能量”。

2)曲线最小二乘拟合的一般形式

因为曲线光顺优化是NURBS的一个重要的研究分支(至今仍然活跃),因此光顺算法也不一而足。其中,广泛应用的一种线性最小二乘拟合的方式一般是这样的:

最小二乘的目标函数由两部分构成,一部分为逼近项,约束曲线上的点与给定的型值点的距离不能太大。一部分为光顺项,或者称优化项,约束曲线达到“光顺要求”。目标函数的形式是:$$f=\alpha(fairing\_component)^2 + \beta(approximating\_component)^2$$两个项都是控制点P的函数,因此通过对目标函数求最小值的方式,可以计算出控制点P。至于系数\(\alpha,\beta\),属于光顺项与逼近项的权重,一般会通过实验事先确定。

在各种方法里,逼近项的形式一般是确定的,给定点\(\{F_1,…,F_x\}\),通过参数化确定他们的参数\(\{u_1,…,u_x\}\),对应曲线上的点为\(MP\),M为对应的基函数矩阵,P为控制点向量。那么逼近项则是$$MP-F=
\begin{bmatrix}
N_0(u_1)&\cdots&N_n(u_1)\\
\vdots&\ddots&\vdots\\
N_0(u_x)&\cdots&N_n(u_x)
\end{bmatrix}
\begin{bmatrix}
P_1\\\vdots\\P_n
\end{bmatrix}-
\begin{bmatrix}
F_1\\\vdots\\F_x
\end{bmatrix}$$

3)光顺项:B样条二次导数的积分的计算

至于光顺项,则有很多不同的实现方法,这里介绍一种 常见的方法:令曲线在给定点\(C(u_i)\)上二阶导数的平方和的积分最小的方法:$$\int_{t=0}^m(\sum_{j=0}^n N_{j,p}^{\prime\prime}
(u)P_j)^2du\rightarrow min$$为简化表达,令\(
W=\int_{t=0}^m(\sum_{j=0}^n N_{j,p}^{\prime\prime}
(u)P_j)^2du\),可以得出曲线拟合的目标函数\(f\)的形式为:$$f=\alpha P^TWP+\beta(MP-F)^2$$按照最小二乘法的解法:$$\frac{\partial f}{\partial P}=
2\alpha P^T(W)+2\beta(MP-F)^TM=0$$整理一下,得到$$P^T(\alpha(W)+\beta(M^TM))=\beta F^TM$$两边同时取转置,可以得到P的解:$$P=\beta(\alpha(W)+\beta(M^TM))^{-1}M^TF$$矩阵M,我已经在B样条的内插一节中涉及,这里不再赘述。比较复杂的,就是拼装矩阵W。

因为B样条的基函数\(N_{i,p}\)在区间\([u_i,u_{i+p+1}\)上非零,并且对于每一个区间\([u_k,u_{k+1})\),基函数都是关于u的三次多项式函数,那么其二次导数的形式为\(au+b\),因此我们需要求出每个区间上的a,b的值。

关于W矩阵的构造,可以构造一个n个基函数\(\times\)n+p个矩阵的辅助矩阵N,其中\(N(i,j)\)代表基函数导数\(N_{i,p}^{\prime\prime}\)在区间\([u_k,u_{k+1})\)上的多项式\(au+b\)的系数a,b。\(W_{i,j}\)的内容就是N第i行与第j行所有系数相乘后积分的结果:

$$W_{i,j}=\sum_{k=0}^{k+p}\int_{u_k}^{u_{k+1}}(N_{i,k}.a u+N_{i,k}.b)(N_{j,k}.a u+N_{j,k}.b)du$$

B样条基函数的二阶导数的递推形式如下:$$N_{i,p}^{\prime\prime}=p\cdot(p-1)\cdot\{
\frac{N_{i,p-2}}{(u_{i+p-1}-u_i)(u_{i+p}-u_i)}-
N_{i+1,p-2}\cdot(\frac{1}{(u_{i+p}-u_i)(u_{i+p}-u_{i+1})}+\frac{1}{(u_{i+p}-u_{i+1})(u_{i+p+1}-u_{i+1})})+\\
\frac{N_{i+2,p-2}}{(u_{i+p+1}-u_{i+1})(u_{i+p+1}-u_{i+2})}
\}$$

$$=c_1\frac{u-u_i}{u_{i+p-2}-u_i}N_{i,p-3}+
[c_1\frac{u_{i+p-1}-u}{u_{i+p-1}-u_{i+1}}+
c_2\frac{u-u_{i+1}}{u_{i+p-1}-u_{i+1}}]N_{i+1,p-3}+\\
[c_2\frac{u_{i+p}-u}{u_{i+p}-u_{i+2}}+
c_3\frac{u-u_{i+2}}{u_{i+p}-u_{i+2}}]N_{i+2,p-3}+
c_3\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+3}}N_{i+3,p-3}
$$

其中c_1,c_2,c_3分别对应上面的常数项。

这篇文章里,我仍然以最常用的三阶B样条曲线的拟合为例。所以对后面的式子里,\(p=3\),如果读者想要更改曲线阶次,需要对我特别针对三次曲线所做的说明进行调整。

当p=3时,\(N_{i,p}\)在区间\([u_i,u_{i+1}),[u_{i+1},u_{i+2}),[u_{i+2},u_{i+3}),[u_{i+3},u_{i+4})\)上对应的a,b的值与二阶导数展开在\(N_{i,p-3},N_{i+1,p-3},N_{i+2,p-3},N_{i+3,p-3}\)的系数一一对应。

当计算出每个区间上的a,b值后,就可以求出W矩阵了。自此,曲线整体光顺逼近算法就介绍完了。通过我的C++实现,可以看到具体的效果。

4)光顺效果

因为实现的代码较长,我们先看效果,把代码实现放到最后。

\(\alpha\)是光顺项的权重,\(\beta\)是逼近项的权重,通过调整两者的比例,可以使曲线获得不同的逼近效果。下图是\(\alpha=0.01,\beta=1\)时逼近曲线与插值曲线的曲率的对比,可见二者值域差别不大,逼近曲线的曲率变化稍微连续一点。二者几何差距不大。

显示曲率软件copyright navinfo

下图是增大光顺项权重\(\alpha=1,beta=1\)是逼近曲线的曲率与几何(红线为逼近曲线,青线为拟合曲线)的效果:

可见,增大光顺项权重后,曲线的曲率变化更加连续了,更符合“光顺”的效果了。下图是继续加大光顺处理效果\(\alpha=1,\beta=0.01\)的曲率与几何效果。可见,逼近曲线自动的平滑了弯曲较大的部分,从而趋近于减小能量。

5)二次型的求导

本节相对于之前的线性的矩阵计算,又引入了矩阵求导,尤其是二次型求导部分,所以,特地将引用的二次型求导公式列举如下:$$\frac{\partial(x^TAx)}{\partial x}=x^T(A+A^T)\\
\frac{\partial(x^TAx)}{\partial z}=x^T(A+A^T)\frac{\partial x}{\partial z}$$

5)我的C++实现

 


参考文献:

罗卫兰. B样条曲线的光顺[D]. 浙江大学, 2003.

兰浩. NURBS曲线整体光顺逼近算法研究与应用[D]. 西安理工大学, 2008.

Randal J. Barnes Matrix Differentiation

发表评论

电子邮件地址不会被公开。