基于區域劃分和三角網格優化的基段選擇算法

吳祿慎，王偉杰，陳華偉，馮 偉

（南昌大學 機電工程學院，江西 南昌330031）

0 引 言

構建四邊形網格的所有算法可以歸類為直接法［1］和間接法［2］。直接法不經過中間環節，把點云直接生成四邊形網格。由于需要處理的點云數據量巨大、模型形狀和邊界復雜，研究人員很難直接獲得分布均勻、拓撲成四邊形的網格模型，因此難以達到B 樣條曲面重構四邊形參數域的要求。間接法先使用成熟Delaunay三角剖化［3］算法構建散亂點云的三角網格，建立空間網格拓撲關系，然后把三角網格轉化為四邊形網格。與直接法相比，間接法程序更易實現，并且能生成質量較好的四邊形網格。

在許多工程分析條件下，沿區域邊界單元尤其重要，理想情況是網格單元齊整地排列在區域邊界處。依據此理論，文獻 ［4］將生成三角形曲面網格的前沿推進算法擴展到四邊形曲面網格生成，新算法從邊界處逐個單元前移前沿，解決了逐排推進鋪路法中出現的魯棒性問題，能夠生成理想的四邊形網格；劉晶等［2］對Q－Morph算法中確定側邊方法進行了改進，雖然使網格轉化過程中殘留三角形單元數目大大減少，但是不能保證生成四邊形網格的質量；陳陽等［5］改進了前沿段生長算法，相比原算法，很好地解決了在幾何質量較差、疏密變化較大的三角網格區域不能生成規則四邊形網格的問題，并通過改變網格頂點度和設置生長約束條件，一定程度上改善了四邊形網格質量，但該算法對疏密度差別較大的三角網格區域使用相同方法生成四邊形，影響網格轉化的效率和質量。

實驗發現在哪一個前沿段上形成四邊形決定四邊形生長方向且嚴重影響最終生成四邊形網格的質量，因此如何選擇基段成為網格轉化的關鍵。本文在生成點云前，使用一種基于多極線約束的重疊點云刪除算法［6］，以避免產生冗余點。對生成的點云數據使用曲率特征簡化［7］后，劃分出平坦區和非平坦區，然后三角化點云，優化形成的三角網格，最后根據提出的基段選擇算法確定待生成四邊形的起始邊。

1 名詞解釋

節點、單元、棱和段：在網格生成中，節點和單元是基本實體。節點是點云的數據點，單元是一系列節點以逆時針方向首尾相連而成。棱是單元的一側，連接兩個端節點。為了區分一般棱，把位于三角網格與四邊形網格邊界上的棱稱為段。如圖1所示。

圖1 名詞解釋

合并前沿、前沿段、額角和基段：合并前沿是待轉換區域中所有邊界段的集合。合并前沿把已轉換區域從待轉換區域中分離開，且必須形成至少一個閉合環。合并前沿上的棱為前沿段；在前沿段兩端節點處與其相鄰前沿段的夾角稱為額角。作為形成四邊形起始邊的前沿段稱為基段。如圖1所示。

2 標記平坦區和非平坦區

定義1 設X＝ ｛x1，x2，…，xn｝是將要建立曲面S上的n 個點集，集合X 中離任意點xi最近的k 個點稱為xi的K－鄰域，記為K（xi）。

本文利用樹的層級結構對空間中每一個點進行K－鄰域搜尋［8］。然后在每一個點如xi的K－鄰域基礎上擬合有約束的最小二乘平面P（xi）作為將要建立曲面在該點的切平面，K（xi）的形心位置是切平面的中心點oi。計算切平面P（xi）的單位法矢ni，使K（xi）內的k 個點到這個切平面距離的平方和達到最小，然后以ni為xi點的單位法矢。為了計算ni的值，需建立K（xi）的協變矩陣

式中：（p－oi）是列向量， （p－oi）T是 （p－oi）的轉置，CV 是一個對稱的半正定3×3 矩陣，記矩陣CV 的最小特征值所對應的單位特征向量為ei，法矢ni取對ei同向化處理后所對應的值。將求得所有數據點的法矢量存在一個一維數組Vector中，順序與原始數據點數組相對應。

如果K（xi）準確表達了將要建立曲面S 在xi點附近的幾何形狀，那么曲面S在xi點上的曲率越大，點xi到切平面P（xi）的垂直距離越大，曲面S在xi點上的曲率越小，點xi到切平面P（xi）的垂直距離越小。可得出下面的定義：

定義2 設X＝ ｛x1，x2，…，xn｝是將要建立曲面S上的n 個點集，對每一個數據點xi用它的K－鄰域擬合將要建立曲面S 在該點處的切平面P（xi），D＝ ｛d1，d2，…，dn｝為每一個點到它的切平面距離的絕對值

