劉凱正， 劉利剛

(中國科學技術大學數(shù)學學院，安徽 合肥 230026)

在計算幾何、計算機輔助設計及計算機圖形學等領域中，布爾運算是一種構造實體的常用方法。網格的布爾運算在19世紀80年代被提出，后續(xù)的研究主要集中在提升算法的效率與準確性方面。而本文提出的布爾運算算法,在兼顧效率與精度的同時，拓寬了布爾運算的適用范圍。

在布爾運算發(fā)展的早期，由于二叉空間分割(binary space partition，BSP)樹內在的空間劃分性質，有很多基于 BSP樹的算法被提出。但這類算法有2個明顯的缺點：一是有很多無必要的分割，如圖1所示，當圖1(a)的多邊形轉為BSP樹時，邊d被分為d1和d2，若未有其他邊與邊d相交，就沒有必要做分割；二是不能對多個模型同時進行運算。基于 BSP樹的布爾運算算法一般是先將其中一個模型轉為BSP樹，再將另一個模型轉為BSP樹并將其轉為多邊形結構，然后插入最初的 BSP樹中，并根據(jù)布爾運算的類型提取多邊形。所以該類布爾運算一次只能對2個模型進行運算，如果有多個輸入，將會導致嚴重的性能問題。由于這2個問題均源于 BSP樹本身，難于做改進，所以現(xiàn)在已鮮有人基于BSP樹來研究布爾運算。

圖1 將多邊形轉為BSP樹

網格布爾運算是求交運算，所以最近很多算法都未借助復雜的數(shù)據(jù)結構，而是直接求出網格的相交元素。其難點在于正確歸類揉合所有輸入網格所劃分的空間區(qū)域。所以現(xiàn)有的大量布爾運算算法對輸入有嚴格要求，如不能有空洞、不能是非流形、不能有退化的三角形及不能有自交等。近些年隨著布爾運算的發(fā)展及完善，這些限制條件已逐步被移除。但現(xiàn)有的大多數(shù)算法均要求輸入是能作為實體模型邊界表示的封閉網格，而在實際應用中，輸入的大多為有邊界邊的開網格，并希望能對其進行運算。如圖2所示，2個開圓環(huán)面A和B按圖2(c)擺放，并進行布爾并運算，希望能得到如圖2(d)所示的封閉圓環(huán)面。為此本文提出了一種能適用于開網格的布爾運算算法。

基于 BSP樹的算法對于網格每一個面片都是根據(jù)其法向分出內外，這種只做局部考慮的特性允許輸入是開網格。但正是因為沒有從整體上考慮網格所劃分的空間區(qū)域，所以會得到與期望不符的結果。如圖3所示，輸入2個拋物面，對其分別進行布爾交與布爾并運算，基于 BSP樹的布爾運算不僅不能保證輸出的是封閉網格(圖 3(d))，而且從圖3(c)和(d)可以看出不符合A∩B?A∪B的邏輯。而本文算法則不會引起類似的錯誤。

圖2 2個非封閉圓環(huán)面做并運算

圖3 本文算法與基于BSP樹算法的比較

布爾運算中有大量的求交計算，普通的浮點運算極易導致數(shù)值上的錯誤，從而可能引起輸出網格頂點位置錯誤及拓撲混亂等一系列問題。為此，很多關于布爾運算的文獻用大量篇幅去闡述如何避免數(shù)值錯誤的方法。但實際上，求交運算只需要用到加減乘除4種運算，且于有理數(shù)域上封閉。此外，常見的保存網格的OBJ文件及OFF文件均用3個有限位的小數(shù)去表示頂點在R3中的位置，所以其屬于Q3。CGAL的幾何庫實現(xiàn)了一種運算核EPECK，且能在有理數(shù)域上進行無損計算(如果所有參與的運算在有理數(shù)域上封閉，如加減乘除滿足條件，而sin及sqrt等則不滿足)。而本文算法的實現(xiàn)正是采用了此運算核，因此不會引起任何數(shù)值上的問題。

1 相關工作

自文獻[1-3]提出關于網格的布爾運算后，有大量工作對其進行了完善及發(fā)展。如文獻[4]提出基于GPU柵格化的多邊形布爾運算，該方法能快速計算，但在進行柵格化時會有精度丟失。文獻[5-6]提出了基于 BSP樹的布爾運算，后續(xù)有很多研究者對此做了大量完善及擴充。文獻[7]基于BSP樹，用平面作為整個算法構建的基本元素，將點用3個平面相交來表示，這樣不用顯式求出交點，從而避免了浮點運算可能引起的錯誤。而文獻[8]在文獻[7]的基礎上，提出用八叉樹將空間分成多個區(qū)域；然后只在有相交產生的區(qū)域局部建立BSP樹，較文獻[7]算法的效率有顯著提高。但基于 BSP的算法均有一個明顯的缺點：在建 BSP樹時會將多邊形轉為平面，而平面的無界特性會使得在劃分空間區(qū)域的過程中做多余的分割，過度加細。

