優化端點條件的平面二次均勻B 樣條插值曲線

李軍成，劉成志，趙文才

（湖南人文科技學院數學與金融學院，湖南婁底 417000）

B 樣條是目前大多數CAD 系統的基礎模塊，利用B 樣條構造插值曲線是工程中經常遇到的問題。反求法是最為常用的一種利用B 樣條構造插值曲線的方法，其基本思路是利用B 樣條曲線滿足的插值條件構造方程，通過求解方程反求B 樣條曲線的控制頂點。為使方程存在唯一解，在利用反求法構造B 樣條插值曲線時，往往需要補充端點條件。選取的端點條件不同，獲得的插值曲線也不同。在實際應用中，端點條件的選取往往較困難。

在CAD 及其相關研究領域，構造光順的曲線是一項重要的研究課題［1-3］。雖然目前尚無法定量描述曲線的光順性，需通過能量極小構造平面光順曲線［4］。例如，利用能量極小構造B 樣條［5-6］、插值曲線 曲 面［7-8］、Hermite 插值［4，9-13］、Bézier 曲 線［14-16］、曲面網格［17-18］等。既然在反求法中端點條件為構造B樣條插值曲線提供了自由度，那么很自然地想到可對端點條件進行優化，使得構造的B 樣條插值曲線滿足光順性要求。

為此，通過極小化曲線的內能，獲得具有較好效果的C1連續二次均勻B 樣條插值曲線，且所構造的插值曲線是保形的。

1 預備知識

1.1 二次均勻B 樣條曲線

給定平面上一列控制頂點pk(k=0，1，…，n)，二次均勻B 樣條曲線可表示為［19］

其中，i=0，1，…，n?2；0 ≤t≤1；Bj(t)為基函數，滿足

二次均勻B 樣條曲線段在端點處滿足

1.2 曲線的內能

拉伸能、應變能與曲率變化能是平面曲線內能的3 種常見形式，分別對應曲線的弧長、曲率和曲率變化率［6，14］。

給定平面參數曲線r(t)，其拉伸能、應變能與曲率變化能可分別定義為［14］

其中，κ(t)為曲線r(t)的曲率。

為簡化計算，在實際應用中，常將式（3）近似地表示為［14］

2 問題的提出

所要討論的問題是：給定平面上n個數據點di(i=1，2，…，n)，如何構造一整條插值于給定數據點且滿足C1連續的二次均勻B 樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}。

由式（1）知，要構造插值曲線r(t)，只需確定各段曲線ri(t) (i=0，1，…，n?2)的控制頂點pk(k=0，1，…，n)。由式（2）知，要使曲線r(t)插值于數據點di(i=1，2，…，n)，則必有

由式（5）可得

式（6）表明，只要確定了控制頂點p0，便可依次確定其他控制頂點pi(i=1，2，…，n)，或只需確定控制頂點pn，也可依次確定其他控制頂點pn?i(i=1，2，…，n)，這種確定插值曲線控制頂點的方法也稱為反求法，p0或pn稱為端點條件。由于二次均勻B樣條曲線具有對稱性［19］，故將p0或pn作為端點條件確定的二次均勻B 樣條插值曲線是完全相同的。下面僅討論將p0作為端點條件的情形。

事實上，當控制頂點p0給定時，控制頂點pi(i=1，2，…，n)可由以下定理確定。

定理1給定控制頂點p0，控制頂點pi(i=1，2，…，n)可表示為

證明用數學歸納法證明。

當i=1 時，由式（6），可得式（7）。 假設當i=m時，有

則當i=m+1 時，由式（6）知

即式（7）成立。

證畢。

由定理1 可知，插值于數據點di(i=1，2，…，n)的C1連續二次均勻B 樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}僅由控制頂點p0（即端點條件）決定。理論上，可通過選取不同的控制頂點p0使插值曲線滿足不同需求。當然，也可優化選取控制頂點p0，使構造的C1連續二次均勻B 樣條插值曲線滿足特定要求。

下面討論如何利用內能極小優化選取控制頂點p0，構造光順的C1連續二次均勻B 樣條插值曲線。

