姚 斌， 何立風(fēng)，2， 康世英， 趙 曉

(1.陜西科技大學(xué) 電子信息與人工智能學(xué)院， 陜西 西安 710021； 2.日本愛知縣立大學(xué) 信息科學(xué)學(xué)院， 日本愛知縣長(zhǎng)久手市 4801198； 3.咸陽師范學(xué)院 計(jì)算機(jī)學(xué)院， 陜西 咸陽 712000)

0 引言

圖像的歐拉特征是圖像重要的拓?fù)湫再|(zhì)之一，用來描述圖像的結(jié)構(gòu)，它與圖像的幾何形狀無關(guān)，并且不會(huì)因?yàn)閳D像的擴(kuò)大、縮小、旋轉(zhuǎn)及變形而改變，有很強(qiáng)的魯棒性.二維圖像的歐拉特征稱為歐拉數(shù)(Euler Number)，它被定義為該圖像中連通域數(shù)量與孔洞數(shù)量之差[1]，被廣泛應(yīng)用于醫(yī)學(xué)診斷、字符識(shí)別以及地質(zhì)勘探中土壤內(nèi)部孔洞連通性分析等方面[2-4].三維圖像的歐拉特征也被稱為三維圖像的歐拉-龐加萊特征(Euler-Poincare Characteristic)或簡(jiǎn)稱三維歐拉數(shù).三維歐拉數(shù)分為三維表面歐拉數(shù)和三維體歐拉數(shù)，對(duì)于確定三維圖像的內(nèi)部結(jié)構(gòu)，三維體歐拉數(shù)會(huì)更有效(在后續(xù)篇幅中，三維圖像歐拉數(shù)特指體歐拉數(shù)).

三維圖像歐拉數(shù)經(jīng)常被用來描述物體的基本特征，例如用來描述結(jié)構(gòu)的連接性等[5]，被廣泛用于解決實(shí)際問題.比如，H.J.Vogel等[6]利用三維圖像歐拉數(shù)確定土壤孔隙結(jié)構(gòu)；A.Pierret等[7]利用三維圖像歐拉數(shù)來重建和量化macropores.P.Lehmann等[8]分析了幾何形態(tài)對(duì)水流和流體分布的影響，他們利用三維圖像歐拉數(shù)來預(yù)測(cè)多孔介質(zhì)中的流體相分布.在工業(yè)生產(chǎn)方面，A.Velichko等[9]利用三維圖像歐拉數(shù)區(qū)分鑄鐵中的石墨形態(tài)，使人們更容易理解鑄鐵增長(zhǎng)機(jī)制和鑄鐵性能.在醫(yī)學(xué)方面，L.Kubínová等[10]用三維圖像歐拉數(shù)來測(cè)量毛細(xì)血管的長(zhǎng)度.由于歐拉數(shù)多用于物體的快速識(shí)別與分類，歐拉數(shù)的計(jì)算速度越快，就越有利于提高相關(guān)實(shí)時(shí)圖像處理系統(tǒng)的整體性能.

1 三維圖像像素的鄰接性及三維圖像歐拉數(shù)

1.1 三維圖像像素的鄰接性

在本文中，圖像特指二值圖像，按照慣例，用“1”表示二值圖像中的前景像素(目標(biāo)像素)，而用“0”表示背景像素(非目標(biāo)像素).

在對(duì)三維圖像的研究中，像素之間的鄰接性可以分為6-鄰域連接和26-鄰域連接[11].對(duì)于兩個(gè)不同的像素X(x1,x2,x3)和Y(y1,y2,y3)，當(dāng)|x1-y1|+|x2-y2|+|x3-y3|=1時(shí)，稱像素X與像素Y為6-鄰域連接.如圖1所示，像素p的6-鄰域連接像素有p1、p3、p5、p7、p17和p26.而當(dāng)max(|x1-y1|,|x2-y2|,|x3-y3|)=1時(shí)，稱像素X與像素Y為26-鄰域連接.在圖1中，像素p的26-鄰域連接像素為p1，p2，…，p26.

1.2 三維圖像歐拉數(shù)

