任利娟,張廣鵬,王 元,王妮娜,黃玉美
(西安理工大學機械與精密儀器工程學院,陜西西安710048)
逆向工程是機械、快速原型、醫學等領域廣泛應用的技術,是目前CAD/CAM研究領域的熱點問題,而復雜曲面的重建是逆向工程研究的核心問題。逆向工程中廣泛應用的重構方法是B樣條曲線算法[1-2],B樣條算法又可以分為B樣條逼近曲線算法和B樣條插值曲線算法。B樣條逼近曲線算法,簡單易操作,擬合曲線不通過控制頂點,一般通過增加控制點數和調整節點向量來提高擬合精度;B樣條插值算法能使插值曲線通過所有的型值點,具有一定的實際意義,但其擬合過程需要求解方程組反算控制點的位置,要對數據點和節點向量做特殊的規定,計算過程較為復雜,也可以通過增加插值點數和調整節點向量來提高擬合精度。3次均勻B樣條可以集逼近插值于一體且形狀可調[3],但均勻B樣條曲線在首末端點處的形狀較難控制。
逆向工程通常要求使用盡可能少的控制點數量來實現高精度的曲線曲面重構。在壓縮控制點方面,Piegl等[4]構造了具有非退化特性的參數區間,增加了節點選擇的靈活性,使控制點數目明顯減少;Park等[5-6]在二分法的基礎上,提出一種基于誤差自適應控制的逼近算法,可以得到滿足給定誤差閾值的控制點較少的逼近曲線;魏棟等[7]提出將離散點的曲率特征點作為初始逼近曲線,在誤差較大處增加新的插值點,通過粒子群算法優化控制點的位置,進一步提高逼近精度;程仙國等[8]提出基于B樣條曲線段復雜度節點配置算法,在保證相同逼近精度的條件下使用更少的控制點;張毓華等[9]提出基于幾何信息均布的B樣條曲線節點設置方法,在控制點數量相同的情況下實現較高精度的逼近;江本赤等[10]提出將曲率優勢點作為輪廓約束點構造初始逼近曲線,在需要改善擬合精度的區域增加約束點,直到獲得滿足精度要求的B樣條曲線。
通過合理的參數化方法可以提高B樣條曲線的逼近精度,潘日晶[11]提出利用數據點的參數化和節點向量的自由度,構造在各數據點滿足切向約束的二次B樣條插值曲線,直觀地控制插值曲線達到預期形狀;于謙等[12]提出利用二次B樣條本身的幾何性質進行參數化,使曲線在每個插值點上都滿足指定的切向,但節點矢量的計算和控制點的反算過程較為復雜。
本文結合B樣條逼近曲線和B樣條插值曲線的優點,提出一種插值于控制點的局部形狀可調整的B樣條逼近曲線算法。所提出的算法集逼近、插值于一體,且能夠實現控制點處的切向特性和曲率特性的控制和調整;所提出的算法能使用更少的控制點數量和最簡單的節點矢量設置實現對離散數據點的高精度逼近。
設控制點集D=di{xi,yi,i=0,1,2,…,n},B樣條逼近曲線方程可寫為[13]:
(1)
其中di為控制點,k為逼近曲線的次數,Ni,k(u)稱為k次B樣條基函數,它是由一個稱為節點向量的非遞減的參數u的序列U=(u0≤u1≤…≤un+k+1)所決定的k次分段多項式,其表達式為:
(2)
式中,Ni,k(u)的雙下標中第二下標k表示曲線的次數,第一下標i表示控制點的序號。為了對曲線在端點的行為有較好的控制, 節點矢量的兩端點取重復度k+1, 即u0=u1=…=uk=0,un+1=un+2=…=un+k+1=1。剩下的n-k個內節點(uk+1,uk+2,…,un)可根據不同的方法進行設置。
性質對3次均勻B樣條來說,若節點ui處對應的控制頂點是di,則曲線在參數ui處的點用德布爾算法的遞推公式可表示為:
(3)
所以,只要取合適控制頂點,B樣條曲線可插值其部分控制頂點。另外,B樣條曲線在c(ui)處的切線平行于向量di+1-di-1。
因為3次均勻B樣條曲線在首末端點處的行為很難控制,而3次準均勻B樣條曲線是最常用的,因此把上述的性質推廣到3次準均勻B樣條曲線中,并對得到的推論進行證明。
推論對于3次準均勻B樣條逼近曲線,除首末端點外,三個順序控制點共線且等間距分布時, 逼近曲線會通過中間的控制點,B樣條曲線在c(ui)處的切線平行于向量di+1-di-1。
證明計算參數u∈[ui,ui+1],對應的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)規定0/0=0。


