曲率單調Bézier曲線G1插值算法及應用

何 川，趙 罡，王 偉，王愛增

(1.北京航空航天大學機械工程及自動化學院，北京 100191；2.北京航空航天大學虛擬現實技術與系統國家重點實驗室，北京 100191)

在道路設計、路徑規劃等問題中，需要避免車輛行駛過程中的突然轉向，因此設計師致力于構造曲率連續變化的緩和曲線。歐拉螺旋線(Euler spiral)具有曲率隨弧長線性變化的特點，被廣泛應用于路徑規劃[1-3]、曲線過渡[4-5]、形狀修復等問題中[6-8]，然而，歐拉螺旋線采用超越方程的形式定義，求解困難，因此出現了許多利用多項式曲線或有理多項式曲線來近似歐拉螺旋線的方法[9-11]。

除了歐拉螺旋線，多項式曲線也常常用來設計緩和曲線。文獻[12-14]采用了平面三次Bézier 螺旋線構造過渡曲線，并對文獻[12,14]內容進行了推廣，增加了Bézier 螺旋線的自由度[15-16]。

由于螺旋線的曲率隨弧長線性變化，這一約束使得螺旋線的構造十分困難，工程實際中如果能夠保證曲率單調遞增或遞減，則認為得到了符合光順性要求的美學曲線。2006 年，FARIN[17]利用Class A矩陣表示Bézier 曲線控制邊之間的幾何關系，提出了一種構造曲率單調變化(monotone curvature variation，MCV) Bézier 曲線的方法，不過該理論被CAO 和WANG[18]證明不夠完善，WANG 和ZHAO[19]給出了反例，并基于幾何方法證明了一類平面任意次Bézier 曲線曲率單調的一個充分必要準則[20]。針對歐拉螺旋線G1 插值時的效率與精度問題，本文提出了一種可滿足G1 插值條件的MCV Bézier 曲線構造方法，并與基于歐拉螺旋線的插值方法進行了對比，相關的數值實例證明了本文方法的有效性與實用性；且相比于文獻[21]采用曲率單調的B 樣條曲線實現G1 插值，本文構造方法簡單直觀，曲線彎曲度更好，算法的適用性與通用性更強。本文算法與圖例均在Windows 10 系統下由MATLAB 語言實現。

1 歐拉螺旋線與MCV Bézier 曲線

1.1 歐拉螺旋線

歐拉螺旋線，又稱為羊角螺旋線(cornu spiral)或回旋線(clothoid)，其曲率κ可表示為弧長L的線性函數，即

其中，0κ為L=0 時的初始曲率值；σ為歐拉螺旋線的曲率變化率。歐拉螺旋線的切線方向或旋轉角θ關于弧長的表達式為[6]

對于規范化的歐拉螺旋線(圖1)，其參數形式為

圖1 規范化的歐拉螺旋線，(x0,y0)=(0,0)，L=4 Fig.1 Normalized Euler spiral,(x0,y0)=(0,0),L=4

其中，C(L)與S(L)為菲涅爾積分(Fresnel integral)[22]，其表達式為

1.2 MCV Bézier 曲線

定理1.平面k次Bézier 曲線定義式為

且當比例因子s≥1 時，曲線P(t)具有單調遞減的曲率；當0＜s＜1 時，曲率單調遞增。圖2 顯示了s≥1 時的一條曲率單減的三次Bézier 曲線及其對應的曲率梳。

圖2 曲率單調遞減的三次Bézier 曲線及其曲率梳，s=1.5，θ＝，V0=(20,0) ((a)三次MCV Bézier 曲線；(b)曲線的曲率梳)Fig.2 Cubic MCV Bézier curve with decreasing curvature and corresponding curvature comb,s=1.5,θ＝,V0=(20,0) ((a) Cubic MCV Bézier curve;(b) Corresponding curvature comb)

2 MCV Bézier 曲線G1 插值

定理 2.給定 2 個插值端點PA(xA,yA)和PB(xB,yB)，及其相應的方向角θA和θB，存在符合邊界約束的k次MCV Bézier 曲線的充分條件為

具有滿足s…co sθ≥ 1(s≥ 1)或s≤cosθ(0＜s＜1)的實數解，其中，θ=(θB-θA)/(k-1)，方向角θA和θB以實際繞過的圈數計算，可大于2π。和s是未知變量。

證明.對于次數k一定時，θ=(θB-θA)/(k-1)即為構造MCV Bézier 曲線的旋轉角度，若方程組存在滿足s…cosθ≥ 1(s≥ 1)或s≤cosθ(0＜s＜1)的實數解，即滿足定理1 要求的MCV 條件，則由此構造的Bézier 曲線曲率必單調。證畢.

下面分析邊界約束條件對MCV Bézier 曲線解的存在性的影響。

由于Bézier 曲線的仿射不變性，故可將其中一個端點PA固定在原點，其方向角θA作為X軸正方向，并以此建立直角坐標系，則PA=(0,0)，θA=0，記2 端點連線與X軸正方向的夾角為θ? (圖3)，為討論方便，本文僅討論PB位于第一象限的情形，即xB＞0，yB＞0。