切平面距離計算如圖2所示。

圖2 切平面距離計算

為標記平坦點和非平坦點，根據模型幾何形狀設定距離閾值Ωd，一般取經驗值Ωd＝0.1。如果di＜Ωd，則該點為平坦點，否則為非平坦點。確定出所有點的類型后，用長度為N的一維數組FP記錄其類型，平坦點標記為1，非平坦點標記為0。在此基礎上再劃分出平坦區域和非平坦區域。

3 三角網格質量評價與修正

三角網格和網格轉化過程中的局部混合網格的質量參差不齊，而修正網格結構 （連通性）會相應提高網格的質量，因此在網格轉化前和網格轉化過程中，以及局部后處理和全局后處理時，將會頻繁對基本網格進行修正。

3.1 網格質量評價

為了使背景三角網格的單元大小分布與度量規范相兼容，在三角網格所有節點處定義曲面度量張量［9］M。

在3D 空間，度量張量是一個對稱正定矩陣M3D，如Det（M3D）＞0，同時矩陣M3D的特征對 （gi，λi） （i＝1，2，3）定義網格的主要延伸方向和節點空間。通過利用M3D，規范化空間中的基本矢量dε和3D 空間中對應矢量dx 之間的長度尺度變換為

在參數空間中，通過使用等式 （x，y，z）T＝r（u，v），把節點密度的用戶規范效果M3D和曲面映射結合在一起的2×2曲面度量張量M 可以定義為

假設參數空間中存在兩點p1（u1，v1）和p2（u2，v2），則相對于M 兩點p1和p2間的距離定義為

通過利用度量規范描述曲面映射 （參數化），三角化和四邊形轉化方法都可以被應用到范圍廣泛的曲面，包括那些具有快速變化的曲面偏導數和奇異點。

若存在△P1P2P3，由長度尺度測量和曲面度量張量，并根據式 （3）和式 （4），它的形狀質量評價如下

式中

在式（7）中，Mi是Pi點處的曲面度量，Det（M）是矩陣M 的行列式。注意等式 （6）中的形狀質量僅僅取決于三角形的形狀變形（相對于度量張量），而與邊長偏差沒有關系。由于任何尺寸等邊三角形的最大值都等于1。鑒于此，使用結合形狀和尺寸的參數來表示三角單元總體質量

本文使用結合形狀和尺寸的參數珘μ 對三角形總體質量進行評價，這比使用α 質量因子更能在質量參差不齊的三角單元中分辨出不滿足要求的三角單元。

3.2 修正策略

完成三角網格質量評價后，采用以下4種網格修正策略對不滿足要求的網格進行修正。

（1）棱交換：用連接相對節點的新棱替換相鄰兩單元原有共有棱，如圖3 （a）所示。對于一個有效的交換操作，在度量空間和參數空間中兩對新形成的角度和都應小1800，如圖3 （b）所示。

圖3 棱交換

（2）棱塌陷（單元刪除）：在棱塌陷操作中，把滿足一定條件的兩個單元共有棱刪除，共有棱兩端的節點合并成一個節點。因此與該棱相鄰的兩個單元同時被刪除。如圖4所示。

圖4 棱塌陷 （單元刪除）

（3）棱分割：棱分割過程中，在相鄰兩單元共有棱的中點處裂開，兩個相鄰單元將會被分為4個單元，如圖5所示。

圖5 棱分割

（4）節點刪除：當一個節點被3個或4個三角形單元圍繞時，可以執行節點刪除操作。如果節點被3個三角形包圍，如圖6 （a）所示，它們將合并為一個三角單元；對于節點由4 個三角形單元圍繞的情況，如圖6 （b）所示，將生成一條新棱，且這條棱為四邊形長度最短的對角線。

圖6 節點刪除

為了改善三角網格質量，利用以上4種網格修正方法優化背景三角單元。首先遍歷計算三角形網格中每一個三角形單元的形狀和尺寸參數，若不在規定的閾值Ω 范圍內 （0.3≤Ω≤0.7），可根據每一個不滿足質量要求的三角形單元的具體情況選擇合適修正策略對其進行處理。

4 基段選擇算法

