張龍飛， 檀結(jié)慶

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

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

張龍飛， 檀結(jié)慶

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

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

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

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

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

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

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

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

2 改進(jìn)的光順準(zhǔn)則

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

其中，k為曲率；s為弧長(zhǎng)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 曲線擬合

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

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

（1）對(duì)給定的數(shù)據(jù)點(diǎn)做累加弦長(zhǎng)參數(shù)化：

其中，lj＝｜Pj－Pj－1｜是Pj到Pj－1的弦長(zhǎng)。

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

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

4 分段光順擬合算法

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

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

如果需要曲線在拼接點(diǎn)處達(dá)到G2連續(xù)，就需要計(jì)算其二階導(dǎo)數(shù)，在拼接點(diǎn)處的二階導(dǎo)數(shù)滿足倍數(shù)關(guān)系。

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

也就是說qn－2、qn－1、qn（＝v0）、v1、v25點(diǎn)共面。

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

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

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

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

分段光順擬合后，拼接點(diǎn)處可達(dá)到G1的具體算法如下：

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

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

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

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

5 應(yīng)用實(shí)例

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

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

圖3 文獻(xiàn)［6］方法光順后的曲線

圖4 文獻(xiàn)［6］方法所得曲線曲率

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

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

圖6 光順處理前曲線曲率

從圖3可以看出，由于文獻(xiàn)［6］方法僅與曲線應(yīng)變能、應(yīng)力能等有關(guān)，使得能量極小時(shí)所得曲線為直線。

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

圖7 本文所得曲線與文獻(xiàn)［6］修正方法所得曲線

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

圖9 文獻(xiàn)［6］修正方法所得曲線曲率

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

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

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

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

6 結(jié)束語(yǔ)

本文給出了分段光順擬合給定數(shù)據(jù)點(diǎn)的算法，并給出了應(yīng)用實(shí)例。通過實(shí)例可以看出，本文算法取得了一定效果，提供了一種新的有效的工程設(shè)計(jì)手段。本文算法在達(dá)到光順效果的情況下盡可能地減小了計(jì)算量，且在曲線拼接處達(dá)到G1連續(xù)，可以滿足工程應(yīng)用的需要。對(duì)于要求在拼接處達(dá)到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］鄭興國(guó)，檀結(jié)慶，三次Bézier曲線的光順［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，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］.武漢大學(xué)學(xué)報(bào)：自然科學(xué)版，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曲線的理論及其應(yīng)用研究［D］.西安：西北大學(xué)，2006.

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

［11］檀結(jié)慶，王 燕，李志明.三次 H－Bézier曲線的分割、拼接及其應(yīng)用［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，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］蘇步青，劉鼎元.計(jì)算幾何［M］.上海：上?？茖W(xué)出版社，1981：217－242.

［14］朱心雄.自由曲線曲面造型［M］.北京：科學(xué)出版社，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

國(guó)家自然科學(xué)基金資助項(xiàng)目（61070227；60773043）；高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目（20070359014）

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

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

（責(zé)任編輯 朱華新）