局部插值的H－Bézier曲線分段光順擬合

張龍飛， 檀結(jié)慶

（合肥工業(yè)大學 數(shù)學學院，安徽 合肥 230009）

局部插值的H－Bézier曲線分段光順擬合

張龍飛， 檀結(jié)慶

（合肥工業(yè)大學 數(shù)學學院，安徽 合肥 230009）

文章提出了一種基于局部插值的光順擬合方法，通過對給定的離散數(shù)據(jù)點進行逐段擬合，可以完整地表現(xiàn)出已知數(shù)據(jù)點的分布情形；還給出一種改進的光順準則，針對不同情況確定光順約束方程之后，通過調(diào)整控制頂點及可能存在的曲線控制參數(shù)，使曲線的光順能達到極小，從而得到較為理想的光順擬合曲線。數(shù)值實例說明了該文算法的有效性。

H－Bézier曲線；分段擬合；光順；能量法；曲率

在計算機輔助幾何設計中，常會遇到給定數(shù)據(jù)點列，要求構(gòu)造一條（張）部分通過此點列的曲線（面），使得曲線（面）能夠較好地反應數(shù)據(jù)點列所表現(xiàn)的形狀，有時候還要求構(gòu)造光順的曲線曲面來擬合給定點列，船舶業(yè)及航空工業(yè)中就經(jīng)常使用光順的曲線曲面來擬合給出的數(shù)據(jù)點，這涉及光順和擬合的問題。工業(yè)設計中，判別曲線光順與否的標準有很多種［1］，較為常見的光順方法是能量法。三次樣條曲線的交互式光順技術由文獻［2－3］提出，并在此后得到了快速的發(fā)展，如節(jié)點刪除插入法［2－4］、Kjellander算法及其改進算法［3－5］等。

近年光順擬合的研究較多，如文獻［6－8］。光順準則一般采用改進后的光順方法，比如既考慮到應變能又考慮到應力能等其他因素的光順準則［6］。擬合時一般是在最小二乘法的基礎上使用一條曲線去擬合，但僅用一條曲線擬合給定的數(shù)據(jù)點所得到的曲線并不能夠很好地反映出數(shù)據(jù)點的整體情況，分段擬合可以更為準確地反映出數(shù)據(jù)點整體分布狀況，局部插值可以在需要的地方增進對數(shù)據(jù)點的逼近程度，將分段擬合與光順處理相結(jié)合在工程設計中有較大的用處。綜上所述，分段光順擬合有其研究的必要性。

本文提出了一種逐段光順擬合的算法，給出了一種改進的光順模型，其中考慮了光順前后的控制頂點的偏差以及光順后曲線與所給定的數(shù)據(jù)點的偏差，解決了光順后曲線與原擬合曲線間偏離可能較大的問題，并在光順擬合的同時實現(xiàn)段與段之間的光滑拼接。大家所熟知的H－Bézier曲線（面）不但保留了Bézier曲線（面）的優(yōu)點，還可精確表示懸鏈線（面）、螺旋線（面）等工程中常見的曲線，雖然目前對于H－Bézier曲線的基本理論、分割拼接及形狀分析都有學者做了研究［9－11］，但關于 H－Bézier曲線（面）光順的研究成果還比較少。鑒于低次曲線在工程實踐中應用范圍較廣，本文使用三次H－Bézier曲線為例來驗證算法的有效性。

1 三次H－Bézier曲線的定義

定義1 設hi∈R2（i＝0，1，2，3）為控制頂點，α是大于零的任意實數(shù)，則曲線

2 改進的光順準則

判別一條曲線光順與否的標準有很多種［1］，較為常見的方法是能量法，曲線光順中的能量法一般是相對于曲線的應變能而言的。所謂應變能是指以應力和應變的形式貯存在物體中的勢能，又稱變形能，在用能量法光順曲線p（t）時，人們通常建立的函數(shù)是對曲線按弧長取關于曲線曲率的積分，即應變能：

其中，k為曲率；s為弧長。

由于（3）式計算起來較為復雜，通常以其近似形式來逼近，即