在四邊形網格形成過程中，基段決定新單元的位置以及合并前沿前進方向。四邊形質量優劣以及轉化過程是否順利進行取決于是否選擇最優基段。Q－Morph根據每一個前沿段的H 值和兩個前沿額角選擇最優基段 （優先權給予具有較小的H 值和額角接近理想值θrect＝1350的前沿段），前沿段和額角如圖 （1）所示。利用這種算法選擇基段之前需要知道每個前沿段的H 值和額角，并比較各前沿段H 值的大小，判斷額角是否接近理想值。因為使用該算法確定基段需要計算、比較兩個值，所以選擇程序復雜，花費時間多。并且在四邊形網格形成過程中沒有考慮不同區域幾何形狀的差異，僅使用一種基段選擇方法。除此之外，QMorph使用的額角計算公式僅適用于常規三角網格，而對于各項異性或幾何質量較差的三角網格，計算精度不高。本文算法思想是首先計算出每個前沿段兩端節點處的額角值，然后依據給定的判別條件確定每個前沿段的狀態；在平坦區域，當前沿段狀態相同時，計算每個前沿段的H值，選擇H 值最大的前沿段作為基段；當前沿段狀態不相同時，則選擇狀態級別高的前沿段作為基段；對于不平坦區域，選擇長度小于cHmin（1.5≤c≤2）且狀態級別相對高的前沿段為基段。由于前沿段狀態不同是大概率事件，因此算法減少了計算量，提高了算法運行效率。此外，本文算法采用適用于各向異性三角網格的額角計算公式，能夠明顯提高額角的計算精度。

4.1 前沿段H 值和額角

第i個前沿段的H 值計算公式如式 （10）所示

式中：λAk和λBk，k＝1，2分別是節點A 和節點B處度量張量M 的特征值 （參考3.1）。

根據度量張量M 和等式 （5），額角的計算如下

式中：du，dv——無窮小矢量。

4.2 前沿段狀態和分級

前沿段在兩端節點處的額角值決定這個前沿段的狀態，額角計算參考式（11）。對合并前沿中的每一段根據它的狀態進行分級，得到的分級結果影響選擇哪一個前沿段作為基段。

由于一條前沿段狀態是根據前沿段兩個端節點處的額角值確定的，因此前沿段狀態由兩個要素進行定義：一個表示左節點位，另一個表示右節點位。如果一個節點處的額角值小于特定閾值Φ （通常Φ＝1350），則該節點位設置為1，否則設置為0。每一個前沿段狀態分類存放在4種狀態列表中的一個，前沿段的4種狀態如圖7所示，圖中當前前沿段用粗實線顯示。

圖7 前沿段狀態

根據狀態對前沿段分級有兩個目的：一是明確處理前沿段的優先次序。狀態為1－1的前沿段給予第一優先權，然后依次是狀態為0－1、1－0和0－0的前沿段；二是為定義側邊做準備，因為通常定義至少一端已設置狀態位的棱為待形成四邊形的側邊。

4.3 基段選擇算法流程

（1）通過計算確定合并前沿中每個前沿段兩端節點處的額角值，判定出它們的狀態；

（2）判斷合并前沿是否在平坦區，若在平坦區，繼續向下進行，否則轉到 （7）；

（3）若前沿段狀態相同，計算每個前沿段的H 值，選擇H 值最大的前沿段作為基段，否則轉到 （8）；

（4）把作為基段的前沿段從前沿列表中刪除；

（5）在基段基礎上生成四邊形后，若前沿列表非空，把新形成的前沿段添加到前沿列表中，同時更新與新四邊形相鄰前沿段的狀態以及相鄰三角網格的結構，然后返回（1），否則轉到 （6）；

（6）若前沿段列表為空，算法結束；

（7）若在不平坦區，選擇長度小于cHmin且狀態級別相對高的前沿段為基段，然后轉到 （4）；

（8）若前沿段狀態不同，選擇等級高的前沿段作為基段，并轉到 （4）。

根據算法選擇出基段之后，由以下步驟構建一個新四邊形單元：①確定待生成四邊形的左右側邊；②執行頂邊再生操作，連接兩個側邊的頂節點，生成頂邊；③合并四條邊圍成的閉合區域內所有三角形，形成一個四邊形單元。

把三角網格轉化為四邊形網格后，執行全局后處理。應用一組標準網格改進程序和網格光順程序［10］，進一步提升四邊形網格質量。轉化完成后形成的四邊形網格與三角網格自由的拓撲連接關系不同，具有規則的連接，網格邊沿主曲率方向分布，能自然表現出模具幾何形狀的變化［11］。

5 實 例

以某汽車模具為例進行了實例驗證，實驗結果如圖8、圖9所示。其中圖8 （a）是模具點云經三角化和優化后得到的三角網格模型，它有156886個三角形；圖8 （c）是轉化后生成的四邊形網格模型，四邊形數為70229 個，圖中顯示出模型曲面處的網格密度大，平坦區域網格密度小，并且生成的四邊形單元勻稱，整個模型很好保持了模具的特征，達到了不同區域分別采用相適宜方法轉化的效果。

圖8 模具點云網格化

圖9 模具點云曲面擬合

表1是Q－Morph改進算法與本文算法在轉化速率和網格質量方面的比較。比較結果表明本文算法達到了預期效果。表中質差率為質量差的四邊形數占總數的百分比。

