基于壓縮控制點的B樣條曲線重構(gòu)算法

任利娟,張廣鵬,王 元,王妮娜,黃玉美

(西安理工大學(xué)機械與精密儀器工程學(xué)院，陜西西安710048)

逆向工程是機械、快速原型、醫(yī)學(xué)等領(lǐng)域廣泛應(yīng)用的技術(shù)，是目前CAD/CAM研究領(lǐng)域的熱點問題，而復(fù)雜曲面的重建是逆向工程研究的核心問題。逆向工程中廣泛應(yīng)用的重構(gòu)方法是B樣條曲線算法[1-2]，B樣條算法又可以分為B樣條逼近曲線算法和B樣條插值曲線算法。B樣條逼近曲線算法，簡單易操作，擬合曲線不通過控制頂點，一般通過增加控制點數(shù)和調(diào)整節(jié)點向量來提高擬合精度；B樣條插值算法能使插值曲線通過所有的型值點，具有一定的實際意義，但其擬合過程需要求解方程組反算控制點的位置，要對數(shù)據(jù)點和節(jié)點向量做特殊的規(guī)定，計算過程較為復(fù)雜，也可以通過增加插值點數(shù)和調(diào)整節(jié)點向量來提高擬合精度。3次均勻B樣條可以集逼近插值于一體且形狀可調(diào)[3]，但均勻B樣條曲線在首末端點處的形狀較難控制。

逆向工程通常要求使用盡可能少的控制點數(shù)量來實現(xiàn)高精度的曲線曲面重構(gòu)。在壓縮控制點方面，Piegl等[4]構(gòu)造了具有非退化特性的參數(shù)區(qū)間，增加了節(jié)點選擇的靈活性，使控制點數(shù)目明顯減少；Park等[5-6]在二分法的基礎(chǔ)上，提出一種基于誤差自適應(yīng)控制的逼近算法，可以得到滿足給定誤差閾值的控制點較少的逼近曲線；魏棟等[7]提出將離散點的曲率特征點作為初始逼近曲線，在誤差較大處增加新的插值點，通過粒子群算法優(yōu)化控制點的位置，進(jìn)一步提高逼近精度；程仙國等[8]提出基于B樣條曲線段復(fù)雜度節(jié)點配置算法，在保證相同逼近精度的條件下使用更少的控制點；張毓華等[9]提出基于幾何信息均布的B樣條曲線節(jié)點設(shè)置方法，在控制點數(shù)量相同的情況下實現(xiàn)較高精度的逼近；江本赤等[10]提出將曲率優(yōu)勢點作為輪廓約束點構(gòu)造初始逼近曲線，在需要改善擬合精度的區(qū)域增加約束點，直到獲得滿足精度要求的B樣條曲線。

通過合理的參數(shù)化方法可以提高B樣條曲線的逼近精度，潘日晶[11]提出利用數(shù)據(jù)點的參數(shù)化和節(jié)點向量的自由度，構(gòu)造在各數(shù)據(jù)點滿足切向約束的二次B樣條插值曲線，直觀地控制插值曲線達(dá)到預(yù)期形狀；于謙等[12]提出利用二次B樣條本身的幾何性質(zhì)進(jìn)行參數(shù)化，使曲線在每個插值點上都滿足指定的切向，但節(jié)點矢量的計算和控制點的反算過程較為復(fù)雜。

本文結(jié)合B樣條逼近曲線和B樣條插值曲線的優(yōu)點，提出一種插值于控制點的局部形狀可調(diào)整的B樣條逼近曲線算法。所提出的算法集逼近、插值于一體，且能夠?qū)崿F(xiàn)控制點處的切向特性和曲率特性的控制和調(diào)整；所提出的算法能使用更少的控制點數(shù)量和最簡單的節(jié)點矢量設(shè)置實現(xiàn)對離散數(shù)據(jù)點的高精度逼近。

1 改進(jìn)的3次準(zhǔn)均勻B樣條曲線算法

1.1 基本概念及性質(zhì)

設(shè)控制點集D=di{xi,yi,i=0,1,2,…,n}，B樣條逼近曲線方程可寫為[13]：

