李鐵瑞 吳慧 王奇勝 高博青

摘 要：為實現存在裁剪、孔洞的復雜自由曲面建筑網格劃分，提出了一種基于離散的、以均勻性為目標的劃分方法.將復雜曲面離散并縫合，形成由大量面片組成的離散曲面，作為多個參數曲面的一體化表示.采用改進的誤差擴散算法，在離散曲面上按一定的密度進行初始布點.采用基于空間距離的粒子動力松弛算法對點云進行初步均勻化，并應用基于曲面距離的k均值算法進行再次均勻化.對均勻的點云求曲面距離的Voronoi圖，并獲得相應網格.對網格進行拓撲優(yōu)化和光順優(yōu)化.算例表明，本文算法可有效處理存在裁剪、孔洞的復雜自由曲面，并得到均勻光順的三角網格.

關鍵詞：復雜曲面；離散化；網格劃分；均勻化；松弛

中圖分類號：TU311.41 文獻標志碼：A

文章編號：1674—2974（2018）07—0048—06

Abstract： A grid generation method for complicated multiple surfaces with trimmings and holes is presented. This method is based on the discretization and concentrates on the aim of homogeneity. The multiple surfaces are discretized separately and seamed together to achieve a discrete surface. The points are distributed on the discrete surface according to the density applying improved error-diffusion method. The points are homogenized by particle dynamics method with Euclid distance and then homogenized once more by k-means algorithm with surface distance. The Voronoi diagram with surface distance is delivered on the discrete surface to obtain the grids. The topological and smooth relaxations are applied on the grids. Eventually， the case study indicates that this method can solve the problem of grid generation for complicated multiple surfaces effectively and achieve the homogeneous and smooth grids.

Key words： complex surfaces；discretization；grid generation；homogenization；relaxation

隨著計算機輔助設計技術，特別是建模技術的發(fā)展，建筑師的創(chuàng)造力可以得到更好的發(fā)揮，自由曲面的建筑形式因其強力的視覺效果得到了人們的青睞.但是，對復雜的曲面進行網格劃分十分困難，尚沒有行之有效的方法，是目前空間結構領域研究的一個難點與熱點.

對自由曲面，采用顯式或隱式方程都無法表達，常采用NURBS技術[1]，或結合形態(tài)優(yōu)化方法[2]建模.而對于特別復雜的曲面，用單個NURBS曲面表達也甚為困難，乃至無法表達，此時需采用多個曲面來建模.T-Splines[3]技術產生于2003年，并由T-Splines 公司.商業(yè)化，近年來也有大量的研究，該技術對多曲面建模提供了良好的支撐.由于NURBS曲面是由高階非線性參數方程組表示的，一些經典的算法難以應用于NURBS曲面，而對于裁剪曲面特別是多曲面，則更為困難.

Owen[4]對經典網格劃分算法進行了總結.早期的學者提出如波前推進法[5]、Delaunay法[6]、映射法[7]等.但是這些方法在自由曲面應用有其局限性.近年來，不少學者對該問題進行了深入研究，并提出了一些針對NURBS曲面的解決方案.現有的設計多采用參數域設計并映射到曲面，通過一些技術手段改善映射變形，或采用空間距離.熊英等[8]使用空橢圓準則代替?zhèn)鹘y(tǒng)的空圓準則解決曲面Delaunay問題；江存等[9]應用自定義單元法改善映射變形；蘇亮等[10]采用基于主應力線的波前推進法，實現網格劃分；危大結等[11]、潘煒等[12]采用等面積曲面展開方法，改善從平面到曲面的映射關系.對于多重曲面，現有的方法多是在各個曲面上分別劃分網格，再進行調整.然而，此法曲面交界處可能會存在明顯的彎折現象，且由于曲線分屬多個曲面，難以共同優(yōu)化調整，難以達到設計要求.

