基于Laplace調(diào)和方程的網(wǎng)格重構(gòu)算法

陳華偉　伍權(quán)　徐衛(wèi)平　余澤云

摘要：針對(duì)三角形網(wǎng)格向四邊形網(wǎng)格的轉(zhuǎn)化問(wèn)題，基于調(diào)和方程構(gòu)建模型梯度場(chǎng)，追蹤表面流線，實(shí)現(xiàn)了參數(shù)化網(wǎng)格重構(gòu)。首先，建立了基于離散Laplace方程的三角網(wǎng)格梯度場(chǎng)理論模型、數(shù)據(jù)結(jié)構(gòu)模型和稀疏矩陣求解方案;其次，提出了局部坐標(biāo)變換和參數(shù)方程相結(jié)合求解流線節(jié)點(diǎn)的統(tǒng)一算法，并針對(duì)流線跟蹤無(wú)交點(diǎn)、有多個(gè)交點(diǎn)等特殊情況，提出了梯度收斂、最短距離和參數(shù)極值等優(yōu)選策略;最后，通過(guò)模型實(shí)驗(yàn)驗(yàn)證了算法。結(jié)果表明，流線網(wǎng)格具有等參、閉合特點(diǎn)，復(fù)雜模型網(wǎng)格劃分沒(méi)有歧義，而且網(wǎng)格質(zhì)量隨網(wǎng)格密度增加而提高。因此，相對(duì)于傳統(tǒng)幾何重構(gòu)算法，數(shù)學(xué)方法對(duì)網(wǎng)格重構(gòu)表達(dá)具有魯棒性和唯一性，且應(yīng)用場(chǎng)景更廣泛。

關(guān)鍵詞：算法理論;網(wǎng)格重構(gòu);參數(shù)化網(wǎng)格;調(diào)和方程;流線追蹤

中圖分類號(hào)：TP3016文獻(xiàn)標(biāo)志碼：A

CHEN Huawei，WU Quan，XU Weiping， et al.Mesh reconstruction algorithm based on Laplace harmonic equation[J].Journal of Hebei University of Science and Technology，2019，40（3）：199-207.Mesh reconstruction algorithm based on Laplace harmonic equation

CHEN Huawei1， WU Quan1， XU Weiping1， YU Zeyun2

（1.School of Mechanical and Electrical Engineering， Guizhou Normal University， Guiyang， Guizhou 550025， China;2.Department of Computer Science， University of Wisconsin-Milwaukee， Milwaukee 53206， USA）

Abstract：Upon the problem of transforming triangular mesh into quadrilateral mesh， gradient field via harmonic equation is created， integral flow is tracked， and parameterized mesh is reconstructed. First， gradient field construction theory， data structure model and solution scheme of sparse matrix are constructed. Then， one uniform algorithm for solving flow line node by integration of local coordinate transformation and parameterized equation is advanced， and schemes like gradient convergence， shortest distance and extreme parameter are optimally occupied upon special cases of no intersection or multiple intersections when tracing flow line. Finally， the algorithm is verified via case studies. The result shows that flow lines are characterized as iso-parameterized and closed， mesh reconstruction of complicate models has no bifurcation， and mesh quality is increased with its density. The result proves that mathematical method has robustness and uniqueness for mesh reconstruction representation compared with traditional geometry method， and it will have more application scenarios.

Keywords：theory of algorithms; mesh reconstruction; parameterized mesh; harmonic equation; integral flow tracking

非結(jié)構(gòu)化網(wǎng)格沒(méi)有規(guī)則的拓?fù)浣Y(jié)構(gòu)，網(wǎng)格劃分靈活性強(qiáng)，但是占用內(nèi)存大，后續(xù)應(yīng)用的數(shù)值計(jì)算量大;結(jié)構(gòu)化網(wǎng)格有更好的網(wǎng)格質(zhì)量，邊界擬合性好，數(shù)值計(jì)算性能高、擴(kuò)散性小。網(wǎng)格轉(zhuǎn)化、重構(gòu)及其衍生問(wèn)題是學(xué)者持續(xù)關(guān)注的課題。以Q-Morph算法[1-2]為代表的幾何推演迭代方法存在計(jì)算量大，優(yōu)質(zhì)網(wǎng)格生成還需要相當(dāng)?shù)暮笾锰幚淼葐?wèn)題。其他的如方向場(chǎng)法[3]、邊界簡(jiǎn)化法[4]、最少奇異點(diǎn)模板法[5]、體素法[6]等，可以通過(guò)參數(shù)優(yōu)化生成自適應(yīng)四邊形網(wǎng)格，但都存在單元方向和單元密度難以控制等問(wèn)題。在數(shù)學(xué)上，離散網(wǎng)格模型即為分段流形?；贛orse理論，為流形曲面選擇合適的Morse函數(shù)，就能保證Morse-Smale復(fù)形將該曲面劃分為純四邊形結(jié)構(gòu)[7-8]。該方法很容易實(shí)現(xiàn)模型特征約束和單元密度控制，但是由于模型的復(fù)雜性和標(biāo)量場(chǎng)數(shù)據(jù)的多樣性，很難魯棒地構(gòu)建Morse-Smale復(fù)形。為此，越來(lái)越多的學(xué)者開(kāi)始采用調(diào)和方程進(jìn)行等參數(shù)化網(wǎng)格的重構(gòu)研究[9-11]。本文以三角網(wǎng)格模型為研究對(duì)象，展開(kāi)離散Laplace調(diào)和方程創(chuàng)建和求解、流線跟蹤及流線網(wǎng)格重構(gòu)方面的研究。該研究可擴(kuò)展至有限元網(wǎng)格劃分，用于材料、流體力學(xué)等綜合性能分析;可擴(kuò)展至圖像處理領(lǐng)域，用于醫(yī)學(xué)分層圖像、地理等高圖像重構(gòu)和網(wǎng)格分析;還可擴(kuò)展至仿生結(jié)構(gòu)設(shè)計(jì)，用于多孔結(jié)構(gòu)、微觀分形結(jié)構(gòu)建模。