(1)

其中di為控制點，k為逼近曲線的次數(shù)，Ni,k(u)稱為k次B樣條基函數(shù)，它是由一個稱為節(jié)點向量的非遞減的參數(shù)u的序列U=(u0≤u1≤…≤un+k+1)所決定的k次分段多項式，其表達(dá)式為：

(2)

式中，Ni,k(u)的雙下標(biāo)中第二下標(biāo)k表示曲線的次數(shù)，第一下標(biāo)i表示控制點的序號。為了對曲線在端點的行為有較好的控制, 節(jié)點矢量的兩端點取重復(fù)度k+1, 即u0=u1=…=uk=0，un+1=un+2=…=un+k+1=1。剩下的n-k個內(nèi)節(jié)點(uk+1,uk+2,…,un)可根據(jù)不同的方法進(jìn)行設(shè)置。

性質(zhì)對3次均勻B樣條來說，若節(jié)點ui處對應(yīng)的控制頂點是di，則曲線在參數(shù)ui處的點用德布爾算法的遞推公式可表示為：

(3)

所以，只要取合適控制頂點，B樣條曲線可插值其部分控制頂點。另外，B樣條曲線在c(ui)處的切線平行于向量di+1-di-1。

因為3次均勻B樣條曲線在首末端點處的行為很難控制，而3次準(zhǔn)均勻B樣條曲線是最常用的，因此把上述的性質(zhì)推廣到3次準(zhǔn)均勻B樣條曲線中，并對得到的推論進(jìn)行證明。

推論對于3次準(zhǔn)均勻B樣條逼近曲線，除首末端點外，三個順序控制點共線且等間距分布時， 逼近曲線會通過中間的控制點，B樣條曲線在c(ui)處的切線平行于向量di+1-di-1。

證明計算參數(shù)u∈[ui,ui+1]，對應(yīng)的B樣條曲線上的值可用德布爾算法的遞推公式表示為[13]：

(4)

式中，u∈[ui,ui+1]∈[uk,un+1]。

(5)

(6)

式中，j=i-k,i-k+1,…,i-l，l=1,2,…,k，對式(6)規(guī)定0/0=0。

圖1 3次準(zhǔn)均勻B樣條曲線上點c(u)的德布爾算法示意圖Fig.1 Schematic diagram of the Debord algorithm for point c(u) on the 3rd quasi-uniform B-spline curve

要證明參數(shù)u對應(yīng)的曲線上的一點c(u)位于控制點di-1處，則遞推過程必需滿足以下兩個條件：

(7)

(8)

(9)

(10)

(11)

(12)

(13)

通過式(5)和(6)計算得到：

(14)

把式(14)代入(15)可得：

(15)

由于di-1是di和di-2的中點, 所以：

(16)

把式(16)代入式(15)即可得：

(17)

證明完畢。

參數(shù)u對應(yīng)的曲線上的點是u∈[ui,ui+1]的曲線段c(u)的分段連接點c(ui+1)， 且c(ui+1) 位于di、di-1和di-2所在的直線的中點，即di-1點。

1.2 本文算法的提出

基于上述結(jié)論，本文提出一種通過增加輔助控制點的方法，使B樣條逼近曲線通過控制點，并引入形狀參數(shù)S和L實現(xiàn)逼近曲線局部形狀的動態(tài)調(diào)整。

圖2 (a) 是4個控制點的3次準(zhǔn)均勻B樣條曲線，圖2(b)是本文算法得到的插值于控制點的3次準(zhǔn)均勻B樣條曲線。圖中d1、d2、d3和d4為控制點，d20、d21、d30和d31為控制點d2和d3的輔助控制點。從圖2(b)中可以看出，引入輔助控制點后，逼近曲線通過對應(yīng)的控制點。

圖2 本文算法實現(xiàn)原理示意圖Fig.2 Schematic diagram of the proposed algorithm implementation principle

1.3 方法的實現(xiàn)

1.3.1形狀參數(shù)初始化