面對復雜曲面，先將其離散為多面體的離散化表示，是一種解決方案.得到離散化表達的曲面后，裁剪邊界、洞口、多曲面，以及映射關系的負面影響將不再存在.一些基于曲面距離的算法也更容易應用.另一種情形是，建筑師首先制作實物模型，再采用3-D掃描技術數字化，得到離散化表示的曲面數據，就可以不必進行復雜的曲面擬合重構步驟[13]，而直接進行網格劃分.

本文基于曲面離散，提出針對復雜自由曲面的網格劃分方法.經過多個復雜自由曲面的驗證，表明網格劃分算法的有效性、可靠性.

1 算法概覽

本算法的基本流程如下：

步驟1：將曲面離散為小面片組成的離散表示；

步驟2：設置密度控制函數與曲率控制函數；

步驟3：在離散曲面上按密度要求生成初始點；

步驟4：對初始點進行密度與曲率控制的均勻化；

步驟5：形成網格；

步驟6：對網格進行拓撲與光順優(yōu)化.

1.1 離散方法

對參數曲面，如NURBS曲面，若曲面完整，則采用參數域均勻布點，形成網格，映射至曲面.對裁剪曲面，先形成均勻點，再將邊界外的點舍棄，接著在參數平面上求剖分，并映射至曲面.對于多重曲面或環(huán)面，需要在曲面交界處進行縫合.

1.2 改進誤差擴散算法

誤差擴散算法（Error-Diffusion，簡稱E-D算法）由Floyd和Steinberg[14]于1976年發(fā)明.最初被應用于計算機圖像處理領域的灰度圖打印.一張圖片上每一個像素點都有一個灰度信息，而打印時只存在有墨點和無墨點兩種狀態(tài)，如何分布這些墨點使得圖片打印后看起來與原圖一致即是誤差擴散算法所要解決的問題.

離散曲面布點和灰度圖生成情形類似，每一小三角形面片都具有一個點密度值，需要解決的問題是怎樣分布這些點，以滿足密度分布要求.

本文采用密度修正函數fρ與曲率修正函數fk對離散曲面上點密度的分布進行調整.該函數根據用戶的美學需求，可任意定義.

根據前人的實踐經驗，誤差擴散算法容易出現重復模式（artifacts）缺陷.現有的研究采用以蛇形[15]或空間填充曲線[16]替代逐行掃描，以及精心設計誤差分配的范圍與系數[17]，甚至采用變系數法[18]等方法規(guī)避這一缺陷.然而這些研究的對象都為像素矩陣，而本文的應用情形不存在行列概念，故難以采用現成的改進方法.

本文參考前人的思路對E-D算法進行改進.傳統(tǒng)算法中，當某點處密度累積到達閾值t后，即布點.作者在該步驟中加入了概率因子p，布點前需要進行隨機判斷，p概率布點，解決了重復模式問題.概率因子p建議值受布點數k與離散曲面頂點數 n的比值k/n影響，該比值越大p的取值應越大，反之亦然.

傳統(tǒng)算法中閾值t一般取常量0.5，即四舍五入.而根據作者的實驗結果，t取常量0.5將會有起始點附近出現空洞的現象.本文針對此現象進行改進，采用變化的閾值，t將隨著訪問點數的增加而逐漸增大至0.5，即初始時更容易布點，以解決起始點附近的空洞現象.t的增長率取值根據比值k/n確定，該比值越大，則增長率取值應越高；反之亦然.根據作者實驗，修正后的算法，起始點選擇對最終結果影響不大，一般取0號點即可.

圖1為該布點算法的平面算例，該算法可以有效地滿足密度分布，且并未出現重復模式現象，為下一步均勻性優(yōu)化提供了一個可行的初值.

1.3 基于空間距離的均勻化算法

前文生成的點陣雖然在曲面上基本滿足密度分布要求，但是均勻性不良，需要應用松弛算法進行優(yōu)化.本文采用粒子動力松弛算法[19].

本算法基本思想為，將各點視為具有質量mi、電量qi的粒子，質量與電量一般取1.0.任意兩個粒子Vi，Vj間存在斥力：