1三角網(wǎng)格模型的場(chǎng)求解

在拓?fù)鋷缀紊希蔷W(wǎng)格模型為節(jié)點(diǎn)V和三角面F集合構(gòu)成的分段線性流形，即M=（V， F）。根據(jù)應(yīng)用需求，為每個(gè)節(jié)點(diǎn)vi∈V賦標(biāo)量值ui，即構(gòu)造M上的標(biāo)量場(chǎng)u：V。

河北科技大學(xué)學(xué)報(bào)2019年第3期陳華偉，等：基于Laplace調(diào)和方程的網(wǎng)格重構(gòu)算法1.1理論分析

對(duì)標(biāo)量場(chǎng)u構(gòu)造Laplace調(diào)和方程[12-13]：Δu=0。? ? （1）給定第一類邊界條件，即狄利克雷（Dirichlet）條件：u|C=φ，? ?（2）其中：C為邊界;φ為邊界函數(shù)。

對(duì)式（1）和式（2）構(gòu)成的約束Laplace方程進(jìn)行離散化，有：Δui=∑j∈Niωij（uj-ui）=0， ?????????（3）

ui=ci，ci∈C 。? ?（4）圖1三角形上的標(biāo)量和矢量

Fig.1Scalar and vector on triangle

Ni表示節(jié)點(diǎn)i的一階鄰域，三角網(wǎng)格中Ni為i的一環(huán)（1-ring）鄰域。三角網(wǎng)格離散化，一般使用cotangent余切值計(jì)算Laplace系數(shù)[14]：ωij=[cot（αij）+cot（βij）]/2 ，? ?（5）式（3）—式（5）構(gòu)成線性方程組，通過(guò)方程組的求解，得標(biāo)量場(chǎng)ui∈u。

進(jìn)一步地，可求解其對(duì)應(yīng)的梯度場(chǎng)（矢量場(chǎng)）。對(duì)于分段線性流形，在每個(gè)分段內(nèi)梯度矢量為常量。如圖1所示，已知三角形T0（v0，v1，v2）或（e0，e1，e2） ，對(duì)應(yīng)的標(biāo)量參數(shù)u0，u1，u2，三角形面的單位法矢為N，則T0的梯度G0可通過(guò)以下線性系統(tǒng)求解：[e0，e1，N]T[G0]=[v1-v0，v2-v1，n]T[G0]=[u1-u0，u2-u1，0]T 。 ????（6）觀察式（6）易知，G0是邊e0和e1的線性組合，即：G0=ae0+be1=[e0，e1][a，b]T 。 ?????????????????????（7）將式（7）代入式（6），則有：[e0，e1]T[G0]=[e0，e1]T[e0，e1][a，b]T=[v20，v0v1;v1v0，v21][a，b]T=[u1-u0，u2-u1]T 。（8）至此，系統(tǒng)降為二維線性系統(tǒng)，可直接通過(guò)手工計(jì)算求解參數(shù)a，b。令Δ=v20×v21-（v0v1）2，則有：a=[（u1-u0）×v21-（u2-u1）×v0v1]/Δ ，

b=[（u2-u1）×v20-（u1-u0）×v0v1]/Δ 。把參數(shù)a，b代入式（7）即可得G0，對(duì)應(yīng)地，可求正交梯度：G⊥=n×G0 。 ???????????????????????????????????（9）1.2數(shù)據(jù)準(zhǔn)備