由于切向參數(shù)S和曲率參數(shù)L決定了輔助控制點的位置，從而決定了該處擬合曲線的形狀，對擬合曲線的精度有決定性作用。根據(jù)大量仿真實驗， 本文給出逆向工程中切向參數(shù)S和曲率參數(shù)L的初始化方法。

1) 在逆向工程中一般給出離散數(shù)據(jù)點集P=pq{xq,yq,q=0,1,2,…,h}，h為離散數(shù)據(jù)點的總數(shù)。控制點處的切向參數(shù)S通過控制點所在位置的前后數(shù)據(jù)點計算，假設(shè)控制點di在離散數(shù)據(jù)點集中為點pq，則Si可表示為：

(18)

2) 參數(shù)Li為相鄰控制點之間的距離乘以比例因子λi，可表示為：

(19)

式中，Li0表示位于該控制點前面的輔助控制點與控制點之間的距離；Li1則是位于該控制點后面的輔助控制點與控制點之間的距離；‖·‖表示兩點之間的歐幾里得距離?？刂泣c集一經(jīng)確定，相鄰控制點之間的距離也確定，可以通過調(diào)整比例因子λi來調(diào)整輔助控制點的位置。

因為參數(shù)λi與該控制點處的曲率特性密切相關(guān)，為了降低初始逼近曲線的誤差，本文給出Li初始化方法，如表1所示。

表1 λi與曲率的對應(yīng)關(guān)系

Ki為控制點di在離散點列中通過U弦長[14]計算得到的曲率。當(dāng)然，初始逼近曲線不可能完全達(dá)到精度要求，曲線在兩個控制點之間的形狀也受到切向參數(shù)Li的影響，對形狀參數(shù)進(jìn)行適當(dāng)調(diào)整可大幅度提高逼近精度。

1.3.2輔助控制點位置

控制點列記為D=di{xi,yi,i=0,1,2,…,n}，對于第i個控制點di(xi,yi)，其前輔助控制點di0(xi0,yi0)的位置坐標(biāo)可以通過聯(lián)立幾個方程來求解：

(20)

式中，Si為斜率參數(shù), 在計算機輔助幾何設(shè)計(CAGD)中可根據(jù)設(shè)計曲線在控制點di處的切向特性進(jìn)行設(shè)置；Li為輔助控制點與對應(yīng)的控制點之間的距離。對應(yīng)的后面的輔助控制點di1(xi1,yi1)的坐標(biāo)同理可求。

輔助控制點的位置由兩個參數(shù)進(jìn)行控制和調(diào)整：切向參數(shù)S和輔助控制點到對應(yīng)的控制點之間的距離L,如圖3所示。

圖3 參數(shù)S和L對逼近曲線形狀的影響Fig.3 Effect of parameters S and L on the shape of the approximation curve

參數(shù)S和L對于曲線形狀具有不同的幾何意義。切向參數(shù)S決定了逼近曲線在控制點處的切向特性，L則決定了曲線在控制點處的曲率特性，L越小，曲線越陡峭，L越大，曲線越平緩。

1.3.3節(jié)點向量

節(jié)點向量的確定需要根據(jù)控制點個數(shù)n和逼近曲線的次數(shù)k來確定。本文中增加的輔助控制點個數(shù)計入控制點數(shù)。如控制點集中有n個控制點， 除去首末端點外其余控制點各增加兩個輔助控制點，增加的輔助控制點個數(shù)為2×(n-2)，則增加后總的控制點數(shù)為：

m=3n-4

(21)

增加輔助控制點后的節(jié)點向量根據(jù)m和k進(jìn)行設(shè)置，對于開曲線和首末端點僅位置連續(xù)的閉曲線，準(zhǔn)均勻B樣條的兩端節(jié)點取重復(fù)度k+1，便于對曲線端點行為有較好的控制，內(nèi)節(jié)點均勻分布，即u0=u1=…=uk=0，un+1=un+2=…=un+k+1=1。剩下的m-k個內(nèi)節(jié)點均勻分布。

1.3.4逼近誤差計算

B樣條曲線的逼近誤差指離散點列到逼近曲線的距離的最大值。本文提出一種逼近誤差的計算方法，設(shè)離散點列P=pq{xq,yq,q=0,1,2,…,h}，p(u)為參數(shù)u對應(yīng)的逼近曲線上的點，pq點處的誤差εq為：

