秦文杰，張舉勇

基于安德森加速的快速B樣條擬合算法

秦文杰，張舉勇

(中國科學技術(shù)大學數(shù)學科學學院，安徽 合肥 230026)

曲線擬合技術(shù)已被廣泛地應用于圖像處理、工程實驗等領(lǐng)域。其中，B樣條曲線擬合是曲線擬合中最常見的方法，它具有局部性好、連續(xù)性好等優(yōu)點，但擬合精度一般較低。在實際應用中，B樣條曲線擬合對于精度和速度的要求都較高。為了提升平面B樣條曲線擬合速度，將安德森加速的想法應用到曲線擬合的方法之中，提出一種基于安德森加速的擬牛頓方法。首先設定一個初始形狀，然后根據(jù)初始形狀找到其每個數(shù)據(jù)點的投影點的位置參數(shù)，然后利用安德森加速計算出控制點的相應位置，迭代進行以上2步，直到結(jié)果收斂。實驗結(jié)果表明，該方法在收斂速度和迭代時間上均優(yōu)于其他方法。

B樣條擬合；安德森加速；擬牛頓方法；曲線擬合；樣條逼近

曲線擬合問題一直是計算機圖形學中重要的組成部分，而隨著計算機科學與工程技術(shù)的發(fā)展，B樣條擬合問題更成為曲線擬合中的熱點問題。在實際的工程應用中，一般獲取到的初始數(shù)據(jù)往往是無序的點云數(shù)據(jù)，因此，如何擬合無序數(shù)據(jù)點是一個重要的問題。本文主要討論如何對平面無序數(shù)據(jù)點進行快速的B樣條曲線擬合的問題。

其中，(,t)為為在擬合曲線上距離點X最近的點，也稱作投影點。

1 相關(guān)工作

1.1 經(jīng)典的B樣條曲線擬合算法

對于絕大多數(shù)B樣條曲線擬合方法，迭代的處理方式均采用了分別計算投影點和控制點的想法：

(1) 投影點計算，即固定當前的曲線控制頂點不變，計算出數(shù)據(jù)點在當前曲線上的投影點。在計算投影點時往往會受到正交性的約束。

(2) 控制點更新，在此步驟中投影點的位置參數(shù)是固定的，可用一個二次函數(shù)E來近似數(shù)據(jù)點與擬合曲線之間的距離關(guān)系，則優(yōu)化該二次目標函數(shù)的能量就可以得到對應的新的控制點位置，即求解一個線性方程組。而相對應的，不同的二次函數(shù)E的取法其實也對應了不同的方法。

對應的是SDM方法[6]。其中，為點(t)處的曲率半徑；為X與(t)間的距離值，以正負號來表示其方向。從能量定義的方式可以看出，該方法包含了二階導數(shù)的信息與曲率信息，因此，與前2種方法相比，其更接近于真實的距離平方函數(shù)。同時，在優(yōu)化方面，SDM方法是一種擬牛頓方法，其修正了Hessian矩陣，使之保留了與幾何性質(zhì)相關(guān)的簡單部分，同時在修正的過程中保證了Hessian矩陣的半正定性。文獻[6]證明了SDM方法在收斂速度和穩(wěn)定性方面明顯優(yōu)于PDM和TDM方法。

L-BFGS方法[7]更多地將注意力放在了優(yōu)化求解速度上。其與傳統(tǒng)方法的區(qū)別在于，在迭代求解的過程中，可通過求出近似的逆Hessian矩陣而避免了矩陣計算，并同時優(yōu)化了投影點與控制點的信息。因此，該方法在收斂速度上要遠遠優(yōu)于其他3種方法。

SDM方法和L-BFGS方法均屬擬牛頓法，但存在一些區(qū)別。如L-BFGS方法能夠在每步迭代當中同時優(yōu)化控制點位置與投影點參數(shù)，并且不需要求解任何線性方程組，因此該方法比其他方法效率更高。

在傳統(tǒng)的B樣條曲線擬合方法中，PDM方法收斂速度較慢；TDM方法本質(zhì)上是高斯-牛頓法，因此具有不穩(wěn)定的特點。SDM方法是一種擬牛頓法，利用近似的Hessian矩陣代替真實的Hessian矩陣，在穩(wěn)定性和收斂速度之間找到了很好的平衡；L-BFGS方法也是一種擬牛頓方法，通過每一步推導近似的逆Hessian矩陣來達到自己優(yōu)化加速的目的。

1.2 局部-全局求解器