CGAL庫是一個眾所周知的開源幾何庫[9]，其布爾運算在魯棒性方面仍然代表最高水準。并實現(xiàn)了一種運算核：EPECK，通過該運算核，能在有理數(shù)域上做無損計算。為了保證算法的魯棒性，本文采用了該運算核。但 CGAL中的布爾運算需要將輸入轉為Nef多邊形[10-11]，這是一種復雜的數(shù)據(jù)結構，運行時會占用大量內存，在效率方面也不盡人意。此外，對輸入也有苛刻的要求，例如不能有空洞，只能具有一個連通分支等。所以 CGAL庫中的布爾運算雖然具有最高水平的精度，但性能表現(xiàn)不佳。

文獻[12]提出將網格沿著交線分成不同的流形塊，最后輸出的網格是流形塊拼裝的結果。而文獻[13]在文獻[12]的基礎上提出胞體的概念，并通過計算胞體的環(huán)繞數(shù)來判定胞體相對于輸入網格的內外關系。文獻[13]中的算法已集成在 libigl庫中[14]，但其只適用于封閉網格，無法計算開網格所形成的胞體的環(huán)繞數(shù)。本文通過添加虛擬流形塊用于封閉胞體，進而正確判斷內外關系。該算法適用于更多類型的網格，網格可以是封閉也可以是非封閉的，可以有多個連通分支，可以有空洞。而布爾運算是一種構造實體的方法，所以本文算法可以保證最后輸出的一定是實體或者空集。

2 流形及環(huán)繞數(shù)

首先介紹一些有關拓撲流形網格的概念[15]。對于網格的任意頂點v，可用N(v)表示網格中共享頂點v的所有三角形形成的鄰域。如果N(v)中所有三角形按序排列且使得序列中所有 2個相鄰三角形都共享一條以頂點v為端點的邊，則將頂點v叫做流形頂點。如果網格中的一條邊e恰好被2個三角形共享或者只屬于某個單個三角形，則將邊e叫做流形邊。如果邊e只屬于一個三角形，將邊e叫做邊界邊。如果網格中所有的頂點和各邊都是流形的，則網格也是流形的。

對于任意一個網格，若其中一個子網格恰好為一個流形，則可將其稱為流形塊(patch)[13]，將流形塊圍成的封閉區(qū)域叫胞體(cell)。

為了便于觀看和闡述，圖 4給出了二維的情形，讀者可類比到三維情形。圖4(a)中不同顏色的曲線段表示流形塊，圖4(b)中不同顏色的區(qū)域表示胞體。

圖4 流形塊和胞體

對于二維空間，定義點p關于閉的利普希茨曲線C的環(huán)繞數(shù)(winding number)為當沿逆時針走完曲線C時繞點p的圈數(shù)[16]。不失一般性，可令點p=0。此時環(huán)繞數(shù)w(p)為

如圖5所示，點P0的環(huán)繞數(shù)為-1，點P1的環(huán)繞數(shù)為1，而點P2環(huán)繞數(shù)為0。

圖5 環(huán)繞數(shù)

同樣地，如果點p∈R3，S為閉的利普希茨曲面，則令點p關于S的環(huán)繞數(shù)為

點p關于S的環(huán)繞數(shù)是一種有符號的整數(shù)，反映了S包圍點p的層數(shù)。從離散的觀點看，如果M={fi}為封閉網格，其中fi為三角形；令Ωi(p)為三角形fi關于點p的有向空間角度[16]，則

如果三角形f的3個定點分別為：v0，v1，v2，令a=v0-p，b=v1-p，c=v2-p，a=||a||，b=||b||，c=||c||，則[16]

其后的布爾運算算法會用一種巧妙的方式求胞體的環(huán)繞數(shù)，不會用到復雜的式(4)。

3 本文算法

對于輸入的n個網格 {Mi}i∈{1,…,n}(為方便闡述，令n=2，即輸入為A和B2個網格，為更具一般性，以開網格為例，如圖 6所示)，須滿足如下條件：在經過步驟1解決自交問題后，檢測每條非流形邊ab，令count表示與ab相關的一個整數(shù)，初始時為0。然后找到所有共享ab邊的三角形，若ab邊在此三角形中的順序恰好為“ab”，則count加 1，若順序為“ba”，則count減 1。若所有非流形邊所對應的count均為0，則是有效輸入。