三維圖像的歐拉數(shù)被定義為：E=O-H+C.其中,O、H和C分別為圖像的連通域、孔洞和空腔的數(shù)量.三維圖像中“孔洞”定義具有多義性與不確定性，在文獻(xiàn)[12]中，C.N.Lee將“孔洞”解釋為不產(chǎn)生分離的最大切割數(shù)，也被稱為“tunnel”，而“空腔”則被稱為“cavity”.由于直接利用定義來計(jì)算三維圖像的歐拉數(shù)需要獲得圖像中連通域、孔洞和空腔的數(shù)量，也就需要對(duì)圖像有全面了解，無法由圖像的局部性質(zhì)計(jì)算得到，因此不能將歐拉數(shù)的計(jì)算分解成每個(gè)網(wǎng)格片的子任務(wù)，無法進(jìn)行局部計(jì)算.

圖1 三維圖像中像素的鄰接性示意圖

除了利用三維圖像歐拉數(shù)定義中的公式，三維圖像歐拉數(shù)也可以通過局部計(jì)算的方法得到.由局部計(jì)算三維圖像歐拉數(shù)的算法最早是在1971年由C.M.Park等[11]的論文提出的，但其算法只對(duì)三維圖像像素間的6-鄰域連接關(guān)系有效，無法計(jì)算圖像像素間是26-鄰域連接關(guān)系時(shí)的圖像歐拉數(shù).1980年，D.G.Morgenthaler[13]提出一種對(duì)三維圖像像素間26-鄰域連接關(guān)系有效的歐拉數(shù)局部算法，該算法需要檢查圖像中每個(gè)2×2×2立方體(稱為“三維立方體模式”，有256種可能的情況)的像素組成情況，而且被后續(xù)的研究所沿用[14].這種方法計(jì)算歐拉數(shù)需要的像素讀取次數(shù)為圖像像素總數(shù)的8倍.在后續(xù)篇幅中為了表示方便，把該算法稱為DGM算法.

1992年，國內(nèi)學(xué)者楊敬安等[15]提出了利用微分幾何及代數(shù)拓?fù)湓碛?jì)算三維圖像歐拉數(shù)的方法并給出了相應(yīng)算法.1996年P(guān).K.Saha等[16]提出了計(jì)算三維圖像歐拉數(shù)的新方法，該方法是通過計(jì)算圖像中3×3×3鄰域內(nèi)的連通域、孔洞和空腔在數(shù)量上的變化來實(shí)現(xiàn)的.由于需要同時(shí)考慮到圖像中連通域、孔洞和空腔的數(shù)量，效率比單純計(jì)算歐拉數(shù)低一些.2010年，國內(nèi)學(xué)者林小竹等[17]提出了計(jì)算三維圖像歐拉數(shù)的新公式，利用圖像中圖段數(shù)及相鄰圖段數(shù)實(shí)現(xiàn)了三維圖像歐拉數(shù)的計(jì)算.這種方法將二維圖像中的方法推廣應(yīng)用到三維圖像中，對(duì)于密度較低的圖像效率較高.2013年，H.Snchez Cruz等[18]提出了一種新的計(jì)算三維圖像歐拉數(shù)算法，利用圖像幾何方面的表面連接性，通過計(jì)算2×2×2三維基本立方體和2×2×1、2×1×2、1×2×2立方體的數(shù)量來實(shí)現(xiàn)計(jì)算.與文獻(xiàn)[13]中提出的算法類似，采用這種方法計(jì)算給定三維圖像歐拉數(shù)需要的像素存取次數(shù)也是該圖像中像素總數(shù)的8倍.

2 DGM算法

在DGM算法中，通過統(tǒng)計(jì)給定三維圖像中包含的特定2×2×2的三維立方體模式的數(shù)量來計(jì)算圖像歐拉數(shù)，如圖2所示.

圖2 DGM算法用到的立方體模式

如果用qn表示給定三維圖像中含有n個(gè)前景像素的2×2×2三維立方體的個(gè)數(shù)，例如，q3是含有3個(gè)前景像素的三維基本立方體的個(gè)數(shù)，即圖2中(6)、(7)和(8)所示三種模式數(shù)量總和，則給定三維圖像在26-鄰域連接情況的下的歐拉數(shù)E26的計(jì)算公式為：E26=q1-q2+q3-q4+q5-q6+q7-q8.由于圖2中列出的三維立方體模式是根據(jù)其中前景像素的個(gè)數(shù)及對(duì)稱性分類的，在處理圖像時(shí)需要從各個(gè)方向考慮，實(shí)際計(jì)算中使用起來并不方便.

