籍 冉 冉， 鄭 曉 朋， 雷 娜， 羅 鐘 鉉

( 大連理工大學(xué) 軟件學(xué)院, 遼寧 大連 116620 )

0 引 言

數(shù)字幾何用于描述空間中的幾何物體，在工業(yè)設(shè)計、影視、游戲等領(lǐng)域都有廣泛的應(yīng)用.數(shù)字幾何首先需要在計算機中對空間幾何模型進行描述，從而在此基礎(chǔ)上進一步進行建模、分析、仿真等.三角形網(wǎng)格是表述中最常用的一種網(wǎng)格形式，生成技術(shù)也比較成熟，但在許多應(yīng)用領(lǐng)域，四邊形網(wǎng)格卻具有無可替代的優(yōu)勢.在數(shù)值模擬仿真中，四邊形網(wǎng)格的穩(wěn)定性、精度、收斂速度等都優(yōu)于三角形網(wǎng)格；四邊形網(wǎng)格也可以用于生成樣條曲面；四邊形網(wǎng)格在模型參數(shù)化、紋理貼圖等領(lǐng)域也有著廣泛的應(yīng)用.在四邊形網(wǎng)格生成中保持模型的特征至關(guān)重要.

曲面四邊形網(wǎng)格生成的方法很多，主要有基于方向場的方法[1-2]、參數(shù)化方法[3-4]、“分而治之”的方法[5-6]和基于Morse理論的方法.Dong等[7]基于Morse理論，分析拉普拉斯矩陣的特征向量得到四邊形網(wǎng)格，算法首先計算輸入網(wǎng)格的拉普拉斯矩陣，求解拉普拉斯矩陣的特征向量，選擇合適的特征向量作為特征函數(shù)；然后利用該函數(shù)生成Morse-Smale復(fù)形完成對三角網(wǎng)格的四邊單元剖分；最后在此剖分的基礎(chǔ)上進行全局分片參數(shù)化，并在參數(shù)域上進行規(guī)則采樣獲得最終的四邊形網(wǎng)格.在之后的工作中，Huang等[8]、Zhang等[9]和Ling等[10]加入了方向、特征對齊和各向異性的約束，通過優(yōu)化求解帶有約束的廣義特征問題來求解特征函數(shù).Huang等[8]將譜分析與離散余弦變換的關(guān)系應(yīng)用在特征函數(shù)的優(yōu)化上，建立了特征函數(shù)與其對應(yīng)Morse復(fù)形的內(nèi)在聯(lián)系，通過改變拉普拉斯算子中的面積項來控制單元密度；通過在廣義拉普拉斯特征問題中引入額外的能量項來控制Morse復(fù)形邊的方向與位置，方向由用戶輸入的方向場控制，特征對齊對應(yīng)于網(wǎng)格邊的具體位置.Zhang等[9]提出一種基于駐波構(gòu)造的各向異性四邊形網(wǎng)格生成技術(shù)，通過挖掘四邊形網(wǎng)格和駐波的內(nèi)在關(guān)系，將所有的用戶需求同時轉(zhuǎn)化為駐波的各種物理屬性約束，通過輸入密度場實現(xiàn)了各向異性和密度控制，同樣也通過輸入方向場來實現(xiàn)方向和特征對齊約束.Ling等[10]通過引入邊界條件來實現(xiàn)特征對齊控制，通過求解Helmh?ltz方程來得到特征函數(shù)，在Helmh?ltz方程中加入了密度場的控制.這些算法均需將所有用戶控制和質(zhì)量需求集中在一個全局優(yōu)化問題中統(tǒng)一求解，算法復(fù)雜，隨著模型規(guī)模的增大，計算代價將急速上升，噪聲增多，也需要后期復(fù)雜的優(yōu)化.2018年，F(xiàn)ang等[11]將參數(shù)化和Morse理論結(jié)合起來，算法通過用戶輸入標(biāo)架場來進行參數(shù)化，在參數(shù)化的過程中也輸入特征線來進行特征約束，在參數(shù)化發(fā)生退化的區(qū)域用Morse理論來生成四邊形網(wǎng)格，通過將兩種方法優(yōu)勢互補，解決了參數(shù)化的退化問題和Morse方法在規(guī)模較大時運算量大的缺陷.

本文將曲率信息加入到拉普拉斯矩陣中，通過迭代的方法計算特征函數(shù)，在迭代的過程中將特征線的信息加入，使得臨界點能夠準(zhǔn)確定位到特征線上，并通過臨界點來構(gòu)造保特征的Morse-Smale復(fù)形，進而參數(shù)化生成四邊形網(wǎng)格.

1 技術(shù)路線

基于Morse理論生成滿足特征的四邊形網(wǎng)格，算法步驟概括為5步，圖1是原始算法與本文算法的對比，本文接下來將對每一步進行更加詳細的介紹.

