等距變形體的矩不變量構造方法

周建存，常亮，郭克華

(1. 湖南城市學院 信息科學與工程學院，湖南 益陽，413000；2. 中南大學 信息科學與工程學院，湖南 長沙，410083)

近年來，變形體的識別問題在模式識別領域內(nèi)越來越多地受到關注。變形體識別具有較為廣闊的應用前景，如表情無關的三維人臉識別[1-2]、三維紋理分析[3]和醫(yī)學圖像處理[4-5]中很多技術都與變形體的識別有關。在變形體識別問題中，目標發(fā)生變形，特征難以提取，該問題成為模式識別領域內(nèi)難點。傳統(tǒng)的三維目標識別方法[6]主要針對剛體進行識別，難以處理目標變形情況下的問題。目前，對變形體的研究主要集中在變形過程的計算機仿真，而對變形體匹配的研究很少。由于對任意情況下變形體進行識別比較困難，因此，目前的研究主要集中在具有條件限制下的變形體識別。在變形體識別的研究中，目前熱門的研究是對等距變形體進行識別。等距變形的定義如下：給定為二維歐幾里德流形為這 2點之間的測地距離。對于變換ψ:S→S'，若滿足：

則稱變換ψ為等距變換，此變換下的變形稱為等距變形。可見：等距曲面變形后，表面相應點之間的測地距離保持不變。因此，測地距離是等距變形體的基本特征，測地距離的不變性成為等距變形體識別的關鍵。在已有的等距變形體的研究中，最傳統(tǒng)的方法是將變形體劃分為幾個未變形的部分，然后，基于這些部分分別進行識別。Bloomenthal等[7-8]基于這種思想，構造了變形體的骨架，在具有關節(jié)結構的變形體識別中，取得了較好效果。但是，這種方法無法針對目標的任意等距變形進行識別。一些研究者提出了基于幾何統(tǒng)計方法的變形體識別策略。Aherne等[9-11]首先計算目標表面一些微分幾何量，然后構造直方圖進行統(tǒng)計，將這種方法在三角網(wǎng)格上得到應用。Osada等[12]則通過評估目標表面任意點的某個歐幾里德值(如距離)的分布函數(shù)來計算目標的統(tǒng)計特征。這些方法比較簡單，但是，只對小范圍內(nèi)的變形具有魯棒性。目前比較好的一種等距變形體識別方法是將等距曲面等價映射成1個高維剛體，任意適用于剛體匹配的算法都能起作用，如Elad等[13]將這種方法應用到表情無關的三維人臉識別中，取得了較好效果，但是，該類算法運算復雜度較高。為此，本文作者針對等距變形這種特殊的變形體識別問題進行研究，提出該變形情況下目標識別的新方法。該方法通過構造不變量來描述等距變形體的特征，首先利用測地距離的概念構造出等距變形體的特征矩陣，然后對特征矩陣進行歸一化，最后利用矩不變量對歸一化特征矩陣的特征進行提取，構造出相應的矩不變量，并對該類矩不變量的平移、尺度和旋轉不變性進行證明。

1 等距變形體的特征數(shù)據(jù)存儲

1.1 測地距離的獲取

等距曲面變形之后，表面相應點之間的測地距離保持不變。因此，測地距離是等距變形體的基本特征，可以通過計算曲面上任意2點之間的測地距離來提取不變量。

對于測地距離的獲取，Sethian[14]提出了基于矩形網(wǎng)格的測地距離計算方法即矩形網(wǎng)格上的快速行進法，Kimmel等[15]將此方法進行改進，提出了三角網(wǎng)格下的快速行進法。對于含有n個頂點的曲面，這種方法能夠計算任意頂點和其余各點之間的測地距離，運算復雜度為O(n)。通過計算每個頂點和其他頂點之間的測地距離，可以在O(n2)運算時間內(nèi)構造1個測地距離矩陣。

1.2 特征矩陣的構造

利用三角網(wǎng)格下的快速行進法算法可以計算網(wǎng)格曲面上任意2點之間的測地距離，構造測地距離矩陣。該矩陣具有如下特征：對于具有n個頂點的等距曲面，此矩陣是n×n的對稱矩陣，對角線上的元素為0。為了表述方便，可將此測地距離矩陣，確定為等距曲面的特征矩陣(SM)，滿足：