面向Laplace方程的構(gòu)造，三角網(wǎng)格模型的數(shù)據(jù)準(zhǔn)備工作主要有2個(gè)方面，即節(jié)點(diǎn)數(shù)據(jù)結(jié)構(gòu)的創(chuàng)建和節(jié)點(diǎn)類型劃分。

1）節(jié)點(diǎn)數(shù)據(jù)結(jié)構(gòu)的創(chuàng)建

節(jié)點(diǎn)數(shù)據(jù)結(jié)構(gòu)應(yīng)同時(shí)反映鄰域點(diǎn)和鄰域三角形，因此，最簡(jiǎn)單的方法是記錄節(jié)點(diǎn)的一環(huán)鄰域邊，本文節(jié)點(diǎn)結(jié)構(gòu)體NodeDetail中主要記錄節(jié)點(diǎn)的對(duì)邊：

struct NodeDetail

{

long v[2];//記錄對(duì)邊兩點(diǎn)的序號(hào)

long idxF;//記錄三角形的序號(hào)

long posF;//posF=0， 1， 2，記錄當(dāng)前點(diǎn)在該三角形三點(diǎn)中的順序號(hào)

};

顯然1個(gè)NodeDetail對(duì)象代表了節(jié)點(diǎn)的1個(gè)鄰接三角形，如果節(jié)點(diǎn)i有Ni個(gè)鄰接三角形，則應(yīng)有Ni個(gè)NodeDetail對(duì)象。

2）節(jié)點(diǎn)分類

為了保證Laplace方程有非零解，還需要設(shè)置邊界約束條件。對(duì)應(yīng)于Laplace方程，將三角形網(wǎng)格中事先指定的約束點(diǎn)和網(wǎng)格邊界點(diǎn)都?xì)w為約束節(jié)點(diǎn)，其他節(jié)點(diǎn)（均為內(nèi)部節(jié)點(diǎn)）全部歸為非約束節(jié)點(diǎn)，如圖2所示。如果模型為閉合曲面，則認(rèn)為模型無(wú)邊界點(diǎn)，只有事先指定點(diǎn)為約束節(jié)點(diǎn)。

1.3Laplace方程的構(gòu)造與求解

式（3）中，每個(gè)節(jié)點(diǎn)只與其臨近的一階鄰域節(jié)點(diǎn)存在系數(shù)關(guān)系，與其他非鄰接節(jié)點(diǎn)無(wú)系數(shù)關(guān)系，具有局部影響性。因此，面向全局網(wǎng)格節(jié)點(diǎn)構(gòu)建Laplace方程，其系數(shù)矩陣是稀疏矩陣[15]，對(duì)此本文采用了SuperLU稀疏矩陣求解方案，其時(shí)間和空間復(fù)雜度分別為T（n）=O（n1.5），S（n）=O（nlogn）。

為了減少數(shù)據(jù)存貯量，需在程序設(shè)計(jì)中直接構(gòu)造稀疏矩陣。建立LaplaceData數(shù)據(jù)結(jié)構(gòu)，存貯與當(dāng)前節(jié)點(diǎn)鄰接節(jié)點(diǎn)的索引和權(quán)重系數(shù)，并使用SortLaplaceData函數(shù)按節(jié)點(diǎn)索引對(duì)鄰接節(jié)點(diǎn)進(jìn)行排序。

struct LaplaceData

{

long idx;//節(jié)點(diǎn)索引

float data;//cotangent計(jì)算值

};

static int SortLaplaceData（const void* e1， const void* e2）

{

LaplaceData* e11 = （LaplaceData*）e1;

LaplaceData* e22 = （LaplaceData*）e2;

return e11->idx - e22->idx;

}

直觀上，鄰域節(jié)點(diǎn)的LaplaceData數(shù)據(jù)按行存貯。但SuperLU采用Harwell-Boeing矩陣存貯格式[16]，即列壓縮格式，需要使用dCompRow_to_CompCol函數(shù)將行存貯參數(shù)轉(zhuǎn)換為列存貯參數(shù)，再使用dCreate_CompCol_Matrix函數(shù)構(gòu)造列壓縮存貯的稀疏矩陣。最后，使用dgssv函數(shù)解矩陣方程，求解結(jié)果就是對(duì)應(yīng)于每個(gè)節(jié)點(diǎn)的參數(shù)ui∈u。進(jìn)一步地，可求解三角網(wǎng)格上的梯度場(chǎng)。

2流線跟蹤

交互輸入種子點(diǎn)，即可通過(guò)該點(diǎn)沿梯度方向追蹤一條唯一的流線。設(shè)種子三角形Ts（v0，v1，v2），T0內(nèi)一點(diǎn)ptSeed為種子點(diǎn)，Ts內(nèi)的梯度矢量為Gs，Ts中梯度線交點(diǎn)為ptInt0和ptInt1，對(duì)應(yīng)的共邊三角形為T0和T1，如圖3所示。

