池寶濤 張見明 鞠傳明

1.湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長沙，4100822.湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，長沙，410082

0 引言

直線與參數(shù)空間NURBS(non-uniform rational B-splines)曲線、直線與NURBS曲面求交問題是計(jì)算機(jī)輔助幾何設(shè)計(jì)(computer aided geometric design, CAGD)的一個(gè)基本問題，也是其中非常重要和復(fù)雜的問題之一。曲線、曲面求交算法是CAD/CAE系統(tǒng)中的一個(gè)核心算法，廣泛應(yīng)用于實(shí)體造型、曲面裁剪等問題中。區(qū)間分析[1]是數(shù)值計(jì)算領(lǐng)域一個(gè)新的分支，它是以“區(qū)間數(shù)”為研究對(duì)象的一種新的數(shù)學(xué)分析方法。直線與參數(shù)空間NURBS曲線、直線與NURBS曲面[2]求交的本質(zhì)問題就是對(duì)非線性方程或非線性方程組進(jìn)行求解，求交算法的穩(wěn)定性、精度和效率直接影響算法的實(shí)用性[3]。

當(dāng)前CAD幾何造型系統(tǒng)中，直線與參數(shù)空間曲線或曲面的求交方法主要?dú)w納為兩大類：一類為離散分割近似方法，其特點(diǎn)是易于實(shí)現(xiàn)，計(jì)算精度較低，存儲(chǔ)量較大，效率不高；另一類為數(shù)值計(jì)算方法，其優(yōu)點(diǎn)是計(jì)算精度高，幾何數(shù)據(jù)存儲(chǔ)量較少，缺點(diǎn)是一般只適用于給定曲面方程的曲面，對(duì)于較為復(fù)雜的方程組幾乎無法求解。按照方法來分，離散分割的近似方法主要包括離散細(xì)分法[4]、網(wǎng)格劃分法[5]等。數(shù)值計(jì)算方法主要包括代數(shù)解析法、點(diǎn)迭代法、跟蹤求交法[6]、同倫迭代法[7]等。

本文提出了改進(jìn)的基于仿射算術(shù)[8]和區(qū)間運(yùn)算的直線與參數(shù)空間NURBS曲線、直線與NURBS曲面求交的算法。與傳統(tǒng)的離散分割近似方法相比，該方法由于采用區(qū)間運(yùn)算，故細(xì)分的閾值可以任意給定，且無需考慮細(xì)分的曲線段或曲面片是否近似線性，從而細(xì)分層數(shù)得到很大程度的減少；與傳統(tǒng)的點(diǎn)迭代法相比，該方法由于采用了區(qū)間運(yùn)算，迭代過程不需給定合適的迭代初始值，具有更好的靈活性；與傳統(tǒng)的區(qū)間迭代法相比，該方法放寬了對(duì)初始區(qū)間的要求，采用基于線曲率和面曲率的子域分解方法，可以快速篩選預(yù)迭代區(qū)間，提高迭代效率，另外，通過運(yùn)用仿射算術(shù)，考慮計(jì)算過程中數(shù)據(jù)的相關(guān)性，可以大大減少迭代次數(shù)，有效彌補(bǔ)了區(qū)間算法的局限性，在相同的容差要求下，可提高迭代求交的效率。

1 算法基本理論

1.1 直線、NURBS曲線與NURBS曲面的定義

(1)直線。設(shè)直線的起始點(diǎn)為O，直線的單位方向矢量為D，則直線L的參數(shù)方程可表示為

L(t)=O+tDt∈(-∞,+∞)

(1)

(2)NURBS曲線。一條p(p≥0)次NURBS曲線定義為

(2)

式中，Pi為控制點(diǎn)；wi為權(quán)因子；Ni,p(t)為定義在非均勻節(jié)點(diǎn)矢量T上的p次樣條基函數(shù)。

基于De Boor-Cox遞歸公式，第i個(gè)p次樣條基函數(shù)Ni,p(t)定義如下：

(3)

(3)NURBS曲面。一張?jiān)趗方向p次和v方向q次的NURBS曲面或T-Spline曲面定義為

(4)

0≤u,v≤1