由于局部-全局求解器在近似求解最優(yōu)化問題上的高效表現(xiàn)，其在圖形學的各個領(lǐng)域均有廣泛地應用。文獻[8]提出了ARAP方法，即通過最小化目標能量來進行曲面建模；文獻[9]也應用了該方法來進行保形或保角的三角網(wǎng)格參數(shù)化工作；文獻[10]也使用了局部-全局求解器的框架對彈簧進行了快速并且接近真實的模擬。

與此同時，針對局部-全局求解器的結(jié)果收斂到高精度時耗時過長的問題，許多改進工作被陸續(xù)提出。SORKINE和RABINOVICH[11]采用Cheybyshev方法來實現(xiàn)收斂的加速；文獻[7]提出了L-BFGS方法，速度比普通的局部-全局求解器更快。本文工作重點在于優(yōu)化局部-全局求解器。

1.3 安德森加速

安德森加速被廣泛地應用于各個領(lǐng)域，該方法能夠加快求解迭代問題的收斂速度。其最初由ANDERSON[12]提出，隨后PULAY[13]在量子化學的計算中引入了同樣的技術(shù)。近年來，越來越多的研究者開始關(guān)注該方法在數(shù)值計算中的應用。文獻[14]用安德森加速來加速求解平流擴散問題，文獻[15]用其來求解低秩張量近似的問題，文獻[16]用來加速幾何優(yōu)化與物理模擬問題的求解。

文獻[16]與本文的工作想法類似，但并不完全相同。前者提出了一種統(tǒng)一的、全新的安德森加速的應用思路，即對于固定點迭代問題，用安德森加速進行優(yōu)化，并舉了幾何優(yōu)化與物理模擬 2個例子進行說明。但其求解的問題往往是單變量的問題，本文通過引入輔助變量，使得問題能夠用局部-全局求解器迭代求解。由于輔助變量由人為引入，所以輔助變量與待求的變量之間關(guān)系簡單，能夠很容易將問題改寫為固定點迭代的格式。而在B樣條曲線擬合的問題中，投影點位置參數(shù)與控制點之間具有較復雜的關(guān)系，也是2個工作的區(qū)別所在。

2 算法

2.1 迭代求解的一般方法

已經(jīng)知道B樣條曲線擬合問題可以化為以下最優(yōu)化問題

一般分為2個步驟：

步驟1.投影點計算，即保持當前擬合曲線的控制頂點不變，計算出關(guān)于數(shù)據(jù)點{X}的對應位置參數(shù)={t}，使得點序列{(t)}均為當前擬合曲線上{X}的投影點。該步驟需要滿足

步驟2.控制點更新，位置參數(shù)={t}被固定，即B樣條的系數(shù)矩陣

2.2 安德森加速

為了加速迭代過程的收斂，首先要意識到在投影點計算的步驟中，待求的位置參數(shù)={t}所對應的矩陣可以表示為當前控制頂點{(t)}的函數(shù)，所以將投影點計算與控制點更新結(jié)合起來，可以將該問題看成固定點迭代問題，即

對于此問題而言

對于固定點迭代問題，其殘差為

在第步迭代中，本文得到了以上結(jié)果，那么對于第+1步的結(jié)果，有

其中，為混合參數(shù)，在絕大部分文獻中，通常有=1，本文也依此取值。

需要注意的是，最小二乘問題也可以寫成如下形式

在求解最小二乘問題時，需要在每一步迭代中都重新計算所有的殘差，明顯影響效率，因此可以采用文獻[17]的做法，通過構(gòu)造法向方程使得效率變高。

2.3 穩(wěn)定性保證

若存在迭代起始位置靠近真解的情況時，文獻[18]已經(jīng)證明安德森加速在弱條件下收斂。而與此同時，也存在起始位置距離最后的真解很遠的情況，在該情況下，安德森加速可能會變得不穩(wěn)定，可能會出現(xiàn)收斂速度緩慢或是收斂到局部解的情況。

2.4 m的選取

根據(jù)以上對于算法的解釋可以知道，用于加速的先前的迭代次數(shù)對于算法的性能有著較大的影響。當取值較大時，更多的信息被用于近似逆Jacobian矩陣，通常會使收斂速度變快。但另一方面，過大的也會造成每一步迭代的計算成本增加，并有可能導致過擬合。本文經(jīng)過試驗發(fā)現(xiàn)，超過6以后對收斂速度影響有限，與文獻[16]的經(jīng)驗一致。在本文中，取5。具體試驗方式為：對于若干例子，將值取1~6進行試驗，得到如圖1所示的對比，可以直觀看出的取值對收斂速度的影響。事實上，是一個需試驗得到的相對最佳值。