εq=min‖pqp(u)‖

(22)

其中，參數(shù)u均勻取值，通過不斷細(xì)化u的取值間距，得到較為準(zhǔn)確的誤差為止。誤差的最大值、平均值和均方差可計算為：

(23)

式中，max表示取最大值;mean表示計算平均值;std表示計算均方差。

1.4 粒子群算法優(yōu)化輔助控制點的位置

粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)的思想源于群鳥捕食行為。用數(shù)學(xué)模型模擬鳥群覓食的過程，每一個數(shù)學(xué)問題的解看作是一只鳥，稱作“粒子”，所有的粒子都由一個適應(yīng)度函數(shù)來判斷當(dāng)前位置的好壞，每一個粒子都有記憶功能，能夠記住搜尋過的最佳位置，每個粒子還有一個速度來決定搜尋的距離和方向，這個速度可以根據(jù)它本身的飛行經(jīng)驗和同伴的飛行經(jīng)驗進(jìn)行動態(tài)調(diào)整。

設(shè)有N個粒子在R維搜索空間搜尋，則第t個粒子的位置為[15]：

Xt=(xt1,xt2,…,xtR),t=1,2,…,N

(24)

第t個粒子經(jīng)歷過的最好位置：

pbest,t=(pt1,pt2,…,ptR)

所有粒子經(jīng)歷過的最好位置：

gbest=(p1,p2,…,pR)

粒子群的飛行速度更新公式為：

(25)

式中，vtr∈[-vmax,r,vmax,r]。

位置更新公式為：

xtre+1=xtre+vtre+1, (r=1,2,…,R)

(26)

粒子群優(yōu)化算法在本文中的應(yīng)用說明：本文要求解的問題是尋找一組解Si和λi使得B樣條逼近曲線的逼近誤差最小。Si和λi的調(diào)整范圍根據(jù)離散數(shù)據(jù)及控制點的特征進(jìn)行設(shè)定。適應(yīng)度函數(shù)為逼近誤差的平均值，誤差的計算方法參照1.3.4節(jié)內(nèi)容。

優(yōu)化輔助控制點位置的PSO算法的具體步驟為：

步驟1 初始化粒子，粒子的初始位置按照參數(shù)Si和λi的初始值計算得到；

步驟2 構(gòu)造逼近曲線，計算初始適應(yīng)度值，將粒子當(dāng)前的適應(yīng)度值ξi作為當(dāng)前粒子的最優(yōu)值pbest，i，并記住此時的參數(shù)Si和λi，然后選出總體最優(yōu)值作為gbest；

步驟3 按照式(25)和(26)更新粒子的速度和位置，重新計算復(fù)制控制點的位置，構(gòu)造逼近曲線，計算粒子的適應(yīng)度值，若ξi

步驟4 若最優(yōu)值gbest計算的適應(yīng)度值小于設(shè)定的閾值或者迭代次數(shù)大于設(shè)定的最大迭代次數(shù)，則停止計算，輸出pbest,i、Si、λi和gbest，繪制逼近曲線；否則轉(zhuǎn)步驟3。

1.5 實現(xiàn)步驟

步驟1 提取控制點。對于給定的待擬合的離散數(shù)據(jù)點，提取曲率優(yōu)勢點作為控制點[12]。

步驟2 確定輔助控制點的初始位置。根據(jù)公式(18)，計算輔助控制點所在直線的斜率Si；通過控制點所在離散點處的曲率值及表1，設(shè)定參數(shù)λi，結(jié)合式(20)，確定輔助控制點的初始位置。

步驟3 節(jié)點向量。根據(jù)式(21)確定控制點的個數(shù)，兩端節(jié)點取重復(fù)度k+1, 便于對曲線端點行為有較好的控制，內(nèi)節(jié)點均勻分布。3次準(zhǔn)均勻B樣條，k=3。

步驟4 由式(1)和(2)繪制B樣條曲線。

