基于離散流形調和基的保特征譜幾何處理方法

秦紅星，洪小洋，孫 穎

（1.重慶郵電大學 計算智能重慶市重點實驗室，重慶400065；2.重慶郵電大學 計算機科學與技術學院，重慶400065）

0 引 言

由于傳統的利用拉普拉斯算子進行幾何處理的方法需要點之間的幾何連接信息，所以無法直接被應用到離散點云上。為了解決這個問題，Belkin等［1］提出了一種計算點云曲面的離散拉普拉斯算子的方法，并且通過理論驗證了該離散拉普拉斯算子的收斂性。研究結果表明，該離散拉普拉斯算子隨著點云之間的距離的靠近 （采樣點數目的增加），會變得更加精確，且收斂的速度也會變快。然而，Belkin的方法構造出的拉普拉斯算子無法保持模型的局部尖銳特征，不能較好表達模型的完整信息，不利于后期的幾何處理。

本文在Belkin等［1］的拉普拉斯算子近似方法的基礎上，利用離散化定義在流形上的函數積分的方式來構造拉普拉斯算子。在此基礎上，通過考慮局部切空間的Voronoi圖以及定義空間內積的方式確保構造出的拉普拉斯算子具有可對稱化性與收斂性。不同于文獻 ［2］的方法，這里在做切空間投影構造局部Voronoi圖時，引入法向約束，使得切空間的上的Voronoi圖變為一個各向異性Voronoi圖，這樣得到的Voronoi圖的面積比直接在切空間上構造要更加精確，可以較好地保存點云模型的局部特征。

利用構造出的拉普拉斯算子，本文做了包括幾何濾波、骨架提取等一系列實驗。實驗結果表明，這種方法在模型的局部特征保持上有明顯效果，且算法具有高效、穩定的特點。

1 相關研究

譜幾何處理［3］作為一種精確有效的幾何處理方法，在過去的十幾年間發展很快，目前已被應用到模型匹配與檢索［4］等多個計算機圖形學領域中。作為譜幾何處理的一個重要組成部分，關于拉普拉斯算子研究一直是幾何處理研究的一大熱點。Vallet［5］為了更有效的研究Laplace算子，引入了離散外積分的計算思想，并且在此基礎上將拉普拉斯算子的特征函數定義為Manifold Harmonics。

算子估計過程如圖1所示。

圖1 算子估計過程

目前，尋找到一種有效的離散化網格曲面的拉普拉斯算子方法是幾何處理的一個重要問題［6］。大多數已有的方法都是應用有限元的相關知識，結合不同的估計條件實現對拉普拉斯算子的離散化，事實上這些方法都可以歸結為cotangent類型近似。Reuter等［7］提出的基于3次有限元的離散拉普拉斯算子在處理連續的問題的效果很好，且證明其可以很好地應用到形狀理解中。Peinecke等［8］提出將特征值譜應用到圖像識別的指紋識別中。事實上，這些方法得到的拉普拉斯算子本質上都是利用三角網格的Voronoi圖面積來近似曲面流形的Voronoi圖面積，但是由于三角網格本身存在不規整性等問題，所以估計出的面積并不足夠精確。

在流形上，拉普拉斯算子和熱擴散方程具有非常緊密的聯系。Belkin等利用分段收斂條件來近似網格曲面的拉普拉斯算子，證明了在估計的高斯核曲面的切平面上估計拉普拉斯算子這一方法的可行性，進一步擴展了二者的聯系，但是并沒有給出算法的收斂性證明，這使得該方法沒有良好的理論支持。Dey等［9］在充分分析Belkin的算子構造方法后，利用矩陣擾動理論證明了該方法的特征值的收斂性。

后來這種方法被擴展應用到離散化點云表達的流形曲面中，目前已經證明這種拉普拉斯算子會隨著點云的密集化逐漸分段收斂。然而這種方法構造出的拉普拉斯算子是無法保證可對稱化性質的，所以不可以用來計算曲面的正交基。

由于缺少幾何連接信息，在點云中計算拉普拉斯算子一般都是采用近鄰處理的方式，即構造圖拉普拉斯算子，但是圖拉普拉斯算子存在著對幾何信息不敏感的問題，無法較好表達原模型的幾何信息。由于拉普拉斯算子是一個二階微分算子，所以在點云上計算拉普拉斯算子與點云上的積分計算有著緊密的聯系。Luo等［10］引入基于測地線距離的Voronoi圖，利用測地線Voronoi圖面積近似曲面流形上的積分計算。本文與Luo等的方法不同是，這里采用的距離度量為歐氏距離度量，利用帶法向約束的切空間上的各向異性Voronoi圖的面積近似流形上Voronoi圖的面積，最終利用該Voronoi圖近似拉普拉斯算子的積分函數。這樣得到的拉普拉斯算子，相比Cotangent型算子利用三角網格的Voronoi圖面積近似流形Voronoi圖面積，以及Liu等［2］的切向面積近似方法，由于更加逼近測地線距離，所以得到Voronoi圖面積更加精確，逼近于準確的Voronoi圖面積。