圖1 3次準均勻B樣條曲線上點c(u)的德布爾算法示意圖Fig.1 Schematic diagram of the Debord algorithm for point c(u) on the 3rd quasi-uniform B-spline curve
要證明參數u對應的曲線上的一點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)
證明完畢。
參數u對應的曲線上的點是u∈[ui,ui+1]的曲線段c(u)的分段連接點c(ui+1), 且c(ui+1) 位于di、di-1和di-2所在的直線的中點,即di-1點。
基于上述結論,本文提出一種通過增加輔助控制點的方法,使B樣條逼近曲線通過控制點,并引入形狀參數S和L實現逼近曲線局部形狀的動態調整。
圖2 (a) 是4個控制點的3次準均勻B樣條曲線,圖2(b)是本文算法得到的插值于控制點的3次準均勻B樣條曲線。圖中d1、d2、d3和d4為控制點,d20、d21、d30和d31為控制點d2和d3的輔助控制點。從圖2(b)中可以看出,引入輔助控制點后,逼近曲線通過對應的控制點。

圖2 本文算法實現原理示意圖Fig.2 Schematic diagram of the proposed algorithm implementation principle
1.3.1形狀參數初始化
由于切向參數S和曲率參數L決定了輔助控制點的位置,從而決定了該處擬合曲線的形狀,對擬合曲線的精度有決定性作用。根據大量仿真實驗, 本文給出逆向工程中切向參數S和曲率參數L的初始化方法。
1) 在逆向工程中一般給出離散數據點集P=pq{xq,yq,q=0,1,2,…,h},h為離散數據點的總數??刂泣c處的切向參數S通過控制點所在位置的前后數據點計算,假設控制點di在離散數據點集中為點pq,則Si可表示為:
(18)
2) 參數Li為相鄰控制點之間的距離乘以比例因子λi,可表示為:
(19)
式中,Li0表示位于該控制點前面的輔助控制點與控制點之間的距離;Li1則是位于該控制點后面的輔助控制點與控制點之間的距離;‖·‖表示兩點之間的歐幾里得距離。控制點集一經確定,相鄰控制點之間的距離也確定,可以通過調整比例因子λi來調整輔助控制點的位置。
因為參數λi與該控制點處的曲率特性密切相關,為了降低初始逼近曲線的誤差,本文給出Li初始化方法,如表1所示。

表1 λi與曲率的對應關系
Ki為控制點di在離散點列中通過U弦長[14]計算得到的曲率。當然,初始逼近曲線不可能完全達到精度要求,曲線在兩個控制點之間的形狀也受到切向參數Li的影響,對形狀參數進行適當調整可大幅度提高逼近精度。
1.3.2輔助控制點位置
控制點列記為D=di{xi,yi,i=0,1,2,…,n},對于第i個控制點di(xi,yi),其前輔助控制點di0(xi0,yi0)的位置坐標可以通過聯立幾個方程來求解:
(20)
式中,Si為斜率參數, 在計算機輔助幾何設計(CAGD)中可根據設計曲線在控制點di處的切向特性進行設置;Li為輔助控制點與對應的控制點之間的距離。對應的后面的輔助控制點di1(xi1,yi1)的坐標同理可求。
輔助控制點的位置由兩個參數進行控制和調整:切向參數S和輔助控制點到對應的控制點之間的距離L,如圖3所示。