為了便于計(jì)算，D.G.Morgenthaler通過對(duì)256種2×2×2的三維基本立方體模式的立體特征進(jìn)行研究、分析，對(duì)其中八個(gè)像素a、b、c、d、e、f、g和h(如圖3所示)的取值情況從00000000到11111111進(jìn)行了分析，計(jì)算出256種模式中每種模式對(duì)于圖像歐拉數(shù)的增量，并給出了相應(yīng)的查找表.對(duì)照查找表進(jìn)行分析，可以看出，在256種模式中，有207種模式的歐拉數(shù)增量為0，30種為-1，17種為1，其余兩種為-2[8].對(duì)于圖像歐拉數(shù)的增量不為0的49種三維立方體模式的編號(hào)及歐拉數(shù)增量△E如表1所示.

表1 DGM算法中對(duì)圖像歐拉數(shù)有貢獻(xiàn)的模式

在實(shí)際進(jìn)行三維圖像歐拉數(shù)的計(jì)算時(shí)，只需要對(duì)圖像進(jìn)行光柵掃描，對(duì)于掃描到的像素，只需要檢查、判斷該像素所在的三維立方體模式屬于256種模式中的哪一種，然后查表即可得到該立方體模式對(duì)圖像歐拉數(shù)的增量值.所有像素掃描完成后，三維圖像中所包含的所有三維立方體模式都被檢查過，所有的三維立方體模式的歐拉數(shù)增量的總和即是給定三維圖像的歐拉數(shù).

如圖3所示，用DGM算法計(jì)算三維圖像歐拉數(shù)時(shí)，需要對(duì)整個(gè)圖像的像素進(jìn)行掃描.對(duì)掃描到的每個(gè)像素，需要檢查其所對(duì)應(yīng)的三維立方體模式中的其他七個(gè)像素才能得到當(dāng)前三維基本立方體對(duì)應(yīng)的歐拉數(shù)增量.例如，對(duì)于當(dāng)前掃描到的像素a，需要檢查像素a所在的三維立方體模式{a,b,c,d,e,f,g,h}，這樣，處理每個(gè)像素一共需要檢查八個(gè)像素.

圖3 DGM算法中檢查三維立方體的像素

3 算法的改進(jìn)

在DGM算法中，為了計(jì)算給定三維圖像的歐拉數(shù)，需要檢查圖像中所包含的所有立方體模式.盡管需要統(tǒng)計(jì)的立方體模式只有49種，但對(duì)當(dāng)前檢查的立方體模式，還需要檢查8個(gè)像素才能確定該立方體模式是不是需要統(tǒng)計(jì).對(duì)于一個(gè)大小為X×Y×Z像素的圖像，為了計(jì)算圖像歐拉數(shù)，像素存取次數(shù)為8×X×Y×Z.下面，本文從兩個(gè)方面對(duì)DGM算法進(jìn)行改進(jìn).

3.1 改變像素檢查順序

DGM算法中，盡管只有一部分立方體模式對(duì)圖像歐拉數(shù)有貢獻(xiàn)，但是為了找出對(duì)圖像歐拉數(shù)有貢獻(xiàn)的立方體模式，還是需要對(duì)圖像中包含的所有立方體模式進(jìn)行檢查.如圖3所示，檢查每一個(gè)立方體模式所包含的a、b、c、d、e、f、g和h八個(gè)像素的值.另一方面，既然只有49種立方體模式對(duì)圖像歐拉數(shù)有貢獻(xiàn)，可以把這些立方體模式的每一個(gè)像素值列出來，進(jìn)行分析，尋找其中規(guī)律，盡可能減少不必要的像素存取以提高算法效率.

將表1中所有三維立方體模式編號(hào)用二進(jìn)制數(shù)表示后(從高位到低位依次為a、b、c、d、e、f、g和h)，可以得到表2.

表2 DGM算法中對(duì)圖像歐拉數(shù)有貢獻(xiàn)的模式

在49種對(duì)圖像歐拉數(shù)有影響的立方體模式中，本文對(duì)所有立方體模式編號(hào)的二進(jìn)制編碼進(jìn)行分析，發(fā)現(xiàn)編碼中數(shù)字“1”出現(xiàn)的位置是有差異的.