圖3 θB＜θ?，無解 Fig.3 θB＜θ?,no solution exists

2.1 次數k=2

圖4 θ B＞，無解 Fig.4 θ B＞,no solution exists

圖5 二次MCV Bézier 曲線控制邊Fig.5 The control edge of quadratic MCV Bézier curve

圖6 約束條件：PB=(20,20)，θA=0 Fig.6 Constraints:PB=(20,20),θA=0

圖7 二次MCV Bézier 曲線，s=2.732 051，||V0||=8.452 995，θ＝θB＝ ((a)二次MCV Bézier 曲線；(b)曲線曲率圖)Fig.7 Quadratic MCV Bézier curve,s=2.732 051,||V0||=8.452 995,θ＝θB＝ ((a) Quadratic MCV Bézier curve;(b) Corresponding curvature plot)

圖8 二次MCV Bézier 曲線，s=0.565 685，=50，θ＝θB＝ ((a)二次MCV Bézier 曲線；(b)曲線曲率圖)

Fig.8Quadratic MCV Bézier curve,s=0.565 685,=50,θ＝θB＝((a) Quadratic MCV Bézier curve; (b) Corresponding curvature plot)

2.2 次數k 可變時

命題2.對于任意次數k，若θB＜θ? (圖9)，則不存在滿足約束的MCV Bézier 曲線。

圖9 θB＜θ?，無MCV Bézier 曲線解 Fig.9 θB＜θ?,without solution of MCV Bézier curve

證明.θB＜θ?的情形如圖9 所示，由k次MCV Bézier 曲線的構造過程可知，控制邊Δbi(i=1,2,…,k-1)的方向角θi=(i-1)…θ是單調遞增的序列，且倒數第二個控制點bk-1位于虛線QPB上第一象限部分，其方向角θk-1=θB。

假設下標為i的控制點bi位于虛線PAPB下側(即三角形ΔPAPPB內)，下標為i+1 的控制點bi+1位于虛線PAPB上側，則控制邊Δbi=bi+1-bi將由下至上穿過PAPB，其方向角滿足θi＞θ?＞θB=θk-1，與已知控制邊的方向角的遞增序列矛盾，故θB＜θ?時不存在滿足約束的MCV Bézier 曲線證畢.

命題2 說明MCV Bézier 曲線不存在S 型曲線。

命題3.當θB＞θ?時，總能找到合適的次數k(k≥2)，使得k次MCV Bézier 曲線滿足約束。

證明.單調遞增的序列，當次數k為自由變量時，MCV Bézier 插值曲線的方程組為

推論.當θB＞θ?時，若k=m≥2 時存在MCV Bézier 曲線滿足約束，則當k≥m+1 時，也存在MCV Bézier 曲線滿足約束，即MCV Bézier 曲線存在多解，且解的極限為0θ→，1s→ 。

2.3 算法

基于G1 邊界約束的k次MCV Bézier 曲線插值算法的步驟如下：

3 數值實例

圖10 G1 插值MCV Bézier 曲線，k=6 ((a)本文算法求解的MCV Bézier 曲線;(b)曲線曲率圖) Fig.10 G1 interpolation MCV Bézier curve,k=6 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

4 與基于歐拉螺旋線的G1 插值算法比較

圖11 G1 插值MCV Bézier 曲線，k=18 ((a)本文算法 求解的MCV Bézier 曲線；(b)曲線曲率圖) Fig.11 G1 interpolation MCV Bézier curve,k=18 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

圖12 G1 插值MCV Bézier 曲線，k=26 ((a)本文算法求解的MCV Bézier 曲線；(b)曲線曲率圖) Fig.12 G1 interpolation MCV Bézier curve,k=26 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

歐拉螺旋線的曲率隨弧長線性變化，是光順 性非常好的一類曲線，然而其復雜的數學形式只能用數值方法進行近似求解，因此帶來了求解效率與計算精度等問題。而本文方法采用的Bézier曲線，屬于簡單的多項式函數，計算結果精確，效率高。與傳統基于歐拉螺旋線的G1 插值算法相比，優點如下：

(1) 構造簡單，計算高效；

(2) 計算精確，插值結果無誤差；

(3) 與現有的NURBS方法兼容。

圖13～16 展示了相同約束條件下采用歐拉螺旋線插值與本文MCV Bézier 曲線插值的結果。

給定約束PA=(0,0)，θA=0，PB=(50,50)，歐拉螺旋線與MCV Bézier 曲線求解結果分別如圖13和圖14 所示。歐拉螺旋線初始曲率0κ=0.007406，弧長L=85.919605，曲率變化率σ=0.000395。本文算法MCV Bézier 曲線插值結果如圖14 所示，k=5，s=0.797364，表1 為2 種方法的誤差與時間對比結果。

表1 誤差與時間對比1Table 1 First comparison of position error and time

圖13 歐拉螺旋線插值第一例((a)歐拉螺旋線；(b)歐拉螺旋線曲率圖(斜率放大20 倍))Fig.13 First interpolation example of Euler spiral ((a) Euler spiral;(b) Corresponding curvature plot (The slope is scaled by 20 times))