步驟5 根據(jù)式(22)和(23)計算擬合曲線的誤差，判斷誤差是否滿足要求，若滿足，結(jié)束；若不滿足，執(zhí)行步驟6。

步驟6 對于誤差大于規(guī)定閾值的控制點，使用1.4節(jié)描述的方法對該控制點對應(yīng)的輔助控制點位置進(jìn)行調(diào)整，調(diào)整后，執(zhí)行本節(jié)步驟4(注：節(jié)點向量不用更新)。

2 應(yīng)用實例分析

B樣條曲線廣泛應(yīng)用于曲線曲面重構(gòu)技術(shù)中，通過增加控制點的個數(shù)來提高B樣條曲線的擬合精度。控制點個數(shù)的增加會大大降低擬合的效率。在逆向工程中，待擬合的點基本都是曲線上或者曲面上的點，要求逼近曲線能最大程度通過待擬合的點以降低逼近誤差。

如圖4所示，以某葉片截面的離散數(shù)據(jù)為例進(jìn)行分析。離散點數(shù)據(jù)為234個，分布范圍約50 mm×25 mm。

圖4 某葉片截面離散數(shù)據(jù)Fig.4 Discrete data of a blade section

2.1 本文方法

應(yīng)用本文算法對某葉片234個離散數(shù)據(jù)點進(jìn)行曲線重構(gòu)。首先，提取曲率特征點作為控制點，10個順序控制點坐標(biāo)為P0=[-20.5, 24; -12.644 5, 18.738 9; 2.815 0, 11.961 8; 13.704 1, 10.389 8; 21.101 1, 10.447 5; 15.786 2, 5.594 7; 3.412 7, 3.713 5; -9.062 0, 9.722 9; -19.182 6, 20.437 7; -20.5, 24]；為了使閉曲線首末端點相連，控制點中首末端點取同一個點； 為了對端點處的曲線形狀也能進(jìn)行控制和調(diào)整， 在兩端點處各增加一個內(nèi)輔助控制點d11和d100； 根據(jù)最終的控制點數(shù)確定節(jié)點向量為NodeVector=[0, 0, 0, 0, 0.04, 0.08, 0.12, 0.16, 0.20, 0.24, 0.28, 0.32, 0.36, 0.40, 0.44, 0.48, 0.52, 0.56, 0.60, 0.64, 0.68, 0.72, 0.76, 0.80, 0.84, 0.88, 0.92, 0.96, 1, 1, 1, 1]。執(zhí)行流程如1.5節(jié)所示。

圖5(a)為初始擬合曲線，從圖中可以看出，擬合曲線和數(shù)據(jù)點之間有較大的誤差。圖5(b)為對參數(shù)進(jìn)行調(diào)整后的擬合曲線。表2為調(diào)整后的形狀參數(shù)表。從圖5(b)中可以看出，擬合曲線和數(shù)據(jù)點具有很好的貼合度，擬合曲線光順。通過形狀參數(shù)的調(diào)整，擬合曲線的平均誤差從0.019 6 mm降低為0.007 5 mm，最大誤差從0.561 5 mm降低為0.034 8 mm。

圖5 本文算法擬合曲線Fig.5 Fitting curve using proposed algorithm

控制點Siλi0λi111.000.302-0.651.001.053-0.261.000.804-0.010.801.005-1.870.320.3560.500.901.007-0.160.800.908-0.751.001.009-1.501.000.30101.000.5

本文算法的逼近誤差如表3所示，可以看出，本文使用10個控制點，通過形狀參數(shù)對逼近曲線的局部調(diào)整，實現(xiàn)了對離散數(shù)據(jù)的高精度逼近，且逼近誤差分布較為均勻。

表3 本文算法逼近誤差

圖6為曲率變化較大區(qū)域的放大圖，從圖中可以看出，本文算法對曲率變化較大的區(qū)域的曲線有較好的控制，相對于其他兩種算法，使用的控制點數(shù)較少。

圖6 局部放大圖Fig.6 Enlarged view of regions I and II

2.2 3次準(zhǔn)均勻B樣條逼近算法