能量法與曲線的絕對曲率相關，能量趨小會使得曲線趨于直線化，但曲線整體的曲率變化卻未必均勻。而剪力躍度的平方和與曲線曲率的變化相關，當剪力躍度的平方和較小時，曲線的曲率變化較為均勻。當知道n＋1個點處的相對曲率以及它們之間的弦長時，可以得到剪力躍度的平方和［12－14］。令

其中，ki為p（ti）點的曲率值；｜di｜為剪力躍度；li為p（ti－1）到p（ti）的弦長。

將剪力躍度的平方和整合進能量法，可以得到新的能量函數(shù)：

其中，a＋b＝1，a＞0，b＞0。在使用以上光順準則對曲線光順之后，原曲線與光順后的曲線之間的偏差可能比較大。

本文針對以上問題將光順前后的控制頂點的偏差以及光順后曲線與所給定的數(shù)據(jù)點｛Pi｝ni＝0的偏差考慮進來。因此，在給定偏差最大值ε1（控制點偏差）及ε2（數(shù)據(jù)點偏差）時，就應在

的條件下，使得能量函數(shù)F極小。其中，m為擬合出來的曲線次數(shù)；qj為原曲線的控制頂點；q1j為光順后的控制點。由此得到改進后的光順模型：

其中a，b，c∈［0，1］，0≤a＋b＋c≤1。對于給定m次曲線c，記其控制頂點為q0，…，qm。

在實驗過程中注意到這樣一個問題，給定的采樣數(shù)據(jù)點中，有一些數(shù)據(jù)點與可以精確逼近大部分數(shù)據(jù)點的最終擬合曲線之間的距離較大，若光順時將所有數(shù)據(jù)點的作用看作同等，則光順后得到的曲線與原擬合曲線間偏差會比較大，減小這類數(shù)據(jù)點的作用，可以得到與擬合曲線偏差不大且光順效果較好的光順曲線。

因此，在確定光順函數(shù)時，應該考慮這一因素，使得擬合后距離擬合曲線較近的數(shù)據(jù)點起的作用大一些，距離擬合曲線較遠的數(shù)據(jù)點所起的作用小一些。

實際操作時，可以給每一個數(shù)據(jù)點加上權(quán)系數(shù)｛wi｝ni＝0。首先計算每一個采樣數(shù)據(jù)點｛Pi｝ni＝0與最終擬合曲線的距離，記為｛edi｝ni＝0。對采樣數(shù)據(jù)點按距離擬合曲線的遠近排序，距離曲線最近的為首數(shù)據(jù)點，距離曲線最遠的為末數(shù)據(jù)點，結(jié)果記為｛P′i｝ni＝0，同時將數(shù)據(jù)點｛Pi｝ni＝0與擬合曲線的距離按從大到小的順序排序，所得結(jié)果記為｛ed′i｝ni＝0。 由 此 得 到 排 序 后 的 每 一 個 數(shù) 據(jù) 點｛P′i｝ni＝0的權(quán)系數(shù)｛w′i｝ni＝0，即

由此得到數(shù)據(jù)點Pi的權(quán)系數(shù)wi。將這一權(quán)函數(shù)加入到光順函數(shù)中有：

令F→min，即可得到光順后曲線c′的控制頂點q10，…，q1m，即

解方程組（9）可求得曲線光順后的新控制頂點，從而得到光順后的曲線。當需要光順后的曲線仍通過首末端點時，有q0＝q10，qm＝q1m，（9）式轉(zhuǎn)化為：

使用方程組（10）代換方程組（9）即可。

當需要光順拼接時，則在（7）式的基礎上加入拼接條件，以G1拼接為例，記前一段m次曲線為c1，其控制頂點為q1，0，…，q1，m，后一段n次曲線為c2，其控制頂點為q2，0，…，q2，n，記其在光順后的控制頂點為q12，0，…，q12，n，為使c2在光順后能與c1達到G1拼接，則應有：

將（11）式作為（9）式的約束條件，得到帶有約束條件的光順模型：

