非剛性三維形狀匹配中基于譜分析的形狀描述符綜述?

張 丹 , 武仲科 , 王醒策 , 呂辰雷 , 劉香圓 , 周明全

1(北京師范大學 信息科學與技術學院,北京 100875)

2(北京師范大學 虛擬現實與可視化技術研究所,北京 100875)

非剛性三維形狀匹配是圖形學中的重要問題,是形狀識別[1]、形狀檢索[2]、形狀配準[3]、形狀分割[4]等工作的研究基礎.同時,非剛性形狀匹配也為三維可視化[5]、生物計算[6]、人臉識別[7]、醫學圖像處理[8]等應用領域提供了堅實的理論依據.在上述研究中,非剛性三維形狀檢索與非剛性三維形狀匹配是兩個非常相似卻不相同的研究問題.非剛性三維形狀檢索的主要思想是:首先,將非剛性三維形狀庫中的所有形狀映射到特征空間中,計算所有形狀的特征值并添加索引;其次,根據用戶的需求設置檢索閾值,并選擇合適的相似度計算方法;最后,提取出滿足閾值的形狀,并按照相似度降序輸出形狀[9].而非剛性三維形狀匹配研究的是形狀相似性問題:同樣將待匹配的形狀映射到特征空間中,選擇形狀的局部特征、全局特征或者兩者的結合代替待匹配的形狀;然后選擇某種代價函數或者距離函數度量特征,并將特征之間的度量值作為非剛性三維形狀匹配度.可將其概括為兩個關鍵步驟:(1) 提取形狀上有效的形狀描述符;(2) 選擇合適的相似度度量.

本文綜述了非剛性三維形狀匹配中基于譜分析的形狀描述符.對于剛性三維形狀匹配,目前已有大量的研究成果[10?12],其中,迭代最近點匹配算法(iterative closest point,簡稱ICP)[13]是最常用的三維形狀匹配算法.ICP將形狀上采樣點的空間位置作為形狀描述符,通過多次迭代最小化源形狀和目標形狀采樣點之間的空間距離,實現剛性三維形狀匹配.然而,ICP 采用人為設置的迭代次數作為迭代終止條件,導致算法容易陷入局部最優.此外,對于有拓撲噪聲的形狀,僅用空間位置作為形狀描述符,無法實現非剛性三維形狀的高精度匹配.因此,研究者需提取更加有效的形狀描述符用于三維形狀匹配.形狀描述符是一種描述形狀語義信息和幾何信息的方法,有時候也被稱為某種算子,研究者通過選擇合適的形狀描述符,可以實現非剛性三維形狀的高效匹配.常用的形狀描述符一般包括4 類:基于形狀表面特征的描述符、基于形狀統計特征的描述符、基于形狀拓撲的描述符以及基于譜分析的形狀描述符,本文重點綜述了基于譜分析的形狀描述符.

第1 類形狀描述符致力于描述形狀表面的特征及其在全局歐氏變換下的不變性.尺度不變特征變換(scale invariant feature transform,簡稱SIFT)描述符[14]是其中應用最廣的描述符,SIFT 于1999 年由Lowe 等人提出,被用于檢測和描述圖像中的局部特征,并在2005 年由Mikolajczy 等人證明其具有很強的魯棒性[15].之后,許多研究者也在此研究的基礎上引入隱馬爾可夫模型、核判別分析及測地線圓環等方法對其進行不斷改進,提高了SIFT的實時性及魯棒性[16?18].形狀上下文(shape context)[19]描述了形狀表面上圖像中的線條,同時存儲每個點相對其他點位置的分布,并給出形狀上點的局部上下文信息,在一定的圖像區域內將點云分布轉化為二維自旋圖,執行三維形狀的表面匹配.為了分析有特殊鉸鏈和關節的三維分子形狀,Feng 等人提出了一種基于節點感知的三維形狀描述符[20],該描述符由形狀邊界上任意點的局部形狀半徑變化的信息編碼定義,用于描述有關節的三維形狀.積分不變量(integral invariant)[21,22]描述符通過對稱組對原始形狀進行重建,用積分不變量作為形狀的特征,該描述符可作為定義形狀之間其他距離的基礎.梯度方向直方圖(mesh HOGs)[23]由離散網格上頂點的幾何特征定義,例如曲率、測地線積分等,該描述符描述了形狀上的紋理特征而非幾何特征.與此類似的還有Tombari等人提出的一種局部描述符——CSHOT 描述符[24],該描述符通過匹配形狀的特征點獲得點對點的對應,主要用于三維形狀表面匹配、目標識別等.

第2 類形狀描述符基于對形狀統計特征的描述,主要描述形狀的全局屬性.Osada 等人[25]在2002 年提出了形狀分布(shape distribution,簡稱SD)描述符,其算法步驟為:(1) 在形狀上選擇合適的歐式度量函數,如D1 距離(測量形狀上某個中心點到其他任意點的距離)、D2 距離(測量形狀上任意兩點間的距離)以及D3 距離(測量形狀上任意3 點組成區域面積的平方根)等,圖1 為Osada 的文章中提到的5 種距離;(2) 計算形狀上所有采樣點對間度量函數的分布直方圖;(3) 將該直方圖作為形狀的線性形狀描述符,并用于形狀分析.該方法基于統計分析思想、易于理解且具有很強的普適性,然而SD 方法中的提到的5 種距離只適合描述剛性物體的屬性,當物體發生非剛性變化時,例如等距變換,其值會隨之變化.等距變換是指形狀保持曲面上任意兩點間曲線長度不變的變換,例如一個人彎曲胳膊后,胳膊上任意兩點長度不變,形狀發生等距變換后的不變性稱為等距不變性.基于上述研究,Mahmoudi 等人[26]使用基于D1 距離改進的測地距離分布直方圖作為新的形狀描述符.測地距離是歐式空間中的直線段在黎曼流形上的推廣,具有等距不變性,測地距離定義了曲面上兩點間的最短距離[27].測地距離不僅具有局部最短性,還含有其他豐富的幾何信息,但當兩點間的曲面部分發生缺失或者有縫隙時,測地距離會因為聯通路徑不能通過而發生改變,對拓撲變化魯棒性不足.此后,其他的工作也對此進行了改進,但始終無法克服測地距離對拓撲變化的敏感性[28].基于形狀統計特征的形狀描述符繼承了統計學中統計量的穩定性,在描述性狀特征時魯棒性較高,很適用于分子形狀比較(molecular shape comparison,簡稱MSC)或者三維關節變形形狀.Liu 等人使用內部距離形狀簽名(inner distance shape signature,簡稱IDSS)描述了三維分子形狀,并計算了IDSS 直方圖之間的度量作為三維分子形狀間的相似度[29].在此基礎上,Liu 等人基于可見性圖提出一種新的內積距離計算方法,計算了有關節的三維體模型的內部距離,其對關節變形形狀發生拓撲變化魯棒,能夠很好地描述三維關節變形形狀[30].

