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

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

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

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

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

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

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

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

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

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

Fig.1Scalar and vector on triangle

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

進一步地，可求解其對應(yīng)的梯度場（矢量場）。對于分段線性流形，在每個分段內(nèi)梯度矢量為常量。如圖1所示，已知三角形T0（v0，v1，v2）或（e0，e1，e2） ，對應(yīng)的標(biāo)量參數(shù)u0，u1，u2，三角形面的單位法矢為N，則T0的梯度G0可通過以下線性系統(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)，可直接通過手工計算求解參數(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，對應(yīng)地，可求正交梯度：G⊥=n×G0 。 ???????????????????????????????????（9）1.2數(shù)據(jù)準(zhǔn)備

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

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

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

struct NodeDetail

{

long v[2];//記錄對邊兩點的序號

long idxF;//記錄三角形的序號

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

};

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

2）節(jié)點分類

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

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

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

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

struct LaplaceData

{

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

float data;//cotangent計算值

};

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

{

LaplaceData* e11 = （LaplaceData*）e1;

LaplaceData* e22 = （LaplaceData*）e2;

return e11->idx - e22->idx;

}

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

2流線跟蹤

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

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

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

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

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

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

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

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

2.1初始梯度線的計算

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

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

t0=-ptSeed.y/ptSeed.x 。

相應(yīng)地：

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

接下來，驗證x0的取值范圍，如果滿足0≤x0≤v1.x，則求得一個交點（x0，0） ，否則直線L與三角形邊v0v1無交點。

如果有交點，則需要如下的后置處理：

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

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

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

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

該問題可描述如下：已知三角形T0（v0，v1，v2），T0內(nèi)的梯度矢量為G0，其中一邊（假設(shè)為v0v1）上有一點ptSeed，過ptSeed且方向為G0的直線L與T0的另外兩邊有一個交點ptInt，求該交點，如圖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)立求解，得交點坐標(biāo)：

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

接下來，判斷交點有效性。構(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，則說明交點（x0，y0）在線段v1v2之內(nèi)，為有效交點，否則認(rèn)為直線L與v1v2邊無交點。

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

2）L與v2v0求交

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

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

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

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

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

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

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

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

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

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

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

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

序號條件交點1t>0， s ≤ 0v12t>0， s ≥ 1v23t>0， 0 < s < 1ptInt4t0≤0無交點

如果出現(xiàn)多個有效節(jié)點，則需要進行沖突消解處理。可根據(jù)實驗效果，實施以下2種消解策略。

1）最短距離策略

對所有有效交點，計算至v0的距離，將最短距離的交點作為輸出。

2）s極值策略

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

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

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

1）分段累計種子線lnSeed的長度，等分lnSeed為nu份，得到（nu-1）個等分點ptSegi，i=0，1，…，nu-2;

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

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

4）等分lnSeed⊥為nv份，再從這些等分點出發(fā)沿G向追蹤，得到（nv-1）條G向流線。

4實驗

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

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

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

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

圖9使用人體網(wǎng)格模型展示了不同種子點對流線網(wǎng)格的影響。圖9 a）和圖9 c）分別在人體正面和背面選定不同的種子點ptSeed，可分別在人體正面和背面追蹤種子線;圖9 b）和圖9 d）為2種情況下的流線網(wǎng)格，對比來看流線網(wǎng)格并未發(fā)生變化。理論上，只要極點vmin，vmax一定，所求的梯度場G及正交梯度場G⊥就是確定的，對應(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è)髁部位（圖示右上部）沒有G向流線，圖10 d）中，在網(wǎng)格加密后，同一部位劃分出了G向流線，對該部位及整體模型的表達也更準(zhǔn)確。

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

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

此外，借助稀疏矩陣方案，以上實驗?zāi)Ｐ偷奶荻葓銮蠼饩鶝]有時間和空間壓力;而流線追蹤過程主要是利用線性方程求解交點，也不會存在性能問題。

5結(jié)語

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

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

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

參考文獻/References：

[1]汪攀， 張見明， 韓磊， 等. 基于帶約束前沿推進的四邊形網(wǎng)格生成方法[J]. 湖南大學(xué)學(xué)報（自然科學(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]. 機械設(shè)計與制造， 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]李天華. 基于方向場的四邊形重網(wǎng)格化技術(shù)研究[D]. 長春：吉林大學(xué)， 2017.

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

[4]徐崗， 舒來新， 朱亞光， 等. 邊界簡化與多目標(biāo)優(yōu)化相結(jié)合的高質(zhì)量四邊形網(wǎng)格生成[J]. 中國圖象圖形學(xué)報， 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)柏良， 林淑金， 等. 針對密集點云的快速自適應(yīng)四邊形網(wǎng)格生成算法[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報， 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]袁潔， 周明全， 耿國華， 等. ?基于Morse-Smale拓?fù)涮卣鞯奈奈锼槠唇铀惴╗J]. 自動化學(xué)報， 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é)報， 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ì)國良， 丁勇， 周曼， 等. 工程計算中大型稀疏矩陣存儲方法研究[J]. 數(shù)值計算與計算機應(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]張永杰， 孫秦. 稀疏矩陣存儲技術(shù)[J]. 長春理工大學(xué)學(xué)報， 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]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報， 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é)報Vol.40，No.3

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