圖3 參數S和L對逼近曲線形狀的影響Fig.3 Effect of parameters S and L on the shape of the approximation curve
參數S和L對于曲線形狀具有不同的幾何意義。切向參數S決定了逼近曲線在控制點處的切向特性,L則決定了曲線在控制點處的曲率特性,L越小,曲線越陡峭,L越大,曲線越平緩。
1.3.3節點向量
節點向量的確定需要根據控制點個數n和逼近曲線的次數k來確定。本文中增加的輔助控制點個數計入控制點數。如控制點集中有n個控制點, 除去首末端點外其余控制點各增加兩個輔助控制點,增加的輔助控制點個數為2×(n-2),則增加后總的控制點數為:
m=3n-4
(21)
增加輔助控制點后的節點向量根據m和k進行設置,對于開曲線和首末端點僅位置連續的閉曲線,準均勻B樣條的兩端節點取重復度k+1,便于對曲線端點行為有較好的控制,內節點均勻分布,即u0=u1=…=uk=0,un+1=un+2=…=un+k+1=1。剩下的m-k個內節點均勻分布。
1.3.4逼近誤差計算
B樣條曲線的逼近誤差指離散點列到逼近曲線的距離的最大值。本文提出一種逼近誤差的計算方法,設離散點列P=pq{xq,yq,q=0,1,2,…,h},p(u)為參數u對應的逼近曲線上的點,pq點處的誤差εq為:
εq=min‖pqp(u)‖
(22)
其中,參數u均勻取值,通過不斷細化u的取值間距,得到較為準確的誤差為止。誤差的最大值、平均值和均方差可計算為:
(23)
式中,max表示取最大值;mean表示計算平均值;std表示計算均方差。
粒子群優化算法(Particle Swarm Optimization,PSO)的思想源于群鳥捕食行為。用數學模型模擬鳥群覓食的過程,每一個數學問題的解看作是一只鳥,稱作“粒子”,所有的粒子都由一個適應度函數來判斷當前位置的好壞,每一個粒子都有記憶功能,能夠記住搜尋過的最佳位置,每個粒子還有一個速度來決定搜尋的距離和方向,這個速度可以根據它本身的飛行經驗和同伴的飛行經驗進行動態調整。
設有N個粒子在R維搜索空間搜尋,則第t個粒子的位置為[15]:
Xt=(xt1,xt2,…,xtR),t=1,2,…,N
(24)
第t個粒子經歷過的最好位置:
pbest,t=(pt1,pt2,…,ptR)
所有粒子經歷過的最好位置:
gbest=(p1,p2,…,pR)
粒子群的飛行速度更新公式為:
(25)
式中,vtr∈[-vmax,r,vmax,r]。
位置更新公式為:
xtre+1=xtre+vtre+1, (r=1,2,…,R)
(26)