其中，kq為用于控制斥力強度的系數，該系數對布點結果影響較大，需要試算并調整；ri，j為節(jié)點Vi與Vj的距離，此處采用歐幾里得距離.本文中不考慮引力，距離大于臨界距離rc后，即認為不存在相互作用力，臨界距離rc一般可取為1.5 ～ 2.5倍目標桿長，亦可根據計算結果調整.其中，kq與rc可受用戶自定義函數fρ與fk的調整，以控制不同區(qū)域的點密度.點將在斥力作用下運動并穩(wěn)定在平衡位置附近.

圖2為該松弛算法的平面算例.實驗表明，該算法可以有效地均勻化點陣， 且對初始布點不敏感.

但應用于曲面時，由于動態(tài)更新兩點間的曲面距離復雜度較高，本算法采用歐幾里得距離，存在不足.因此僅作為初步優(yōu)化，為后文步驟提供初值.

1.4 基于曲面距離的均勻化算法

Lloyd算法（又稱Voronoi迭代、Voronoi松弛）由Lloyd[20]于1957年發(fā)明，1982年發(fā)表.Macqueen[21]改進后用于離散對象，也稱k均值算法（k-means）.用于機器學習領域的聚類分析.

對于空間中需要聚類（cluster）的n個d維向量（x1，x2，…，xn），k均值算法的目標是將其分為k個簇 S = {S1，S2，…，Sk}，且簇內各點到其中心的距離之和最小.即：

直接求解優(yōu)化函數（2）較困難.k均值算法采用迭代思想，交替進行以下兩個步驟直至收斂[22]：分配步驟：根據距離函數確定xi屬于哪個簇Si.更新步驟：用各個簇的中心更新其均值.

k均值算法目標與本文均勻分布離散曲面點的目標一致.Du Qiang等[23]采用該算法進行二維點的變密度分布，周炎濤等[24]將其應用于復雜三維點的聚類，取得良好的效果.

本文將該算法應用于離散化曲面.且可受函數fρ與fk的調整.實驗表明該算法可以有效地實現空間曲面上點的變密度分布.

1.5 離散曲面Voronoi圖計算方法

Voronoi圖的基本思想是存在k個站點（site），平面分成k個區(qū)域，使得每個區(qū)域內的點到它所屬區(qū)域站點的距離最近.

本文將曲面離散化后，可以應用Voronoi圖的定義求解.將Voronoi圖定義推廣到離散曲面.離散曲面可視作一個以面片為頂點，相鄰面片之間存在一條邊的圖.當小面片的相對尺度足夠小時，即可認為該圖兩點間的最短路徑即是對應的曲面距離.

求得Voronoi圖后，即可形成網格.

1.6 拓撲優(yōu)化

然而，三角形質量最優(yōu)的網格并非建筑上要求的網格.在1.5節(jié)步驟之后，需要對網格進行拓撲優(yōu)化.Frey等[25]認為，對于三角網格，內部節(jié)點度數為6是較優(yōu)的.基本思想為遍歷每一個四邊形，若交換該四邊形的對角線能使內部頂點度數更接近6，則進行交換.

如圖3所示，左圖圈中四邊形箭頭所指對角線對應著質量更高的兩個三角形，但四個頂點皆為奇異點，其度數分別為5，7，5，7.進行邊交換后，得到右圖，對應的四個頂點度數都為6，規(guī)整性提高.

1.7 光順優(yōu)化

本文采用彈簧質點法[26].其基本思想是將桿件視為具有剛度ks的彈簧，節(jié)點視為具有質量的質點.其中， m一般取常量1.0即可；ks的大小影響計算速度與收斂性，需要根據經驗調整；計算過程中通過調整各桿件的原長控制彈簧力的大小與方向，并影響最終結果.原長可受密度修正函數fρ與曲率修正函數fk的修正，以控制不同區(qū)域的密度.節(jié)點將在彈簧力作用下運動，并穩(wěn)定在平衡位置.實驗表明，最終可得到較為光順的網格.

2 算例分析

2.1 完整單曲面算例