式中，Pi,j為控制點(diǎn)；wi,j為權(quán)因子；Ni,p(u)、Nj,q(v)分別為定義在非均勻節(jié)點(diǎn)矢量U和V上的非有理樣條基函數(shù)。

1.2 區(qū)間分析

區(qū)間分析的主要思想是在整個(gè)數(shù)值計(jì)算過程以區(qū)間為運(yùn)算對(duì)象，用區(qū)間代替取值不精確的浮點(diǎn)數(shù)進(jìn)行計(jì)算，從而自動(dòng)地記錄計(jì)算機(jī)浮點(diǎn)運(yùn)算中所產(chǎn)生的截尾誤差和舍入誤差，以區(qū)間形式輸出計(jì)算結(jié)果，并且結(jié)果區(qū)間能夠保證包含數(shù)據(jù)運(yùn)算的真實(shí)結(jié)果，從而保證計(jì)算的可靠性。

1.2.1Newton算子

Newton算子定義如下：

(5)

其中，X為區(qū)間向量。

假定f:Rn→Rn，F(xiàn)為f的區(qū)間擴(kuò)展[9]，F(xiàn)′為F的一階導(dǎo)數(shù)，則新區(qū)間X(k+1)及Newton算子區(qū)間迭代序列[10]為

(6)

k=0,1,2,…

由區(qū)間除法的定義可知，區(qū)間Newton法的區(qū)間劃分判斷是自適應(yīng)的，這是區(qū)間Newton法的最大優(yōu)勢。

Newton算子具有如下性質(zhì)。

假定映像f:X?Rn→Rn,求解非線性方程組f(x)=0,x∈X?Rn。

(1)若x*為f(x)=0的解，且x*∈X，則x*∈N(X)。

(2)若X∩N(X)=?，則f(x)在X中無解。

(3)若N(X)?X，則f(x)在X中必有解存在，當(dāng)且僅當(dāng)W(N(X))?W(X)時(shí)，f(x)在X中存在唯一解。

1.2.2Krawczyk算子

Krawczyk算子定義如下：

K(X)=y-Yf(y)+[I-YF′(X)](X-y)

(7)

Krawczyk算子[11]區(qū)間迭代規(guī)則如下：

(8)

k=0,1,2,…

其中，y∈X，I為n階單位矩陣，Y為n階非奇異矩陣。特別地，若在式(8)中可取y=m(X),Y=[m(F′(X))]-1。

Krawczyk算子具有如下性質(zhì)。

假定映像f:X?Rn→Rn,求解非線性方程組f(x)=0,x∈X?Rn。

(1)若x*為f(x)=0的解，且x*∈X，則x*∈K(X)。

(2)若X∩K(X)=?，則f(x)在X中無解。

(3)若K(X)?X，則f(x)在X中必有解存在，當(dāng)且僅當(dāng)W(K(X))?W(X)，f(x)在X中存在唯一解。

1.2.3區(qū)間迭代幾何意義

對(duì)一任意函數(shù)f(x)，設(shè)初始迭代區(qū)間為[a,b]，圖1給出了f(x)在不同迭代區(qū)間存在唯一解或多個(gè)解的情況示意圖。

(a)迭代區(qū)間有唯一解(b)迭代區(qū)間有多個(gè)解圖1 區(qū)間迭代的幾何意義Fig.1 Geometrical interpretation of the univariate interval Newton method

1.3 仿射算術(shù)

區(qū)間算術(shù)的缺點(diǎn)在于其保守性，經(jīng)區(qū)間算術(shù)運(yùn)算得到的結(jié)果區(qū)間往往被很大程度地放大，放大倍數(shù)取決于區(qū)間運(yùn)算計(jì)算鏈的長度，理想的實(shí)際區(qū)間為結(jié)果區(qū)間的子區(qū)間，而且結(jié)果區(qū)間往往比理想的實(shí)際區(qū)間大得多，這一問題在長計(jì)算鏈的區(qū)間運(yùn)算中尤為突出，因此會(huì)產(chǎn)生“誤差爆炸”現(xiàn)象，而這樣的長計(jì)算鏈在實(shí)際計(jì)算中經(jīng)常出現(xiàn)[12]。仿射算術(shù)為以解決區(qū)間運(yùn)算中相關(guān)性問題為目的而建立的有效的數(shù)值計(jì)算模型，該模型可以解決區(qū)間算術(shù)過于保守的問題。與區(qū)間算術(shù)一樣，仿射算術(shù)能夠自動(dòng)記錄浮點(diǎn)數(shù)的截尾和舍入誤差，此外，仿射算術(shù)還能保持輸入和計(jì)算中各個(gè)不確定量之間的依賴關(guān)系，并自動(dòng)應(yīng)用于原始運(yùn)算中。基于此優(yōu)點(diǎn)，仿射算術(shù)能得到比區(qū)間算術(shù)區(qū)間范圍更小更合理的區(qū)間，特別是在實(shí)際計(jì)算的長計(jì)算鏈中優(yōu)勢更加明顯。