3 利用內能極小優化端點條件

由于對二次均勻B 樣條曲線求三階導數后歸零，故只討論如何通過拉伸能極小與應變能極小優化選取控制頂點p0，構造光順的C1連續二次均勻B樣條插值曲線。

3.1 利用拉伸能極小優化端點條件

由式（6），式（1）可改寫為

由式（8）可得

由式（4）與式（9），插值曲線r(t)={ri(t)，i=0，1，…，n?2}的拉伸能可表示為

由式（7）知，在式（10）中僅p0為自由變量。因此，為確定最優的p0，插值曲線須具有極小拉伸能，得到無約束優化問題

為求解問題（11），首先給出以下引理。

引 理 1設 向 量a=(ax，ay)，記則有

證明令向量b=(bx，by)，由于a?b=axbx+ayby，故有

從而有，

于是，由式（7）和式（12），可得

證畢。

定理2要使插值于數據平面上的數據點di(i=1，2，…，n)的C1連續二次均勻B 樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}具有極小拉伸能，控制頂點p0應為

證明當n=2 時，由引理1 和式（10），可得

當n＞2 時，由引理1 和式（10），可得

證畢。

3.2 利用應變能極小優化端點條件

式（8）經計算可得

由式（4）和式（16），插值曲線r(t)={ri(t)，i=0，1，…，n?2}的應變能可表示為

由式（7）可知，在式（17）中僅控制頂點p0為自由變量。為確定最優的p0使得插值曲線具有極小應變能，可得無約束優化問題

定理 3給定平面上的數據點di(i=1，2，…，n)，具有極小拉伸能的C1連續二次均勻B樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}也具有極小應變能。

證明當n=2 時，由引理1 和式（17），可得

當n＞2 時，由引理1 和式（17），可得

由定理2 可知，當插值曲線r(t)具有極小拉伸能時，控制頂點p0應滿足式（13）。將式（13）代入式（19）和式（20），可得，故式（13）也為式（18）的解。因此，具有極小拉伸能的C1連續二次均勻B 樣條插值曲線也具有極小應變能。

證畢。

3.3 保形插值

在數據插值中，保形插值一直以來深受重視［20-22］。若給定平面上的數據點di(i=1，2，…，n)，設s(t)：=si(t)=(1?t)di+tdi+1，0 ≤t≤1，顯然線性插值s(t)是一種最簡單的保形插值［23］。要使插值于數據點di(i=1，2，…，n)的C1連續二次均勻B樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}具有良好的保形性，須對控制頂點p0進行優化，使得r(t)盡可能地接近s(t)，于是，得到無約束優化問題

定理4給定平面上的數據點di(i=1，2，…，n)，具有極小拉伸能或極小應變能的C1連續二次均勻B 樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}為保形插值。

證明由式（8），經計算可得

當n=2 時，由引理1 和式（22），可得

當n＞2 時，由引理1 和式（22），可得

由定理2 可知，當插值曲線r(t)具有極小拉伸能時，控制頂點p0應滿足式（13）。將式（13）代入式（23）和式（24），可得=0，故式（13）也為問題（21）的解。因此，具有極小拉伸能的C1連續二次均勻B 樣條插值曲線為保形插值。由定理3 可得，具有極小應變能的C1連續二次均勻B 樣條插值曲線也為保形插值。

證畢。

綜上，利用拉伸能極小或應變能極小構造光順的C1連續二次均勻B 樣條保形插值曲線算法的步驟為：

算法1利用拉伸能極小或應變能極小構造C1連續二次均勻B 樣條保形插值曲線。

步驟1輸入數據點di(i=1，2，…，n)。

步驟2利用式（13）計算控制頂點p0。

步驟3利用式（7）計算控制頂點pi(i=1，2，…，n)。

步驟4由式（1）生成C1連續二次均勻B 樣條插值曲線r(t)={ri(t)，i=0，1，…，n?2}。

4 數值算例與應用

下面通過幾個數值算例說明本文算法的有效性。

例1給定數據點