Fig.1 Five distances defined in shape distribution圖1 形狀分布中定義的5 種距離

第3 類描述符基于對形狀拓撲結構的特征提取,該類描述符將三維形狀匹配問題轉換成其拓撲結構匹配問題,應用兩個形狀拓撲結構的匹配結果作為形狀的匹配結果.三維形狀的拓撲結構精確地描述了形狀的全局和局部幾何形態特征,并且保留了形狀的層次結構.具有代表性的兩類描述符分別是基于Reeb 圖理論的描述符和基于形狀的骨架線理論的描述符,圖2 為兩個不同形狀的Reeb 圖和骨架線圖示.

Fig.2 Diagram of Reeb graph and skeleton with different shapes圖2 不同形狀Reeb 圖與骨架線圖示

Reeb 在1946 年基于形狀的拓撲結構提出了Reeb 圖的概念,其具體步驟為:首先,在三維形狀的頂點上定義連續光滑函數f:M→R;其次,根據形狀的頂點坐標計算頂點處的函數值,并將頂點進行分類;最后,根據頂點分類結果將形狀M映射為Reeb 圖R,函數值相同且位于同一連通區域的頂點在Reeb 圖中映射為一個節點.Reeb 圖能夠很好地刻畫形狀的拓撲結構,且能起到降維的作用.定義合適的函數f是Reeb 圖算法的關鍵,常用的函數f有高度函數和Morse 函數.Shinagawa 等人[31]采用高度函數、權重函數和形狀上孔的數量等先驗知識自動地生成三維形狀的Reeb 圖.Hilaga 等人[32]通過測地距離、測地線定義了Morse 函數,并基于此提出了多分辨率Reeb算法(MRG),Bespalov 等人[33]在此研究上改進了MRG 算法,并用于CAD 模型匹配.Biaso 等人[34]基于Morse 函數,提出了擴展Reeb 圖算法(extend reeb graph,簡稱ERG),該算法刻畫了形狀上臨界點之間的拓撲關系,但該算法對形狀發生拓撲變化時魯棒性較差.Tierny J 等人[35]基于特征點的思想構造Reeb 圖,其通過特征提取算法提取形狀的特征點,并通過圖構造方法生成Reeb 圖,并應用Reeb 圖進行部分形狀檢索.骨架線,也被稱為三維形狀的中軸,是刻畫三維形狀拓撲結構的另一個重要方法,不但能用線段很好地描述模型的結構信息,而且能高效地保存形狀,提高形狀空間存儲率和形狀檢索率.Sundar 等人[36]使用拓撲細化算法提取了形狀的骨架線.Cao 等人[37]提出了一種基于拉普拉斯壓縮的骨架線提取算法,該算法可快速提取點云模型的骨架線.Sfikas 等人[38]基于形狀的拓撲信息和共形幾何特征,提出了一種非剛性三維形狀檢索方法,該算法具有穩健又高效的檢索精度和計算速度.

以上3 類形狀描述符大多應用于描述剛性形狀,對于形狀發生等距、拓撲等非剛性變化不魯棒,因此不適用于非剛性三維形狀匹配.近年來,應用基于譜分析的形狀描述符進行非剛性三維形狀匹配成為了一個新的研究熱點,部分研究工作[39?44]表明,基于譜分析的形狀描述符對非剛性形狀的拓撲噪聲有較好的魯棒性,同時具有等距不變性.譜分析源于黎曼流形表面上的拉普拉斯-貝爾塔拉米(Laplace-Beltrami,簡稱LB)算子[45,46],LB 算子是一個著名的內蘊算子,它被分解為特征值和特征向量的組合,研究者常常將LB 算子的特征值稱為“譜”.為了研究方便,研究者將三維非剛性形狀定義為黎曼流形,通過LB 算子描述形狀特征,將形狀用“譜”的方法表示.這種在譜空間中進行形狀分析的方法被稱為譜分析[47].譜分析源于形狀的內蘊屬性,該方法提供了一種自然的方式進行非剛性三維形狀匹配.

譜分析包括譜形狀描述符及誘導出的譜距離,常用的譜形狀描述符包括形狀DNA(shapeDNA)[48]、全局點簽名(global point signature,簡稱GPS)[49]、熱核簽名(heat kernel signature,簡稱HKS)[50]、雙調合簽名(biharmonic signature,簡稱BS)[51]、波動核簽名(wave kernel signature,簡稱WKS)[52]等.在一個形狀表面上,由譜形狀描述符誘導出的譜距離[53]是一種較好的度量結構,具有等距不變性以及對拓撲變化的魯棒性.常用的譜距離包括交換時間距離(commute-time distance,簡稱CD)[54]、熱擴散距離(heat diffusion distance,簡稱HDD)[55]、雙調和距離(biharmonic distance,簡稱BD)[56]及波核距離(wave kernel distance,簡稱WKD)[57].使用譜形狀描述符進行三維非剛性形狀匹配時,需要選擇數量一致的采樣點,為了避免非剛性形狀匹配時采樣點的選擇問題,基于SD 方法,Bronstein 等人提出使用譜距離分布函數作為一種新的形狀描述符[58,59].譜距離分布函數將一對形狀的相似性轉化為其形狀上譜距離分布函數的相似性,同時兼具SD 方法和譜形狀描述符的優點,能描述形狀的全局統計屬性.Cao 等人也應用譜距離分布函數進行了三維非剛性形狀分類,并比較了幾種譜距離分布函數的性能[60].因此,基于現有研究方法,本文對基于譜分析的譜形狀描述符和譜距離分布函數用于三維非剛性形狀匹配進行了詳細的研究.