(a) 原始算法步驟

(b) 本文算法步驟

2 Morse理論與Morse-Smale復(fù)形

2.1 Morse理論

(a) 極小值點

(b) 鞍點

(c) 極大值點

Banchoff[13]用分段線性函數(shù)替代光滑函數(shù)f將Morse理論推廣到三角網(wǎng)格上，f的值定義在三角網(wǎng)格頂點上，通過判斷頂點的一環(huán)領(lǐng)域內(nèi)的點來判斷臨界點.三角網(wǎng)格上光滑的Morse函數(shù)通常是求解一個拉普拉斯方程得到，即

(1)

式中：Ni為頂點i的一環(huán)領(lǐng)域的所有頂點的集合，wij表示與頂點i相連接的邊(i，j)的權(quán)重，在三維流形曲面網(wǎng)格上，通過用拉普拉斯矩陣的形式來表示拉普拉斯算子，即

(2)

式中：αij、βij為邊(i，j)的兩個三角形的兩個對角.-Lf=λf，特征函數(shù)即為計算特征值所對應(yīng)的特征向量，特征值0=λ1≤λ2≤…≤λn所對應(yīng)的特征向量e1，e2，…，en，隨著特征值的增大，頻率越來越高，臨界點的個數(shù)越來越多，在模型上的可視化如圖3所示.拉普拉斯方程的解稱為調(diào)和函數(shù)，因此拉普拉斯矩陣的特征向量也是離散情形下球面調(diào)和函數(shù)的基.

2.2 Morse-Smale 復(fù)形

Morse函數(shù)f，如果它的穩(wěn)定流形和不穩(wěn)定流形之間只有橫向相交，則f是Morse-Smale函數(shù)[12]，穩(wěn)定流形與不穩(wěn)定流形相交形成的胞腔稱為Morse-Smale復(fù)形.Morse函數(shù)f定義在三角

網(wǎng)格上，f的值定義在三角網(wǎng)格頂點上，通過比較頂點的值與一環(huán)領(lǐng)域內(nèi)點的值來判斷臨界點的類型，如圖4所示.

Morse-Smale復(fù)形是在一個流形上標(biāo)量函數(shù)的胞腔剖分.如圖5所示，從鞍點出發(fā)沿梯度上升或者下降最快的方向，直到達到極大值點或者極小值點，一般來說會從鞍點引出兩條上升梯度線和兩條下降梯度線，這些路徑把流形分成若干四邊形區(qū)域，所有的四邊形區(qū)域都有兩個鞍點、一個極大值點、一個極小值點.

Morse-Smale復(fù)形的生成直接關(guān)乎到最終生成的四邊形網(wǎng)格的質(zhì)量，而生成Morse-Smale復(fù)形的關(guān)鍵在于求解Morse函數(shù)，在計算Morse函數(shù)時通常會充分考慮特征對齊約束的要求，但由于離散化中的精度損失以及噪聲點的影響，臨界點不能完全準(zhǔn)確地定位在特征線上，導(dǎo)致Morse-Smale復(fù)形中的邊沒能與特征對齊，因此需要對臨界點進行局部擾動，使其準(zhǔn)確定位到特征線上.本文針對這一點提出了相應(yīng)的改進算法，下面將詳細說明.

3 加入特征約束的Morse-Smale復(fù)形

生成滿足模型特征的網(wǎng)格在曲面四邊形化中是非常重要的一個約束要求.而生成滿足模型特征的Morse-Smale復(fù)形，關(guān)鍵在于求解拉普拉斯矩陣的特征函數(shù)，即Morse函數(shù).本文給出如何求解特征函數(shù)，使得臨界點準(zhǔn)確地定位在特征線上.

3.1 加入曲率信息的拉普拉斯矩陣

在計算三角網(wǎng)格上標(biāo)量函數(shù)時，通過拉普拉斯矩陣的某個特征值的特征向量來定義，即

-Lf=λf

(3)

通過借鑒之前的研究[8，14]，對拉普拉斯矩陣的對角線進行擾動.在拉普拉斯方程中加入了模型曲率的信息，使用三角網(wǎng)格曲面模型頂點上的高斯曲率，并歸一化：

(4)

其中θi為點v的一環(huán)領(lǐng)域的三角形中所對應(yīng)的頂點為v的角的大小.曲率反映了模型表面的彎曲程度，刻畫了模型的幾何特征信息，求解特征函數(shù)的方程定義如下：

-Lf=λ(I+K)f?-(I+K)-1Lf=λf

(5)

其中矩陣K為對角矩陣，對角線上的元素ki等于第i個頂點vi處的高斯曲率.

拉普拉斯矩陣在加入曲率信息后，臨界點會因為加入的曲率值而進行相應(yīng)的位置擾動，臨界點能夠更加精確地定位到模型特征的地方，生成的Morse-Smale復(fù)形更加貼合模型的特征，如圖6所示.