設(shè)最小極點(diǎn)和最大極點(diǎn)分別為vmin和vmax，則有如下流線跟蹤的總體過(guò)程：

1）在Ts內(nèi)，過(guò)點(diǎn)ptSeed和方向Gs做直線，與Ts的2個(gè)交點(diǎn)ptInt0和ptInt1，即為初始節(jié)點(diǎn);

2）在共邊三角形T0內(nèi)，從ptInt0出發(fā)，沿T0的負(fù)梯度方向搜索下一個(gè)交點(diǎn)，記為ptInt2;

3）重復(fù)步驟2），直到ptInt2=vmin;

4）在共邊三角形T1內(nèi)，從ptInt1出發(fā)，沿T1的梯度方向搜索下一個(gè)交點(diǎn)，記為ptInt2;

5）重復(fù)步驟4），直到ptInt2=vmax;

6）以上交點(diǎn)按vmin至vmax的正梯度方向整理輸出為流線。

2.1初始梯度線的計(jì)算

該問(wèn)題可描述為求過(guò)ptSeed且方向?yàn)镚s的直線L與Ts各邊的2個(gè)交點(diǎn)ptInt0和ptInt1。為了便于敘述，假設(shè)網(wǎng)格內(nèi)三角形方向DT為逆時(shí)針，法向?yàn)镹。令直線L與Ts的三條邊分別求交，最終解出2個(gè)交點(diǎn)。

本文使用局部坐標(biāo)變換和參數(shù)方程法統(tǒng)一求解。首先，以三角形一個(gè)角點(diǎn)為坐標(biāo)原點(diǎn)，以該點(diǎn)的下一條有序邊（DT方向）為X軸，使用Y=N×X，即X軸沿DT方向旋轉(zhuǎn)90°構(gòu)造Y軸，形成三角面上的局部坐標(biāo)系XOY。假設(shè)當(dāng)前角點(diǎn)為v0，則X軸與v0v1邊重合，方向一致。接下來(lái)，將原三維坐標(biāo)系參數(shù)轉(zhuǎn)換至XOY二維坐標(biāo)系，分別構(gòu)造直線L和v0v1的參數(shù)方程，并聯(lián)立求解即可求得直線L和v0v1的交點(diǎn)，如圖4所示。

v0v1直線方程為

y=0，0≤x≤v1.x。

直線L的參數(shù)方程為

x=ptSeed.x+t×Gsx，y=ptSeed.y+t×Gsy，0≤t≤1。

聯(lián)立求解，可得交點(diǎn)參數(shù)t，記為t0：

t0=-ptSeed.y/ptSeed.x 。

相應(yīng)地：

x0=ptSeed.x+t0×Gs.x 。

接下來(lái)，驗(yàn)證x0的取值范圍，如果滿足0≤x0≤v1.x，則求得一個(gè)交點(diǎn)（x0，0） ，否則直線L與三角形邊v0v1無(wú)交點(diǎn)。

如果有交點(diǎn)，則需要如下的后置處理：

1）通過(guò)逆變換將二維點(diǎn) （x0， 0） 轉(zhuǎn)換至原三維坐標(biāo)系，求得交點(diǎn)的三維坐標(biāo);

2）搜索與Ts共v0v1邊的三角形T0或T1，用于梯度線的連續(xù)搜索。

依次將原點(diǎn)輪換至另外2個(gè)三角形角點(diǎn)，按上述過(guò)程繼續(xù)求解直線L與三角形另外兩條邊v1v2，v2v0的交點(diǎn)。三次求解結(jié)束后，最終可得2個(gè)交點(diǎn)ptInt0和ptInt1。

2.2從邊出發(fā)追蹤交點(diǎn)

該問(wèn)題可描述如下：已知三角形T0（v0，v1，v2），T0內(nèi)的梯度矢量為G0，其中一邊（假設(shè)為v0v1）上有一點(diǎn)ptSeed，過(guò)ptSeed且方向?yàn)镚0的直線L與T0的另外兩邊有一個(gè)交點(diǎn)ptInt，求該交點(diǎn)，如圖5所示。

1）L與v1v2求交

構(gòu)造圖5所示的局部坐標(biāo)系XOY。顯然，XOY坐標(biāo)系中，ptSeed.y=0。直線L的參數(shù)方程同2.1節(jié)，直線v1v2采用斜截式方程：

x=ky+b 。

不難得出，斜率k=（v2.x-v1.x）/v2.y，X軸截距b=v1.x。以上各式聯(lián)立求解，得交點(diǎn)坐標(biāo)：

（x0，y0）=（ptSeed.x+t0×G0.x，t0×G0.y）。

接下來(lái)，判斷交點(diǎn)有效性。構(gòu)造直線v1v2的參數(shù)方程：

