李軍成，劉成志，趙文才

（湖南人文科技學(xué)院數(shù)學(xué)與金融學(xué)院，湖南婁底 417000）

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

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

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

1 預(yù)備知識

1.1 二次均勻B 樣條曲線

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

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

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

1.2 曲線的內(nèi)能

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

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

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

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

2 問題的提出

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

由式（1）知，要構(gòu)造插值曲線r(t)，只需確定各段曲線ri(t) (i=0，1，…，n?2)的控制頂點pk(k=0，1，…，n)。由式（2）知，要使曲線r(t)插值于數(shù)據(jù)點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作為端點條件的情形。

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

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

證明用數(shù)學(xué)歸納法證明。

當(dāng)i=1 時，由式（6），可得式（7）。 假設(shè)當(dāng)i=m時，有

則當(dāng)i=m+1 時，由式（6）知

即式（7）成立。

證畢。

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

下面討論如何利用內(nèi)能極小優(yōu)化選取控制頂點p0，構(gòu)造光順的C1連續(xù)二次均勻B 樣條插值曲線。

3 利用內(nèi)能極小優(yōu)化端點條件

由于對二次均勻B 樣條曲線求三階導(dǎo)數(shù)后歸零，故只討論如何通過拉伸能極小與應(yīng)變能極小優(yōu)化選取控制頂點p0，構(gòu)造光順的C1連續(xù)二次均勻B樣條插值曲線。

3.1 利用拉伸能極小優(yōu)化端點條件

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

由式（8）可得

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

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

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

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

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

從而有，

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

證畢。

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

證明當(dāng)n=2 時，由引理1 和式（10），可得

當(dāng)n＞2 時，由引理1 和式（10），可得

證畢。

3.2 利用應(yīng)變能極小優(yōu)化端點條件

式（8）經(jīng)計算可得

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

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

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

證明當(dāng)n=2 時，由引理1 和式（17），可得

當(dāng)n＞2 時，由引理1 和式（17），可得

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

證畢。

3.3 保形插值

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

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

證明由式（8），經(jīng)計算可得

當(dāng)n=2 時，由引理1 和式（22），可得

當(dāng)n＞2 時，由引理1 和式（22），可得

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

證畢。

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

算法1利用拉伸能極小或應(yīng)變能極小構(gòu)造C1連續(xù)二次均勻B 樣條保形插值曲線。

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

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

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

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

4 數(shù)值算例與應(yīng)用

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

例1給定數(shù)據(jù)點

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

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

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

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

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

例2設(shè)數(shù)據(jù)di=(xi，yi)取自函數(shù)

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

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

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

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

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

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

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

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

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

圖6 GUIFig.6 The graphical user interface

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

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

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

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

5 結(jié) 語

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