解這一含有約束條件的優(yōu)化問題，可求得滿足條件的解。

3 曲線擬合

本文在做曲線擬合的時候，采用最小二乘法來擬合給定的數(shù)據(jù)點列。設給定數(shù)據(jù)點列為擬合出來的曲線為p（t），其控制頂點為q0，…，qm。

局部插值的最小二乘擬合方法如下：

（1）對給定的數(shù)據(jù)點做累加弦長參數(shù)化：

其中，lj＝｜Pj－Pj－1｜是Pj到Pj－1的弦長。

（2）設若希望擬合出來的曲線通過n＋1個數(shù)據(jù)點中的N（N＜m＋1）個，記這N個點參數(shù)化后的點為tN，0，…，tN，N－1（將這 N 個點按順序表示，僅為方便記）。將tN，0，…，tN，N－1代入p（t），則有p（tN，i）＝0，i＝0，…，N－1。以上N 個方程，共包含m＋1個未知量，將q0，…，qN－1用余下的m＋1－N個自由未知量表示出來。再求p（t）對這m＋1－N個自由未知量的偏導并使偏導等于零，從而得如下方程組：

（3）解（13）式可得q0，…，qm，求出m 次擬合曲線。如果希望曲線通過首末數(shù)據(jù)點則有q0＝P0，…，qm＝Pn，代入（13）式可簡化計算量。

4 分段光順擬合算法

本文提出的局部插值分段光順擬合算法，可以逐段擬合光順，通過上一段的光順對下一段的光順做相應調(diào)整，使得前后段擬合曲線在連接點處能夠達到G1連續(xù)。設拼接點的前一段曲線控制頂點為q0，q1，q2，…，qn，后一段曲線的控制頂點為v0，v1，v2，…，vn，qn＝v0為拼接點。若要連接處達到G1連續(xù)，則需2段曲線的控制頂點滿足如下關系：

從而保持它們切向連續(xù)，可以允許它們模長不相等［14－15］。

如果需要曲線在拼接點處達到G2連續(xù)，就需要計算其二階導數(shù)，在拼接點處的二階導數(shù)滿足倍數(shù)關系。

拼接點處的前后2段曲線sfront（u1）和sback（u2）的一階導數(shù)和二階導數(shù)滿足如下關系：將（15）式展開化簡并結(jié)合其在拼接點處二階導數(shù)滿足倍數(shù)關系可以得到：

也就是說qn－2、qn－1、qn（＝v0）、v1、v25點共面。

分段光順擬合，首先要考慮的是對給定的點列進行分組，以分段處理。設給定數(shù)據(jù)點數(shù)為n，需擬合為m次曲線（m＜n），則每組的數(shù)據(jù)點數(shù)應不少于m＋2個，不然將得到一條完全插值于給定點的曲線。首先從數(shù)據(jù)點列中取出m＋2個點作為第1組，并將第1組的最后一個數(shù)據(jù)點作為下一段的起始點，記第i個數(shù)據(jù)點分組內(nèi)的第j個點為pi，j，第i組內(nèi)的數(shù)據(jù)點數(shù)為si，則j＝0，1，…，si。計算分組連接點與前后數(shù)據(jù)點所構(gòu)成向量之間的夾角及連接點到其前后數(shù)據(jù)點所構(gòu)成的直線段的距離，保證在連接點處不會出現(xiàn)曲線振動較大的情形，如圖1所示。給定夾角與距離的范圍（本文給定為夾角0°～90°，距離選定為0到前后數(shù)據(jù)點所構(gòu)成直線段長度的1／3），夾角與距離均在給定范圍則以此點為連接點，否則判斷下一點是否符合連接點標準，如圖2所示。若剩下的點數(shù)小于m＋2則并入上一組，若剩下的點數(shù)大于等于m＋2且小于2m＋4，則全部作為第2組。依次類推，對數(shù)據(jù)點列進行分組。

圖1 連接點處曲線出現(xiàn)較大振動的情形

圖2 連接點與前后數(shù)據(jù)點的向量構(gòu)成