不同控制頂點p0所對應的C1連續二次均勻B 樣條插值曲線如圖1 所示。

由圖1 可知，選取的控制頂點p0不同，所構造的插值曲線也不同。當p0選取不恰當時，構造的插值曲線也不理想，如圖1（a）所示。在實際應用中，控制頂點的選取往往比較困難，此時可利用本文算法構造光順的保形插值曲線，見圖2。由式（13），可得p0=(?1，2)。

由圖2 可知，本文算法構造的插值曲線較圖1（a）中的插值曲線效果更好，較圖1（b）中的插值曲線保形性更好。

圖1 不同控制頂點p0 所對應的插值曲線Fig.1 The interpolation curves with different p0

圖2 光順的保形插值曲線Fig.2 The fairing shape-preserving interpolation curve

例2設數據di=(xi，yi)取自函數

由式（13），經計算可得p0=(0.130 0，152.760 1)，由本文算法構造的插值曲線（實線）及原函數曲線（虛線）如圖3 所示。由圖3 知，利用本文算法構造的插值曲線能較好地逼近函數，且具有較好的保形性。

圖3 插值于函數的光順保形曲線Fig.3 The fairing shape-preserving curve interpolating the functional

例3考慮封閉插值曲線情形。設數據di=(xi，yi)取自單位圓

當n=4 時，由式（13）可得p0=(1，?1)；當n=8 時，由式（13）可得p0=(1，?0.414 2)。由本文算法構造的插值曲線（實線）及單位圓曲線（虛線）如圖4所示。

圖4 插值于單位圓的光順保形曲線Fig.4 The fairing shape-preserving curve interpolating the unit circle

由圖4 可知，即使給定的數據點只有5 個，利用本文算法構造的插值曲線也能較好地逼近單位圓；當給定的數據點為9 時，利用本文算法構造的插值曲線與單位圓幾乎重合。

例4考慮特定曲線情形。圖5 為利用本文算法生成的α型插值曲線和β型插值曲線，其中β型插值曲線由2 段插值曲線合成。

圖5 特定曲線造型Fig.5 The specific curve modeling

為便于實際應用，在MATLAB 平臺上設計了算法所對應的圖形用戶界面（graphical user interface，GUI）。該GUI 主要由生成插值曲線、保存插值曲線、清除數據等按鈕，以及輸入插值數據點坐標、插值曲線等顯示框組成，如圖6 所示。

圖6 GUIFig.6 The graphical user interface

用戶只需在界面指定位置輸入插值數據點坐標，點擊相應按鈕即可生成、保存光順的C1連續二次均勻B 樣條保形插值曲線，如圖7 所示。

圖7 GUI 演示1Fig.7 The graphical user interface demo 1

用戶也可在界面指定位置輸入插值數據點的橫坐標，利用給定函數計算相應縱坐標，點擊相應按鈕即可生成、保存逼近于給定函數的光順C1連續二次均勻B 樣條保形插值曲線，如圖8 所示。

圖8 GUI 演示2Fig.8 The graphical user interface demo 2

5 結 語

為構造C1連續二次均勻B 樣條插值曲線，首先給出了二次均勻B 樣條插值曲線的分控制頂點與首個控制頂點（即端點條件）的遞推關系式，然后分別利用拉伸能極小和應變能極小對首個控制頂點進行優化選取，從而構造光順的插值曲線。事實表明，具有極小拉伸能的C1連續二次均勻B 樣條插值曲線也具有極小應變能，且所構造的插值曲線為保形插值。進一步，給出了基于MATLAB 平臺設計的GUI，用戶只需輸入插值數據點坐標并點擊按鈕即可獲得光順的C1連續二次均勻B 樣條保形插值曲線。本文算法及給出的GUI 為B 樣條插值曲線的構造提供了新選擇。同時注意到，本文只給出了適用于構造C1連續二次均勻B 樣條插值曲線的算法，若要構造更高次數或滿足更高連續性要求的均勻B樣條插值曲線，則需要推導分控制頂點與端點條件之間的遞推關系，有待下一步研究。