作者創(chuàng)建了一個月牙形曲面模型，如圖4（a）所示.長約100 m，跨度約40 m，矢高約20 m.該模型建模時采用一個完整曲面，無裁剪及孔洞，作為基準算例.將其離散為320 000個小三角形面片.點布置以及Voronoi圖如圖4（c）所示，網格效果如圖4（d）所示.本算法生成的網格均勻性較好，但由于曲面變化較大，為保證均勻，網格流線存在發(fā)散和收縮現象，存在少量奇異節(jié)點，但過渡較為流暢.

作為對照的映射法網格如圖4（b）所示.由于該模型為單曲面，且無裁剪及孔洞，適合應用映射法，由于曲面形態(tài)、建模方式等原因，映射變形較大，生成的網格流暢性較好，但是存在明顯均勻性問題.

均勻性指標如表1所示，桿長相近的情況下，本算法網格均勻性顯著優(yōu)于映射法.

2.2 裁剪孔洞曲面算例

作者對2.1節(jié)曲面進行了裁剪，如圖5（a）所示.模型尺寸同2.1節(jié)，洞口直徑為12 m，遠大于網格尺寸，因此將會影響網格生成.

算例表明，本算法可有效處理裁剪和孔洞，優(yōu)化后的點布置及Voronoi圖如圖5（c）所示，生成的網格如圖5（d）所示.洞口附近通過少量奇異點過渡，無明顯不流暢現象.洞口網格質量較好，整體流暢性尚可.均勻性指標如表2所示，桿長相近的情況下，本算法網格均勻性顯著優(yōu)于映射法.

作為對照的映射法網格如圖5（b）所示.從圖中，可以發(fā)現，映射法受孔洞影響較大，孔洞和裁剪邊界周圍網格難以處理，出現明顯的不流暢、不均勻現象.且由于過于不均勻，1.7節(jié)的優(yōu)化算法難以應用，優(yōu)化過程中易出現網格折疊現象.

2.3 復雜曲面算例

復雜的建筑模型，采用多重曲面建?？梢院喕５碾y度，提高建模質量，提高建模效率[3].

然而對于復雜曲面，網格劃分困難.現有的方法多是在各個曲面上分別劃分網格，再對多個網格進行縫合，最后進行調整和優(yōu)化.由于一條曲線分屬兩個曲面，難以共同調整優(yōu)化.流暢性與均勻性都難以保證，曲面交界處往往存在明顯彎折.

圖6是一個船形曲面的局部.由圖6（a）可見，由于模型形式復雜，建筑師建模時采用了兩個曲面組合來表達.如圖6（a）的曲面結構線所示，上部曲面是裁剪曲面，且一端存在急劇收縮的現象，這都給傳統(tǒng)網格劃分方法造成了很大困難.

對于該算例，采用1.1節(jié)中所述的方法，先分別將兩個曲面離散為約150 000、210 000個小三角形面片，再對兩個網格進行縫合，最終形成共360k個面片.

實驗證明本算法可以處理裁剪曲面、多重曲面等形式復雜的模型.如圖6（b）所示，生成的網格均勻性較好，桿長均值為2.313 m，桿長方差僅為0.044 m2.但由于模型變化較大，為保證上下網格密度相同，則網格流線不可避免存在發(fā)散和收縮現象，個別內部節(jié)點所連桿件數量不為6，優(yōu)化后的網格整體流暢性尚可.

作者對上下兩個曲面分別采用映射法劃分網格，如圖6（c）所示，下部曲面是完整曲面，映射法效果較好；而上部曲面由于是裁剪曲面，且一端急劇收縮，映射法僅保證了拓撲流暢，實際空間效果較差.由圖6（c）可見，映射法存在嚴重的映射變形問題.對上下兩曲面分別應用1.7節(jié)的均勻化后，如圖6（d）所示，上部曲面的映射變形問題得到一定改善.但整體仍存在下部桿件過于密集，上部桿件過長的問題.而且網格流線在兩曲面分界線處存在明顯的彎折現象，而流線分屬兩個曲面，難以調整.而圖6（b）所示的本算法網格并沒有受到曲面分界線的影響.