在仿射算術(shù)中，一個(gè)不確定量x的仿射形式為

給定兩個(gè)仿射形式進(jìn)行運(yùn)算：

2 直線與NURBS曲線快速求交算法

2.1 構(gòu)建參數(shù)空間直線與NURBS曲線求交目標(biāo)函數(shù)

在二維參數(shù)空間中，直線與NURBS曲線交點(diǎn)(u,v)應(yīng)滿足

(9)

和

(10)

其中，Pui、Pvi分別為NURBS曲線控制點(diǎn)的坐標(biāo)，(u0,v0)及(u1,v1)為直線上兩個(gè)給定點(diǎn)的參數(shù)坐標(biāo)。在二維參數(shù)空間求解直線與NURBS曲線交點(diǎn)的目標(biāo)函數(shù)如下。

(1) 當(dāng)u1≠u0且v1≠v0時(shí)

(11)

(2)當(dāng)u1≠u0且v1=v0時(shí)

(12)

(3)當(dāng)u1=u0且v1≠v0時(shí)

(13)

2.2 基于仿射算術(shù)的Newton算子運(yùn)算

對(duì)于Newton區(qū)間迭代法，利用區(qū)間除法的性質(zhì)，對(duì)迭代區(qū)間進(jìn)行自適應(yīng)劃分，這是Newton區(qū)間迭代法的優(yōu)勢，而在樣條基函數(shù)Ni,p(t)的區(qū)間運(yùn)算中需要多次調(diào)用遞歸運(yùn)算，屬于長計(jì)算鏈的區(qū)間運(yùn)算，在計(jì)算過程中如果迭代區(qū)間多次累計(jì)放大，會(huì)導(dǎo)致迭代次數(shù)呈非線性增長，區(qū)間迭代效率降低或出現(xiàn)不收斂的情況；此外，對(duì)不包含真實(shí)解的無效區(qū)間進(jìn)行判斷是影響區(qū)間迭代效率的另一主要原因，如何快速準(zhǔn)確地定位包含真實(shí)解的有效區(qū)間直接影響算法的效率。

(a)NURBS曲線 (b)細(xì)分步驟(1)

(c)細(xì)分步驟(2) (d)構(gòu)建包圍盒

(e)直線與參數(shù)空間NURBS曲線求交圖2 曲線離散和直線與包圍盒求交Fig.2 Curve discretization and intersection of line and boxes

算法具體步驟描述如下：

(2)對(duì)離散的曲線段創(chuàng)建包圍盒(圖2d、圖2e)。如圖2d、圖2e所示，判斷直線與包圍盒是否相交，若不相交則去除這些包圍盒對(duì)應(yīng)曲線段的區(qū)間，避免對(duì)這些不包含交點(diǎn)的區(qū)間進(jìn)行無效迭代，同時(shí)篩選出與直線相交的包圍盒作為下一步進(jìn)行區(qū)間迭代的預(yù)迭代區(qū)間。

(3)對(duì)篩選的區(qū)間進(jìn)行區(qū)間Newton迭代。若存在多個(gè)解，迭代過程中迭代區(qū)間會(huì)進(jìn)行自適應(yīng)劃分，進(jìn)一步對(duì)細(xì)分區(qū)間進(jìn)行判斷；若無解，則停止迭代,終止計(jì)算；若存在唯一解，則在這一區(qū)間內(nèi)給定一個(gè)初值，通過牛頓迭代法或擬牛頓迭代法求出真解。

