池寶濤 張見明 鞠傳明

1.湖南大學機械與運載工程學院，長沙，4100822.湖南大學汽車車身先進設(shè)計制造國家重點實驗室，長沙，410082

0 引言

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

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

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

1 算法基本理論

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

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

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

(1)

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

(2)

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

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

(3)

(3)NURBS曲面。一張在u方向p次和v方向q次的NURBS曲面或T-Spline曲面定義為

(4)

0≤u,v≤1

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

1.2 區(qū)間分析

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

1.2.1Newton算子

Newton算子定義如下：

(5)

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

假定f:Rn→Rn，F(xiàn)為f的區(qū)間擴展[9]，F(xiàn)′為F的一階導數(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中必有解存在，當且僅當W(N(X))?W(X)時，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中必有解存在，當且僅當W(K(X))?W(X)，f(x)在X中存在唯一解。

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

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

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

1.3 仿射算術(shù)

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

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

給定兩個仿射形式進行運算：

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

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

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

(9)

和

(10)

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

(1) 當u1≠u0且v1≠v0時

(11)

(2)當u1≠u0且v1=v0時

(12)

(3)當u1=u0且v1≠v0時

(13)

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

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

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

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

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

算法具體步驟描述如下：

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

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

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

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

(14)

0≤u,v≤1

圖3 絕對坐標系下的直線與NURBS曲面Fig.3 Line and NURBS surface in the absolute coordinate system

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

(15)

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

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

算法具體步驟描述如下：

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

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

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

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

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ù)計算G(X(k))=Y(k)F′(X(k)),將仿射矩陣G(X(k))轉(zhuǎn)化為新的區(qū)間矩陣H(X(k))。代入Krawczyk算子公式，應(yīng)用區(qū)間算術(shù)計算

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

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

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

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

4 算法實例驗證

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

(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) 計算時間(ms)M1M211.208×10-1427.517.126.758×10-1448.226.437.933×10-1452.628.3

(a)直線與NURBS曲面求交存在一個交點情況

(b)直線與NURBS曲面求交存在兩個交點情況

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

直線誤差(mm) 計算時間(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)直線與彎管曲面求交存在兩個交點情況

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

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

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

5 結(jié)論

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

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

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