圖14 本文算法MCV Bézier 曲線插值第一例((a) MCV Bézier 曲線；(b)曲線曲率圖) Fig.14 First interpolation example of MCV Bézier curve based on our algorithm ((a) MCV Bézier curve; (b) Corresponding curvature plot)

給定約束PA=(0,0)，θA=0，PB=(40,30)，θB=π-arctan歐拉螺旋線與MCV Bézier 曲線求解結果分別如圖15 和圖16 所示。歐拉螺旋線初始曲率κ0=-0.010878，弧長L=67.290185，曲率變化率σ=0.001427。本文算法MCV Bézier 曲線插值結果如圖16 所示，k=4，s=0.5876780，對比結果見表2。

圖15 歐拉螺旋線插值第二例((a)歐拉螺旋線；(b)歐拉螺旋線曲率圖(斜率放大20 倍)) Fig.15 Second interpolation example of Euler spiral ((a) Euler spiral;(b) Corresponding curvature plot (The slope is scaled by 20 times))

圖16 本文算法MCV Bézier 曲線插值第二例 ((a) MCV Bézier 曲線；(b)曲線曲率圖) Fig.16 Second interpolation example of MCV Bézier curve based on our algorithm ((a) MCV Bézier curve; (b) Corresponding curvature plot)

表2 誤差與時間對比2Table 2 Second comparison of position error and time

5 分析與討論

5.1 算法推廣

本文給出了基于第一象限2 端點的G1 約束MCV Bézier 樣條曲線插值算法，該算法能夠推廣至更加復雜的端點約束情形，如非第一象限的端點插值及多數據點插值。

針對尾點位于非第一象限的情況，可以通過補充中間插值點的方式進行分段處理，使得每一段的尾點位于局部坐標系的第一象限，將原問題轉化為多個第一象限的G1 插值問題。

對于給定一個平面數據點序列的插值問題，同樣需要采取分段處理的方式，除首尾2 端點外，中間數據點沒有切線方向的限制。由于減少了每一段插值曲線尾點的切向約束(最后一段除外)，因此相鄰2 點的插值自由度增加，曲線設計更加靈活。

分段處理引入的一個問題是，分段點處曲率可能會發生突變，且多段拼接而成的整條曲線，只能保證分段曲線的曲率單調性。要解決以上問題，需要引入G2 插值，這是將來的研究方向。

給定平面數據點序列PA=P0=(0,0)，P1=(14,3)，P2=(18,16)，P3=(12,22)，PB=P4=(-2,26) (圖17，粉色數據點)，以及首尾切線方向θA=0，θB=π。一種分段插值結果如圖17 藍色曲線所示，該曲線是由4段二次MCV Bézier 曲線組合而成，每個中間數據點即為分段連接點。組合曲線的曲率圖分段單調(圖18)，在中間3 個數據點處曲率發生突變。

圖17 給定平面數據點列的MCV Bézier 曲線插值 Fig.17 MCV Bézier curve interpolation of given planar data point

圖18 組合曲線的曲率圖Fig.18 Curvature plot of the whole composite curve

5.2 與基于平面三次B 樣條曲線G1 插值算法比較

文獻[21]提出了一種曲率單調的平面三次B 樣條曲線構造算法，且將其應用于形狀修復，屬于一類G1 插值問題。該算法是通過前2 條控制邊矢量的夾角及模長，確定后一條控制邊矢量的范圍以保證曲率的單調性。但未具體給出如何選取控制邊矢量的方向及大小，算法隨機性較強；且相鄰邊矢量的夾角逐漸減小，隨著控制點的增多，曲線趨于直線，曲線的角度變化范圍有限。

對于5.1 節所示的復雜平面點列，文獻[21]算法受限于首尾切線方向的角度變化，無法生成通過各點的插值曲線，相關理論有待進一步研究。而本文方法基于變換矩陣，構造更加快速、直觀，曲線整體的彎曲度更大，能夠處理的切向約束范圍更廣，因此算法的適用性與通用性更強。

6 結束語

本文提出了一種基于G1 約束的MCV Bézier 樣條曲線插值算法。首先闡述了給定邊界約束下MCV Bézier 曲線解的存在性條件并進行證明；接著，理論與實踐表明，該算法可以靈活構造滿足G1 插值條件的曲率單調Bézier 曲線；對比基于歐拉螺旋線的插值算法，本文算法更加精確且高效；作為處理第一象限2 端點G1 約束問題的基礎，該算法可以推廣至更加復雜的約束問題求解，相比于曲率單調的平面三次B 樣條曲線，MCV Bézier 曲線的構造更加快速、直觀，算法的適用性與通用性更強。

未來工作是尋求當PB位置與方向角θB任意取值時，構造滿足約束的MCV Bézier 樣條曲線的有效解法，以解決復雜情形下的路徑規劃等工程問題；分段處理將引入曲率的突變，G2 插值問題相比G1 更加重要且復雜，這是未來的研究重點。