3 直線與NURBS曲面快速求交算法

3.1 構(gòu)建直線與NURBS曲面求交目標(biāo)函數(shù)

(14)

0≤u,v≤1

圖3 絕對(duì)坐標(biāo)系下的直線與NURBS曲面Fig.3 Line and NURBS surface in the absolute coordinate system

圖4 局部坐標(biāo)系下的直線與NURBS曲面Fig.4 Line and NURBS surface in the local coordinate system

(15)

3.2 基于仿射算術(shù)的Krawczyk算子運(yùn)算

直線與NURBS曲面求交中同樣存在直線與NURBS曲線求交遇到的問題，另外，在Krawczyk算子計(jì)算中也需要考慮相關(guān)性問題。

算法具體步驟描述如下：

(1)對(duì)NURBS曲面進(jìn)行區(qū)域劃分。在二維參數(shù)空間考慮黎曼度量和面曲率生成二叉樹網(wǎng)格，對(duì)于裁剪曲面無需將二叉樹網(wǎng)格與邊界進(jìn)行求交處理。為簡化計(jì)算，采用弦長比作為衡量曲面曲率的標(biāo)準(zhǔn)。如圖5所示，曲線ab為細(xì)分單元內(nèi)一條等參線所對(duì)應(yīng)曲線的剖視圖，根據(jù)面曲率對(duì)目標(biāo)面進(jìn)行細(xì)分，然后判斷各子域是否滿足給定細(xì)分閾值，若滿足則將此子域細(xì)分，否則停止細(xì)分，直到單元內(nèi)的面曲率滿足預(yù)給定閾值，然后將細(xì)分結(jié)果映射到三維空間。

(2)對(duì)目標(biāo)面的剖分區(qū)域創(chuàng)建體包圍盒。將直線與體包圍盒進(jìn)行相交判斷，篩選出與直線相交的包圍盒作為下一步進(jìn)行區(qū)間迭代的預(yù)迭代區(qū)間。

(3)對(duì)篩選的區(qū)間進(jìn)行區(qū)間迭代。若在u,v的某一區(qū)域有多個(gè)解，對(duì)U、V采用四叉樹方法進(jìn)行細(xì)分，分成四個(gè)小區(qū)域，對(duì)各個(gè)子區(qū)域計(jì)算Krawczyk算子，并對(duì)這些子區(qū)域進(jìn)行判斷。若無解則終止計(jì)算；若存在唯一解，則在這一區(qū)域內(nèi)給定一初值，通過牛頓迭代法或擬牛頓迭代法迭代求出真解。

(4)在計(jì)算過程中對(duì)Krawczyk算子進(jìn)行改進(jìn)：

X(0)=X=M+Rε= (m1+r1ε1,m2+r2ε2,…,mn+rnεn)′
X(k+1)=X(k)∩K(X(k))

將X(k+1)轉(zhuǎn)化為仿射形式，應(yīng)用仿射算術(shù)計(jì)算G(X(k))=Y(k)F′(X(k)),將仿射矩陣G(X(k))轉(zhuǎn)化為新的區(qū)間矩陣H(X(k))。代入Krawczyk算子公式，應(yīng)用區(qū)間算術(shù)計(jì)算

K(X(k))=m(X(k))-Y(k)f(m(X(k)))+
[I-H(X(k))](Rε)(k)

(a)NURBS曲線 (b)細(xì)分步驟(1)

(c)細(xì)分步驟(2) (d)細(xì)分步驟(3)

(e)NURBS曲面細(xì)分實(shí)例 (f)曲面細(xì)分俯視圖圖5 基于曲率的NURBS曲面細(xì)分Fig.5 NURBS surface subdivision based on curvature

4 算法實(shí)例驗(yàn)證

