周建存，常亮，郭克華

(1. 湖南城市學(xué)院 信息科學(xué)與工程學(xué)院，湖南 益陽(yáng)，413000；2. 中南大學(xué) 信息科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙，410083)

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

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

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

1.1 測(cè)地距離的獲取

等距曲面變形之后，表面相應(yīng)點(diǎn)之間的測(cè)地距離保持不變。因此，測(cè)地距離是等距變形體的基本特征，可以通過(guò)計(jì)算曲面上任意2點(diǎn)之間的測(cè)地距離來(lái)提取不變量。

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

1.2 特征矩陣的構(gòu)造

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

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

1.3 特征矩陣的歸一化

構(gòu)造出特征矩陣后，目標(biāo)之間的匹配轉(zhuǎn)化為特征矩陣之間的匹配。但是，由于曲面上像素點(diǎn)提取的順序是任意的，所以，同一個(gè)曲面可能構(gòu)造出不同的特征矩陣。這些特征矩陣對(duì)應(yīng)著一些同構(gòu)的帶權(quán)圖，若直接針對(duì)特征矩陣進(jìn)行匹配，則其復(fù)雜度等價(jià)于同構(gòu)圖的匹配，運(yùn)算時(shí)間為O(n!)。所以，必須對(duì)特征矩陣進(jìn)行歸一化，以保證同一個(gè)目標(biāo)得到的是同一個(gè)特征矩陣。

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

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

首先，計(jì)算每個(gè)頂點(diǎn)的均值和方差：

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

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

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

對(duì)于K(i)的構(gòu)造，實(shí)際上基于每個(gè)頂點(diǎn)的均值和方差，將頂點(diǎn)按升序進(jìn)行排列。而對(duì)特征矩陣中的頂點(diǎn)進(jìn)行的排序，可以保證特征矩陣的唯一性。

然后，構(gòu)建變換矩陣T來(lái)存儲(chǔ)頂點(diǎn)順序的調(diào)整情況，T為1個(gè)n×n矩陣，滿足：

T中其他元素為0。

最后，構(gòu)造歸一化后的特征矩陣NS。該矩陣可通過(guò)下式計(jì)算：

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

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

Step 1: 利用式(3)計(jì)算每個(gè)頂點(diǎn)的均值和方差。

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

Step 3：基于Step 2中計(jì)算的個(gè)頂點(diǎn)對(duì)應(yīng)的K(i)，根據(jù)式(6)構(gòu)建變換矩陣T，確定對(duì)特征矩陣進(jìn)行歸一化的變換矩陣。

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

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

2 等距變形體的變形矩及其性質(zhì)

歸一化特征矩陣的特征可以由許多方法來(lái)提取，如傅里葉描述子、矩特征等。本文將對(duì)矩不變量進(jìn)行擴(kuò)展，以提取NS的特征。

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

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

顯然，由于對(duì)特征矩陣進(jìn)行了歸一化，該矩陣中頂點(diǎn)的順序進(jìn)行了調(diào)整，中心矩vpq具備了平移不變性。

為保證尺度變換下的不變性，需要對(duì)中心矩進(jìn)行歸一化，為此，構(gòu)造出歸一化的中心矩μpq：

根據(jù)已經(jīng)構(gòu)造出來(lái)的歸一化的中心矩μpq，進(jìn)一步構(gòu)造2個(gè)針對(duì)等距變形體的變形矩：

要使全速變形矩得到應(yīng)用，必須具有不變性。對(duì)于前面秘構(gòu)造的變形矩，下面從平移、尺度和旋轉(zhuǎn) 3個(gè)方面闡述其不變性質(zhì)。

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

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

需指出的是：旋轉(zhuǎn)變換不會(huì)改變測(cè)地距離的值。由于歸一化特征矩陣中頂點(diǎn)的順序進(jìn)行了調(diào)整，因此，旋轉(zhuǎn)變換只會(huì)改變特征矩陣中頂點(diǎn)的順序，但不會(huì)改變歸一化特征矩陣。因此，式(11)構(gòu)造的變形矩具有旋轉(zhuǎn)不變性。

3 實(shí)驗(yàn)結(jié)果和復(fù)雜度分析

因?yàn)樯鲜隼碚撌腔谌S變形曲面的，所以，必須首先獲取圖像，然后構(gòu)造三角網(wǎng)格。目前，對(duì)于三維目標(biāo)圖像的數(shù)據(jù)獲取一般采用 2種方法：(1)由于二維圖像的存儲(chǔ)比較方便，可以通過(guò)從二維圖像進(jìn)行三維建模來(lái)獲取三維數(shù)據(jù)；(2)直接使用可以獲取三維位置信息的工具，如三維掃描儀等。在實(shí)驗(yàn)中，選取幾個(gè)比較典型的三維物體來(lái)進(jìn)行驗(yàn)證。

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

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

圖1 數(shù)據(jù)庫(kù)中的20個(gè)目標(biāo)Fig.1 20 Objects in Database

為了采用更多的樣本，對(duì)于每個(gè)目標(biāo)，對(duì)其實(shí)施任意的等距變形、旋轉(zhuǎn)和縮放變換，獲得10幅圖像?？倲?shù)據(jù)庫(kù)中有200幅三維目標(biāo)圖像。

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

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

表1 2種算法的計(jì)算時(shí)間比較Table 1 Computational time of two algorithms ms

從表1可看出：本文算法在運(yùn)算方面具有較低的運(yùn)算時(shí)間，因而復(fù)雜度低。

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

表2 2種算法的計(jì)算復(fù)雜度Table 2 Computational complexity of two algorithms

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

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

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

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

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

從表3可看出：本文算法在預(yù)處理方面具有較低的運(yùn)算復(fù)雜度，但是，并沒(méi)有降低識(shí)別率。

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

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

從圖2可以看出：本文算法和傳統(tǒng)算法對(duì)噪聲都具有較好的魯棒性。實(shí)際上，因?yàn)檫@2種方法都是基于距離進(jìn)行計(jì)算的，因此，噪聲對(duì)計(jì)算結(jié)果的影響不是太大。

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

4 結(jié)論

(1)針對(duì)等距變形這種特殊的變形體識(shí)別問(wèn)題進(jìn)行研究，提出了該變形情況下目標(biāo)識(shí)別的新方法。該方法通過(guò)構(gòu)造不變量來(lái)描述等距變形體的特征，首先利用測(cè)地距離的概念構(gòu)造出等距變形體的特征矩陣，然后對(duì)特征矩陣進(jìn)行歸一化，最后利用矩不變量對(duì)歸一化特征矩陣的特征進(jìn)行提取，構(gòu)造出相應(yīng)的矩不變量，并對(duì)該類矩不變量的平移、尺度和旋轉(zhuǎn)不變性進(jìn)行了證明。

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

(3)下一步工作是對(duì)有遮擋情況下的等距變形體識(shí)別問(wèn)題以及任意拓?fù)渥冃吻闆r下的目標(biāo)識(shí)別問(wèn)題進(jìn)行研究。

[1]文沁, 汪增福. 基于三維數(shù)據(jù)的人臉表情識(shí)別[J]. 計(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.