分块高斯-塞德尔迭代的曲线曲面拟合

常清俊，邓重阳

(杭州电子科技大学理学院，浙江 杭州 310018)

0 引 言

曲线曲面拟合技术在许多领域扮演着重要的角色，广泛应用于计算机辅助设计、计算机图形学、数据可视化、虚拟现实、数字几何处理、数据挖掘等领域。在几何建模等应用中，曲线曲面的具体表达式可以构造出曲线曲面的具体模型。但是，在工业应用中，通常是先有具体的曲线曲面模型，对这些模型进行数字化储存或者显示时，再构造一个解析式来表达这样的曲线曲面模型。一般情况下，选择相对简单的B样条多项式作为解析式，因为B样条具有几何不变性、分段光滑和局部支撑等优点[1]。用B样条进行曲线曲面拟合时，需要确定节点向量、控制顶点和参数，通常的做法是先确定节点向量和参数，然后求解控制顶点。一旦节点向量和参数确定，问题就变成一个求解控制顶点的线性问题，即求解一个线性方程组[2]。通常情况下，拟合数据点的数量多于控制顶点，故只能对数据点进行逼近，所以采用最小二乘方法。最小二乘拟合是求解B样条拟合的经典方法，可以直接求解一个线性系统[3]。但是，当拟合数据点的规模较大时，系数矩阵维数较高，需要求解一个大型的稀疏线性方程组，用直接求逆矩阵的方式进行求解消耗大量计算资源，且效率较低，常用的做法是用迭代法求解线性方程组，例如雅可比迭代、高斯-塞德尔迭代，除此之外，还有Lin H.W.等[4-5]提出的基于渐进逼近迭代法的B样条曲线曲面拟合算法(Progressive and Iterative Approximation，PIA)和Deng C.Y.等[6]提出的最小二乘渐进逼近迭代法(Least Square Progressive and Iterative Approximation，LSPIA)。由于传统的雅可比迭代法和高斯-塞德尔迭代法每一次迭代过程只更新一个未知数的值，在一定程度上影响收敛速度。为此，本文通过使用分块高斯-塞德尔迭代法求解方程组的解来加快迭代的收敛速度。

1 曲线曲面拟合问题

1.1 曲线拟合

(1)

式中，Bi(t)为B样条基函数。令P=[P0P1…Pn]T为控制顶点，

将数据点列代入式(1)，写成矩阵形式：

AP=Q

(2)

1.2 曲面拟合

(3)

式中，n1和n2分别表示曲面u方向和v方向上的控制顶点个数，令

将数据点列代入式(3)，写成矩阵形式：

AP=Q

(4)

1.3 求解控制顶点

在进行B样条曲线曲面拟合时，由于数据点数量通常比控制顶点的数量多，式(2)和式(4)的系数矩阵A不一定为方阵，导致该线性方程组通常没有精确解，故需采用最小二乘法求解一个较优解。求解式(2)和式(4)的较优解，转换为求解线性方程组

ATAP=ATQ

(5)

的解。

由于B样条基函数具有局部支撑性，所以当配置矩阵维数较大时，A通常是稀疏矩阵。在求解大型稀疏线性系统时，常用的方法是迭代法，如高斯-塞德尔迭代法。

2 快速分块迭代求解控制顶点

2.1 高斯-塞德尔迭代法

迭代法是求解大型线性方程组的常用方法。迭代法从初始估计开始，然后在每次迭代步骤中不断细化估计，最后达到收敛，得到方程的解向量[7]。

(6)

从式(6)的迭代形式可以看出：在高斯-塞德尔迭代法中，解向量的更新是单独进行，每次只更新解向量中的一个元素，在一定程度上影响了收敛速度，为此，本文基于分块的思想提出一种分块迭代的方法——分块高斯-塞德尔迭代法，使得每次更新解向量中的多个元素。

2.2 分块高斯-塞德尔迭代法

线性方程组Ax=b，其中A为大型稀疏矩阵，将A,x,b分别分块为：

(7)

式中，分块矩阵的对角元素Aii为ni维非奇异方阵，xi∈Rni，bi∈Rni。则分块高斯-塞德尔迭代法的分量形式为：