x=v1x+s×（v2.x-v1.x），y=s×v2.y，0≤s≤1。

將y0代入，可解參數(shù)s：

s=y0/v2.y 。

如果滿足0≤s≤1，則說(shuō)明交點(diǎn)（x0，y0）在線段v1v2之內(nèi)，為有效交點(diǎn)，否則認(rèn)為直線L與v1v2邊無(wú)交點(diǎn)。

有交點(diǎn)的情況下，按前述方法做后置處理，提前結(jié)束搜索;否則，繼續(xù)求解L與v2v0的交點(diǎn)。

2）L與v2v0求交

參數(shù)方程可表示為

x=s×v2.x，y=s×v2.y，0≤s≤1。

同理求解交點(diǎn)坐標(biāo)，并判斷參數(shù)s的有效性。

上述計(jì)算中，需要針對(duì)以下2種情況做特殊處理：

1）當(dāng)s=0或1，說(shuō)明交點(diǎn)與節(jié)點(diǎn)重合。此時(shí)，無(wú)法輸出共邊三角形，只能按2.3節(jié)的方法繼續(xù)搜索。

2）如果均不滿足0 ≤s≤1條件，即有向直線L與線段v1v2和v2v0均無(wú)交點(diǎn)，此時(shí)，將會(huì)出現(xiàn)梯度收斂現(xiàn)象。如圖6所示，當(dāng)前三角形T0內(nèi)，過(guò)點(diǎn)ptSeed沿G0方向，無(wú)交點(diǎn)。但是，綜合T0及其共邊三角形T2可知，梯度矢量G0和G2是向節(jié)點(diǎn)v1收斂的?？梢?jiàn)，該情況下節(jié)點(diǎn)v1應(yīng)為所求點(diǎn)。

2.3從節(jié)點(diǎn)出發(fā)追蹤交點(diǎn)

當(dāng)搜索起點(diǎn)為三角形節(jié)點(diǎn)時(shí)，本文按以下方法處理：遍歷當(dāng)前節(jié)點(diǎn)的所有鄰接三角形，在每個(gè)三角形內(nèi)分別求交點(diǎn)，并判斷交點(diǎn)有效性。

該問(wèn)題可描述如下：已知節(jié)點(diǎn)v0及其一階鄰域三角形Ti，i=0，1，…，m，m為鄰域三角形數(shù)，對(duì)應(yīng)梯度矢量為Gi，對(duì)每個(gè)三角形Ti，求過(guò)v0且方向?yàn)镚i的直線L與v0對(duì)邊（v1v2）的交點(diǎn)，最終在所有有效交點(diǎn)中輸出一個(gè)最優(yōu)點(diǎn)，如圖7所示。

圖7Case3：從節(jié)點(diǎn)出發(fā)追蹤

Fig.7Case3： Track from node如2.2節(jié)所述，只需將搜索起點(diǎn)置于原點(diǎn)，即ptSeed=v0，即可同理求解。

接下來(lái)根據(jù)參數(shù)t0和s判斷交點(diǎn)有效性，如表1所示。

序號(hào)條件交點(diǎn)1t>0， s ≤ 0v12t>0， s ≥ 1v23t>0， 0 < s < 1ptInt4t0≤0無(wú)交點(diǎn)

如果出現(xiàn)多個(gè)有效節(jié)點(diǎn)，則需要進(jìn)行沖突消解處理?？筛鶕?jù)實(shí)驗(yàn)效果，實(shí)施以下2種消解策略。

1）最短距離策略

對(duì)所有有效交點(diǎn)，計(jì)算至v0的距離，將最短距離的交點(diǎn)作為輸出。

2）s極值策略

對(duì)所有有效交點(diǎn)，比較對(duì)應(yīng)的參數(shù)s，s<0時(shí)，令s0=|s|，s>1時(shí)，令s0=|s-1|，將s0極小值對(duì)應(yīng)的交點(diǎn)作為輸出。

3流線網(wǎng)格生成

按上述方法，可在三角網(wǎng)格表面追蹤一條過(guò)種子點(diǎn)ptSeed，且兩端分別抵達(dá)極點(diǎn)vmin，vmax的流線，記為種子線lnSeed。接下來(lái)基于種子線，重建模型表面的四邊形網(wǎng)格。設(shè)梯度方向及其正交方向網(wǎng)格密度為nu×nv，則有以下重建流程：

1）分段累計(jì)種子線lnSeed的長(zhǎng)度，等分lnSeed為nu份，得到（nu-1）個(gè)等分點(diǎn)ptSegi，i=0，1，…，nu-2;

2）以等分點(diǎn)ptSegi為新的種子點(diǎn)ptSeedi，從該點(diǎn)出發(fā)，沿G⊥方向跟蹤正交流線（等參線）;