2 拉普拉斯算子

2.1 定義

拉普拉斯算子是一個簡單的二階微分算子，傳統的拉普拉斯算子可以定義為維歐氏空間的梯度的散度

與歐氏空間相似，具有緊致邊界條件n 維的流形上的拉普拉斯算子同樣可以定義為

2.2 特征值問題

由上述定義及矩陣論的相關知識，拉普拉斯算子的特征值問題可以定義為

式中：λ和H ——拉普拉斯算子的特征值與該特征值對應的特征函數 （離散的形式中特征函數就是特征向量）。這里的負號用來保證所有的特征值λ大于0。這里用λi表示拉普拉斯算子的第i個特征值，用Hi表示第i 個特征值對應特征函數。

2.3 對稱性

具有緊致邊界條件的流形，其對應的拉普拉斯算子滿足對稱性，即其內積滿足

式中：f，g——定義在n維流形M 上的函數，Δ——梯度子。由此可得，拉普拉斯算子的不同特征值對應的特征函數是相互正交的。

2.4 流形調和基

通過解特征值問題 （1），可以得到定義在流形上的特征值集合｛λi｝以及其對應的特征向量集合｛Hi｝。通過特征值｛λi｝，可以獲得流形的一些基本數據，例如體積、邊界長度以及虧格數等，另一方面，｛Hi｝作為拉普拉斯算子對應的特征向量，其反映了曲面的本征屬性，可以用來特征提取、曲面分析等，其離散形式稱為流形調和基（MHB）。

2.4.1 構造離散拉普拉斯算子

為了獲得MHB，需要定義離散的拉普拉斯算子用來近似流形的拉普拉斯算子。為了避免直接離散化拉普拉斯算子，這里采用文獻 ［2］離散化函數積分的方式構造拉普拉斯算子。這種方法需要引入以下幾個定義：

定義1 對任意ε＞0，如果流形M 上的一個有限采樣集合S（S M）滿足

則稱S 為流形M 的一個ε采樣。

如果該采樣除了滿足ε采樣條件外還滿足

則稱S 為流形M 的一個（ε，ζ）采樣。本文的給定的點云模型P 都默認為流形M 的一個（ε，ζ）采樣。

定義2 如果存在一個球B 除了流形M 的兩個邊界點之外，不包含任何其它位于M 內的點，則稱球B 為M 的一個中球，稱M 的所有的中球的球心組成的集合為中軸。

定義3 每個采樣自流形M 的點集S，S 中的任意點s的Voronoi圖可以定義為流形M 的子集，即

這里的 · 表示任意兩點之間的歐氏距離。

定義4 設ρ為流形M 上所有點到其中軸的距離函數，則稱該距離函數ρ為局部特征尺度。

2.4.2 計算近似的拉普拉斯算子ΔM

在2.4.1中已經介紹這里通過離散化積分的方法來近似拉普拉斯算子，該積分可以表示為

其中，f 為定義在流形上的函數。為了構造一個較好的對拉普拉斯算子離散化近似，所以需要對每個點單獨進行運算，這里我們引入Voronoi圖面積來離散化近似上述積分，近似的函數可以通過此公式表達

具體的計算過程為：

（1）估計點s處的切平面：對模型進行ε采樣，得到采樣點集合S，取S 中一采樣點s，考慮距離s長度為r 的點集Sr，即Sr＝S ∩B（s，r），B（s，r）為中心位于點s半徑為r 的球。設P＊為經過點s的最適應的平面，故d（Sr，P＊）可以取到最小值。利用Har－Peled和Varadarajan的方法［11］，可以構造一個P＊的二階近似，是一個滿足dH（Sr）≤2dH（Sr，P＊）經過點s的平面，這里的dH（·，·）是兩點間的Hausdorff距離。

（2）估計點S 的Voronoi圖：將Sr中的點投影到切平面上，當r足夠小的時候，投影映射是一個雙射，記為，然后我們在切平面上構建Voronoi圖，將切平面上的Voronoi圖（s）面積作為流形上的Voronoi圖VrM（s）的面積的一個近似，為了保證估計的Voronoi圖更加精確，這里考慮利用法向約束切空間的投影過程，法向約束可以描述為

其中，Voronoi為約束后的Voronoi圖，Voronoioriginal為約束前的Voronoi圖。隨著采樣點變得密集，這種近似會變得更加精確，后者收斂到前者的效率變得更高。