本文第1 節給出基于譜分析的三維非剛性形狀匹配的一般框架、LB 算子的詳細介紹及離散化計算.第2節首先詳細介紹了幾種譜形狀描述符:shapeDNA,GPS,HKS,BS,WKS,給出了譜形狀描述符的推導過程及其在離散網格上的計算方法;其次,總結與分析了這幾種形狀描述符在非剛性三維形狀匹配中的表現和特性.第3 節給出譜距離的定義和形式化表達,同時給出不同譜距離在三角網格上的離散計算方法以及譜距離分布函數的計算方式.第4 節是實驗驗證部分,實驗中使用不同譜形狀描述符和譜距離分布函數進行非剛性三維形狀匹配,觀察不同譜形狀描述符參數變化對匹配效果的影響,同時驗證了第2 節提出的預測的有效性,并做出合理分析.

現有的形狀描述符研究工作[9,61?66]對于譜分析的總結和介紹相對較少,大多數研究工作只是從理論上研究了譜形狀描述符或譜距離分布函數,很少有文章系統地對于兩類形狀描述符進行理論分析和實驗驗證.本文基于現有的研究成果,彌補了表1 中工作的不足,主要貢獻如下.

(1) 提供應用基于譜分析的形狀描述符進行三維非剛性形狀匹配的框架,并給出該方法的原理分析和數值計算;

(2) 系統地對比不同形狀描述符的數學定義及算法特性;從計算精度、魯棒性、時間復雜度等多方面比較其各自優缺點;并且在非剛性三維形狀標準庫中進行了兩類描述符的實驗比較;

(3) 給出不同譜形狀描述符和譜距離分布函數的最優使用場景,討論了譜分析應用于非剛性三維形狀匹配中存在的問題以及未來的發展趨勢,對譜分析進行推廣,為研究者選擇基于譜分析的形狀描述符提供參考.

Table 1 Several studies about summary and analysis of spectral analysis表1 譜分析的主要文章綜述與分析

1 基于譜分析的非剛性三維形狀匹配框架

本文首先對基于譜分析的非剛性三維形狀匹配的框架進行介紹.在數學中,譜分析是一個廣義的方法,它將一個矩陣的特征向量和特征值理論擴展到一個含有更廣泛運算符結構的譜空間中.在形狀匹配中,譜分析是指將形狀上的LB 算子離散化表示為LB 矩陣,對LB 矩陣進行譜分解得到LB 算子的特征向量和特征值.利用LB算子的特征值和特征向量,可以定義出不同的譜形狀描述符及其誘導出的不同的譜距離.通過計算一對形狀上的譜形狀描述符離散值或譜距離分布函數值,研究者可以對比一對形狀的局部或整體對應關系,得到一對形狀的非剛性匹配結果.本節首先給出非剛性匹配譜分析的一般框架,其次給出LB 算子的定義及離散化計算及譜分解形式.

1.1 非剛性匹配譜分析框架

LB 算子的特征值與特征向量常常被用來描述模型的形狀特性.利用譜形狀描述符和譜距離可以很好地進行非剛性形狀匹配,本文給出基于譜分析的非剛性三維形狀匹配的一般框架如圖3 所示.三維非剛性形狀匹配的具體流程如下.

· 第1 步:輸入一對3D 非剛性形狀(點云模型、三角片模型等).

· 第2 步:計算形狀上每個采樣點的LB 算子值,并將其進行譜分解,由LB 算子特征值和特征向量定義不同的形狀描述符,譜形狀描述符可以誘導出譜距離.

· 第3 步:對譜形狀描述符和譜距離分布函數進行離散化求值,得到譜形狀描述符矩陣和譜距離分布函數矩陣.

· 第4 步:應用方差或其他度量方法計算一對形狀間譜形狀描述符或譜距離分布函數數值,并選擇合適的度量函數進行形狀匹配,形狀匹配結果可以應用于形狀檢索、形狀分類、形狀對應等.

整個匹配過程中,形狀描述符的選擇是重要步驟,通過選擇合適的形狀描述符,研究者就可找到一對形狀間的局部或整體匹配關系.

Fig.3 Non-rigid shape matching framework using spectral analysis圖3 非剛性形狀匹配譜分析框架

1.2 拉普拉斯-貝爾塔拉米算子

作為譜分析中的重要算子,拉普拉斯-貝爾塔拉米算子是Laplace 算子黎曼流形上的推廣.Laplace 算子是歐氏空間中作用于光滑函數f的二階微分算子,描述為f的梯度的散度[67].任意二階可微函數的Laplace 算子定義為

根據黎曼流形梯度和散度的定義,若g為流形上的度量張量,G為矩陣{gij}的行列式,則函數f的LB 算子在局部坐標系中定義為[68]

在非剛性三維形狀匹配中,研究者需要計算離散網格上每個頂點的LB 算子值.網格上某頂點vi周圍三角形面片示意圖如圖4 所示.

在離散數學中,有限維離散LB 算子通常稱為離散LB 矩陣,是對連續LB 算子的一種逼近.在頂點數為n的三角網格上定義函數f,則該函數在網格頂點vi處的離散LB 算子可以定義為[69]

等式(3)中,當計算點vi的LB 算子時,考慮網格中的所有頂點.對于網格上某頂點vi,若僅對其周圍三角形面片能量求和,然后計算其偏導數并合并同類項,得到該點對應離散LB 算子的值:

其中,αj和βj分別表示連接vi和vj的邊eij兩側的對角,Neigh(vi)表示與頂點vi相鄰的頂點的集合.在完備有界的緊致流形上定義的LB 算子,具有對稱性和非負性.若將LB 算子進行譜分解(或稱特征分解)為特征值和特征向量的矩陣乘積形式,可得到流形上的LB 譜,即LB 算子的特征值和特征向量表達式:

λ0≥λ1≥…≥λn是LB 算子的非負特征值序列,λi是第i個特征值,φi是第i個特征值對應的特征向量.如果在封閉區域使用Neumann 邊界條件[70],第1 個特征值為0,一般將最小的非零特征值定義為λ1.由于LB 算子是半正定算子,所以λn≥0.LB 算子可以解析地計算一些形狀幾何量(例如矩形、圓柱、圓盤或球面等).如果一些形狀,例如動物、植物等,變換其形體姿態,例如在其關節處只有輕微的拉伸,這種情況被稱為近似等距變化,LB 算子同樣對近似等距變化魯棒.

Fig.4 Diagram of a vertex vion a triangular mesh圖4 網格上某頂點vi 的三角面片圖示

2 譜形狀描述符

由LB 算子的特征值和特征向量以定義不同的譜形狀描述符,例如上文提到的shapeDNA,GPS,HKS,BS 和WKS,本節對幾種譜形狀描述符的詳細定義以及離散計算方法進行介紹.

2.1 ShapeDNA

ShapeDNA 是Reuter 等人在2005 年通過提取黎曼流形表面的LB 算子的特征值序列進行非剛性形狀檢索,它的主要優點是易于表示形狀,計算簡單[48,71].ShapeDNA 是全局描述符,不能用于局部形狀分析.Reuter 等人對LB 算子特征值序列的數目進行了研究,當特征值序列數目等于500 時算法收斂,在工程應用中,數目取值一般為50～100.基于LB 算子的定義,形狀M上某點x的shapeDNAx可被計算為

g是被定義在形狀M上的度量,n為特征值序列的編號,shapeDNA 為非負值.ShapeDNA 很好地描述了形狀的內蘊屬性,不依賴形狀的參數化表示.形狀上,所有采樣點的shapeDNA 組成了shapeDNA 矩陣,確定n的取值后,可以唯一確定該形狀的shapeDNA 矩陣,但相似形狀的shapeDNA 矩陣非常近似,因此,其對于相似形狀的區分度較差.

2.2 全局點簽名算子

由于同類相似形狀的shapeDNA 值很近似,為了提高shapeDNA 對同類形狀的區分度,Rustamov 等人在其基礎上定義了一種新的譜形狀描述符——全局點簽名.如果將形狀內在的對稱性轉化為特征空間,將非剛性形狀的特征空間映射到一個無限維空間——全局點嵌入域(global point signature embedding dominant),那么在該無限維空間中,可以定義M上的每點x的GPS(x)可定義為

和shapeDNAx一樣,GPS(x)描述形狀的全局特征.形狀上每個點的GPS(x)都表示一個向量,一個形狀上所有采樣點的GPSM表示為一個m×n的矩陣,其中,m為形狀上采樣點的數量,n為LB 算子特征值數量,如等式(8)所示.

由上述定義可知,一個形狀的GPS 矩陣維數很高.研究者需要根據應用選擇合適的特征值數量n,以避免較高的計算量.

2.3 熱核簽名算子

根據熱擴散理論,假定在形狀上每點有初始熱源μ0(x),并隨時間t在形狀M表面上進行熱量擴散.在一定時刻,形狀表面上達到熱平衡狀態.在這個過程中,定義熱核ht(x,y)為t時刻從x點到y點轉移所需的熱量,表示熱量從一個點傳遞到另一個點的可能性.等式(9)為形狀上的熱擴散偏微分方程,描述了形狀表面上溫度隨時間變化狀態.

其中,μ(x,t)是形狀M上時間t對應的熱量分布函數,Δ是LB 算子.該方程的解為

同樣,對熱核進行譜分解:

熱核能完全表征一個形狀表面的幾何信息,如果將熱核限制在時間域內,可得到一個簡潔的形狀描述符——熱核簽名:

HKS 具有多尺度特性,能通過調節時間t改變其描述的是形狀的局部特征或者全局特征.形狀M在不同時間尺度下HKS 值分布可表示為

其中,每一列代表形狀在不同時間t下的HKS 值分布.圖5 給出了在較小t時刻下,3 個形狀的熱核簽名示意圖.可以從圖中看出,當t很小時,熱核簽名描述了形狀的局部特征[44].

Fig.5 Diagram of heat kernel function for a small fixed t on the hand,Homer,and trim-star圖5 較小t 值下手掌、人偶及五角星的熱核簽名圖示

基于HKS,Bronstein 等人對HKS 進行改進,提出了具有比例不變性熱核簽名(scale-invariant heat kernel signature,簡稱SIHKS)[72].該描述符采用對數采樣以及傅里葉變換,消除了一對形狀縮放前后的縮放倍數,在原有的HKS 上增加了縮放比例不變的特性.其具體過程如下.

· 首先,設縮放系數為β,對于形狀M,其發生縮放后的形狀為M′=βM.參照HKS 定義,縮放后的特征值和特征向量滿足λ′=β2λ,φ′=βφ,則縮放后,形狀M′上某點x處HKS(x)的譜分解形式可寫為

· 最后,對h取對數,消除β2的影響:,則.接著對hτ˙進行傅里葉變換,使時域的平移變換轉移到復數域:

對等式(15)兩邊取傅里葉模后得到等式(16):

文獻[72]從理論上證明了形狀縮放前后的的熱核簽名僅有時間軸上的偏移,SIHKS 具有尺度變換不變性.圖6 為原始馬和縮放變化后,馬的縮放不變熱核簽名圖示.

Fig.6 Scale-invariant heat kernel signature for the initial and scaled shape圖6 原始形狀和縮放變化后形狀的縮放不變熱核簽名圖示

2.4 雙調和算子

為了同時兼顧形狀的局部特性和全局特性,在HKS 和GPS 的基礎上,將LB 算子的特征值和特征向量進行另一種組合,在形狀M上的某點x定義另一種譜形狀描述符,即雙調和簽名:

與GPS 類似,形狀上的每個點的BS(x)都表示一個n維向量.一個形狀的BSM表示為一個m×n的矩陣,其中,m為形狀上點的數量,n為LB 算子譜分解數量.