3）在（nu-1）條G⊥向流線中，選擇最長(zhǎng)的一條作為G⊥向種子線lnSeed⊥;

4）等分lnSeed⊥為nv份，再?gòu)倪@些等分點(diǎn)出發(fā)沿G向追蹤，得到（nv-1）條G向流線。

4實(shí)驗(yàn)

在Windows系統(tǒng)下，以VC和OpenGL為開(kāi)發(fā)平臺(tái)，集成SuperLU稀疏矩陣求解方案，實(shí)現(xiàn)了上述算法。

圖8使用球體網(wǎng)格模型（節(jié)點(diǎn)數(shù)N=2 562，三角形數(shù)T=5 120）展示了算法各主要步驟的結(jié)果。圖8 b）中，極點(diǎn)附近梯度矢量從四面八方指向極點(diǎn)，這保證了所有G向流線都能最終收斂至極點(diǎn)，如圖8 f）圖8球體模型（N=2 562， T=5 120）

Fig.8Sphere model（N=2 562， T=5 120）

所示;圖8 d）中，G⊥向流線即等參線均為閉合曲線，從lnSeed上的每個(gè)等分點(diǎn)均可追蹤出一條等參線，該等分點(diǎn)是對(duì)應(yīng)等參線的起止點(diǎn)。

圖9使用人體網(wǎng)格模型展示了不同種子點(diǎn)對(duì)流線網(wǎng)格的影響。圖9 a）和圖9 c）分別在人體正面和背面選定不同的種子點(diǎn)ptSeed，可分別在人體正面和背面追蹤種子線;圖9 b）和圖9 d）為2種情況下的流線網(wǎng)格，對(duì)比來(lái)看流線網(wǎng)格并未發(fā)生變化。理論上，只要極點(diǎn)vmin，vmax一定，所求的梯度場(chǎng)G及正交梯度場(chǎng)G⊥就是確定的，對(duì)應(yīng)地，模型表面流線和等參線也是確定的。

圖10使用股骨網(wǎng)格模型展示了不同網(wǎng)格密度下的流線網(wǎng)格。圖10 b）—d）分別是網(wǎng)格密度nu×nv為10×10，100×20，400×100的結(jié)果。圖10 c）中外側(cè)髁部位（圖示右上部）沒(méi)有G向流線，圖10 d）中，在網(wǎng)格加密后，同一部位劃分出了G向流線，對(duì)該部位及整體模型的表達(dá)也更準(zhǔn)確。

以角度偏斜度（Skewness）為網(wǎng)格質(zhì)量評(píng)價(jià)參數(shù)，Skewness=max[（θmax-θe）/（180-θe），（θe-θmin）/θe]，其中θmax和θmin為單元中的最大角和最小角，θe為單元均衡角，本文網(wǎng)格在極點(diǎn)處為三角形網(wǎng)格（圖8 f）），因此θe有60（三角形）和90（四邊形）2個(gè)取值。對(duì)圖10的網(wǎng)格質(zhì)量進(jìn)行分析，10×10，100×20，400×100網(wǎng)格的平均偏斜度分別為0.68，0.33，0.25，分別達(dá)到良、優(yōu)、優(yōu)的網(wǎng)格質(zhì)量等級(jí)，3種情況下，均未出現(xiàn)偏斜度大于090的網(wǎng)格。

與傳統(tǒng)幾何方法相比，如文獻(xiàn)＼[2＼]的Q-Morph方法，調(diào)和方程法是以底層網(wǎng)格的梯度場(chǎng)為指導(dǎo)重新劃分網(wǎng)格，這一點(diǎn)對(duì)復(fù)雜網(wǎng)格同樣實(shí)用，而且，網(wǎng)格精度和質(zhì)量很容易通過(guò)網(wǎng)格密度進(jìn)行調(diào)整。這些特點(diǎn)也體現(xiàn)在實(shí)驗(yàn)結(jié)果之中：1）分段流形的梯度場(chǎng)是流線追蹤的依據(jù)，相對(duì)于傳統(tǒng)幾何方法，數(shù)學(xué)方法對(duì)網(wǎng)格重構(gòu)表達(dá)的魯棒性更高;2）給定模型、約束條件及網(wǎng)格密度等條件，可確定唯一的梯度場(chǎng)和流線網(wǎng)格;3）網(wǎng)格密度越高，流線網(wǎng)格越能準(zhǔn)確表達(dá)原始模型，但網(wǎng)格密度過(guò)高，將延長(zhǎng)流線追蹤時(shí)間，增加存貯空間。

此外，借助稀疏矩陣方案，以上實(shí)驗(yàn)?zāi)Ｐ偷奶荻葓?chǎng)求解均沒(méi)有時(shí)間和空間壓力;而流線追蹤過(guò)程主要是利用線性方程求解交點(diǎn)，也不會(huì)存在性能問(wèn)題。