其中：dS(pi，pj)表示曲面S上頂點pi和pj之間的測地距離，此距離由三角網(wǎng)格下的快速行進法算法求出。

1.3 特征矩陣的歸一化

構造出特征矩陣后，目標之間的匹配轉化為特征矩陣之間的匹配。但是，由于曲面上像素點提取的順序是任意的，所以，同一個曲面可能構造出不同的特征矩陣。這些特征矩陣對應著一些同構的帶權圖，若直接針對特征矩陣進行匹配，則其復雜度等價于同構圖的匹配，運算時間為O(n!)。所以，必須對特征矩陣進行歸一化，以保證同一個目標得到的是同一個特征矩陣。

歸一化的目的是調整頂點的順序。實際上，頂點vi和vj的調換將導致特征矩陣中第i行與第j行、第i列與第j列之間的調換。這里需要對特征矩陣中的頂點進行排序，以保證特征矩陣的唯一性。

該算法的基本思想是：基于每個頂點的可量化特征，將頂點按升序排列，本算法中選取頂點的均值和方差進行計算。步驟如下。

首先，計算每個頂點的均值和方差：

其中：1≤i≤n；n為頂點個數(shù)；SM為特征矩陣。

其次，構造向量K來存儲特征矩陣中每一行在歸一化特征矩陣中的位置。K初始化為零向量，K(i)通過以下準則確定：

其中：n′為示集合的元素個數(shù)；(mj，vj)須滿足如下條件：

對于K(i)的構造，實際上基于每個頂點的均值和方差，將頂點按升序進行排列。而對特征矩陣中的頂點進行的排序，可以保證特征矩陣的唯一性。

然后，構建變換矩陣T來存儲頂點順序的調整情況，T為1個n×n矩陣，滿足：

T中其他元素為0。

最后，構造歸一化后的特征矩陣NS。該矩陣可通過下式計算：

其中：SM為等距曲面的特征矩陣；T′為矩陣T的轉置。

綜上所述，歸一化特征矩陣的構造算法可以描述如下。

Step 1: 利用式(3)計算每個頂點的均值和方差。

Step 2：基于Step 1中計算的均值和方差，根據(jù)式(4)創(chuàng)建向量K(i)，計算每個頂點在歸一化矩陣中的位置。

Step 3：基于Step 2中計算的個頂點對應的K(i)，根據(jù)式(6)構建變換矩陣T，確定對特征矩陣進行歸一化的變換矩陣。

Step 4：利用式(7)，根據(jù)Step 3中計算的變換矩陣，對特征矩陣進行變換，計算出歸一化特征矩陣NS。

在以上的算法中，計算每個頂點的均值和方差，運算時間為O(n2)；利用快速排序算法，運算時間為O(nlgn)；在Step 3中，變換矩陣T的行向量和列向量都是單位向量，所以，NS的每行能夠根據(jù)變換矩陣中1的位置從SM中提取，其運算復雜度為O(n2)。綜上所述，歸一化步驟的運算復雜度為O(n2)。

2 等距變形體的變形矩及其性質

歸一化特征矩陣的特征可以由許多方法來提取，如傅里葉描述子、矩特征等。本文將對矩不變量進行擴展，以提取NS的特征。

對于歸一化特征矩陣，其p+q階中心矩定義如下：

其中：N為自然數(shù)集合。x0和y0為歸一化特征矩陣的質心坐標，滿足：

顯然，由于對特征矩陣進行了歸一化，該矩陣中頂點的順序進行了調整，中心矩vpq具備了平移不變性。

為保證尺度變換下的不變性，需要對中心矩進行歸一化，為此，構造出歸一化的中心矩μpq：

根據(jù)已經(jīng)構造出來的歸一化的中心矩μpq，進一步構造2個針對等距變形體的變形矩：

要使全速變形矩得到應用，必須具有不變性。對于前面秘構造的變形矩，下面從平移、尺度和旋轉 3個方面闡述其不變性質。

(1)平移不變性。由于歸一化特征矩陣中頂點的順序進行了調整，中心矩vpq具備了平移不變性，因此，根據(jù)式(10)構造的歸一化的中心矩μpq也具有平移不變性。由此可見：式(11)構造出來的變形矩也具有平移不變性。