圖1 m不同取值對比圖

2.5 算法流程

步驟2.對于每個數(shù)據(jù)點X找到其投影點的位置參數(shù)，從而確定出初始的位置參數(shù)。

步驟3.計算P+1，其中=1,2,...。

步驟3.2.比較目標能量，根據(jù)結(jié)果修改相應P，值，修改原則可見2.3節(jié)。

王鶴鵬[4]等在基于VR技術(shù)的汽車拆裝實訓課程教學改革探索中認為，VR技術(shù)的應用可以有效解決實訓場地要求較大，試驗材料及教學成本較高且存在一定事故隱患等問題。蔣斌[5]在發(fā)動機拆裝實訓課程中引入翻轉(zhuǎn)課堂的教學方法。李軍功[6]在發(fā)動機構(gòu)造與拆裝實訓課程的教學中采用項目教學的方法。

步驟3.3. 計算P，求解相應的最小二乘問題，求出P+1。

步驟4.重復步驟3直到能量收斂。

在實際應用中，初始曲線一般由人為確定，本文中采用圓形、橢圓等形狀。對于步驟3中計算投影點的–步驟，參考了文獻[7, 19]。需要注意的是，在本文的模型中，并未引入相應的光滑項。

3 實驗結(jié)果及分析

3.1 實驗的相關(guān)定義與參數(shù)

首先對每一步的擬合誤差給出定義，對于第步，擬合誤差為

該函數(shù)是在穩(wěn)定性保證過程中需要比較的目標函數(shù)。

在本文的實例中，B樣條曲線的參數(shù)域為=[0,1]，所有曲線擬合數(shù)據(jù)點都會被縮放至包圍盒[0,1]×[0,1]中。

同時，由于文獻[7]中已經(jīng)有了L-BFGS方法與其他經(jīng)典方法的完整比較，因此在本文中，只與L-BFGS方法作對比。L-BFGS方法的值取5，與文獻[16]保持一致，同時，不添加光滑項。本文方法與L-BFGS方法均采用C++編寫，L-BFGS方法的具體實現(xiàn)可參考文獻[16]。

同時，所有操作均在如下配置下進行：Intel? Core? i5-8050U @ 1.60 GHz，1.80 GHz，8.00 GB。

3.2 實驗對比

對比3個實例，并展示最終結(jié)果、誤差-時間關(guān)系圖、誤差-迭代步數(shù)關(guān)系圖(圖2~4)，及統(tǒng)計最終收斂的總時間與每一步迭代需要花費的時間 (表1~3)。

圖2 擬合曲線實例1 (8個控制頂點，100個數(shù)據(jù)點)

表1 實例1收斂總時間及迭代平均時間(s)

圖3 擬合曲線實例2 (20個控制點，90個數(shù)據(jù)點，包含尖銳特征)

表2 實例2收斂總時間及迭代平均時間(s)

由3個例子可以看出，在最初的幾步迭代過程中，L-BFGS與本文方法的效果接近，原因是這 2種方法都需要前步的結(jié)果用來優(yōu)化近似之后的結(jié)果，在最初的迭代中2種方法都沒有使用各自的優(yōu)化思路。在之后的迭代中，安德森加速的收斂速度要快于L-BFGS方法，兩二者在每一步迭代上所花費的時間基本保持在同一量級。

圖4 擬合曲線實例3 (8個控制頂點，200個數(shù)據(jù)點，存在噪音數(shù)據(jù)點)

表3 實例3收斂總時間及迭代平均時間(s)

同時，從收斂效果上來看，擬合誤差是衡量收斂效果的指標，在相同的收斂時間中，本文方法的擬合誤差結(jié)果總是小于L-BFGS方法。

3.3 單步迭代時間穩(wěn)定性

3個實例的數(shù)據(jù)點的個數(shù)都在102量級，為了更好地比較算法的穩(wěn)定性，更希望看到在較大數(shù)據(jù)點個數(shù)下L-BFGS算法與本文算法的表現(xiàn)。因此，對于實例1進行了加點處理。加點原則為：對于每一個已有的數(shù)據(jù)點，在以其為中心、大小為[0,]×[0,]的包圍盒內(nèi)，隨機取個點作為新的數(shù)據(jù)點；取=5×10–2。