本文提出一種改進(jìn)的基于仿射算術(shù)和區(qū)間運(yùn)算的直線與參數(shù)空間NURBS曲線及直線與NURBS曲面快速求交算法，分別驗(yàn)證了三種情況下的直線與5次NURBS參數(shù)空間曲線求交，并將本文改進(jìn)的方法與傳統(tǒng)的區(qū)間迭代求交方法進(jìn)行比較，如圖6及表1所示。其中，M1為傳統(tǒng)區(qū)間迭代求交方法，M2為本文改進(jìn)的區(qū)間迭代求交方法。為驗(yàn)證本文算法的可行性與一般性，將本文算法與傳統(tǒng)的區(qū)間迭代求交算法進(jìn)行對(duì)比，驗(yàn)證了直線與控制點(diǎn)為7×3的雙三次NURBS曲面、直線與控制點(diǎn)為12×3的雙三次NURBS曲面，在相應(yīng)容差要求下的計(jì)算結(jié)果和時(shí)間的對(duì)比情況，如圖7及表2所示。將本文改進(jìn)方法與直線與復(fù)雜曲面求交或?qū)嶋H工程算例進(jìn)行對(duì)比，在相應(yīng)容差要求下的計(jì)算結(jié)果和時(shí)間對(duì)比情況如圖8、圖9及表3、表4所示。其中，M3為文獻(xiàn)[13]中直線與復(fù)雜曲面求交的方法，M4為文獻(xiàn)[14]中直線與彎管曲面求交的方法。圖8為直線與控制點(diǎn)為13×3的u方向四次、v方向三次的復(fù)雜NURBS曲面求交算例[14]。圖9為直線與控制點(diǎn)為14×43的u方向四次、v方向三次的復(fù)雜彎管曲面求交的工程算例[15]。計(jì)算結(jié)果表明采用本文方法計(jì)算多個(gè)交點(diǎn)的求解時(shí)間與其他算例中的平均時(shí)間相當(dāng)或優(yōu)于其他算例求解時(shí)間，證明了本文算法的可行性。

(a)直線u=u0與參數(shù)空間NURBS曲線求交

(b)直線v=v0與參數(shù)空間NURBS曲線求交

(c)一般直線與參數(shù)空間NURBS曲線求交圖6 直線與參數(shù)空間NURBS曲線求交Fig.6 Intersection between the line and the NURBS curve

直線誤差(mm) 計(jì)算時(shí)間(ms)M1M211.208×10-1427.517.126.758×10-1448.226.437.933×10-1452.628.3

(a)直線與NURBS曲面求交存在一個(gè)交點(diǎn)情況

(b)直線與NURBS曲面求交存在兩個(gè)交點(diǎn)情況

(c)直線與NURBS曲面求交存在多個(gè)交點(diǎn)情況圖7 直線與NURBS曲面求交Fig.7 Intersection between the line and the NURBS surface

直線誤差(mm) 計(jì)算時(shí)間(ms)M1M213.238×10-1457.837.126.712×10-1469.745.934.693×10-1371.449.6

圖8 直線與復(fù)雜NURBS曲面求交Fig.8 Intersection between the line and the complex NURBS surface

(a)直線與彎管曲面求交存在兩個(gè)交點(diǎn)情況

(b)直線與彎管曲面求交存在多個(gè)交點(diǎn)情況圖9 直線與彎管曲面求交Fig.9 Intersection between the line and the NURBS surface of canal

誤差(mm) 計(jì)算時(shí)間(ms) M1M2M37.114×10-1394.361.7平均時(shí)間65～70

表4 不同直線與彎管曲面求交時(shí)的誤差與計(jì)算時(shí)間比較(迭代容差為10-9 mm)

5 結(jié)論

(1)本文將仿射算術(shù)與區(qū)間運(yùn)算相結(jié)合，一定程度上解決了傳統(tǒng)區(qū)間運(yùn)算的“保守性”，減少了迭代過程中不必要的求交判斷，并放寬了對(duì)初始區(qū)間的要求。

(2)將基于邊曲率或面曲率的子域分解方法應(yīng)用到求交算法中，快速定位預(yù)迭代區(qū)間，在相同的容差要求下，比傳統(tǒng)區(qū)間迭代求交算法更快，提高了迭代求交算法的效率。

(3)通過計(jì)算區(qū)間算子判斷給定直線與參數(shù)空間NURBS曲線或直線與NURBS曲面有無交點(diǎn)和存在交點(diǎn)時(shí)的交點(diǎn)數(shù)目，可保證求解交點(diǎn)精度，且不會(huì)遺漏交點(diǎn)。