3.2 特征函數(shù)的迭代算法

通過計算拉普拉斯矩陣L的某個特征值λk的特征向量來作為特征函數(shù)f，在這里一般使用的是ARPACK稀疏矩陣特征問題求解器來求解，求解出來的只能是特定的數(shù)值λk所對應(yīng)的特征函數(shù)，特征函數(shù)的選取不具有靈活性，具體體現(xiàn)在以下兩個方面：

(a) 加入曲率信息之前的結(jié)果

(b) 加入曲率信息之后的結(jié)果

(1)2.1節(jié)分析的臨界點個數(shù)隨著特征值的數(shù)值增大而增多，如果在求解中λk所對應(yīng)的特征函數(shù)fk的臨界點個數(shù)偏少，而λk+1所對應(yīng)的特征函數(shù)fk+1的臨界點個數(shù)偏多，根據(jù)原始的方法不能求得fk和fk+1之間的特征函數(shù)；

(2)想要對特征值λk所對應(yīng)的特征函數(shù)fk做少許的擾動，或者是對局部頂點上的特征函數(shù)做擾動，使得臨界點移動到模型特征的地方，根據(jù)原始的方法也沒有辦法實現(xiàn).

從以上兩點可以看出，特征函數(shù)的選取受到局限，因此本文提出了一種迭代算法來計算特征函數(shù).

(6)

對于任意的頂點vi，式(6)的第i行可以寫為

(7)

這里λ不僅限于特征值，可以是任意數(shù)值，則相對應(yīng)求解出來的f就更具有靈活性，則迭代公式即為式(7).根據(jù)拉普拉斯算子的離散形式(1)，定義全局M上的調(diào)和能量[15]

(8)

當(dāng)?shù)蠼廒呌诜€(wěn)定時，調(diào)和能量逐漸減小并趨于0.想要對特征值λk所對應(yīng)的特征函數(shù)fk做相應(yīng)擾動，就將fk作為迭代初值.算法流程如下：

輸入： 拉普拉斯矩陣L.

輸出： 特征函數(shù)f.

Step 1計算矩陣L的特征值λk所對應(yīng)的特征函數(shù)fk，將fk作為迭代的初值f0.

Step 2計算調(diào)和能量E0.

Step 3對所有網(wǎng)格上的點進行更新：

Step 4計算調(diào)和能量Em+1，當(dāng)|Em+1-Em|<ε時，停止迭代，否則返回Step3.

3.3 加入特征信息

在模型表面生成四邊形網(wǎng)格時，希望生成的網(wǎng)格的邊能與模型的特征線對齊，在基于Morse理論生成四邊形網(wǎng)格的過程中，Morse-Smale復(fù)形的邊一定是四邊形網(wǎng)格的邊，所以生成的Morse-Smale復(fù)形的邊與模型的特征線對齊.在本文算法中，關(guān)鍵是讓Morse-Smale復(fù)形的臨界點準(zhǔn)確定位到模型的特征線上，步驟如下：

(1)首先在模型的表面畫出特征線，在本文中，為了驗證算法的效果，只選取了部分特征線，特征線的集合為Se，定義這些邊的特征屬性作為輸入，如圖7所示

圖7 加入特征線信息

(2)在3.2節(jié)的迭代算法流程中，針對Step 3，如果[vi，vj]∈Se，μij>1，否則μij=1.針對不同的模型μij的取值不同，在本文的模型中μij∈[1，2].特征函數(shù)的迭代公式為

(9)

即增大特征線的權(quán)重.通過增大特征線的權(quán)重，就可以使得原本在特征線附近的臨界點，準(zhǔn)確定位到特征線上，并且采用迭代的算法，對不是特征線上的點沒有產(chǎn)生太大的影響，如圖8給出了增加特征線權(quán)重前后的結(jié)果對比.

圖8 加入特征線的臨界點前后對比圖

從圖中可以看出，臨界點都準(zhǔn)確地定位到了特征線上，并且其他位置的特征點沒有發(fā)生太大的變化.從之前的分析也可以知道，Morse-Smale復(fù)形的邊都是從鞍點出發(fā)，連向極大值點或極小值點，因此，只要臨界點在特征線上，Morse-Smale復(fù)形的邊大多就會沿著特征線，少數(shù)情況下在特征線上存在兩個相鄰的臨界點為極值點，兩個點之間沒有復(fù)形的邊，此時可以通過對偶操作或者補齊操作來優(yōu)化，補齊操作的優(yōu)化過程在下一節(jié)會詳細介紹，優(yōu)化使得生成的復(fù)形的邊可以沿著模型的特征線，如圖9所示，藍線是補齊算法補齊的特征線.

圖9 沿特征線的復(fù)形結(jié)果圖

