基于jerk能量最小化的空間五次G2插值曲線構造

陳 娟，何 歆，張起航，孫義環，陸利正

(浙江工商大學統計與數學學院，浙江 杭州 310018)

光順曲線具有簡單的曲率輪廓(例如螺線的曲率是單調的),它的構造方法是計算機輔助設計和相關領域的核心問題[1-2].在路徑設計和軌道規劃時,使用光順曲線以滿足幾何造型的需求就顯得尤為重要[3].通常來說,曲線光順由某些內在的能量函數來衡量;除了滿足一定端點條件外,曲線表示的其余自由度往往通過優化某些能量函數來確定[4].在曲線擬合時,常見做法也是采用能量作為目標函數的構成部分,以改善擬合曲線的光順程度.

光順曲線的構造方法大致上可分為兩大類.第一類方法使用螺線分段擬合自由曲線,因此具有分段單調的曲率分布.Walton等[5]借助歐拉螺線提出G1插值方法.但更多研究側重于使用多項式曲線構造G2插值的螺線[6-9],以適應計算機輔助設計(CAD)系統的要求.然而,這類方法須滿足一定端點條件才能保證得到的插值曲線是螺線,這意味著在某些端點條件下不存在G2插值的螺線,因此不適用于任意的端點數據,從而給外形設計造成不便.

然而,已有研究工作主要集中在平面曲線情形,很少專門針對空間曲線的系統研究.部分G1插值構造方法[10-11,13-14]可直接推廣用于空間曲線的G1插值問題.由于空間曲線必須考慮在插值點處的曲率和Frenet標架,已有平面曲線G2插值構造方法[6-9]還不能應用于空間曲線的G2插值問題.

本文提出基于jerk能量最小化的空間五次G2插值曲線構造方法.對于任意給定的G2數據,提出端點滿足特定曲率和Frenet標架且含4個自由參數的空間五次Bézier曲線模型;然后,借助能量優化確定這4個自由參數,以生成光順曲線.目前這方面研究工作仍不多,Lu[12]的工作由于采用最小化應變能優化曲線,導致目標曲線的曲率輪廓在許多情況下不夠理想,并且只適用于平面曲線的插值問題.本文主要的目的是在Lu[12]的工作基礎上做出改進,使其更好地滿足實際需要.

借鑒文獻[14-17]的做法,本文使用jerk能量近似表示曲線的曲率變化.首先,定義目標函數為jerk能量和曲線長度的加權組合,將曲線長度作為修正項以抑制曲線長度；然后,化簡后的目標函數是二元四次多項式,并且很容易求得它的梯度和Hessian矩陣；最后,使用投影牛頓法確定參數的最優取值.相比于Lu的工作[12],本文方法生成曲線的曲率變化更加簡單平緩,因此更符合光順要求.另外,本文方法也適用于空間曲線情形,滿足幾何造型的實際需求.

1 空間五次G2插值曲線模型

假設兩端點處的G2數據為{P0,T0,κ0N0;P1,T1,κ1N1}，其中κi,Ti和Ni分別表示曲線在端點Pi的曲率、單位切向量和單位法向量.那么曲線在端點處的Frenet標架分別表示為{T0,N0,B0}和{T1,N1,B1},Bi=Ti×Ni.

空間五次G2插值曲線的構造問題是指針對給定的G2數據,構造一條滿足插值條件的空間五次Bézier曲線

(1)

圖1 空間五次G2插值方法Fig.1 Spatial quintic G2 interpolating method

假設五次Bézier曲線b(t)的長度為L以及弧長參數表示為b(s),s=s(t)∈[0,L]，那么曲線b(t)的導數計算為：

b′(t)=s′(t)b′(s),

b″(t)=s″(t)b′(s)+(s′(t))2b″(s).

由于曲線在兩端點處須滿足給定的曲率和Frenet標架,再借助Frenet公式[19],所以：

其中自由參數αi=s′(i),βi=s″(i),i=0,1.再根據Bézier曲線的性質[1],控制頂點計算為：

(2)

(3)

(4)

(5)

只要確定αi和βi的參數值,就可代入計算得到五次Bézier曲線的控制頂點.對于任意給定的G2數據,該五次模型不僅能生成相應的插值曲線,而且為后續形狀優化提供4個可修改參數.同時,該五次模型只取決于相鄰兩個插值點的端點條件,適合自由曲線的分段構造.

2 jerk能量和曲線優化構造

對于曲線b(t),t∈[0,1]， Farin[13]提出曲率變化e是較好衡量光順的觀點,其定義為

(6)

雖然b(t)是五次多項式曲線,但κ(t)是非線性的,[κ′(t)]2是次數很高的有理分式,這不利于后續優化求解.尤其當碰到自由曲線由許多段構成時,計算量就非常大且很難做到實時求解.因此,使用jerk能量

(7)

近似表示曲率變化.當表示曲線的參數是弧長參數時,該曲線的jerk能量和曲率變化剛好相等；對于一般的多項式曲線而言,jerk能量是曲率變化一個較好的近似.

另一方面,外形設計時曲線長度也是需要考慮的重要因素.對于一系列插值點來說,構造的插值曲線應盡量接近插值點表示的曲線輪廓,同時避免過多振蕩.這要求兩個插值點之間的曲線形狀不能變化太大,因此在構造光順曲線時就要求曲線長度也保持適當小.基于上述考慮,構造目標函數

(8)