(2)尺度不變性。設目標邊界Σ為光滑的空間曲面，坐標x，y和z縮放同一因子k(k＞0)，此時，曲面上各點的坐標也將縮放因子k，特征矩陣SM和歸一化特征矩陣NS中的元素也將縮放因子k；但是，由于式(10)中對vpq進行了處理，其分母也是以因子k進行縮放，因此，upq具有尺度不變性，由式(11)構造的變形矩也具有平移不變性。

需指出的是：旋轉變換不會改變測地距離的值。由于歸一化特征矩陣中頂點的順序進行了調整，因此，旋轉變換只會改變特征矩陣中頂點的順序，但不會改變歸一化特征矩陣。因此，式(11)構造的變形矩具有旋轉不變性。

3 實驗結果和復雜度分析

因為上述理論是基于三維變形曲面的，所以，必須首先獲取圖像，然后構造三角網(wǎng)格。目前，對于三維目標圖像的數(shù)據(jù)獲取一般采用 2種方法：(1)由于二維圖像的存儲比較方便，可以通過從二維圖像進行三維建模來獲取三維數(shù)據(jù)；(2)直接使用可以獲取三維位置信息的工具，如三維掃描儀等。在實驗中，選取幾個比較典型的三維物體來進行驗證。

3.1 本文方法與傳統(tǒng)方法計算復雜度的比較

利用一些三維圖像測試本文算法的分類效果。變形體研究目前還沒有公認的數(shù)據(jù)庫，為此建立1個包含20個目標的數(shù)據(jù)庫(見圖1)。在這20個目標中，有些目標具有相似之處。為了對效果進行對比，將一些傳統(tǒng)方法應用到分類過程之中。場景像素為 450。對于數(shù)據(jù)庫中的每個目標，首先計算它們的特征矩陣，然后歸一化，最后計算變形矩。

圖1 數(shù)據(jù)庫中的20個目標Fig.1 20 Objects in Database

為了采用更多的樣本，對于每個目標，對其實施任意的等距變形、旋轉和縮放變換，獲得10幅圖像。總數(shù)據(jù)庫中有200幅三維目標圖像。

首先，利用本文算法對這200個不同的樣本進行分類，然后，與傳統(tǒng)的高維映射方法進行對比。

此外，傳統(tǒng)的高維映射方法也用于本文數(shù)據(jù)庫中。用MATLAB在P IV計算機上對圖1所示的三維圖像，分別采用本文算法和傳統(tǒng)的高維映射算法進行計算。這2種算法的運算時間如表1所示。

表1 2種算法的計算時間比較Table 1 Computational time of two algorithms ms

從表1可看出：本文算法在運算方面具有較低的運算時間，因而復雜度低。

實際上，高維映射方法是將本文的特征矩陣映射到高維空間的剛體，其過程主要利用了 MDS算法。在高維映射算法中，對傳統(tǒng)MDS算法、LSMDS算法和快速MDS算法進行比較。由于快速MDS誤差率較高，最終采用LSMDS算法，其運算復雜度為O(n2×N)(N為迭代次數(shù))，其精確度依賴于迭代次數(shù)。本文算法中，歸一化特征矩陣的構造的運算時間為O(n2)，小于LSMDS算法的運算時間。此外，預處理后的目標之間的匹配還需要花額外的時間，在這一點上，這2種算法所花費時間區(qū)別不大。2種算法的計算復雜度如表2所示。

表2 2種算法的計算復雜度Table 2 Computational complexity of two algorithms

因此，從理論上分析，本文的算法也具有較低的復雜度。

3.2 本文方法與傳統(tǒng)方法識別率的比較

本實驗中，基于樣本數(shù)據(jù)庫，選取圖1中的3個目標進行分類。首先，對圖像庫中的樣本進行了旋轉或尺度變化，選取J1和J2矩不變量來進行計算。此實驗中，利用一階Minkowski距離來衡量目標之間的差別。令S和S′為2個等距變形曲面。其一階Minkowski距離定義為：

本實驗中閾值取0.1。 當物體之間的距離小于或等于0.1時，就將它們歸入同一類。2種算法分類效果(識別率)如表2所示。

表3 2種算法的分類效果比較Table 3 Performance comparison of two algorithms

從表3可看出：本文算法在預處理方面具有較低的運算復雜度，但是，并沒有降低識別率。

3.3 噪聲情況下的魯棒性比較