經(jīng)過仔細(xì)分析后，發(fā)現(xiàn)數(shù)字“1”在所有立方體模式編號(hào)的二進(jìn)制編碼中“b”的位置都沒有出現(xiàn).根據(jù)立方體模式二進(jìn)制編碼的這個(gè)特點(diǎn)，在檢查當(dāng)前立方體模式的八個(gè)像素值時(shí)，可以先檢查像素“b”的值，如果“b”的值為1，由表2得知，當(dāng)前立方體模式對(duì)圖像歐拉數(shù)沒有影響，因此可以跳過該立方體模式的檢查，直接檢查下一個(gè)立方體模式.這樣，對(duì)于前景像素較多的圖像，在特定情況下(像素b為前景像素)只需要檢查一個(gè)像素就可以確定該立方體模式是不是對(duì)圖像歐拉數(shù)有影響，這樣可以節(jié)省很大一部分時(shí)間，算法效率也隨之提高.

因此，在算法實(shí)現(xiàn)時(shí)，本文將檢查像素的順序由原來的“abcdefgh”改變?yōu)椤癰acdefgh”.盡管只改變了兩個(gè)像素的檢查順序，在檢查圖像中包含的立方體模式時(shí)，部分立方體模式只需要檢查一個(gè)像素就可以確定該模式是否對(duì)圖像歐拉數(shù)有影響.

3.2 合并部分立方體模式

對(duì)表1進(jìn)行進(jìn)一步分析，發(fā)現(xiàn)在對(duì)圖像歐拉數(shù)有影響的49種立方體模式中，部分立方體模式的編號(hào)是連續(xù)的，如編號(hào)24-27，36-39，44-47，56-59，148-151，156-159，180-183和188-191的立方體模式，這些編號(hào)連續(xù)的立方體模式還有個(gè)共同點(diǎn)，就是對(duì)圖像的歐拉數(shù)增量相同，如表3所示.

表3 可以合并的立方體模式

由表3可以得到，對(duì)于編號(hào)為24-27的四種立方體模式，這些立方體模式中a、b、c、d、e、f、g和h八個(gè)像素形成的二進(jìn)制編碼的值為00011000-00011011.不難發(fā)現(xiàn)，這四個(gè)編碼的八位數(shù)值中只有最后兩位不同(00-11)，前六位數(shù)值是相同的(000110).與此同時(shí)，可以發(fā)現(xiàn)這些立方體模式對(duì)圖像歐拉數(shù)的增量相同，都是“-1”.因此，在算法實(shí)現(xiàn)時(shí)，可以把這四種立方體模式合并，檢查立方體模式的像素時(shí)，只需要檢查前六個(gè)像素的值就可以確定該模式對(duì)圖像歐拉數(shù)的影響.

對(duì)于當(dāng)前正在處理的立方體模式，當(dāng)檢查其像素串“abcdef”的值為“000110”時(shí)，就可以確定當(dāng)前立方體模式是編號(hào)為24-27的模式中的某一個(gè)，因此可以直接將圖像歐拉數(shù)增加“-1”，不需要再檢查最后兩個(gè)像素的值.同理，當(dāng)檢查立方體模式時(shí)，如果二進(jìn)制編碼“abcdef”的值為“001001”、“001011”或“001110”時(shí)，圖像歐拉數(shù)增加“-1”.

而對(duì)于正在處理的立方體模式，當(dāng)檢查像素串“abcdef”的值為“100101”、“100111”、“101101”或“101111”時(shí)，同樣可以不用檢查最后兩個(gè)像素的值，直接可以將圖像歐拉數(shù)增加1.這樣，同樣可以檢查六個(gè)像素就可以判斷一個(gè)立方體模式的類型，與DGM算法相比，在處理部分立方體模式時(shí)可以減少像素存取次數(shù)，同樣可以提高算法效率.

算法在實(shí)際實(shí)現(xiàn)時(shí)，將改變像素檢查順序和合并立方體模式兩種改進(jìn)方法相結(jié)合，對(duì)于當(dāng)前正在處理的立方體模式，首先將像素檢查順序更改為“bacdefgh”，如果像素“b”的值為“1”，直接跳過該立方體模式，檢查下一個(gè)；如果正在處理的立方體模式中像素“b”的值為“0”；再檢查其他像素的值.如果前六個(gè)像素的值是表3所列的八組之一，則由前六個(gè)像素的值就能確定正在處理的立方體模式對(duì)圖像歐拉數(shù)的影響，盡可能減少像素存取次數(shù)，提高算法效率.

4 實(shí)驗(yàn)結(jié)果及分析

為了測(cè)試提出的算法的正確性和有效性，本文采用結(jié)構(gòu)較為復(fù)雜的噪聲類型的圖像對(duì)算法進(jìn)行測(cè)試，并把本文提出算法(Ours)的運(yùn)行結(jié)果與DGM算法進(jìn)行對(duì)比.

