常清俊，鄧重陽

(杭州電子科技大學(xué)理學(xué)院，浙江 杭州 310018)

0 引 言

曲線曲面擬合技術(shù)在許多領(lǐng)域扮演著重要的角色，廣泛應(yīng)用于計算機(jī)輔助設(shè)計、計算機(jī)圖形學(xué)、數(shù)據(jù)可視化、虛擬現(xiàn)實(shí)、數(shù)字幾何處理、數(shù)據(jù)挖掘等領(lǐng)域。在幾何建模等應(yīng)用中，曲線曲面的具體表達(dá)式可以構(gòu)造出曲線曲面的具體模型。但是，在工業(yè)應(yīng)用中，通常是先有具體的曲線曲面模型，對這些模型進(jìn)行數(shù)字化儲存或者顯示時，再構(gòu)造一個解析式來表達(dá)這樣的曲線曲面模型。一般情況下，選擇相對簡單的B樣條多項式作為解析式，因為B樣條具有幾何不變性、分段光滑和局部支撐等優(yōu)點(diǎn)[1]。用B樣條進(jìn)行曲線曲面擬合時，需要確定節(jié)點(diǎn)向量、控制頂點(diǎn)和參數(shù)，通常的做法是先確定節(jié)點(diǎn)向量和參數(shù)，然后求解控制頂點(diǎn)。一旦節(jié)點(diǎn)向量和參數(shù)確定，問題就變成一個求解控制頂點(diǎn)的線性問題，即求解一個線性方程組[2]。通常情況下，擬合數(shù)據(jù)點(diǎn)的數(shù)量多于控制頂點(diǎn)，故只能對數(shù)據(jù)點(diǎn)進(jìn)行逼近，所以采用最小二乘方法。最小二乘擬合是求解B樣條擬合的經(jīng)典方法，可以直接求解一個線性系統(tǒng)[3]。但是，當(dāng)擬合數(shù)據(jù)點(diǎn)的規(guī)模較大時，系數(shù)矩陣維數(shù)較高，需要求解一個大型的稀疏線性方程組，用直接求逆矩陣的方式進(jìn)行求解消耗大量計算資源，且效率較低，常用的做法是用迭代法求解線性方程組，例如雅可比迭代、高斯-塞德爾迭代，除此之外，還有Lin H.W.等[4-5]提出的基于漸進(jìn)逼近迭代法的B樣條曲線曲面擬合算法(Progressive and Iterative Approximation，PIA)和Deng C.Y.等[6]提出的最小二乘漸進(jìn)逼近迭代法(Least Square Progressive and Iterative Approximation，LSPIA)。由于傳統(tǒng)的雅可比迭代法和高斯-塞德爾迭代法每一次迭代過程只更新一個未知數(shù)的值，在一定程度上影響收斂速度。為此，本文通過使用分塊高斯-塞德爾迭代法求解方程組的解來加快迭代的收斂速度。

1 曲線曲面擬合問題

1.1 曲線擬合

(1)

式中，Bi(t)為B樣條基函數(shù)。令P=[P0P1…Pn]T為控制頂點(diǎn)，

將數(shù)據(jù)點(diǎn)列代入式(1)，寫成矩陣形式：

AP=Q

(2)

1.2 曲面擬合

(3)

式中，n1和n2分別表示曲面u方向和v方向上的控制頂點(diǎn)個數(shù)，令

將數(shù)據(jù)點(diǎn)列代入式(3)，寫成矩陣形式：

AP=Q

(4)

1.3 求解控制頂點(diǎn)

在進(jìn)行B樣條曲線曲面擬合時，由于數(shù)據(jù)點(diǎn)數(shù)量通常比控制頂點(diǎn)的數(shù)量多，式(2)和式(4)的系數(shù)矩陣A不一定為方陣，導(dǎo)致該線性方程組通常沒有精確解，故需采用最小二乘法求解一個較優(yōu)解。求解式(2)和式(4)的較優(yōu)解，轉(zhuǎn)換為求解線性方程組

ATAP=ATQ

(5)

的解。

由于B樣條基函數(shù)具有局部支撐性，所以當(dāng)配置矩陣維數(shù)較大時，A通常是稀疏矩陣。在求解大型稀疏線性系統(tǒng)時，常用的方法是迭代法，如高斯-塞德爾迭代法。

2 快速分塊迭代求解控制頂點(diǎn)

2.1 高斯-塞德爾迭代法

迭代法是求解大型線性方程組的常用方法。迭代法從初始估計開始，然后在每次迭代步驟中不斷細(xì)化估計，最后達(dá)到收斂，得到方程的解向量[7]。

(6)