作者根據標高設置了分段密度控制函數，上部疏下部密，并得到了圖6（e）和（f）的最終結果，網格流暢，且滿足密度分布，僅存在少量奇異點.實驗表明本算法可以根據用戶的美學需求進行變密度網格生成.

3 結 論

針對復雜曲面，本文提出了一種以密度控制為目標的網格劃分方法.首先，對復雜曲面離散并縫合，一體化為由大量三角面片組成的離散曲面，可消除裁剪邊界、孔洞、多曲面拼接、閉合等負面影響.其次，采用改進的誤差擴散算法進行初始點的分布，該算法布點結果雖不夠均勻，但滿足密度分布，為之后的步驟提供了可靠的初值.然后，采用基于空間距離的粒子動力松弛算法，對點云進行初步均勻化，應用基于曲面距離的k均值算法對點云進行再次均勻化，得到曲面上滿足密度要求的均勻點陣.接著求離散曲面Voronoi圖，由此得到曲面Delaunay圖.之后，進行拓撲優(yōu)化，提高規(guī)整性.最后，采用彈簧質點松弛算法對該三角網格進行光順優(yōu)化，得到均勻流暢的最終網格.

算例表明，本文算法可以處理裁剪、開洞、變化劇烈、環(huán)面、多曲面拼接等復雜情況，適應性好.最終形成的網格能滿足密度要求，但網格流線可能存在發(fā)散和收縮情況，存在內部節(jié)點連有5根或7根桿件的現象.

奇異節(jié)點不可避免，本文算法的缺陷在于，奇異節(jié)點的位置并不能控制，本文采用的拓撲優(yōu)化算法并無全局優(yōu)化能力，可能需要在形成網格后交互地修改拓撲，消除奇異點的負面影響.對此仍需要進行更多的研究探索，在形成網格前就交互確定奇異點，或在算法中全局地優(yōu)化奇異點的數量與位置，以獲得更規(guī)整的網格.

參考文獻

[1] PIEGLL L，TILLER W. The NURBS Book[M] .2nd ed. Berlin Heidelberg： Springer-Verlag，1997.

[2] 王磊，楊彬，張其林. 非均勻有理B樣條曲面形狀優(yōu)化方法[J]. 湖南大學學報（自然科學版）， 2012，39（7）：14—19.

WANG L， YANG B，ZHANG Q L. Shape optimization of non-uniform rational B-spline surface[J].Journal of Hunan University（Natural Sciences），2012，39（7）：14—19. （In Chinese）

[3] SEDERBERG T W， ZHENG J， BAKENOV A， et al. T-splines and T-NURCCs[J]. ACM Transactions on Graphics，2003，22（3）： 477 —484.

[4] OWEN S. A survey of unstructured mesh generation[C]// International Meshing Roundtable. Dearborn， USA： IMR， 1998： 239—267.

[5] LHNER R. Progress in grid generation via the advancing front technique[J]. Engineering with Computers，1996，12（3）：186—210.

[6] BOROUCHAKI H，HECHT F，SALTEL E，et al.Reasonably efficient Delaunay based mesh generator in 3 dimensions[C]// International Meshing Roundtable. South Lake Tahoe， USA： IMR， 1999： 3—14.

[7] COOK W A，OAKES W R. Mapping method for generating three-dimensional meshes： past and present[C]// International Computer Engineering Conference.San Diego，USA：ASME，1982.

[8] 熊英，胡于進，趙建軍. 基于映射法和Delaunay方法的曲面三角網格劃分算法[J]. 計算機輔助設計與圖形學學報， 2002， 14（1）： 56—60.

XIONG Y， HU Y J， ZHAO J J. An algorithm of surface triangulation based on mapping and Delaunay method[J]. Journal of Computer-Aided Design & Computer Graphics，2002，14（1）： 56—60.（In Chinese）

[9] 江存，高博青. 基于自定義單元法的自由曲面網格劃分研究[J]. 建筑結構， 2015， 45（5）： 44—48.