5結(jié)語(yǔ)

將調(diào)和理論應(yīng)用至分段流形網(wǎng)格重構(gòu)可有效解決幾何重構(gòu)法中存在的算法不穩(wěn)定、迭代時(shí)間長(zhǎng)等問(wèn)題。論文針對(duì)三角化網(wǎng)格重構(gòu)問(wèn)題，使用離散Laplace調(diào)和方程求解模型梯度場(chǎng)，進(jìn)而通過(guò)流線追蹤重構(gòu)等參數(shù)化網(wǎng)格，并通過(guò)算法設(shè)計(jì)、開(kāi)發(fā)和模型實(shí)驗(yàn)驗(yàn)證了該過(guò)程。相對(duì)傳統(tǒng)幾何算法，調(diào)和理論算法更具魯棒性和網(wǎng)格表達(dá)的唯一性，在參數(shù)化網(wǎng)格重構(gòu)研究領(lǐng)域具有擴(kuò)展應(yīng)用的能力，具體如下。

1）雙極點(diǎn)約束模型很好地驗(yàn)證了前文的理論和方法，而且每個(gè)網(wǎng)格內(nèi)都有確定的梯度方向，因此能有效規(guī)避零曲率曲面或平面區(qū)域網(wǎng)格再劃分問(wèn)題。但是，雙極點(diǎn)法只能表達(dá)一個(gè)方向的等參線（G⊥向的閉合流線），流線網(wǎng)格并非等參網(wǎng)格。而且，流線網(wǎng)格也非純四邊形網(wǎng)格，在極點(diǎn)處網(wǎng)格退化為三角形。通過(guò)改變約束條件可有效解決以上問(wèn)題，例如，文獻(xiàn)＼[17—18＼]使用8點(diǎn)約束在模型表面構(gòu)造六面體模型，從而很好地重構(gòu)出純四邊形的等參數(shù)化網(wǎng)格。

2）對(duì)多虧格（genius-n）模型＼[19＼]，本文的等參線追蹤方法需要改進(jìn)。如圖9，人體模型兩腿之間有一個(gè)孔洞，為genius-1模型，在雙腿部位，只在一側(cè)追蹤了等參線。genius-n模型中，在同一高度參數(shù)下可能追蹤出等參線，這一點(diǎn)可通過(guò)完善流線追蹤方法實(shí)現(xiàn)。此外，也可以采用區(qū)域分割法將genius-n模型“降維”為無(wú)虧格（genius-0）模型進(jìn)行處理。在分支特征（圖9手臂部位）或尖銳特征（圖10外側(cè)髁部位）處，即使加密了整體網(wǎng)格，流線還是難以到達(dá)、局部網(wǎng)格加密困難等問(wèn)題，此時(shí)可以采用區(qū)域分割法＼[20＼]對(duì)這些特征區(qū)域單獨(dú)劃分網(wǎng)格。

參考文獻(xiàn)/References：

[1]汪攀， 張見(jiàn)明， 韓磊， 等. 基于帶約束前沿推進(jìn)的四邊形網(wǎng)格生成方法[J]. 湖南大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2017， 44（8）： 29-34.

WANG Pan， ZHANG Jianming， HAN Lei， et al. An advancing front quadrilateral mesh generation method with constraint[J]. Journal of Hunan University（Natural Sciences）， 2017， 44（8）： 29-34.

[2]吳祿慎， 王偉杰， 陳華偉， 等. 基于區(qū)域劃分的各向異性四邊形網(wǎng)格重建方法[J]. 機(jī)械設(shè)計(jì)與制造， 2015（3）： 212-216.

WU Lushen， WANG Weijie， CHEN Huawei， et al. Anisotropic quadrilateral remeshing based on region division[J]. Machinery Design & Manufacture， 2015（3）： 212-216.

[3]李天華. 基于方向場(chǎng)的四邊形重網(wǎng)格化技術(shù)研究[D]. 長(zhǎng)春：吉林大學(xué)， 2017.

LI Tianhua. Research on Quad Remeshing based on Vector Field [D]. Changchun： Jilin University， 2017.

[4]徐崗， 舒來(lái)新， 朱亞光， 等. 邊界簡(jiǎn)化與多目標(biāo)優(yōu)化相結(jié)合的高質(zhì)量四邊形網(wǎng)格生成[J]. 中國(guó)圖象圖形學(xué)報(bào)， 2018， 23（1）： 61-73.

XU Gang， SHU Laixin， ZHU Yaguang， et al. High-quality quadrilateral mesh generation by combining boundary simplification and multi-objective optimization[J]. Journal of Image and Graphics， 2018， 23（1）：61-73.

