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

1）光顺的定义

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

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

假设样条曲线是一条钢丝（就像最早的造船业使用的样条一样），被上面的四个黑色的“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\)时逼近曲线与插值曲线的曲率的对比，可见二者值域差别不大，逼近曲线的曲率变化稍微连续一点。二者几何差距不大。

下图是增大光顺项权重\(\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