JIANG C，GAO B Q. Research on free-form surface meshing based on self-defined element method[J]. Building Structure， 2015， 45（5）： 44—48. （In Chinese）

[10] SU L， ZHU S， XIAO N， et al. An automatic grid generation approach over free-form surface for architectural design[J]. Journal of Central South University， 2014， 21（6）： 2444—2453.

[11] 危大結，舒贛平. 自由曲面網格的劃分與優(yōu)化方法[J]. 建筑結構， 2013， 43（19）： 48—53.

WEI D J， SHU G P. Mesh generation and optimization method for free-form surface grid[J]. Building Structure， 2013， 43（19）： 48—53. （In Chinese）

[12] 潘煒，吳慧，李鐵瑞，等.基于曲面展開的自由曲面網格劃分[J]. 浙江大學學報（工學版）， 2016， 50（10）： 1973—1979.

PAN W， WU H， LI T R，et al. Grid generation on free-form surface based on surface flattening[J]. Journal of Zhejiang University（Engineering Science），2016，50（10）：1973—1979.（In Chinese）

[13] SCHALL O， SAMOZINO M， FALCIDIENO B， et al. Surface from scattered points： a brief survey of recent developments[C]// Proccedings of the first international workshop towards Semantic Virtual Environments. Villars， Switzerland： SVE， 2005： 138—147.

[14] FLOYD R W， STEINBERG L. Adaptive algorithm for spatial greyscale[C]// Proceeding SID. New York， USA： SID， 1976：75—77.

[15] ULICHNEY R. Digital halftoning[M].Cambrige，MA：MIT Press， 1987：12—13.

[16] VELHO L， GOMES J D M. Digital halftoning with space filling curves[J]. ACM SIGGRAPH Computer Graphics， 1991， 25（4）：81—90.

[17] FAN Z. Set of easily implementable coefficients in error diffusion with reduced worm artifacts[C]// Proceeding SPIE. San Jose， USA： SPIE，1996：222—225.

[18] OSTROMOUKHOV V. A simple and efficient error-diffusion algorithm[C]// Proceeding SIGGRAPH. Los Angeles， USA： ACM， 2001： 567—572.

[19] BRONSON J R， LEVINE J A， WHITAKER R T. Particle systems for adaptive， isotropic meshing of CAD models[J]. Engineering with Computers， 2012， 28（4）： 331—344.

[20] LLOYD S. Least squares quantization in PCM[J]. IEEE Transactions on Information Theory， 1982， 28（2）： 129—137.

[21] MACQUEEN J. Some METHODS for classification and analysis of multi variate observations[C]//Proceedings of Berkeley Symposium

on Mathematical Statistics and Probability. Oakland， USA： University of California Press， 1967：281—297.

[22] MACKAY D.Information theory，inference，and learning algorithms[M]. Cambridge， UK： Cambridge University Press，2003： 640.

[23] DU Q， FABER V， GUNZBURGER M. Centroidal Voronoi tessellations： applications and algorithms[J]. Siam Review， 1999， 41（4）： 637—676.

[24] 周炎濤，吳正國，易興東. 基于網格帶有參考參數的擴展聚類算法[J]. 湖南大學學報（自然科學版），2009，36（2）：48—52.

ZHOU Y T， WU Z G， YI X D. Extended grid-based clustering algorithm with referential parameters[J]. Journal of Hunan University（Natural Sciences），2009，36（2）：48—52.（In Chinese）

[25] FREY W H， FIELD D A. Mesh relaxation： a new technique for improving triangulations[J]. International Journal for Numerical Methods in Engineering， 1991， 31（6）： 1121—1133.

[26] 李基拓，陸國棟. 基于邊折疊和質點-彈簧模型的網格簡化優(yōu)化算法[J].計算機輔助設計與圖形學學報，2006，18（3）：426—432.

LI J T， LU G D. Mesh simplification and optimization with edge collapse and mass-spring model[J]. Journal of Computer-Aided Design & Computer Graphics，2006，18（3）：426—432.（In Chinese）