BS 通過正則化拉普拉斯算子的特征值,很好地平衡了形狀的局部特征和全局特征.BS 算子來源于雙調和微分方程,該算子在形式上與GPS 非常相像,但是性能卻有很大提升,分母由LB 算子的特征值的平方根變為LB算子的特征值,大大加快了描述符的歸一化.與GPS 一樣,當我們選用BS 表示形狀時,需要根據應用場景選擇合適的譜分解數量n,以避免較高的計算量.

2.5 波核簽名

對于形狀上的每個點,通過測量不同能量級的量子粒子的平均概率分布,文獻[52]定義了一種新的形狀描述符——波核簽名,形狀表面上的量子粒子的演化由波函數控制.通過求解Schr?dinger 方程可得:

波函數的形式表達類似于熱核函數,但意義卻截然不同:熱核函數表示是熱量耗散,波動函數表示了能量的振蕩.其中,x是形狀上任意一點,Δ是LB 算子,i 是虛數,LB 算子和i 的乘積確保能量經過不同頻率振蕩后不會衰減.通過及譜分解理論可得,波函數φ(x,t)的譜分解形式為

其中,fk(t)為t時刻粒子的第k個頻率,φk(x)為第k個頻率對應的特征向量,可計算如下:

當t=0 時,表示期望值為E,頻率λk的概率分布.Laplace 譜沒有重復值,結合公式(20)及公式(21)可得:

|φ(x,t)|2為點x處粒子的概率分布.由于時間參數t對概率分布沒有直接影響,若只考慮能量參數,將WKS 算子定義為點x處能量為E的一個粒子可被測量到平均概率:

由上述可知,WKS 采用帶通濾波器,因此可以很好地分離形狀,如圖7 所示.

公式(24)中,WKS 的表達式具有一般性,可以通過選擇不同能量概率分布函數fE(λk)得到不同的WKS.選擇對數正態分布函數作為能量概率分布函數,則WKS 可寫為

eN為能量規模參數,eN=log(E),λk為LB 算子第k個特征向量,σ為正態分布的方差,Ce為正則化WKS 函數.在實驗中,本文采用與文獻[52]一樣的參數設置,則形狀的WKS 在不同能量規模下分布為

其中,WKSM(eN)形狀第m個頂點在能量規模為N下的WKS 值,每列代表不同能量規模下形狀M每點的WKS值.當我們選用WKS 表示形狀時,需要根據應用選擇合適的能量規模eN,以突出WKS 的優勢.

Fig.7 Wave kernel signature on a dog圖7 狗的波動核簽名圖示

2.6 譜形狀描述符比較

在譜形狀描述符中,shapeDNA 的研究時間最早,因此整體性性能相對較差,但其為之后譜形狀描述符的發展奠定了基礎.每點的shapeDNA 由LB 算子的前n個特征值確定,忽略了LB 算子的特征向量的作用,shapeDNA最大的優點是易于理解,計算簡單.但對局部特征的描述能力較弱,不適合相似形狀的快速區分.

GPS 由LB 算子的前n個特征向量比上特征值得到.從定義上看,GPS 更加注重低頻相關的信息,但對于形狀發生非剛性形變(變化較小)時,形狀上一點的GPS 值也會完全改變,增加額外的算法復雜度,性能并不好.總體來說,GPS 能夠很好地反映形狀上所有采樣點的上下文信息,但對局部特征的描述能力較弱,適合非剛性形狀的粗糙匹配,不適用于局部匹配.

HKS 定義了點x處的局部和全局屬性.由于形狀上的熱擴散本質近似布朗運動,因此有較強的魯棒性以弱化局部噪聲的影響.作為低通濾波器的集合,HKS 主要由低頻傳輸,能夠很好地描述形狀的局部幾何信息.當時間t比較短時,形狀上每個點的HKS 值與該點的高斯曲率直接相關,具有很強的信息存儲性.但HKS 過于強調低頻信息,會過濾掉一些高頻率信息,損害精確定位特征的能力.因此相比較其他3 種算子,HKS 表征形狀時不能很好地分離不同頻率區域.此外,由于HKS 對時間參數敏感,所以在某個時刻下,不能同時表征形狀的局部屬性和全局屬性.SIHKS 擁有HKS 所有的優點,還具有縮放不變性,但是其理解相對較難且計算復雜.

BS 平衡了大尺度距離(反映全局特性)和小尺度距離(反映局部特性),具有多尺度特性.它不依賴于時間參數,克服了HKS 依賴時間參數重的缺點,同時克服了GPS 沒有多尺度特性的缺點.然而,由于BS 具有調和性質,單獨表征形狀的局部及全局屬性性能較差.

WKS 同樣對時間參數自由,其最大優點是采用帶通濾波器能清楚地分離形狀上的不同頻率集合,且允許訪問高頻率信息,從而增加算子的精確匹配能力.此外,WKS 通過選擇不同的能量規模而具有多尺度特性,若選擇能量級別較高的量子粒子,波長越短,其分布越靠近形狀上的點,此時反映形狀的局部特性;反之,能量級別較低的量子粒子反映形狀的全局特性.所以在匹配時,研究者應該根據應用場景選擇一個合適的能量規模.

幾種譜形狀描述符在不同變換下魯棒性等級見表2.

Table 2 Robustness levels of several spectral shape descriptors under different transformations表2 幾種譜形狀描述符在不同變換下魯棒性等級

3 譜距離分布函數

譜距離(shape spectral distance distribution)源于譜分析,譜距離由形狀表面上定義的譜形狀描述符誘導得到,包括熱擴散距離、交換時間距離、雙調和距離等.若在形狀M上定義度量空間,則由譜形狀描述符誘導出的譜距離可定義為[62]

其中,φ(λi)為不同譜形狀描述符使用的濾波器.在三角面片上,譜距離的離散化形式為[63]

其中,p,q為三角面片上的頂點,分別代表頂點p和q上LB 算子第i個特征向量.前文中提到的另一類基于譜分析的形狀描述符就是譜距離分布函數,基于SD 方法的研究,Brostein 等人通過統計形狀上任意兩點間的譜距離,構造了譜距離分布函數最為一種新的形狀描述符.假定形狀M上任意兩點的譜距離為d(x,y),若δ是距離閾值,μ是定義在M中的范數矩陣,χ是指示函數,則譜距離頻率直方圖可以計算為