本文算法流程如下：

輸入：網格A和B。

輸出：布爾運算得到的網格。

步驟1.合并A和B，并解決自交問題；

步驟2.劃分流形塊并刪除含邊界的流形塊；

步驟3.標記胞體；

步驟4.建立虛擬流形塊；

步驟5.計算胞體的環(huán)繞數(shù)；

步驟6.根據(jù)環(huán)繞數(shù)輸出網格。

圖6 輸入網格

對于任意網格M，可用V(M)={v1,v2,···,vr}表示M的r個頂點，用F(M)={f1,f2,···,fs}表示M的s個面。

3.1 標記流形塊及胞體

步驟1.將A和B合并到一張網格M0上。此過程只是純粹地將V(A)和V(B)以及F(A)和F(B)進行集合并操作。此時M0可能會有自交。M0中每個三角形對應一個平面，P為所有平面形成的集合。那么，對于任意平面p∈P，可找出所有與此平面相交的元素，這些元素可能是三角形、線段或點。將所有這些元素作為限制條件，在平面P上進行二維的帶約束的 Delaunay三角化(constrained Delaunay triangulation，CDT)。記C(M0)為CDT之后的所有三角形形成的集合。對于三角形f∈C(M0)，如果f被F(M0)中某個三角形所包含，則保留，否則舍棄。將所有保留的三角形所形成的網格記作M1。M1與M0相比，V(M0)?V(M1)，且任意f∈F(M1)一定被某個f′∈F(M0)所包含。所以M1是對M0的加細，M1不包含任何相交的三角對，并且沒有丟失或改變M0的任何形狀信息。如圖7所示。

圖7 解決自交問題

步驟 2.為了確認M1中的每個三角形所在的流形塊，可采用寬度優(yōu)先搜索算法，通過遍歷每個三角形f∈F(M1)，然后檢查f的每條邊e，如果e為流形邊，則將包含e的另一個三角形f′與f歸在同一個流形塊。如圖8所示，用不同顏色來區(qū)分流形塊。

圖8 劃分流形塊

本文算法允許輸入是開網格，而開網格并不具有決定一個唯一實體的完備信息。修補開網格的方法，即使其成為一個封閉的網格。采取不同的修補算法可得到不同的被修補的網格，從而影響最終布爾運算的結果。本文算法則是通過所有輸入去決定開網格所代表的實體模型。其并不是一成不變的，而是隨著輸入的變化而變化。如圖9所示，一個拋物面與不同的網格做布爾運算，其潛在所表達的實體模型是不同的。至于具體用哪些流形塊來修補,則是通過步驟 5建立虛擬流形塊來計算胞體環(huán)繞數(shù)決定的，而這一步僅僅是將包含邊界邊的流形塊舍去，如圖10所示。

圖9 根據(jù)所有輸入決定開網格所代表的實體

圖10 去掉包含邊界邊的流形塊

步驟3.為了標記胞體，首先遍歷網格M2中的每條邊，如果這條邊為非流形邊(共享此條邊的三角形個數(shù)大于2)，此時共享此邊的每個三角形必定屬于不同的流形塊。這樣，就可知流形塊之間的相鄰關系。三角形3個頂點的排列順序指定了其法向的方向，因為流形塊中所有三角形的法向是一致的，所以可將其中任意一個三角形的方向記作此流形塊的正方向。每個流形塊p與2個胞體相鄰，將p正方向對應的胞體記作C+(p)，負方向記作C-(p)。以每條非流形邊e為軸心，按逆時針方向排列共享該條邊的三角形，得到循環(huán)序列Sf(e)。而每個三角形fi∈Sf(e)對應一個流形塊pi，從而得到一個對應的流形塊序列Sp(e)={···,pi,pi+1,···}。對于任意一對相鄰的流形塊pi，pi+1∈Sp(e)，令胞體pi與pi+1等同。如圖11(a)箭頭所示2個流形塊，左邊記為p0，右邊記為p1，那么按照先前的敘述則應C_(p0)=C_(p1)。如果M2是連通的，則此時所有胞體已被標記。