粒子群優化算法在本文中的應用說明:本文要求解的問題是尋找一組解Si和λi使得B樣條逼近曲線的逼近誤差最小。Si和λi的調整范圍根據離散數據及控制點的特征進行設定。適應度函數為逼近誤差的平均值,誤差的計算方法參照1.3.4節內容。
優化輔助控制點位置的PSO算法的具體步驟為:
步驟1 初始化粒子,粒子的初始位置按照參數Si和λi的初始值計算得到;
步驟2 構造逼近曲線,計算初始適應度值,將粒子當前的適應度值ξi作為當前粒子的最優值pbest,i,并記住此時的參數Si和λi,然后選出總體最優值作為gbest;
步驟3 按照式(25)和(26)更新粒子的速度和位置,重新計算復制控制點的位置,構造逼近曲線,計算粒子的適應度值,若ξi 步驟4 若最優值gbest計算的適應度值小于設定的閾值或者迭代次數大于設定的最大迭代次數,則停止計算,輸出pbest,i、Si、λi和gbest,繪制逼近曲線;否則轉步驟3。 步驟1 提取控制點。對于給定的待擬合的離散數據點,提取曲率優勢點作為控制點[12]。 步驟2 確定輔助控制點的初始位置。根據公式(18),計算輔助控制點所在直線的斜率Si;通過控制點所在離散點處的曲率值及表1,設定參數λi,結合式(20),確定輔助控制點的初始位置。 步驟3 節點向量。根據式(21)確定控制點的個數,兩端節點取重復度k+1, 便于對曲線端點行為有較好的控制,內節點均勻分布。3次準均勻B樣條,k=3。 步驟4 由式(1)和(2)繪制B樣條曲線。 步驟5 根據式(22)和(23)計算擬合曲線的誤差,判斷誤差是否滿足要求,若滿足,結束;若不滿足,執行步驟6。 步驟6 對于誤差大于規定閾值的控制點,使用1.4節描述的方法對該控制點對應的輔助控制點位置進行調整,調整后,執行本節步驟4(注:節點向量不用更新)。 B樣條曲線廣泛應用于曲線曲面重構技術中,通過增加控制點的個數來提高B樣條曲線的擬合精度。控制點個數的增加會大大降低擬合的效率。在逆向工程中,待擬合的點基本都是曲線上或者曲面上的點,要求逼近曲線能最大程度通過待擬合的點以降低逼近誤差。 如圖4所示,以某葉片截面的離散數據為例進行分析。離散點數據為234個,分布范圍約50 mm×25 mm。 圖4 某葉片截面離散數據Fig.4 Discrete data of a blade section 應用本文算法對某葉片234個離散數據點進行曲線重構。首先,提取曲率特征點作為控制點,10個順序控制點坐標為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];為了使閉曲線首末端點相連,控制點中首末端點取同一個點; 為了對端點處的曲線形狀也能進行控制和調整, 在兩端點處各增加一個內輔助控制點d11和d100; 根據最終的控制點數確定節點向量為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]。執行流程如1.5節所示。 圖5(a)為初始擬合曲線,從圖中可以看出,擬合曲線和數據點之間有較大的誤差。圖5(b)為對參數進行調整后的擬合曲線。表2為調整后的形狀參數表。從圖5(b)中可以看出,擬合曲線和數據點具有很好的貼合度,擬合曲線光順。通過形狀參數的調整,擬合曲線的平均誤差從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個控制點,通過形狀參數對逼近曲線的局部調整,實現了對離散數據的高精度逼近,且逼近誤差分布較為均勻。 表3 本文算法逼近誤差 圖6為曲率變化較大區域的放大圖,從圖中可以看出,本文算法對曲率變化較大的區域的曲線有較好的控制,相對于其他兩種算法,使用的控制點數較少。 圖6 局部放大圖Fig.6 Enlarged view of regions I and II 應用傳統的3次準均勻B樣條逼近曲線算法進行曲線重構,使用與2.1節同樣的曲率優勢點作為初始控制點,通過逐個增加控制點的方法來提高逼近曲線的精度,直到平均誤差小于0.05。 仿真結果如圖7所示。圖7(a)為初始逼近曲線,從圖中可以看出,逼近曲線比較光順,但逼近曲線與離散點之間有明顯的偏差;圖7(b)為控制點增加到68個時,平均誤差小于0.05的最終逼近曲線,在曲率變化較大的地方控制點較多。 擬合曲線的誤差如表4所示,從表中可以看出,隨著控制點數的增加,逼近精度得到明顯改善,最終的逼近曲線使用了68個控制點。 表4 逼近誤差 B樣條逼近曲線的優點在于算法簡單,擬合曲線較為光順,如圖7所示,可以通過增加控制點提高擬合曲線精度。但即使將所有的點都作為控制點,控制點與擬合曲線之間存在的誤差也無法消除,如圖8所示。 圖7 B樣條逼近算法仿真結果Fig.7 Simulation results of B-spline approximation algorithm 圖8 圖7中III區域放大圖Fig.8 Enlarged view of region III in Fig.7 根據文獻[12]構造的二次B樣條插值曲線算法對離散點列進行擬合,通過特殊的節點矢量設置,增加插值點處的切向約束,使曲線在每個插值點上都滿足指定的切向。通過將誤差最大處的離散點添加到控制點列來提高插值曲線的精度,直到平均誤差小于0.05。 仿真結果如圖9所示。從圖中可以看出,圖9(a)中逼近曲線與離散數據點之間有微小偏差;圖9(b)逼近誤差明顯降低,同樣在曲率變化較大的地方控制點較多。 表5為文獻[12]方法的擬合誤差,在擬合誤差最大值小于0.05時,插值點數為39,最后的平均擬合誤差較小。 對三種方法的控制點數量、初始逼近誤差及最終逼近誤差進行對比分析,結果如表6和表7所示。從表6中可以看出,在使用相同的控制點時,本文方法的初始最大逼近誤差、平均誤差和均方差遠小于其他兩種方法。從表7可以看出,本文方法使用最少的控制點數,實現了與其他兩種方法相近的逼近精度。 表6 初始逼近曲線誤差對比 表7 最終逼近曲線誤差對比 圖10為三種方法最終擬合誤差分布圖,從圖中可以看出,在兩段曲率變化平緩的區域,文獻[12]的逼近精度最好,在中間曲率變化較快的位置,本文算法的誤差得到了較好的控制,且在該區域使用的控制點數遠小于其他兩種算法。 圖10 逼近誤差分布Fig.10 Approximation error distribution 本文對逆向工程中的B樣條曲線重構關鍵技術進行了研究,提出了一種集逼近、插值于一體的B樣條曲線算法,且能夠實現控制點處的切向特性和曲率特性的控制和調整。所提出的曲線重構方法大大壓縮了曲線重構的控制點數量,且算法簡單易操作,其主要結論為: 1) 使用本文方法,通過10個控制點確定的初始B樣條曲線的逼近誤差遠小于傳統的3次準均勻B樣條逼近曲線算法和2次B樣條插值算法; 2) 使用本文算法獲得的最終逼近曲線可達到與其他兩種方法同等的逼近誤差,且使用的控制點數遠少于其他兩種方法; 3) 本文算法對曲率變化較大區域的逼近曲線誤差有較好的控制。1.5 實現步驟
2 應用實例分析

2.1 本文方法




2.2 3次準均勻B樣條逼近算法



2.3 B樣條插值算法
2.4 結果分析



3 結 論