2.4.3 離散化拉普拉斯算子

由于要構造的離散拉普拉斯算子具有可對稱性，故根據空間內積的一般定義

式中：S——實正定對稱矩陣，那么構造需要的目標就是

這里的L是離散拉普拉斯算子的矩陣形式，那么問題就轉換為確保SL＝LTS＝（SL）T為對稱矩陣，令R＝SL，則可得拉普拉斯算子對應的矩陣L為S－1R，即若L滿足條件L＝S－1R，那么L 一定是對稱矩陣，這樣就可以確保生成的調和基是正交的。我們采用方程 （4）來構造離散拉普拉斯算子，那么上述的情形就可重新變形為：＝，其中是一個N 維向量，f ＝ ［f（s1），f（s2），…f（sN）］T表示在采樣點處的連續函數的f 對應的函數值，N 為采樣點的個數。所以我們可以得到定義在模型上的拉普拉斯算子的矩陣表達形式為

由于上述矩陣S為對角矩陣，R 為對稱矩陣，綜合式（5）的內積定義，可以證明得到的矩陣是一個可對稱化的矩陣。

2.4.4 譜濾波器

MHT 和逆MHT：MHT 是一種將函數投影到MHB下的變換，通過找到適當的系數來確定最終的變換，且是一種正交變換。對應于MHT，逆MHT 是將經過變換的對象變回原對象的變換。

利用矩陣論的相關知識，MHT 可以定義為如下形式

式中：f——采樣點的函數值向量，Hi——矩陣的第i個特征向量。同理，逆MHT 可以定義為

利用上述變換與信號處理中的相關方法可以實現對原有對象的變形，濾波等操作。

設F（ωk）為一個頻率空間的濾波器，那么濾波后的函數值就變為了

即對譜坐標應用濾波器，得到新的模型，最終通過逆MHT 將譜坐標變換回原有的坐標，從而獲得最終的模型，這樣也得到了濾波后的模型，可以通過改變F（·）來獲取不同的濾波效果。

3 算法對比與分析

如圖2所示，文獻 ［2］的方法核心就是直接用垂直投影在切空間上的兩點距離來近似曲面流形的上測地線距離，實際上當角度θ變得很大的時候，這種近似精度急速下降，導致近似出的距離有相對大的偏差。本文采用的近似通過融合切向距離和法向距離兩個方面，由于這兩個角度本身近似互余的關系，所以當角度θ變得很大的時候，角度φ變小，則后一項變大，整體距離更加趨近于實際距離，所以在特征尖銳的部分保持局部特征具有很好的效果。反之，當角度θ變小時，角度φ變大，則后一項變小，這樣使得特征平滑的部分變得更加平滑。由誤差分析圖可知本文的算法的誤差比文獻 ［6］的方法的誤差要小很多，在特征保持方面效果更好。

圖2 切空間投影

4 實驗

本文算法采用C＋＋與Matlab 實現，實驗結果利用OpenGL渲染顯示。為了方便與其它方法進行對比，這里使用幾個常用的三維點云模型 （Kitty以及Lion等）。所有實驗在一臺Intel CORE 2Duo 2.00GHz CPU，2GB DDR2 RAM 的Windows 7平臺上完成。以下幾組實驗分別證明了本文算子的有效性及較已有方法在性能上的可靠性。為了證明本文拉普拉斯算子的有效性，這里利用譜幾何處理中的基于拉普拉斯算子的濾波方法以及基于拉普拉斯算子的特征向量重建模型的方法，結合本文的拉普拉斯算子做了幾組實驗。其中，圖3 為譜濾波實驗結果，圖3 （a）為Rocker－arm，Kitty以及Lion的原始模型，圖3 （b）為應用文獻 ［2］的拉普拉斯算子對Rocker－arm，Kitty 以及Lion等模型的濾波結果，圖3 （c）為應用本文的拉普拉斯算子對Rocker－arm，Kitty 以及Lion 等模型的濾波結果。從Rocker的尾部，Kitty的眼睛和尾巴，Lion的脖子和鈴鐺可以發現，本文算子可以有效地應用于譜濾波，且在局部尖銳特征保持方法具有更好的效率，這直接證明了近似算子的準確性和有效性。

圖3 譜濾波結果