本文采用41幅噪聲圖像作兩種算法的圖像密度與程序執(zhí)行時(shí)間的關(guān)系對(duì)比實(shí)驗(yàn).噪聲圖像的密度指圖像中前景像素的比例，密度越大，前景像素所占比例越大.當(dāng)圖像密度為1時(shí)，所有像素均為前景像素；反之，當(dāng)圖像密度為0時(shí)，所有像素均為背景像素.不同密度的三維圖像截面如圖4所示.

41幅噪聲圖像采用門限法隨機(jī)生成，閾值從0到1 000，步長(zhǎng)值為25，大小為128×128×128像素.測(cè)試平臺(tái)為臺(tái)式機(jī)(Intel Core i7-6700 CPU@3.40 GHz，8 GB內(nèi)存，操作系統(tǒng)Windows 7)，編譯器為GNU C compiler(Version 4.6.1).為保證實(shí)驗(yàn)結(jié)果的準(zhǔn)確性，所有計(jì)算圖像的歐拉數(shù)的程序均被運(yùn)行1 000次，取平均值作為實(shí)驗(yàn)結(jié)果.

圖4 不同密度的三維圖像截面

為了驗(yàn)證本文改進(jìn)算法的正確性，分別采用兩種算法計(jì)算給定圖像的歐拉數(shù).實(shí)驗(yàn)結(jié)果表明用兩種算法計(jì)算給定圖像的歐拉數(shù)結(jié)果是相同的.表4給出了用兩種算法計(jì)算時(shí)部分圖像的歐拉數(shù)值.

表4 兩種算法計(jì)算圖像歐拉數(shù)的結(jié)果

算法效率實(shí)驗(yàn)結(jié)果如圖5所示.從圖5可以看出，幾乎對(duì)于所有圖像，本文提出的算法的效率都高于DGM算法，特別是在圖像密度值在0.4～0.8之間的圖像，本文提出算法的優(yōu)勢(shì)更為明顯.在本文提出的算法中，對(duì)于當(dāng)前正在檢查的立方體模式，首先檢查像素“b”的值，如果該像素的值為“0”，然后再檢查其他像素的值；如果像素“b”的值為“1”，則可以跳過當(dāng)前立方體模式，直接檢查下一個(gè).在密度比較小的圖像中，圖像中包含的“1”的個(gè)數(shù)較少，因此，本文提出算法的效率和DGM算法相比優(yōu)勢(shì)不是很明顯.從圖5可以看出，對(duì)于密度接近0的圖像，兩種算法計(jì)算歐拉數(shù)的時(shí)間幾乎相等.隨著圖像中前景像素的增加，本文提出算法的優(yōu)勢(shì)逐漸增大.另一方面，計(jì)算圖像歐拉數(shù)所需要的時(shí)間除了對(duì)像素進(jìn)行讀取、判斷之外，還需要進(jìn)行一些計(jì)算，前景像素越多，需要的計(jì)算量也越多，因此當(dāng)圖像中前景像素的數(shù)量增加到一定程度后兩種算法效率之間的差異反而逐漸縮小了.

圖5 兩種算法的運(yùn)行時(shí)間對(duì)比圖

5 結(jié)論

為了解決三維二值圖像歐拉數(shù)的快速計(jì)算問題，本文在現(xiàn)有算法的基礎(chǔ)上，通過改變立方體模式的像素檢查順序、合并部分立方體模式兩種方法，盡可能減少處理一個(gè)立方體模式所需要檢查的像素個(gè)數(shù)，提高算法效率.實(shí)驗(yàn)結(jié)果表明，本文提出的算法的計(jì)算結(jié)果與現(xiàn)有算法得到的結(jié)果相同，對(duì)于絕大多數(shù)被測(cè)圖像來說，本文提出算法的性能優(yōu)于現(xiàn)有算法.

在算法實(shí)現(xiàn)過程中，發(fā)現(xiàn)在對(duì)圖像進(jìn)行逐行掃描、檢查每個(gè)三維立方體模式來計(jì)算圖像歐拉數(shù)時(shí)，仍然存在像素重復(fù)檢查的現(xiàn)象.所以在接下來的工作中，筆者將對(duì)算法進(jìn)行優(yōu)化和改進(jìn)，以期得到性能更好的二值圖像歐拉數(shù)算法.