從圖5和圖6中不難看出，對于L-BFGS算法和本文方法，在誤差允許的范圍內(nèi)，每一步的迭代平均時間與數(shù)據(jù)點的個數(shù)成線性關(guān)系，且每一步的迭代時間與控制點的個數(shù)也成線性關(guān)系。由文獻[7]可知，傳統(tǒng)的擬合方法如PDM，TDM和SDM方法，其平均迭代時間隨著控制點和數(shù)據(jù)點數(shù)目的增大而快速增長。

圖5 擬合曲線實例4 (8個控制點，數(shù)據(jù)點個數(shù)依次為100，200，400，800，1 600)

圖6 擬合曲線實例5 (100個數(shù)據(jù)點，控制點個數(shù)依次為8，28，48，68，88)

因此，L-BFGS方法與安德森加速方法都有著更適用于大型擬合問題的優(yōu)點，尤其是控制點較多并且數(shù)據(jù)點個數(shù)較大時，相比其他方法，這2種方法有著較大的優(yōu)勢。

3.4 初值穩(wěn)定性

算法所取的初始曲線對于算法的效率所造成的影響也是本文所關(guān)心的問題。因此，對于實例1，取3種不同的初始曲線進行測試，如圖7和圖8所示。

圖7 3種不同初始曲線

圖8 3種不同初始曲線收斂情況

從圖8可以看出，不同的初始曲線對于最終收斂情況的影響較小。

3.5 算法分析

從上文實例可以看出，和L-BFGS方法一樣，安德森加速方法的計算速度要明顯優(yōu)于一階的方法，同時，也比L-BFGS方法快一些。

文獻[16]和文獻[19]已證明了安德森加速是一種用于求解以下非線性方程的擬牛頓方法，即

等價于

其中，G為在P處的逆Jacobian矩陣的近似。

類似的優(yōu)化策略在L-BFGS方法中也可以看到，因此L-BFGS方法與安德森加速方法都明顯優(yōu)于一階方法。在L-BFGS方法中，先為L-BFGS求解器構(gòu)建初始的下降方向。其求解器的目標是搜索方程()=0的根，其中，為目標能量函數(shù)的梯度。在每一個迭代過程中，均需先計算下降方向并進行線搜索來確定下一次迭代。

為了計算下降方向，L-BFGS方法從初始估計-(P)開始計算，其中，是目標函數(shù)的逆Hessian矩陣的近似，然后用two-loop recursion算法的步驟來更新出新的逆Hessian矩陣的近似與下降方向。

從以上解釋可以看出，L-BFGS方法構(gòu)造初始逆Hessian矩陣的思路與本文在求解過程中構(gòu)造逆Jacobian矩陣的思路類似。2種方法的區(qū)別就在于構(gòu)造逆矩陣的方式的不同。在L-BFGS方法中，是基于two-loop recursion算法通過求解次割線方程來更新相應逆矩陣，而本文方法中并沒有2步的循環(huán)求解的操作。因此，L-BFGS方法在每一步中會有3次點積的操作，而本文方法只需要次，在計算效率上略優(yōu)于L-BFGS 方法。

4 結(jié)束語

本文提出了一種利用安德森加速來進行曲線擬合的方法。該方法本質(zhì)上是一種求解非線性系統(tǒng)的擬牛頓方法，與傳統(tǒng)曲線擬合方法(如TDM，PDM等)相比，不需要在每一步迭代過程中花費大量時間重新進行投影點計算，因而節(jié)約了大量的時間。與同為擬牛頓方法的L-BFGS相比，由于在每一步迭代過程中加入的目標能量判斷，使得該方法擁有更快的收斂速度。

本文的方法也有一些不足，比如沒有加入類似SDM方法中的光滑項，原因是積分形式的光滑項約束很難寫成穩(wěn)定的矩陣形式參與迭代過程。因此，在今后的工作中，期望可以添加合適的光滑項，使得結(jié)果更加魯棒，同時，也考慮將該方法應用到曲面擬合的優(yōu)化當中。

[1] HOSCHEK J. Intrinsic parametrization for approximation[J]. Computer Aided Geometric Design, 1988, 5(1): 27-31.

[2] PAVLIDIS T. Curve fitting with conic splines[J]. ACM Transactions on Graphics, 1983, 2(1): 1-31.

[3] PLASS M, STONE M. Curve-fitting with piecewise parametric cubics[C]//Proceedings of the 10th Annual Conference on Computer Graphics and Interactive Techniques-SIGGRAPH’83. New York: ACM Press, 1983: 229-239.

[4] SAUX E, DANIEL M. An improved Hoschek intrinsic parametrization[J]. Computer Aided Geometric Design, 2003, 20(8/9): 513-521.