其中，圖2中的θ為向量間的夾角；h為連接點到前后數(shù)據(jù)點所構(gòu)成直線段的距離。

分段光順擬合后，拼接點處可達到G1的具體算法如下：

（1）將數(shù)據(jù)點列按本節(jié)所述方法分組，記所分組列為ri（i＝1，…，k），k為總組數(shù)。

（2）利用介紹的局部插值擬合方法擬合r1內(nèi)的數(shù)據(jù)點列，使擬合曲線通過r1的首末數(shù)據(jù)點及其間所選定的若干（少于m－1個）數(shù)據(jù)點，得m次曲線c1。使用改進的光順模型（8）式處理曲線c1，得光順后的曲線c′1。

（3）擬合ri（i＝2，…，k）內(nèi)的數(shù)據(jù)點列，使得擬合的曲線通過ri的首末數(shù)據(jù)點及其間所選定的若干（少于m－1個）數(shù)據(jù)點，得到m次曲線ci。運用改進的光順模型（12）式，對曲線ci做光順處理，得到與前一段已光順曲線光滑拼接的光順曲線c′i。重復步驟（3），直到i＝k。

按照上述步驟，邊擬合邊光順，逐段進行，可以得到一條整體較為光滑且與所給數(shù)據(jù)點整體走勢較為相符的曲線。

5 應用實例

以三次H－Bézier為例，對給定數(shù)據(jù)點分2段

表1 數(shù)據(jù)點Pi（xi，yi）

圖3 文獻［6］方法光順后的曲線

圖4 文獻［6］方法所得曲線曲率

曲線及曲率如圖3～圖9所示。處理，簡單實現(xiàn)本文所提出的分段光順擬合方法。光順方程系數(shù)取a＝2×10－3，b＝3.9×10－1，c＝2×10－3，光順模型（12）式中的λ取1。給定一組數(shù)據(jù)點Pi（xi，yi），i＝0，1，…，9，見表1所列。

圖5 未做處理前的擬合曲線拼接

圖6 光順處理前曲線曲率

從圖3可以看出，由于文獻［6］方法僅與曲線應變能、應力能等有關，使得能量極小時所得曲線為直線。

盡管這時的曲線曲率極小（如圖4所示），但所得結(jié)果已不再是一條包含形狀變化的曲線，失去了應用價值，偏離了光順的初衷。在文獻［6］方法的基礎上加入控制頂點變化范圍的約束條件可以避免出現(xiàn)曲線直化的問題。在此，使用這種修正的方法與本文方法做對比。

圖7 本文所得曲線與文獻［6］修正方法所得曲線

圖8 本文方法所得曲線曲率

圖9 文獻［6］修正方法所得曲線曲率

由圖5～圖9可知，按照本文算法處理后的曲線整體走勢與擬合拼接曲線基本相同，但處理后的曲線曲率比未處理曲線的曲率變化更小，曲率分布區(qū)間相對光順處理前較小，達到了光順的目的。

本文方法相對于修正的文獻［6］方法更能準確反映給定數(shù)據(jù)點的走勢，處理之后得到的曲線與給定數(shù)據(jù)點的偏差遠小于文獻［6］方法所得的曲線與給定數(shù)據(jù)點的偏差，見表2所列，且能插值于指定點，更能適應工業(yè)上對復雜器件外形設計時的需求。由于文獻［6］所采用的光順準則是能量法準則，光順能愈小，曲率愈小，曲率變化區(qū)間變小，但由此帶來曲線變化愈是趨于直線化的結(jié)果，失去了設計曲線應有的多樣性。修正后的文獻［6］方法雖然可以避免曲線直化的問題，但曲線仍不能夠更準確地逼近數(shù)據(jù)點。

表2 光順處理后曲線與給定數(shù)據(jù)點的偏差

注意到本文方法光順后的曲線曲率圖中有間斷點，這是經(jīng)過本文算法處理后，在拼接點處達到了G1連續(xù)，雖然這已經(jīng)可以達到工業(yè)設計的要求，滿足了人們對設計光滑曲線的需求，但并不能保證在拼接點處達到G2連續(xù)，因此會有中間斷點的出現(xiàn)，這也是下一步要做的工作。

