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

秦文杰，張舉勇

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

秦文杰，張舉勇

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

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

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

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

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

1 相關工作

1.1 經典的B樣條曲線擬合算法

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

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

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

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

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

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

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

1.2 局部-全局求解器

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

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

1.3 安德森加速

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

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

2 算法

2.1 迭代求解的一般方法

已經知道B樣條曲線擬合問題可以化為以下最優化問題

一般分為2個步驟：

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

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

2.2 安德森加速

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

對于此問題而言

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

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

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

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

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

2.3 穩定性保證

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

2.4 m的選取

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

圖1 m不同取值對比圖

2.5 算法流程

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

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

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

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

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

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

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

3 實驗結果及分析

3.1 實驗的相關定義與參數

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

該函數是在穩定性保證過程中需要比較的目標函數。

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

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

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

3.2 實驗對比

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

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

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

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

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

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

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

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

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

3.3 單步迭代時間穩定性

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

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

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

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

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

3.4 初值穩定性

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

圖7 3種不同初始曲線

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

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

3.5 算法分析

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

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

等價于

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

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

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

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

4 結束語

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

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

[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–)，男，江蘇常州人，碩士研究生。主要研究方向為計算機圖形學、最優化算法。E-mail：darkqin@mail.ustc.edu.cn

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