如果M2有多個連通分支，則對于M2的任意一個連通分支T，可將所有與T相鄰的胞體中包含無窮遠處的胞體記作該連通分支的外胞體O(T)。對于任意異于T的連通分支T′，找到T′中到T距離最短的點v(T,T′)。如果v(T,T′)不位于T′的外胞體且位于T的外胞體，則可將T′加入到候選集合A(T)中，找到A(T)中距離T最短的分支T0，同時將T0中包含v(T,T0)的胞體與O(T)等價。最后再將所有未做等價的外胞體等價，并標注所有胞體，如圖11所示。

圖11 標記胞體

3.2 建立虛擬流形塊

步驟 4.因為輸入的網格可能是非封閉的，文獻[13]中計算環(huán)繞數(shù)的方法并不適用。如圖 12所示，當從外胞體穿過流形塊p進入胞體C時，因為p∈A卻?B，所以按照文獻[13]胞體C關于B的環(huán)繞數(shù)應與外胞體的環(huán)繞數(shù)相同，而外胞體的環(huán)繞數(shù)為0，但胞體C的關于B的環(huán)繞數(shù)應為1?；谏鲜鲈?，本文建立虛擬流形塊來封閉胞體的具體規(guī)則如下：

圖12 穿過p計算C的環(huán)繞數(shù)

對于任意一個胞體C，用P(C)表示所有與C相鄰的流形塊的集合。對于任意流形塊p∈P(C)，用(p)表示p所在的初始輸入網格，并令(C) ={(p)|p∈P(C)}。那么對于任意M0∈(C)且M0≠(p)，總存在流形塊p0使得(p0) =M0。如果p0的負方向剛好指向C，則可新建一個虛擬的流形塊p′=p。如果p的正方向指向C，則翻轉p′使得p′的負方向指向C。最后，令p′所在的輸入網格為M0，即(p′)=M0。如圖13所示。建立虛擬流形塊的算法可用如下偽代碼表示：

圖13 建立虛擬的流形塊

ForpinP(C):

Ifp0的負方向指向C:

p′=p

Ifp的正方向指向C:

翻轉p′

3.3 計算環(huán)繞數(shù)

步驟5.首先，將連接無窮遠處的胞體關于任何網格的環(huán)繞數(shù)記為0。假定胞體C的環(huán)繞數(shù)W(C)已知，用P(C,C′)表示從胞體C到胞體C′需要穿過的流形塊的集合。對于任意p∈P(C,C′)，如果穿過p且從正方向到負方向，則W(C′)加上1，否則減去1。

如圖14所示，假定已知胞體C0關于網格A的環(huán)繞數(shù)WA(C0)=0，關于B的環(huán)繞數(shù)為WB(C0)=1。需要計算C1關于A和B的環(huán)繞數(shù)WA(C1)和WB(C1)。從C0到C1需要穿過3個流形塊，其中一個屬于網格A(紅色，記為p0)，2個屬于網格B(綠色，記為p1和p2)。穿過p0為從正方向到負方向，所以WA(C1) =WA(C0)+1=1。而穿過p1為從負方向到正方向，穿過p2為從正方向到負方向，所以WB(C1)=WB(C0)+(- 1 )+ 1 =1。

圖14 計算環(huán)繞數(shù)

3.4 輸出結果

步驟 6.計算完所有胞體的環(huán)繞數(shù)后，可以根據(jù)得到的環(huán)繞數(shù)及布爾運算的類型決定保留、翻轉或舍棄f∈M2。最后保留的三角形形成最終的輸出網格。

當一個胞體C關于輸入網格M的環(huán)繞數(shù)大于0，則表示胞體C在M的內部，否則在外部。當進行A-B布爾運算時，希望保留那些在A內部但在B外部的胞體。所以對每個胞體有如下保留規(guī)則：

當進行布爾并時：

當進行布爾交時

當進行布爾減時

每一個流形塊p與2個胞體相鄰，當2個胞體同時保留或舍棄時，則流形塊p不能作為輸出實體的邊界，所以應舍棄。因此，當一個流形塊若被保留，則其相鄰的2個胞體必須有互斥的保留屬性。流形塊保留的具體規(guī)則為舍棄(s,s)、翻轉(s,n)、保留(n,s)、舍棄(n,n)，其中s為保留胞體，n為舍棄胞體，括號內第一個值為p正方向胞體是否保留，第二個值為p負方向胞體是否保留。

4 實驗結果

實驗在Ubuntu19.04下進行的。實驗所采用的CPU是Intel i5 6500，算法由C++實現(xiàn)，使用的C++編譯器為g++8.3。

對于有邊界的非封閉網格，也能像實體一樣進行布爾運算。圖15為對2個有邊界(已用綠色標出)的網格進行交并差后的結果。