從式(6)的迭代形式可以看出：在高斯-塞德爾迭代法中，解向量的更新是單獨(dú)進(jìn)行，每次只更新解向量中的一個元素，在一定程度上影響了收斂速度，為此，本文基于分塊的思想提出一種分塊迭代的方法——分塊高斯-塞德爾迭代法，使得每次更新解向量中的多個元素。

2.2 分塊高斯-塞德爾迭代法

線性方程組Ax=b，其中A為大型稀疏矩陣，將A,x,b分別分塊為：

(7)

式中，分塊矩陣的對角元素Aii為ni維非奇異方陣，xi∈Rni，bi∈Rni。則分塊高斯-塞德爾迭代法的分量形式為：

(8)

從式(8)可以看出：在分塊高斯-塞德爾迭代法中，每一步迭代更新多個解向量的元素，即x中的一個分塊。式(8)的迭代形式每次迭代需要多計算q+1個非奇異方陣的逆矩陣，為了克服反復(fù)計算逆矩陣的困難，在迭代開始之前預(yù)先計算出q+1個逆矩陣。分塊高斯-塞德爾迭代算法如下。

輸入:線性方程組Ax=b中的A和b,分塊矩陣的維數(shù)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 實(shí)驗結(jié)果與分析

實(shí)驗平臺的操作系統(tǒng)為Windows 10，測試工具為MATLAB R2019b，分別對10個測試數(shù)據(jù)集進(jìn)行擬合。分別使用高斯-塞德爾迭代法和本文提出的分塊高斯-塞德爾迭代法進(jìn)行計算得到控制頂點(diǎn)矩陣P，對于曲線擬合的情況，將P=[P0P1…Pn]T代入式(1)得到最終擬合曲線；對于曲面擬合情況，將P=[P00P01…Pij…Pn1n2]T代入式(3)中得到最終的擬合曲面。

高斯-塞德爾與分塊高斯塞德爾迭代法在迭代次數(shù)和迭代時間的比較結(jié)果如表1所示。迭代時間是多次測試的平均值，平均誤差表示為：

表1 2種迭代法的迭代次數(shù)和迭代時間比較

圖1 分塊高斯-塞德爾迭代法用于“天鵝”數(shù)據(jù)點(diǎn)擬合

通過表1可以看出：無論在時間上還是迭代次數(shù)上，分塊高斯-塞德爾迭代都比直接高斯-塞德爾迭代法的效率提高近1倍。

圖1展示了使用分塊高斯-塞德爾迭代法對“天鵝”原始數(shù)據(jù)點(diǎn)進(jìn)行三次B樣條曲線擬合的擬合結(jié)果，迭代40次后可達(dá)到給定的精度ε=1.0×10-13。

圖2展示了對“天鵝”數(shù)據(jù)使用2種迭代法的迭代過程，此處展示了前3次迭代的結(jié)果。實(shí)驗中，2種迭代法的控制頂點(diǎn)初始值設(shè)為零矩陣，參數(shù)化方法為弦長參數(shù)化，可以看出：分塊高斯-塞德爾可以很快達(dá)到收斂條件。在第1次迭代時，高斯-塞德爾迭代結(jié)果的幾何特征比較明顯，而分塊高斯-塞德爾迭代幾何特征不明顯，這是因為高斯-塞德爾每次迭代時只更新解向量的其中一個元素，而分塊高斯-塞德爾迭代每次更新解向量中的一塊元素。從圖2最后1列圖可以看出：分塊高斯-塞德爾迭代法在3次迭代后的效果明顯優(yōu)于未分塊的高斯-塞德爾迭代法。

圖2 2種迭代法擬合“天鵝”數(shù)據(jù)點(diǎn)的迭代過程

圖3展示了使用分塊高斯-塞德爾迭代法得到的其他曲線例子的擬合結(jié)果，其中黑色實(shí)心點(diǎn)為原始數(shù)據(jù)點(diǎn)，黑色曲線是擬合結(jié)果，可以看出：分塊高斯-塞德爾迭代法用于曲線擬合時具有較好的擬合結(jié)果。圖4展示了使用分塊高斯-塞德爾迭代法對Nefertiti人臉模型的數(shù)據(jù)點(diǎn)進(jìn)行曲面擬合的結(jié)果。在曲面擬合時，選用的參數(shù)化方法為Floater的均值坐標(biāo)參數(shù)化[8]。

例2

例3

例4

例5

例6

例7

例8

例9

圖4 分塊高斯-塞德爾迭代法曲面擬合結(jié)果

3 結(jié)束語

在高斯-塞德爾迭代法的基礎(chǔ)上，本文引入分塊的思想，將分塊高斯-塞德爾迭代法應(yīng)用到三次B樣條曲線曲面擬合中。與高斯-塞德爾迭代法相比，在收斂速度上有一定的提升。今后，將研究其他迭代法的分塊形式，并改善分塊方式使得迭代效率達(dá)到最優(yōu)。