對數(shù)據(jù)庫中的每幅圖像增加高斯噪聲(N(0，σ))，σ在0～2.0 mm之間變化。與傳統(tǒng)的高維映射算法的分類效果對比結果如圖2所示。

從圖2可以看出：本文算法和傳統(tǒng)算法對噪聲都具有較好的魯棒性。實際上，因為這2種方法都是基于距離進行計算的，因此，噪聲對計算結果的影響不是太大。

圖2 高斯噪聲變化時的識別效果Fig.2 Performance comparison in gaussian noise

4 結論

(1)針對等距變形這種特殊的變形體識別問題進行研究，提出了該變形情況下目標識別的新方法。該方法通過構造不變量來描述等距變形體的特征，首先利用測地距離的概念構造出等距變形體的特征矩陣，然后對特征矩陣進行歸一化，最后利用矩不變量對歸一化特征矩陣的特征進行提取，構造出相應的矩不變量，并對該類矩不變量的平移、尺度和旋轉不變性進行了證明。

(2)在不降低識別效果的前提下，本文構造的方法和傳統(tǒng)的高維映射算法相比，不需要將三維目標映射到高維空間，具有較低的運算復雜度，并具有較強的魯棒性。

(3)下一步工作是對有遮擋情況下的等距變形體識別問題以及任意拓撲變形情況下的目標識別問題進行研究。

[1]文沁, 汪增福. 基于三維數(shù)據(jù)的人臉表情識別[J]. 計算機仿真, 2005, 22(7): 99-103.WEN Qin, WANG Zeng-fu. Recognition of human facial expression based on 3D data[J]. Computer Simulation, 2005,22(7): 99-103.

[2]PAN Gang, HAN Song, WU Zhao-hui, et al. Removal of 3D facial expressions: A learning-based approach[C]//2010 IEEE Conference on Computer Vision and Pattern Recognition, San Francisco, USA, IEEE Computer Society Press, 2010:2614-2621.

[3]Burner A, Donner R, Mayerhoefer M, et al. Texture bags:Anomaly retrieval in medical images based on local 3D-texture similarity[J]. Lecture Notes in Computer Science, 2012, 7075:116-127.

[4]Heckela F, Konradb O, Hahna H K, et al. Interactive 3D medical image segmentation with energy-minimizing implicit functions[J]. Computers & Graphics, 2011, 35(2): 275-287.

[5]DU Xin-wei, CHE Xiang-jiu. A hierarchical approach to 3D scattered data interpolation with radial basis functions[C]//12th International Conference on Computer-Aided Design and Computer Graphics. Jinan: IEEE Computer Society Press, 2011:262-267.

[6]Besl P J, Jain R C. Three-dimensional object recognition[J].ACM Computing Surveys, 1985, 17(1): 75-145.

[7]Bloomenthal J, Lim C. Skeletal methods of shape manipulation[C]//Proceedings of the International Conference on Shape Modeling and Applications. Washington: IEEE Computer Society Press, 1999: 44-47.

[8]TIAN Guang-feng. Three-dimension rheological model and parameters solving for unsaturated soil[J]. Advanced Materials Research, 2012, 368: 2698-2701.

[9]Aherne F, Thacker N, Rockett P. Optimal pairwise geometric histograms[C]//Proc British Machine Vision Conference Colchester: UK, 1997: 480-490.

[10]Park S. 2DPCA-based method for place classification using range scan[J]. Electronics Letters, 2011, 47(25): 1364-1366.

[11]ZHANG Yu, QI Mei-xing, LI Shang. Palmprint recognition based on two-dimensional gabor wavelet transform and two-dimensional principal component analysis[J]. Lecture Notes in Computer Science, 2012, 6838: 405-411.

[12]Osada R, Funkhouser T,Chazelle B, et al. Shape distribution[J].ACM Transactions on Graphics, 2000, 21(4): 807-832.

[13]Elad A, Kimmel R. On bending invariant signatures for surfaces[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2003, 25(10): 1285-1295.

[14]Sethian J A. Theory, algorithms and applications of level set method for propagating interfaces[J]. Acta Numerica, 1996, 5:309-395.

[15]Kimmel R, Sethian J A. Computing geodesic paths on manifolds[J]. Proc of National Academy of Science, 1998,95(15): 8431-8435.