表1 兩種算法的數據比較

6 結束語

本文在優化三角網格后，使用基于區域劃分的基段選擇算法確定待形成四邊形基段，其優勢有以下3點：

（1）采用長度尺度測量和曲面度量張量相結合的三角形質量評價準則，可以更準確檢測出質量差的三角形單元。

（2）使用新額角計算公式，提高了計算精度，從而更準確對前沿段進行分級；

（3）簡化選擇基段的判斷條件，并能夠根據不同區域采用與之相適宜的基段選擇方法，因此提高了網格轉化效率以及增大了選擇出最優基段的概率。

實例驗證在優化三角網格的基礎上，使用本文算法提高了四邊形網格生成效率，明顯改善了最終生成的四邊形網格質量，保證了隨后進行的曲面重構效果，并為下一步模型改進、缺損模具修復以及數控加工創造了有利條件。

此外，算法還存在不足之處：平坦區網格密度要求盡可能小，但是在平坦區由于沒有結合前沿段的級別和長度選擇最優基段，會造成平坦區網格沒有達到最理想的密度。如果兩者都考慮，會大幅增大計算量，也會對網格質量產生一定影響。

［1］Park C，Noh JS，Jang IS，et al.A new automated scheme of quadrilateral mesh generation for randomly distributed line constraints ［J］.Computer－Aided Design，2007，39 （4）：258－267.

［2］LIU Jing，NIE Yufeng，SU Shaopu.Indirect method of quadrilateral mesh generation ［J］.Computer Engineering and Application，2010，46 （2）：44－47 （in Chinese）. ［劉晶，聶玉峰，蘇少普.四邊形網格間接生成方法 ［J］.計算機工程與應用，2010，46 （2）：44－47.］

［3］Lo SH.Delaunay triangulation of nonuniform point distributions by means of multigrid insertion ［J］.Finite Elements in Analysis and Design，2013，63：8－22.

［4］CHEN Jianjun，ZHEN Jianjing，JI Tingwei，et al.Advancing front quadrilateral surface mesh generation ［J］.Chinese Journal of Computational Mechanics，2011，28 （5）：779－784 （in Chinese）.［陳建軍，鄭建靖，季廷煒，等.前沿推進曲面四邊形網 格 生 成 算 法 ［J］.計 算 力 學 學 報，2011，28 （5）：779－784.］

［5］CHEN Yang，CUI Hanguo，LIU Jianxin，et al.Improved algorithm for advancing front growth in quadrilateral mesh generation ［J］.Computer Engineering，2011，37 （9）：291－293 （in Chinese）.［陳陽，崔漢國，劉健鑫，等.四邊形網格生成中的前沿邊生長改進算法 ［J］.計算機工程，2011，37（9）：291－293.］

［6］GUO Jin，YUAN Jianying，CHEN Xiaoning.Algorithm for removing redundancy points based on multiview geometric［J］.Computer Engineering ＆Design，2014，35 （3）：958－962 （in Chinese）.［郭進，袁建英，陳小寧.基于多視幾何的重疊點云 刪 除 算 法 ［J］.計 算 機 工 程 與 設 計，2014，35 （3）：958－962.］

［7］DAI Xing，CUI Hanguo，HU Huaiyu.Fast data point simplification algorithm based on curvature character［J］.Journal of Computer Applications，2009，29 （11）：3030－3032 （in Chinese）.［代星，崔漢國，胡懷宇.基于曲率特征的點云快速簡化算法 ［J］.計算機應用，2009，29 （11）：3030－3032.］

［8］Connor M，Kumar P.Fast construction of K－nearest neighbor graphs for point clouds ［J］.IEEE Transactions on Visualization ＆Computer Graphics，2010，16 （4）：599－608.

［9］Xie Hehu，Yin Xiaobo.Metric tensors for the interpolation error and its gradient in Lp norm ［J］.Journal of Computational Physics，2014，256：543－562.

［10］HU Shimin，LAI Yukun，YANG Yongliang.A curvature flow based fairing algorithm of quad－dominant meshes ［J］.Chinese Journal of Computers，2008，31 （9）：1622－1628 （in Chinese）.［胡事民，來煜坤，楊永亮.基于曲率流的四邊形主導網格的光順方法 ［J］.計算機學報，2008，31 （9）：1622－1628.］

［11］ZHU Weipeng，GAO Chengying，LUO Xiaonan.Global anisotropic quad－dominant remeshing ［J］.Journal of Soft，2012，23 （5）：1305－1314 （in Chinese）.［朱為鵬，高成英，羅笑男.全局各向異性四邊形主導網格重建方法 ［J］.軟件學報，2012，23 （5）：1305－1314.］