在概率論與統計學中,概率密度函數(probability density function,簡稱PDF)是一個實值隨機變量,用于描述多隨機變量的分布,再由公式(29)可得譜距離分布函數可計算為

譜距離分布函數作為一種線性形狀描述符,繼承了譜距離的特征,具有以下特點.

(1) 采樣不變性:對于形狀M,如果對M的三角網格模型的頂點進行采樣,包括上采樣和下采樣,則采樣前后的形狀的譜距離分布函數保持不變.

(2) 等距不變性:由于LB 算子是形狀的內蘊算子,所以等距形狀中任意兩點的譜距離具有等距不變形.因此,等距形變前后,形狀的譜距離分布函數理論上保持不變.但是在下文中,我們給出了不同譜距離的譜分解計算形式,在實際應用中,一般取前100 個特征值和特征向量.因此在實際的實驗中,等距形狀的譜距離分布函數與原始形狀的譜距離分布函數值相似.

(3) 拓撲魯棒性:相對測地距離分布,譜距離分布對拓撲變化的敏感性較低,譜距離分布函數具有較強的拓撲魯棒性.

(4) 無需預處理:相對于譜形狀描述符,應用譜距離分布進行三維非剛性形狀匹配時,不需要尋找數量相同的采樣點,也不需要配準采樣點.

下文就詳細對這4 種譜距離進行介紹.

3.1 交換時間距離

根據第2.3 節,在GPS 域中的內積可定義交換時間距離:

G(x,y)是兩個無限維向量的內積,交換時間距離可寫為

交換時間距離反映了連接一對點之間隨機游走的平均時間.通過譜分解,交換時間距離可以表達為

其離散化形式為

熱擴散距離和交換時間距離的關系為

熱擴散距離反映了形狀表面上兩個點在時間t內的路徑連通性,而交換時間距離是一對點之間在平均時間t內隨機游走的擴散長度之和.

3.2 熱擴散距離

熱擴散距離由擴散核導出,并應用于降維和數據參數化等問題.擴散距離描述了形狀M上點x與點y之間在時刻t時的連通率.將形狀M上的隨機運動看作是布朗運動,在這種情況下,擴散距離是對形狀M上時間t時兩點間的布朗運動的平均,更直觀上的理解為兩個熱核之間的信息交互.所以形狀上的兩點x和y之間的擴散距離可以由下面的等式定義[39,51,73].

根據熱核的譜分解形式以及熱擴散理論,熱擴散距離(也稱熱核距離)表示為

為了不失一般性,特征值從1 開始,離散化形式為

熱擴散距離反映了擴散時間t內兩點間的熱流連通性.由于該距離是定義在形狀表面的距離,所以是一個內蘊距離,當擴散時間t的取值很小時,此時熱量擴散范圍較小,該距離只能描述形狀的局部特性;當t的取值較大時,該距離可以描述形狀的全局屬性,最后熱量擴散直至達到熱平衡狀態.但t取值并不是越大越好,合適的t取值[60]為1/λj,λj為第LB 算子的第1 個非零特征值.

3.3 雙調和距離

類似GPS 映射,雙調和映射定義了一個無限維的雙調和空間.

雙調和域中的內積可定義雙調和距離[40]:

B(x,y)是兩個無限維向量的內積,雙調和距離可表示為

通過譜分解,雙調和距離可寫為

離散形式可寫為

3.4 波核距離

文獻[66]用L2范式定義波核距離,類似于熱核距離,波核距離的譜分解形式為

離散化形式為

4 實驗與分析