通常取λ為較小的正數,因此稱第2項為修正項,用于抑制目標曲線的長度.在目標函數中加入修正項有兩個好處:第一,可進一步增強目標函數的凸性,提高優化求解時的穩定性;第二,通過抑制‖b′(t)‖的整體變化使曲線參數更接近弧長參數,從而使jerk能量更接近曲率變化.

利用Bézier曲線的導數公式[1]

其中Δbi=bi+1-bi和Δ3bi=bi+3-3bi+2+3bi+1-bi是控制頂點的差分算子,以及Bernstein多項式的積分公式[20]

式(8)展開為關于控制頂點的二次函數:

(9)

其中

[gij]=

[hij]=

將控制頂點的表達式(2)～(5)代入目標函數(9)后,f本質上是α0,α1的四次函數和β0,β1的二次函數.顯然,曲線b(t)在兩端點處切向量的長度分別等于α0和α1.為使曲線盡量光順,α0和α1的取值不能太大也不能太小.因此,限定(α0,α1)的可行域為

D={(α0,α1):d0≤α0≤d1,d2≤α1≤d3},

其中di>0是給定的界.最后通過求解如下約束優化問題

(10)

確定4個未知數α0,α1,β0,β1， 從而最終得到五次Bézier曲線的控制頂點.

3 算法實現

目標函數(9)的梯度計算如下.首先,

根據式(2)～(5),bi的偏導數中只有6個不為零,羅列如下:

?b1/?α0=(1/5)T0,

?b2/?α0=(2/5)T0+(κ0α0/10)N0,

?b4/?α1=(-1/5)T1,

?b3/?α1=(-2/5)T1+(κ1α1/10)N1,

?b2/?β0=(1/20)T0,

?b3/?β1=(1/20)T1.

應用鏈式法則,得

令?f/?β0=?f/?β1=0， 得到關于β0,β1的線性方程組(其中cij是依賴于兩端點處G2數據的常數):

從該方程組可知,β0,β1是α0,α1的二次函數,即β0=β0(α0,α1),β1=β1(α0,α1).這意味著當f取到極值時,β0,β1總是二次依賴于α0,α1.

最后,將β0=β0(α0,α1)和β1=β1(α0,α1)代入目標函數(9)后,函數f是α0,α1的二元四次多項式.于是,問題(10)最終化簡為

(11)

使用投影牛頓法[21]對上述問題進行數值求解,通常取初值為(α0,α1)=(1,1).函數f關于α0,α1的梯度和Hessian矩陣根據鏈式法則重新計算;雖然具體公式表達式有些復雜,但借助數學軟件容易計算.

4 實 例

4.1 平面曲線

圖中圓點為控制頂點，折線表示控制多邊形，下同.圖2 本文方法(虛線)與Lu[12]方法(實線)的比較Fig.2 Comparison of the method of this article(dashed) and Lu’s[12](solid)

下面給出一些平面曲線例子來驗證本文方法的可行性,并與文獻[12]進行比較.在所有例子中,取權重λ=0.01和可行域D=[0.1,5]2.

考慮6個代表性的G2插值問題,不失一般性,令插值點為P0=(0,0),P1=(1,0)， 以及相應的單位切向量為T0=(cosφ0,-sinφ0),T1=(cosφ1,sinφ1).圖2所示為本文方法與Lu[12]的方法的對比結果:曲線圖列在左邊,而對應的曲率圖列在右邊.對比表明,本文方法得到曲線具有曲率變化更加平緩的特點,因此更符合光順要求.表1列出這兩種方法得到曲線的曲率變化和jerk能量.對比數據可以看出,本文方法得到曲線的曲率變化值更小.

表1 圖2中兩種方法得到曲線的曲率變化和jerk能量Tab.1 Curvature variation and jerk energy of thecurve obtained by the two methods in figure 2

歐拉螺線具有線性曲率，是一類廣泛使用的光順曲線,常用于路徑設計.考慮一段歐拉螺線(圖3(a)中實線):

首先將該螺線分成3段,然后分別用文獻[6,12]中的方法和本文方法得到五次G2插值曲線.從圖3可知,本文方法得到的曲線不僅更接近歐拉螺線,而且具有更完美的曲率分布.

圖3 歐拉螺線的五次G2分段插值Fig.3 Quintic G2 piecewise interpolation of Euler spiral

4.2 空間曲線

本文方法適用于空間曲線的分段插值構造.例如,對于Viviani曲線和Eightknot曲線[19],其定義分別為：

圖4 Viviani曲線的五次G2分段插值Fig.4 Quintic G2 piecewise interpolation of Viviani curve

圖5 Eightknot曲線的五次G2分段插值Fig.5 Quintic G2 piecewise interpolation of Eightknot curve

5 結 論

本文提出基于jerk能量最小化的空間五次G2插值曲線構造方法.衡量曲線光順的標準有幾個,而曲率變化是較好的標準.用jerk能量近似表示曲率變化以降低約束優化求解時的難度,因此可滿足大規模插值問題快速求解的需要.許多對比實例表明本文方法能夠生成更光順的曲線,也可用于解決一系列空間數據點的插值曲線構造問題.下一步研究重點是開展在路徑設計等方面的應用[3].

由于目標函數是jerk能量和曲線長度的加權組合,所以權重λ的取值對最終結果有一定影響.如果λ越小,那么生成曲線的jerk能量就越小.對于大多數例子建議取默認值λ=0.01，而少部分例子需要交互修改.

此外,本文關于空間曲線高質量插值構造方法可用于曲面的插值與擬合問題.未來將進一步開展曲面問題的研究.