(8)

从式(8)可以看出：在分块高斯-塞德尔迭代法中，每一步迭代更新多个解向量的元素，即x中的一个分块。式(8)的迭代形式每次迭代需要多计算q+1个非奇异方阵的逆矩阵，为了克服反复计算逆矩阵的困难，在迭代开始之前预先计算出q+1个逆矩阵。分块高斯-塞德尔迭代算法如下。

输入:线性方程组Ax=b中的A和b,分块矩阵的维数q,以及精度ε输出:线性方程组的解向量x将矩阵A和b按照式(7)要求分块,得到Aij和bi(i,j=0,1,…,q)计算分块矩阵A每个对角矩阵Aii的逆矩阵,记为Aii(i=0,1,…,q)初始化解向量x(0)← x(0)0,x(0)1,…,x(0)q T,k=0repeat计算x(k+1)i←Aii bi-∑i-1j=0Aijx(k+1)j-∑qj=i+1Aijx(k)j ,i=0,1,…,qk←k+1x(k)← x(k)0,x(k)1,…,x(k)q TuntilAx(k)-b2≤εx←x(k)

2.3 实验结果与分析

实验平台的操作系统为Windows 10，测试工具为MATLAB R2019b，分别对10个测试数据集进行拟合。分别使用高斯-塞德尔迭代法和本文提出的分块高斯-塞德尔迭代法进行计算得到控制顶点矩阵P，对于曲线拟合的情况，将P=[P0P1…Pn]T代入式(1)得到最终拟合曲线；对于曲面拟合情况，将P=[P00P01…Pij…Pn1n2]T代入式(3)中得到最终的拟合曲面。

高斯-塞德尔与分块高斯塞德尔迭代法在迭代次数和迭代时间的比较结果如表1所示。迭代时间是多次测试的平均值，平均误差表示为：

表1 2种迭代法的迭代次数和迭代时间比较

图1 分块高斯-塞德尔迭代法用于“天鹅”数据点拟合

通过表1可以看出：无论在时间上还是迭代次数上，分块高斯-塞德尔迭代都比直接高斯-塞德尔迭代法的效率提高近1倍。

图1展示了使用分块高斯-塞德尔迭代法对“天鹅”原始数据点进行三次B样条曲线拟合的拟合结果，迭代40次后可达到给定的精度ε=1.0×10-13。

图2展示了对“天鹅”数据使用2种迭代法的迭代过程，此处展示了前3次迭代的结果。实验中，2种迭代法的控制顶点初始值设为零矩阵，参数化方法为弦长参数化，可以看出：分块高斯-塞德尔可以很快达到收敛条件。在第1次迭代时，高斯-塞德尔迭代结果的几何特征比较明显，而分块高斯-塞德尔迭代几何特征不明显，这是因为高斯-塞德尔每次迭代时只更新解向量的其中一个元素，而分块高斯-塞德尔迭代每次更新解向量中的一块元素。从图2最后1列图可以看出：分块高斯-塞德尔迭代法在3次迭代后的效果明显优于未分块的高斯-塞德尔迭代法。

图2 2种迭代法拟合“天鹅”数据点的迭代过程

图3展示了使用分块高斯-塞德尔迭代法得到的其他曲线例子的拟合结果，其中黑色实心点为原始数据点，黑色曲线是拟合结果，可以看出：分块高斯-塞德尔迭代法用于曲线拟合时具有较好的拟合结果。图4展示了使用分块高斯-塞德尔迭代法对Nefertiti人脸模型的数据点进行曲面拟合的结果。在曲面拟合时，选用的参数化方法为Floater的均值坐标参数化[8]。

例2

例3

例4

例5

例6

例7

例8

例9

图4 分块高斯-塞德尔迭代法曲面拟合结果

3 结束语

在高斯-塞德尔迭代法的基础上，本文引入分块的思想，将分块高斯-塞德尔迭代法应用到三次B样条曲线曲面拟合中。与高斯-塞德尔迭代法相比，在收敛速度上有一定的提升。今后，将研究其他迭代法的分块形式，并改善分块方式使得迭代效率达到最优。