應(yīng)用傳統(tǒng)的3次準(zhǔn)均勻B樣條逼近曲線算法進(jìn)行曲線重構(gòu)，使用與2.1節(jié)同樣的曲率優(yōu)勢點作為初始控制點，通過逐個增加控制點的方法來提高逼近曲線的精度，直到平均誤差小于0.05。

仿真結(jié)果如圖7所示。圖7(a)為初始逼近曲線，從圖中可以看出，逼近曲線比較光順，但逼近曲線與離散點之間有明顯的偏差；圖7(b)為控制點增加到68個時，平均誤差小于0.05的最終逼近曲線，在曲率變化較大的地方控制點較多。

擬合曲線的誤差如表4所示，從表中可以看出，隨著控制點數(shù)的增加，逼近精度得到明顯改善，最終的逼近曲線使用了68個控制點。

表4 逼近誤差

B樣條逼近曲線的優(yōu)點在于算法簡單，擬合曲線較為光順，如圖7所示，可以通過增加控制點提高擬合曲線精度。但即使將所有的點都作為控制點，控制點與擬合曲線之間存在的誤差也無法消除，如圖8所示。

圖7 B樣條逼近算法仿真結(jié)果Fig.7 Simulation results of B-spline approximation algorithm

圖8 圖7中III區(qū)域放大圖Fig.8 Enlarged view of region III in Fig.7

2.3 B樣條插值算法

根據(jù)文獻(xiàn)[12]構(gòu)造的二次B樣條插值曲線算法對離散點列進(jìn)行擬合，通過特殊的節(jié)點矢量設(shè)置，增加插值點處的切向約束，使曲線在每個插值點上都滿足指定的切向。通過將誤差最大處的離散點添加到控制點列來提高插值曲線的精度，直到平均誤差小于0.05。

仿真結(jié)果如圖9所示。從圖中可以看出，圖9(a)中逼近曲線與離散數(shù)據(jù)點之間有微小偏差；圖9(b)逼近誤差明顯降低，同樣在曲率變化較大的地方控制點較多。

表5為文獻(xiàn)[12]方法的擬合誤差，在擬合誤差最大值小于0.05時，插值點數(shù)為39，最后的平均擬合誤差較小。

2.4 結(jié)果分析

對三種方法的控制點數(shù)量、初始逼近誤差及最終逼近誤差進(jìn)行對比分析，結(jié)果如表6和表7所示。從表6中可以看出，在使用相同的控制點時，本文方法的初始最大逼近誤差、平均誤差和均方差遠(yuǎn)小于其他兩種方法。從表7可以看出，本文方法使用最少的控制點數(shù)，實現(xiàn)了與其他兩種方法相近的逼近精度。

表6 初始逼近曲線誤差對比

表7 最終逼近曲線誤差對比

圖10為三種方法最終擬合誤差分布圖，從圖中可以看出，在兩段曲率變化平緩的區(qū)域，文獻(xiàn)[12]的逼近精度最好，在中間曲率變化較快的位置，本文算法的誤差得到了較好的控制，且在該區(qū)域使用的控制點數(shù)遠(yuǎn)小于其他兩種算法。

圖10 逼近誤差分布Fig.10 Approximation error distribution

3 結(jié) 論

本文對逆向工程中的B樣條曲線重構(gòu)關(guān)鍵技術(shù)進(jìn)行了研究，提出了一種集逼近、插值于一體的B樣條曲線算法，且能夠?qū)崿F(xiàn)控制點處的切向特性和曲率特性的控制和調(diào)整。所提出的曲線重構(gòu)方法大大壓縮了曲線重構(gòu)的控制點數(shù)量，且算法簡單易操作，其主要結(jié)論為：

1) 使用本文方法，通過10個控制點確定的初始B樣條曲線的逼近誤差遠(yuǎn)小于傳統(tǒng)的3次準(zhǔn)均勻B樣條逼近曲線算法和2次B樣條插值算法；

2) 使用本文算法獲得的最終逼近曲線可達(dá)到與其他兩種方法同等的逼近誤差，且使用的控制點數(shù)遠(yuǎn)少于其他兩種方法；

3) 本文算法對曲率變化較大區(qū)域的逼近曲線誤差有較好的控制。