為了對幾種譜形狀描述符和譜距離分布函數進行對比,本文在64 位32G 內存,win10 系統的Matlab2015 上進行了實驗上進行了實驗(注:由于shapeDNA 最早被研究,性能較差,因此在本文不加入比較).評估這種方法所使用的是TOSCA 2010(tools for surface comparison and analysis)數據集(高分辨率,http://tosca.cs.technion.ac.il/data/toscahires-mat.zip)[74]、SHREC 2011 數據集(魯棒性,http://tosca.cs.technion.ac.il/book/shrec_robustness2010.html)[75]和 SHREC 2015 數據集(標準型,http://www.cs.cf.ac.uk/shaperetrieval/shrec15/SHREC15.zip).TOSCA 2010 數據集為非剛性形變的形狀匹配提供了大量高清三維形狀,圖8 為部分TOS CA2010 庫中形狀.TOSCA 2010 數據庫共包含80 個對象,包括11 只貓、9 只狗、3 只狼、8 匹馬、6 人馬、4 只大猩猩、12 個女性人物、2 個不同的男性形象,典型的頂點計數大約是50 000,數據庫適用于非剛性形狀分析.SHREC 2011 包含13 類形狀的的各種變化形狀,變化分為12 類,包括等距變化及在等距變化上加入洞、縮放、噪聲、下采樣等變化,每種變化共有5 個強度,圖9 為部分SHREC 2011 庫中形狀.

Fig.8 Part of the shapes of TOSCA 2010 database圖8 TOSCA 2010 數據庫中的部分形狀

Fig.9 Part of the shapes of SHREC 2011 database圖9 SHREC 2011 數據庫中的部分形狀

SHREC 2015 數據庫由SHREC 2011[76]和SHREC 2014[77]中的部分形狀組合而成,包含10 類形狀,每類形狀包含了基礎形狀的等距變化、近似等距變化、拓撲變化及加洞等非剛性變化,共計100 個形狀,為非剛性三維形狀檢索提供標準形式,圖10 為SHREC 2015 中的部分形狀.

Fig.10 Part of the shapes of SHREC 2015 database圖10 SHREC 2015 數據庫中的部分形狀

4.1 譜形狀描述符參數比較

HKS 具有多尺度特性,由時間參數t決定該點描述的是形狀的的局部或全局屬性.圖11 中選取不同時刻應用HKS 進行一組形狀等距對應,顏色由黃到藍代表HKS 值由大到小,當t=0.5 時,此時熱量剛開始擴散,只能描述行david 足部的局部幾何信息,此時,HKS 值與足部的高斯曲率值直接相關;當t=1 時,熱量擴散到形狀的大部分區域;當t=3.0 時,熱量擴散到形狀整個表面,此時,HKS 描述形狀的全局幾何信息.

Fig.11 Non-rigid shape matching using heat kernel signature under different time t (shape:david 1 & david 10)圖11 在不同時間t 下應用熱核簽名進行一組非剛性形狀匹配(shape:david 1 & david 10)

WKS 對時間參數自由,圖12 中選取不同能量級下進行一組非剛性形狀匹配.當能量級下的數量增大時,WKS 能更精確地表達形狀的局部特征,但其數量并不是越大越好,過大的能量級會增加導致WKS 無法刻畫形狀的全局特性,間接增加更多計算誤差;如圖12 所示,本文選取e100作為合適的能量規模.

Fig.12 Non-rigid shape matching using wave kernel signature under different energy scale eN(shape:david 1 & david 10)圖12 在不同能量級(eN)下應用波動核簽名進行一組非剛性形狀對應(shape:david 1 & david 10)

圖13 基于庫SHREC 2010 驗證了表1 每種譜形狀描述符在不同等距變化下的魯棒性.可以看出,GPS 和WKS 總體魯棒性較高.

Fig.13 Compared with GPS,HKS,BS,WKS robustness in isometric,sampling,cave,noise and topological changes (t=3.0,λNum=100,e100)圖13 對比GPS、HKS、BS、WKS 在等距、采樣、加洞、噪聲及拓撲變化下魯棒性(t=3.0,λNum=100,e100)

4.2 譜形狀描述符用于三維非剛性形狀匹配比較

從圖14 中可以看出,GPS 作為一個全局形狀描述符,對cat 0 和cat 1 足部和腿部的細節描述不夠,不能應用到局部匹配中,但是能夠分清cat 0 和cat 1 前足和后足.

Fig.14 Non-rigid shape matching using GPS,HKS,BS,WKS (shape:cat 0 & cat 1,t=3.0,λNum=100,e100)圖14 應用GPS、HKS、BS、WKS 進行非剛性形狀匹配(shape:cat 0 & cat 1,t=3.0,λNum=100,e100)

而當時間參數足夠大時,HKS 能表征cat 0 與cat 1 的局部幾何信息和全局幾何信息,但由于HKS 使用的都是一些低通濾波器,形狀的高頻率信息被抑制,不能精確地表示形狀,相比GPS,HKS 能分清cat 的腿部和身體,但沒有辦法區分貓的前足和后足;BS 表現最佳,cat 1 相對cat 0 尾巴發生較大扭曲,此時,cat 0 和cat 1 為近似等距變化,BS 能夠明確地描述尾部的近似等距變化且匹配度高;同時,WKS 在貓的尾部匹配度同樣較高,且相比HKS,WKS 使用帶通濾波器,減少低頻的影響,在圖中能夠清楚地分離出形狀的頻帶區域,具有優越的特征定位,且能區分貓的四足,適合高精度的匹配,但算法時間復雜度較高.

通過10 次實驗求取平均值,幾個形狀描述符耗費時間如表3 所示,應用4 種譜形狀描述符進行非剛性三維形狀匹配的空間復雜度為O(n),n為三維形狀上采樣點的個數.由于GPS 和BS 都是高維向量,因此其耗費時間要多于HKS;同時,WKS 采用帶通濾波器分離形狀上的不同頻率集合,對于形狀的細節刻畫較多,因此時間耗費相對最高.為了比較應用不同譜形狀描述符進行非剛性三維形狀匹配的匹配度,實驗中,我們計算一對形狀上采樣點的譜形狀描述符的相關系數R=corr2(A,B)作為三維非剛性形狀匹配的匹配度(即一對形狀的相似度),結果如表4 所示,WKS 和BS 都對參數自由,BS 能夠調和地描述形狀的全局和局部屬性,WKS 能夠在描述形狀全局屬性的同時刻畫形狀的細節.因此,應用WKS 和BS 的形狀匹配相對較高.

Table 3 Time-consuming of non-rigid shape matching using spectral shape descriptors表3 應用譜形狀描述符進行非剛性形狀匹配所耗費時間

Table 4 Non-rigid shape matching using spectral shape descriptors表4 應用譜形狀描述符進行非剛性形狀匹配

4.3 譜距離分布用于三維非剛性形狀匹配比較

有效的譜距離分布函數可以區分不同類別的形狀,且對于形狀發生非剛性變化,譜距離概率分布趨勢差別較小,故可以通過匹配形狀的譜距離分布,進行三維非剛性形狀匹配[78].圖 15 是計算一對形狀(選用最TOSCA 2010 中最復雜的centuar0 和centuar1)的4 種譜距離分布:CD,HDD,BD,WKD,給出譜距離概率分布(注:由于整體分布趨勢接近,無法看出差別,故同時給出分布概率小于等于0.1 的分布圖作比較),可以很直觀地從圖15中看出,centuar0和centuar1的WKD概率分布趨勢一致.說明WKD具有良好的精確匹配性,通過圖中centuar0和centua1 的WKS 值對應,可以發現:應用WKD 概率分布能夠區分半人馬的胳膊、4 條腿;與上述相同,BD 概率分布同樣趨勢一致,但在區分本人馬的前腿和后腿時區分度不大,只能分清前腿和后腿;CD 概率分布略有差別,由于GPS 對局部細節表述不夠,導致CD 值有差異,圖中的半人馬僅能區分胳膊和腿部;HDD 概率分布總體走勢一致,但差異較大,無法區分半人馬的胳膊和腿,但同樣的,對于局部細節刻畫清楚.

通過10 次實驗求取平均值,幾個譜距離分布函數的耗費時間見表5.應用4 種譜形狀描述符進行非剛性三維形狀匹配的空間復雜度為O(n2),n為三維形狀上采樣點的個數.為了比較非剛性形狀應用不同譜形狀分布函數進行匹配的匹配度,實驗中,為了直接得到一對形狀的匹配度,我們同樣計算一對形狀的譜形狀距離分布函數的相關系數作為三維非剛性形狀匹配的匹配度(即一對形狀的相似度),結果見表6.相比應用譜形狀描述符進行形狀匹配,應用譜距離分布進行形狀匹配時,所有形狀的匹配度都有提升.原因在于:應用譜距離分布進行形狀匹配時不考慮形狀的局部細節匹配,得到的是一對形狀的全局匹配結果.因此,相比于譜形狀描述符,譜距離分布函數在進行三維非剛性形狀匹配時無需尋找一對形狀的對應點,穩定性更強,更適用于非剛性形狀整體匹配.

Fig.15 Distribution function of four spectral distances for non-rigid shapes using matching(shape:centaur 0 & centua 1)圖15 4 種譜距離分布函數進行非剛性形狀匹配(shape:centaur 0 & centua 1)

Table 5 Time-consuming of non-rigid shape matching using the distribution function of four spectral distances表5 應用4 種譜距離分布函數進行非剛性形狀匹配所耗費時間

Table 6 Non-rigid shape matching using the distribution function of four spectral distances表6 應用4 種譜距離分布函數進行非剛性形狀匹配

結合表3 和表5 我們可以看出,應用譜形狀描述符進行形狀匹配時的時間復雜度和空間復雜度要比應用譜距離分布函數的時間復雜度低,因為計算形狀上任意兩點的譜距離會耗費較多的時間和占用較多的空間.同時,結合表4 和表6 我們可以看出,應用譜形狀描述符進行形狀匹配時,隨著形狀大小的增加,形狀的匹配度會略微降低;反之,應用譜距離分布函數進行形狀匹配時,隨著形狀大小的增加,形狀的匹配度會略微增高.因為當形狀增大時,形狀的三角網格面片數也會增加,導致采樣點的譜形狀描述符的計算量大大增加,降低形狀匹配度;反之,形狀的三角網格面片數增加,會增大譜距離分布函數的樣本量,更能反映形狀的全局屬性,進一步提高形狀的匹配度.

圖16 為基于SHREC 2015,應用4 種譜距離分布函數進行非剛性匹配的形狀相似性熱力圖.為了區分不同類形狀的差異性,本文采用最常用的歐氏距離計算一對形狀的相似性.

Fig.16 Thermodynamic diagram of non-rigid 3D shape matching using four spectral distances based on SHREC 2015 database圖16 基于SHREC 2015 數據庫,應用4 種譜距離進行非剛性三維形狀匹配的熱力圖

如圖16 所示(其中,每連續10 個編號表示一類形狀,1～10 為centaur;11～20 為ants;21～30 為gorilla;31～40 為Male 0;41～50 為female-thin;51～60 為male 13;61～70 為gdog;71～80 為male16;81～90 為plies;91～100 為malebodybuilder).在4 種譜距離分布函數中:CD 描述了形狀的全局屬性,因此匹配結果較為集中;HDD 描述了形狀的局部屬性,且相對于其他3 種譜距離分布函數的匹配性能較差,原因在于HDD 對于時間參數和噪聲非常敏感,在實際實驗中很難尋找到最優的時間參數;BD 調和地描述了形狀的局部和全局屬性,且不依賴于時間參數t;WKD 既能描述形狀的局部屬性,又能描述形狀的全局屬性,同時,相對于其他3 種譜距離分布函數,WKD 能夠精準地描述不同類的形狀,對于male0,male13 和male16 的區分度更大.

表7 為應用圖16 的結果進行不同類非剛性三維形狀檢索的查準率,表中結果與圖16 結果相一致.

Table 7 Precision ratio of 3D Non-rigid shape retrieval using the distribution function of four spectral distances based on SHREC 2015 database表7 基于SHREC 2015 數據庫,應用4 種譜距離分布函數進行非剛性三維形狀檢索查準率

4.4 問題與展望

通過實驗比較了不同譜形狀描述符和譜距離分布函數進行非剛性三維形狀匹配的性能,可以發現利用基于譜分析的形狀描述符進行非剛性形狀匹配效果較好.在4 類譜形狀描述符中,WKS 和WKD 整體匹配表現性能最優,適用于精細匹配,但時間復雜度較高,不適用于大規模形狀的快速匹配;BS 和BD 性能次之,能同時調和地表示形狀的局部及全局信息,但過于強調函數同時描述形狀的局部和全局性質,也會弱化形狀的真實全局和局部特性;HKS 和HDD 對時間參數敏感,所以僅憑某時刻t下HKS 值進行形狀的非剛性匹配性能較差,若能比較形狀上每個點或部分特征點的一段時間序列下的HKS 值,則能提高匹配性能,且HKS 適用于部分匹配;GPS和CDD 整體匹配度較低,適用于快速匹配和粗匹配,但不適用部分匹配.同時,使用4 種譜距離分布函數進行三維非剛性形狀匹配時,無法得到形狀部分匹配結果,只能得到一對形狀的全局匹配結果,其時間復雜度和空間復雜度也比譜形狀描述符復雜度高.在未來的研究工作中,針對不同變化類型的形狀(如一些近似等距變化或大變形形狀)及不同的應用場景,應結合幾種描述符的優點,考慮同時使用多個描述符的權重組合或其他改進,提升非剛性三維形狀匹配度.

5 總結

本文給出基于譜分析的形狀描述符進行非剛性三維形狀匹配的方法流程,詳細介紹了幾種譜形狀描述符和譜距離分布函數,并在以下幾方面對比了不同形狀描述符的性能:(1) 局部及全局屬性;(2) 有無參數及參數選擇;(3) 時間復雜度;(4) 最優匹配度;(5) 適用匹配場景;(6) 整體表現性能.通過實驗,證明了本文預估的正確性.譜分析是一個易于理解、普適且魯棒的分析方法,基于譜分析的形狀描述符在非剛性三維形狀整體匹配中表現出了優異的性能,我們希望提升基于譜分析的形狀描述符在非剛性形狀匹配中重要的理論意義及并推動其在工程應用價值的發展.同時,在未來的工作中,我們會對基于譜分析的形狀描述符進行非剛性三維形狀局部匹配進行進一步研究.