[5] BLAKE A, ISARD M. Active contours[EB/OL]. [2020-03-29]. https://link.springer.com/book/10.1007/ 978-1-4471-1555-7.

[6] WANG W P, POTTMANN H, LIU Y. Fitting B-spline curves to point clouds by curvature-based squared distance minimization[J]. ACM Transactions on Graphics, 2006, 25(2): 214-238.

[7] ZHENG W N, BO P B, LIU Y, et al. Fast B-spline curve fitting by L-BFGS[J]. Computer Aided Geometric Design, 2012, 29(7): 448-462.

[8] SORKINE O, MARC A. As-rigid-as-possible surface modeling[C]//Proceedings of Eurographics/ACM Siggraph Symposium on Geometry Processing. New York: ACM Press, 2007: 109-116.

[9] LIU L G, ZHANG L, XU Y, et al. A local/global approach to mesh parameterization[J]. Computer Graphics Forum, 2008, 27(5): 1495-1504.

[10] LIU T T, BARGTEIL A W, O’BRIEN J F, et al. Fast simulation of mass-spring systems[J]. ACM Transactions on Graphics, 2013, 32(6): 214:1-214:7.

[11] SORKINE O, RABINOVICH M. Least-squares rigid motion using SVD[EB/OL]. [2020-03-29]. https://www.mendeley.com/catalogue/498963b1-7979-3736-915b-85ff25716ed4/.

[12] ANDERSON D G. Iterative procedures for nonlinear integral equations[J]. Journal of the ACM, 1965, 12(4): 547-560.

[13] PULAY P. Convergence acceleration of iterative sequences the case of SCF iteration[J]. Chemical Physics Letters, 1980, 73(2): 393-398.

[14] LIPNIKOV K, SVYATSKIY D, VASSILEVSKI Y. Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems[J]. SIAM Journal on Scientific Computing, 2013, 35(2): A1120-A1136.

[15] DE STERCK H. A nonlinear GMRES optimization algorithm for canonical tensor decomposition[J]. SIAM Journal on Scientific Computing, 2011, 34(3): A1351-A1379.

[16] PENG Y, DENG B L, ZHANG J Y, et al. Anderson acceleration for geometry optimization and physics simulation[J]. ACM Transactions on Graphics, 2018, 37(4): 1-14.

[17] FANG H R, SAAD Y. Two classes of multisecant methods for nonlinear acceleration[J]. Numerical Linear Algebra with Applications, 2009, 16(3): 197-221.

[18] TOTH A, ELLIS J A, EVANS T, et al. Local improvement results for Anderson acceleration with inaccurate function evaluations[J]. SIAM Journal on Scientific Computing, 2017, 39(5): S47-S65.

[19] HU S M, WALLNER J. A second order algorithm for orthogonal projection onto curves and surfaces[J]. Computer Aided Geometric Design, 2005, 22(3): 251-260.

Anderson acceleration for B-spline curve fitting

QIN Wen-jie, ZHANG Ju-yong

(School of Mathematical Sciences of University of Science and Technology of China, Hefei Anhui 230026, China)

In recent years, curve fitting technology has been widely used in image processing, engineering experiments and other fields. Among them, B-spline curve fitting is the most common method in curve fitting, the method of B-spline curve fitting has the advantages of locality, continuity but the fitting precision is relatively low. In practical application, B-spline curve fitting requires higher accuracy and speed. In order to increase the speed of planar B-spline curve fitting, Anderson acceleration is applied to the method of planar B-spline curve fitting. And then a quasi-Newton method based on Anderson acceleration is proposed. Firstly, an initial shape is set, and then the position parameters of the projection point of each data point are found according to the initial shape. Then, the corresponding position of control points is calculated by Anderson acceleration, and the above two steps are iterated until the result converges. The experimental results show that the proposed method in this paper outperforms other methods with respect to convergence speed and iteration time.

B-spline fitting; Anderson acceleration; quasi-Newton method; curve fitting; spline approach

TP 391

10.11996/JG.j.2095-302X.2020020246

A

2095-302X(2020)02-0246-08

2019-08-19；

2019-09-28

秦文杰(1994–)，男，江蘇常州人，碩士研究生。主要研究方向為計算機圖形學、最優(yōu)化算法。E-mail：darkqin@mail.ustc.edu.cn

張舉勇(1984–)，男，重慶人，副教授，博士。主要研究方向為計算機圖形學、計算機視覺、數(shù)字圖像處理、數(shù)值最優(yōu)化。E-mail：juyong@ustc.edu.cn