圖15 非封閉網格的布爾運算

通過布爾運算可做一些有趣的事情，如圖 16所示，讓塑像A和拋物面B做布爾并后，就像是給塑像加上了一頂帽子。

圖17同Campen的算法[8]，將不同分辨率的塑像網格旋轉90°并與原網格做布爾并。本文算法能得到拓撲良好的輸出，且未產生任何非流形邊。由于文獻[8]未能公布源碼，所以只能在不同的實驗配置下進行效率上的比較。其所用的 CPU為 i7 2.67 GHz，本文所用的CPU為i5 3.2 GHz。本文選取了3種不同分辨率的網格，三角面片的個數(shù)分別為200 K，400 K及800 K，旋轉90°后與自身做布爾并，最終各自的算法的所需時間見表1。從表1中可以看出，本文算法相較于Campen的算法，節(jié)省了約一半時間。

圖16 給塑像加上帽子

圖17 2個塑像的并

表1 本文算法與Campen算法的比較(s)

為了進一步測試本文算法的效率性，與3個著名的布爾運算開源庫CGAL，Cork[17]及CARVE[18]做比較(表2)。為此，構造了面數(shù)分別為4 K，24 K，240 K及950 K的球面做布爾并運算。QuickCSG能對多個輸入快速進行布爾運算[19]，由于只是發(fā)布了一個預編譯好的二進制文件，并未公布源碼，在使用該文件進行測試時，發(fā)現(xiàn)極易崩潰，所以此次未能將其列入到比較中。從表2中可以看出CGAL庫在效率上存在嚴重不足，當網格面數(shù)過多時，程序會直接崩潰；Cork庫在面數(shù)較少時，性能最佳，當面數(shù)較多時，效率上有明顯下降；而CARVE庫則始終能高效計算；本文算法雖然在面數(shù)較少時，效率上遜于 CAEVE及 Cork，而當面數(shù)較多時，幾乎與CARVE一樣高效。雖然從實驗數(shù)據(jù)上看，本文算法在效率上與CARVE相比并無優(yōu)勢，但有一個不可忽略的事實——本文算法采用了CGAL的EPECK運算核。而根據(jù)文獻[20]表明，這將會導致比浮點運算慢 4～5倍。所以實際上，本文算法在架構上是優(yōu)于 3種算法的。另外，只有本文算法允許輸入開網格，而其他算法只能被應用到封閉網格。

表2 本文算法與3個開源庫進行比較(s)

為了進一步測試本文算法的魯棒性，構造了一張復雜的球面。如圖18(a)所示，其不僅非封閉(邊界已用綠色標出)，而且有很多褶皺和溝壑；在圖18(b)中畫出了特征邊，且非常不光滑。3DSMAX是一款專業(yè)的建模軟件，所集成的布爾運算允許輸入是開網格。采用3DS MAX對2張球面做布爾并時，其無法處理如此復雜的情形，拋出了運行錯誤：Invalid Boolean。而本文算法能得出正確的布爾并的結果，如圖 18(c)所示。本文已將代碼公布在 GitHub上(https://github.com/liukaizheng/GeneralBoolean.git)。

5 結 束 語

本文提出一種基于三角網格的布爾運算方法。通過建立虛擬流形塊，可準確計算每個胞體的環(huán)繞數(shù)，進而合理判斷一般網格的內外關系。因此，相比許多其他算法，本文算法能應用到更多類型的網格上，允許輸入是開網格、允許有空洞、多個連通分支及自交等。此外，采用CGAL庫中的EPECK運算核，能無損計算布爾運算所產生的相交元素，保證了算法的魯棒性。經過多次實驗表明，本文算法是一種魯棒的且適用廣的布爾運算算法。

圖18 非封閉球面的布爾并

圖19 包含相交元素的三角形

但本文算法也有一些不足之處。在算法的實現(xiàn)過程中，對于步驟1中三角化前的輸入未做預處理。如圖19所示，對這樣一個包含相交元素的三角形進行三角化時，本文輸入的是邊AB，BC，CA，DE及AF。沒有將AB分成AE及AC，也沒有預先算出AF及DE的交點，而是交由 CGAL的CDT函數(shù)去做的。所以三角化函數(shù)不僅需要計算交點，而且還會生成多余的退化三角形AEC及三角形BFC等，隨后還需將這些三角形過濾掉。如果進行預處理，會有效提升效率，所以在將來的工作，將考慮解決該問題。另外，在每個平面上進行三角化的工作是相互獨立的，這為并行計算提供了可能，所以，后續(xù)也會考慮用GPU進行加速。