[5]VERMA C S， SURESH K. Alpha-MST： A robust unified algorithm for quadrilateral mesh adaptation[J]. Computer-Aided Design， 2018， 103（SI）： 47-60.

[6]鄭永川， 關(guān)柏良， 林淑金， 等. 針對(duì)密集點(diǎn)云的快速自適應(yīng)四邊形網(wǎng)格生成算法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 2019， 31（1）： 39-46.

ZHENG Yongchuan， GUAN Boliang， LIN Shujin， et al. Rapid and adaptive quadrilateral mesh generation algorithm from dense point cloud [J]. Journal of Computer-Aided Design & Computer Graphics， 2019， 31（1）： 39-46.

[7]PATANE G. Mesh-based and meshless design and approximation of scalar functions[J]. Computer Aided Geometric Design， 2017， 57： 23-43.

[8]袁潔， 周明全， 耿國(guó)華， 等. ?基于Morse-Smale拓?fù)涮卣鞯奈奈锼槠唇铀惴╗J]. 自動(dòng)化學(xué)報(bào)， 2018， 44（8）： 1486-1495.

YUAN Jie， ZHOU Mingquan， GENG Guohua， et al. Automatic reassembly of fractured fragments using Morse topological features[J]. ACTA Automatica Sinica， 2018， 44（8）： 1486-1495.

[9]齊學(xué)義， 楊帆. 雙調(diào)和方程生成正交曲線坐標(biāo)網(wǎng)格技術(shù)的研究[J]. 水力發(fā)電學(xué)報(bào)， 2000（3）： 84-89.

QI Xueyi， YANG Fan. Study on technique for generating boundary-orthogonal curvilinear coordinate grid using the biharmonic equations [J]. Journal of Hydroelectric Engineering， 2000（3）： 84-89.

[10]LING R， HUANG J， SUN F， et al. Spectral quadrangulation with feature curve alignment and element size control[J]. ACM Transactions on Graphics， 2014， 34（1）： 1-11.

[11]MARCHANDISE E， WIART C， VOS W G， et al. High-quality surface remeshing using harmonic maps（part II）： Surfaces with high genus and of large aspect ratio[J]. International Journal for Numerical Methods in Engineering， 2011， 86（11）： 1303-1321.

[12]ABDULHADI Z， MUHANNA Y， PONNUSAMY S. Dirichlet problem， univalency and schwarz lemma for biharmonic mappings[J]. Mediterranean Journal of Mathematics， 2018， 15（4）： 187.

[13]BEZRODNYKH S， VALSOV V. On a problem of the constructive theory of harmonic mappings[J]. Journal of Mathematical Sciences， 2014， 201（6）： 705-732.

[14]GLICKENSTEIN D. Geometric triangulations and discrete Laplacians on manifolds[J]. Mathematics， 2005：1-43.

[15]紀(jì)國(guó)良， 丁勇， 周曼， 等. 工程計(jì)算中大型稀疏矩陣存儲(chǔ)方法研究[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用， 2018， 39（3）： 217-230.

JI Guoliang， DING Yong， ZHOU Man， et al. Research on the storage method of large-scale sparse matrix in engineering calculation[J]. Journal on Numerical Methods and Computer Applications， 2018， 39（3）： 217-230.

[16]張永杰， 孫秦. 稀疏矩陣存儲(chǔ)技術(shù)[J]. 長(zhǎng)春理工大學(xué)學(xué)報(bào)， 2006， 29（3）： 38-41.

ZHANG Yongjie， SUN Qin. Sparse storage technique for sparse matrix[J]. Journal of Changchun University of Science and Technology， 2006， 29（3）： 38-41.

[17]CHEN Huawei， GUO Ye， ROSTAMI R， et al. Porous structure design using parameterized hexahedral meshes and triply periodic minimal surfaces[C]// Computer Graphics International 2018. Bintan Island：[s.n.]， 2018.

[18]MARTIN T， COHEN E， KIRBY R M. Volumetric parameterization and trivariate B-spline fitting using harmonic functions[J]. Computer Aided Geometric Design， 2009（26）： 648-664.

[19]錢坤， 張家玲， 李映華， 等. 高虧格曲面共形參數(shù)化方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 2017， 29（12）： 2225-2234.

QIAN Kun， ZHANG Jialing， LI Yinghua， et al. Conformal parameterization for high genus surfaces[J]. Journal of Computer-Aided Design & Computer Graphics， 2017， 29（12）： 2225-2234.

[20]LI Bo， QIN Hong. Feature-aware reconstruction of volume data via trivariate splines [C]// Pacific Graphics 2011. Kaohsiung：[s.n.]， 2011.第40卷第3期河北科技大學(xué)學(xué)報(bào)Vol.40，No.3

2019年6月Journal of Hebei University of Science and TechnologyJune 2019