4 Morse-Smale復(fù)形的優(yōu)化

4.1 噪聲點的消除

初步生成的Morse-Smale復(fù)形會產(chǎn)生很多的噪聲點，如圖10所示.

通過持久性簡化[16]來消除不穩(wěn)定的點對.在生成的Morse-Smale復(fù)形中，對于兩個相鄰的臨界點vi、vj，它們對應(yīng)于特征函數(shù)所對應(yīng)的標(biāo)量場的一組拓撲特征，持久性值定義為兩個點的函數(shù)值之差的絕對值，即|f(vi)-f(vj)|，其是描述移除這個拓撲特征時，特征函數(shù)的變化量，持久性值較小的點對屬于不穩(wěn)定點對，要進行消除.一對相鄰的臨界點，必然一個是鞍點，另一個是極值點，定義在它們上的消除操作是指將鞍點和極值點消除，同時將所有連向該極值點的梯度線連向鞍點的另一個同屬性的極值點.如圖11所示.

(a) 存在大量噪聲的結(jié)果圖

(b) 消除噪聲的結(jié)果圖

圖11 持久性簡化操作

4.2 幾何優(yōu)化

希望生成的Morse-Smale復(fù)形是近似于正方形或者矩形，但是由于根據(jù)特征函數(shù)計算梯度線過程中的不可控制性，初始復(fù)形的形狀可能會出現(xiàn)菱形、梯形等，如果不對其優(yōu)化，會對后續(xù)參數(shù)化生成四邊形網(wǎng)格的質(zhì)量產(chǎn)生影響，因此對全局的復(fù)形進行形狀的優(yōu)化.

(1)文獻[14]中提到通過參數(shù)化和松弛迭代的方法來重新定位復(fù)形的邊緣點，文獻[8]中由于加入了特征線，為了不使復(fù)形的邊緣偏離特征線，對文獻[14]的方法進行了擴充，固定特征線上的參數(shù)化值，借鑒文獻[8]的算法，可以對復(fù)形的邊緣進行優(yōu)化；

(2)本文迭代算法中，需要對特征線上相鄰的兩個極值點之間進行補齊操作，由此一個復(fù)形區(qū)域被分割成了兩個三角形，針對因為這種情況生成的三角形區(qū)域，可以對這些區(qū)域采用“分而治之”的思想，將1個三角形剖分成3個四邊形.

5 四邊形網(wǎng)格生成

優(yōu)化完成后的Morse-Smale復(fù)形可以保證生成高質(zhì)量的四邊形網(wǎng)格.將每一個復(fù)形區(qū)域參數(shù)化到平面上，這一步可以通過調(diào)和映射完成[15]，根據(jù)密度要求d，在二維平面上生成規(guī)則的d×d維的四邊形網(wǎng)格，再將其映回到模型表面，輸出四邊形網(wǎng)格.

本文給出了幾種模型迭代結(jié)果，如圖12(c)，臨界點可以很準(zhǔn)確地與特征線對齊，通過補齊操作可以實現(xiàn)最終的保特征四邊形網(wǎng)格生成；圖12(d)給出了補齊操作的結(jié)果.幾何優(yōu)化部分和最后的分區(qū)域參數(shù)化四邊形網(wǎng)格生成都是非常成熟的算法，在本文中沒有顯示，可以在以后的優(yōu)化中實現(xiàn).

圖12 4個模型的結(jié)果

6 結(jié) 語

本文在Morse理論的基礎(chǔ)上，提出了一種特征函數(shù)的新計算方法.通過迭代算法計算特征函數(shù)，使得特征函數(shù)的選取更加具有靈活性；在迭代算法中加入特征約束，可以將臨界點準(zhǔn)確定位在特征線上，生成Morse-Smale復(fù)形；最后對初始形成的Morse-Smale復(fù)形進行拓撲優(yōu)化和幾何優(yōu)化，進而生成保特征的四邊形網(wǎng)格.本文所展示的結(jié)果與文獻[8-11]的結(jié)果相比，存在不足，但是本文提出了計算特征函數(shù)的一種特殊的新算法，并且這個算法實現(xiàn)了臨界點對齊這一好的結(jié)果，對后續(xù)的保特征的四邊形網(wǎng)格生成起到了關(guān)鍵作用，在優(yōu)化完成后也可以生成質(zhì)量較好的四邊形網(wǎng)格，并且本文提出的算法簡單，易于實現(xiàn)，輸入信息較少，計算效率高.

后續(xù)的研究工作包括以下內(nèi)容：(1)實現(xiàn)對目前四邊形網(wǎng)格的優(yōu)化，并將本文算法應(yīng)用到實際問題中，驗證其有效性；(2)將本文的算法推廣到高維，生成三維的Morse-Smale復(fù)形，進而生成六面體網(wǎng)格.