6 結(jié)束語

本文給出了分段光順擬合給定數(shù)據(jù)點的算法，并給出了應用實例。通過實例可以看出，本文算法取得了一定效果，提供了一種新的有效的工程設計手段。本文算法在達到光順效果的情況下盡可能地減小了計算量，且在曲線拼接處達到G1連續(xù)，可以滿足工程應用的需要。對于要求在拼接處達到G2連續(xù)的需求，是下一步的研究方向。

［1］Renz W.Interactive smoothing of digitized point data［J］.Computer Aided Geometric Design，1982，14 （5）：267－269.

［2］Farin G，Rein G，Spapidis N，et al.Fairing cubic B－spline curves［J］.Computer Aided Geometric Design，1987，4（1／2）：91－103.

［3］Kjellander J A P.Smoothing of cubic parametric spline［J］.Computer－Aided Design，1983，15（3）：175－179.

［4］Spapidis N，F(xiàn)arin G.Automatic fairing algorithm for B－spline curves［J］.Computer－Aided Design，1990，22（2）：121－129.

［5］鄭興國，檀結(jié)慶，三次Bézier曲線的光順［J］.合肥工業(yè)大學學報：自然科學版，2005，28（1）：109－112.

［6］Yang Yadi，Qin Xinqiang，Hu Gang，et al.Fairing and approximation algorithm of C－Bézier curves［J］.Journal of Computer Applications，2008，28（12）：3132－3134.

［7］沈海鷗，陳淑珍，孫曉安.曲線的二次有理Bézier曲線擬合［J］.武漢大學學報：自然科學版，1997，43（1）：119－123.

［8］Fl¨ory S，Hofer M.Constrained curve fitting on manifolds［J］.Computer Aided Geometric Design，2006，23（1）：17－44.

［9］王 媛.H－Bézier曲線的理論及其應用研究［D］.西安：西北大學，2006.

［10］吳榮軍.平面三次 H－Bézier曲線的形狀分析［J］.應用數(shù)學學報，2007，30（5）：813－821.

［11］檀結(jié)慶，王 燕，李志明.三次 H－Bézier曲線的分割、拼接及其應用［J］.計算機輔助設計與圖形學學報，2009，21（5）：584－588.

［12］Zhang Caiming，Zhang Pifu，Cheng Fuhua.Fairing spline curves and surfaces by minimizing energy［J］.Computer－Aided Design，2001，33（13）：913－923.

［13］蘇步青，劉鼎元.計算幾何［M］.上海：上?？茖W出版社，1981：217－242.

［14］朱心雄.自由曲線曲面造型［M］.北京：科學出版社，2000：206－365.

［15］Farin G，Sapidis N.Curvature and the fairness of curves and surfaces［J］.Computer Graphics and Applications，1989，9（2）：52－57.

Piecewise fairing and fitting of H－Bézier curves by partial interpolations

ZHANG Long－fei， TAN Jie－qing

（School of Mathematics，Hefei University of Technology，Hefei 230009，China）

This paper presents a new method of fairing and fitting based on partial interpolations，which can completely describe the distribution of the given data by piecewise fitting.An improved fairing method is also put forward.Once the fairing constraint equations are determined case by case，the fairing energy can reach the minimum level through adjusting the control points and choosing the suitable control parameters to get an ideal approximation of the given points.Numerical examples are given to show the effectiveness of the presented methods.

H－Bézier curve；piecewise fitting；fairing；energy method；curvature

TP391

A

1003－5060（2011）10－1570－06

10.3969／j.issn.1003－5060.2011.10.029

2011－04－22；

2011－06－20

國家自然科學基金資助項目（61070227；60773043）；高等學校博士學科點專項科研基金資助項目（20070359014）

張龍飛（1986－），男，安徽蕭縣人，合肥工業(yè)大學碩士生；

檀結(jié)慶（1962－），男，安徽望江人，博士，合肥工業(yè)大學教授，博士生導師.

（責任編輯 朱華新）