為了證明本文方法構造出的算子在性能方面的可靠性，我們給出3個模型的多分辨率重建結果。通過實驗來分析和突出本文算子在局部特征的保持方面的優勢。本文在求解拉普拉斯算子特征值和特征向量的基礎上，對特征值進行從小到大排序，再在此基礎上依次挑選特征向量。圖4為利用近似的拉普拉斯算子對Vase－lion、Armadillo 以及Julius等模型的結果。其中圖4 （a）為原始模型，圖4 （b）為利用前1%的特征向量重建的結果，圖4 （c）為利用前10%的特征向量重建的結果，圖4 （d）為50%的特征向量重建的結果，圖4 （e）為100%的特征向量重建的結果。由Vase－lion底部和嘴巴、Armadillo的肘部和腿部以及Julius的眼睛可以看出本文的方法可以較好重建原始模型尖銳特征，從側面顯示本文的方法構造出的拉普拉斯保存了更多局部特征信息。

如圖5所示本文方法在Maxplanck模型的眼部、嘴部，Block的邊緣部分以及Hand 的掌紋保護方面效果較文獻［2］的方法更加顯著，說明本文的近似的拉普拉斯算子更加逼近于原始流形的拉普拉斯算子。為了證明本文算法的精確性，這里對比文獻 ［2］的方法與本文方法的最大重建誤差，即用本文的利用各向異性Voronoi圖的構造出的拉普拉斯算子的特征向量和文獻 ［2］的利用切空間Voronoi圖構造出的拉普拉斯算子特征向量，對原始模型進行重建。將重建出的模型與原始模型在歸一化空間中進行誤差分析。圖6為6個模型的重建最大誤差圖，可以清楚的看到本文方法重建結果的最大誤差比文獻 ［2］方法重建出的最大誤差更加小，效果更加好。

綜合以上實驗結果，可以發現本文估計出的拉普拉斯算子更加逼近于流形上的拉普拉斯算子，可以很好的保持原始模型的尖銳特征信息。

圖4 多分辨率重建結果

圖5 重建實驗對比

5 結束語

本文提出了一種有效的計算拉普拉斯算子的近似方法，在切空間投影生成Voronoi圖時引入法向約束，這樣在切空間中生成的各向異性Voronoi圖相對于直接在切空間中構造Voronoi圖在近似精度上更加逼近于曲面流形上的Voronoi圖。在此基礎上利用譜幾何處理的相關知識，對幾何模型進行高通、低通濾波等操作。本文實驗通過模型圖以及平均曲率圖很好地顯示了算法在保持特征方法的優異，且算法易于實現，本文貢獻主要有以下幾點：

（1）提出了一種新的近似離散拉普拉斯算子的方法；

（2）將其直接應用到點云中，實現離散數字幾何處理包括骨架提取、幾何濾波等；

圖6 重建最大誤差

（3）結合譜濾波器，實現對噪聲的保特征處理，改進了傳統方法在特征保持等方面的缺陷。

［1］Belkin M，Sun J，Wang Y.Constructing Laplace operator from point clouds in R d ［C］／／Proceedings of the twentieth Annual ACM－SIAM Symposium on Discrete Algorithms.Society for Industrial and Applied Mathematics，2009：1031－1040.

［2］Liu Y，Prabhakaran B，Guo X.Point－based manifold harmonics［J］.IEEE Transactions on Visualization and Computer Graphics，2012，18 （10）：1693－1703.

［3］Zhang H，Van Kaick O，Dyer R.Spectral methods for mesh processing and analysis ［C］／／Proceedings of Eurographics State－of－the－Art Report，2007：1－22.

［4］Hu J，Hua J.Salient spectral geometric features for shape matching and retrieval［J］.The Visual Computer，2009，25（5－7）：667－675.

［5］Vallet B，Lévy B.Spectral geometry processing with manifold harmonics［C］／／Computer Graphics Forum.Blackwell Publishing Ltd，2008，27 （2）：251－260.

［6］Qin H，Yang J，Zhu Y.Nonuniform bilateral filtering for point sets and surface attributes ［J］.The Visual Computer，2008，24 （12）：1067－1074.

［7］Reuter M，Biasotti S，Giorgi D，et al.Discrete Laplace－Beltrami operators for shape analysis and segmentation ［J］.Computers ＆Graphics，2009，33 （3）：381－390.

［8］Peinecke N，Wolter FE，Reuter M.Laplace spectra as fingerprints for image recognition ［J］.Computer－Aided Design，2007，39 （6）：460－476.

［9］Dey TK，Ranjan P，Wang Y.Convergence，stability，and discrete approximation of Laplace spectra ［C］／／Symposium on Discrete Algorithms.Society for Industrial and Applied Mathematics，2010：650－663.

［10］Luo C，Sun J，Wang Y.Integral estimation from point cloud in d－dimensional space：A geometric view ［C］／／Proceedings of the 25th Annual Symposium on Computational Geometry.ACM，2009：116－124.

［11］Har－Peled S，Varadarajan K.Projective clustering in high dimensions using core－sets ［C］／／Proceedings of the 18th Annual Symposium on Computational Geometry.ACM，2002：312－318.