王忠勇, 賈　萌, 侯中新, 陳在旭

(1.鄭州大學(xué) 信息工程學(xué)院,河南 鄭州 450001； 2.河南高遠公路養(yǎng)護設(shè)備股份有限公司, 河南 新鄉(xiāng) 453000)

基于形態(tài)學(xué)特征的顆粒圖像分割和計數(shù)

王忠勇1, 賈萌1, 侯中新2, 陳在旭2

(1.鄭州大學(xué) 信息工程學(xué)院,河南 鄭州 450001； 2.河南高遠公路養(yǎng)護設(shè)備股份有限公司, 河南 新鄉(xiāng) 453000)

摘要：針對黏連嚴重、分辨率和清晰度較低的顆粒圖像分割和計數(shù)問題，提出了一種基于分水嶺和形態(tài)學(xué)特征的新方法.首先使用分水嶺算法分割圖像，得到過分割的結(jié)果，接著通過本文定義的區(qū)域形態(tài)學(xué)特征，再根據(jù)加權(quán)馬氏距離和區(qū)域連接圖(Region Adjacency Graph)進行區(qū)域合并.采用實際圖像進行的實驗表明，該算法效果良好，對原始圖像的灰度、對比度和噪聲變化具有不變性，在準確率-查全率(Precision-recall)曲線的表現(xiàn)上優(yōu)于現(xiàn)有方法.

關(guān)鍵詞：顆粒分割計數(shù)；分水嶺算法；形態(tài)學(xué)特征；區(qū)域連接圖

0引言

基于圖像分割的物體計數(shù)有很多實際應(yīng)用，例如谷物顆粒計數(shù)，微生物和細胞計數(shù)，基于截面圖像的棒材計數(shù)等[1].大部分是依靠分水嶺對計數(shù)目標(biāo)進行分割.例如結(jié)合分水嶺算法和歐拉數(shù)的計數(shù)方法[2]，基于距離變換[3]和分水嶺算法分割二值化的黏連顆粒[4]是目前應(yīng)用最廣泛的方法.

實際工程中的成像條件往往十分復(fù)雜，在同步碎石封層的公路施工中，常常需要檢測石子撒布覆蓋率，目前僅能依靠人工計數(shù)樣本方格內(nèi)的石子數(shù)目，效率十分低下.筆者針對清晰度差，分辨率不足的圖像，提出了一種魯棒性較強的分割和計數(shù)方法.

由于分水嶺算法能夠響應(yīng)非常微弱的邊緣，形成的邊緣連續(xù)且閉合，因此該算法針對成像質(zhì)量差的圖像有效.基于不同的實現(xiàn)思想，分水嶺算法可分為模擬漲水法[5]、模擬降雨法[6]和模擬洪水消退法[7].分水嶺分割算法也存在過分割問題，抑制過分割的方法一是預(yù)處理，采用濾波抑制原始圖像中的偽極小值，同時保護邊緣[8]；二是基于標(biāo)記抑制過分割[9]，只有標(biāo)記之間最強的分水線會被返回，三是后處理，即先得到過分割的結(jié)果，再進行區(qū)域合并[8].

本方法是一種后處理的技術(shù).筆者充分發(fā)揮了分水嶺算法將圖像看作三維地形的思想，定義了若干形態(tài)學(xué)的特征來描述區(qū)域，并提出了根據(jù)特征合并區(qū)域的方法，以及停止規(guī)則.

1算法流程概述

本算法的流程如圖1所示.首先對原始的顆粒圖像進行高斯平滑預(yù)處理，這是為了減輕過分割的程度(而不是消除過分割)，以減小區(qū)域合并的運算量.

其次采用分水嶺算法對其負圖像進行分割.圖像分割操作可以看作對像素分配標(biāo)簽，因此，稱分水嶺分割的結(jié)果為label圖，它的尺寸與原始圖像img相同，存儲著對相應(yīng)像素分配的標(biāo)簽(即區(qū)域標(biāo)號).區(qū)域標(biāo)號從1開始，每次發(fā)現(xiàn)新的局部極小值，就把標(biāo)號加1，賦給該極小值像素，漲水過程中，對屬于同一集水盆地的像素分配號碼n，并稱該區(qū)域為區(qū)域n.當(dāng)漲水造成兩個區(qū)域?qū)⒁噙B時,對此處的像素標(biāo)簽賦值0,以標(biāo)志該像素是邊緣.初始化形態(tài)學(xué)特征和區(qū)域連接圖的工作可以在漲水過程中同步地完成，以減少運算量.

最后進行區(qū)域合并.以分割結(jié)果、形態(tài)學(xué)特征和區(qū)域連接圖為輸入數(shù)據(jù)，循環(huán)執(zhí)行，每一次循環(huán)可完成一對相鄰區(qū)域的合并，直至算法認定所有錯誤區(qū)域合并殆盡.

2形態(tài)學(xué)特征和區(qū)域連接性圖

過分割會將圖像分為大量細小的區(qū)域，如果一個區(qū)域里包括且僅包含一個顆粒，筆者定義其為正確分割的區(qū)域.如果一個區(qū)域不包括顆?；騼H僅是顆粒的一部分，筆者定義其為過分割的錯誤區(qū)域.圖2為上述兩種情形的例子，圖2(a)表現(xiàn)了第一類不包括顆粒錯誤區(qū)域，圖2(b)表現(xiàn)了第二類僅僅是顆粒的一部分.顆粒稀疏時傾向于發(fā)生第一類錯誤，顆粒密集時發(fā)生第二類.

本方法引入深度、面積和容量的幾何概念，定義了6種形態(tài)學(xué)特征描述區(qū)域，并且在區(qū)域鏈接性圖中使用區(qū)域之間分水嶺的絕對高度衡量兩區(qū)域的連接強度，以指導(dǎo)區(qū)域合并.

2.1區(qū)域形態(tài)學(xué)特征

如圖3所示，該地形圖的截面上有兩個局部極小值，高度分別是hi和hj，因此初始的分割將在這里產(chǎn)生區(qū)域i和區(qū)域j. 點A和點B分別是區(qū)域i邊緣上的最高點和最低點，記為hmin(i)和hmax(i)，區(qū)域j可類推.

定義1定義在漲水過程中，水面上升至hmin時的面積為盆地面積，單位為平方像素，記為Sbasin(如圖3(a)中穿過B點的水平線所示)，

(1)

式中:δ和Sgn分別是狄拉克函數(shù)和符號函數(shù).

定義2定義漲水至hmax時，水面面積為區(qū)域面積，記為Sregion(如圖3(a)中穿過A點和C點的水平線所示)，

(2)

定義3定義漲水至hmin時，水池的深度為盆地深度，單位為像素，記為Hbasin(即水面面積為Sbasin時的深度)，

Hbasin(i)?hmin(i)-hi.

(3)

定義4定義漲水至hmax時，水池的深度為區(qū)域深度，記為Hregion(即水面面積為Sregion時的深度)，

Hregion(i)?hregion(i)-hi.

(4)

定義5定義漲水至hmin時，蓄水的體積為區(qū)盆地蓄水量，單位為立方像素，記為Vbasin(即Sbasin以下的體積)，

(5)

定義6定義漲水至hmax時，蓄水的體積為區(qū)域蓄水量，記為Vregion(Sregion以下的體積)，

(6)

2.2區(qū)域連接圖

區(qū)域連接圖[8,10]用于記錄區(qū)域之間的連接關(guān)系和連接強度.它是一個無向圖，記為G=(N,E).其中N={1,2，…，K},是節(jié)點的集合，代表區(qū)域，K是初始分割得到區(qū)域的數(shù)目，E∈RK×K是節(jié)點間連接的集合，代表區(qū)域間的邊緣.

圖4為一個區(qū)域連接圖的例子，其中(a)模擬了一個包括6個子區(qū)域的分割;(b)用節(jié)點和連接的形式將它們表示出來;(c)是用于保存區(qū)域連接圖的矩陣E，它的元素是布爾值，表示邊緣，元素在矩陣中的位置表示所連接的節(jié)點.

除此之外，合并過程中還需要一個K×1的矩陣Validity，元素也是布爾值.在區(qū)域合并過程中，矩陣E的尺寸不會隨著合并而動態(tài)改變，如果一個區(qū)域被合并而消失掉了，它會在Validity中被標(biāo)記為無效，E中的對應(yīng)行列也同樣失效.

Strength是與E同尺寸的矩陣，其元素代表E中對應(yīng)連接的連接強度.筆者采用分水嶺的高度衡量2個區(qū)域之間的連接強度，分水嶺的高度就是2個區(qū)域之間的邊緣平均高度，

(7)

區(qū)域連接圖節(jié)點的合并如圖4(d)所示，所有與區(qū)域a和區(qū)域b有連接的區(qū)域，在合并后都將與新區(qū)域有連接.以Ai表示與區(qū)域i連接的區(qū)域集合

(8)

3合并算法

合并過程是循環(huán)執(zhí)行的，每次循環(huán)使2個區(qū)域的合并，每次循環(huán)中需要完成：①搜索出一個錯誤區(qū)域i；②從與i相連的區(qū)域中選出連接強度最大的區(qū)域j；③執(zhí)行i與j的合并；④更新區(qū)域特征和區(qū)域連接圖；⑤檢測停止規(guī)則.

3.1基于加權(quán)馬氏距離搜索錯誤區(qū)域

如前所述，錯誤區(qū)域的面積、深度和蓄水量較小，因此首先統(tǒng)計所有區(qū)域形態(tài)學(xué)特征的均值向量M，然后找出距離均值最遠(在負方向上)的一個樣本，即為待合并的區(qū)域i.

顯然，筆者定義的6種形態(tài)學(xué)特征不但量綱和數(shù)值變化范圍不同，而且不同的形態(tài)學(xué)特征在描述區(qū)域時的重要性也不同.采用馬氏距離可以消除量綱的影響，但不能反應(yīng)重要性差別；加權(quán)歐氏距離沒有考慮特征之間的相關(guān)性和量綱，因此，一種直觀的想法是對馬氏距離進行加權(quán).

馬氏距離的計數(shù)公式為

(9)

式中：Σ是樣本協(xié)方差矩陣.而所謂的加權(quán)馬氏距離[11]計算公式為

(10)

式(10)所示的加權(quán)馬氏距離滿足距離定義的三條基本公理，且對一切非奇異線性變換具有不變性(說明可以排除量綱影響),但是加權(quán)馬氏距離不能嚴格地反映特征之間的相關(guān)性.

如果樣本在xy平面上呈現(xiàn)二維高斯分布，歐氏距離的等距離曲線如圖5(a)所示，歐氏空間中到某點等距離的點集是圓(高維時，為球面或超球面).加權(quán)歐氏距離的等距離曲線如圖5(b)、(c)所示，它們是以均值為中心的橢圓簇.圖5(d)中顯示與樣本數(shù)據(jù)的分布有關(guān)的馬氏距離的等距離曲線.圖5(e)、(f)是加權(quán)馬氏距離的等距離曲線，這些曲線長軸的方向與樣本分布有關(guān)，同時也受到加權(quán)的影響.

加權(quán)馬氏距離是對馬氏距離等距離曲線進行拉伸或壓縮，可以看作是加權(quán)距離和馬氏距離的結(jié)合.筆者采用線性搜索的方法優(yōu)化權(quán)值.

3.2基于近鄰函數(shù)確定合并入的區(qū)域

搜索得到要合并的區(qū)域i之后，需要從區(qū)域i的相鄰區(qū)域中找到一個區(qū)域j，使區(qū)域i與區(qū)域j合并.確定區(qū)域j的依據(jù)是區(qū)域連接圖的連接強度(式7)，采用近鄰函數(shù)的方法.

間隔分水嶺越低的區(qū)域連接強度越大，合并后的損失越小.對于相連的任意2個區(qū)域i和j，如果i是j的第I個近鄰的區(qū)域，則定義區(qū)域i對區(qū)域j的近鄰系數(shù)為I，記為d(i，j)=I.于是區(qū)域i和區(qū)域j之間的近鄰函數(shù)值

Dij=d(i，j)+d(j，i)-2.

(11)

逐個檢查與區(qū)域i相連的區(qū)域的近鄰函數(shù)，找出最小值(最近鄰)，然后使其與區(qū)域i合并.

3.3停止規(guī)則

如果算法的任務(wù)僅僅是對顆粒圖像進行分割，停止規(guī)則可選取為：平均區(qū)域大小、區(qū)域數(shù)目、區(qū)域形狀因子等等，其參數(shù)可以根據(jù)分割的目的進行設(shè)置.如果算法的任務(wù)是自動地統(tǒng)計顆粒的數(shù)目，則需要選擇一種客觀的衡量標(biāo)準.

在顆粒圖像中，局部的顆粒分布大致是均勻的，據(jù)此，筆者采用了區(qū)域面積方差的停止規(guī)則.方差可以度量數(shù)據(jù)和其均值的偏離程度，初始的過分割區(qū)域雜亂無章，大小不一，其方差較大，隨著區(qū)域合并的進行，區(qū)域逐漸正確地與顆粒相一致，這時區(qū)域的面積差異應(yīng)該最小.

由于合并過程中，區(qū)域面積的均值不斷變化，為了使不同合并進程中的方差具有可比性，采取如下方式標(biāo)準化數(shù)據(jù)，

(12)

4實驗結(jié)果

圖6是采用實際的同步碎石封層施工形成的石子圖像實驗的結(jié)果，可見本算法在實施區(qū)域合并中，能夠在減少錯誤區(qū)域的同時避免將兩顆相鄰石子合并入同一區(qū)域，因而能對石子進行計數(shù).

為了與一些現(xiàn)有方法[2，4]對比，人工標(biāo)記出Groundtruth圖像，用色點標(biāo)示出了顆粒的位置.如果一個區(qū)域里包含色點，則該區(qū)域是一個正確分割，否則為錯誤分割，進而可以畫出準確率(Precision)隨返回區(qū)域數(shù)目(Recall)變化的曲線，如圖7(a)所示.直接分水嶺算法[2]僅通過預(yù)處理控制過分割，而基于距離變換方法[4]的準確性依賴于圖像二值化的結(jié)果，在圖像模糊、分辨率較低的數(shù)據(jù)上實驗，它們的性能相對較差.

為了證明該分割方法對灰度、對比度和噪聲變化的魯棒性，以及計數(shù)方法的不變性，分別使原始圖像經(jīng)歷Gamma變換、圖像平滑、圖像增強和加白噪聲，再應(yīng)用本算法，結(jié)果如表1所示.

對于不同的變換，表1第二列的參數(shù)意義不同：對圖像平滑，是高斯模板的標(biāo)準差；對Gamma變換，是γ值；對拉普拉斯模板圖像增強，其為增強系數(shù)；對加白噪聲，其為標(biāo)準差.

進一步將表1中8個實驗的準確率-查全率曲線繪制出來，如圖7(b)所示，該結(jié)果說明本方法對圖像降質(zhì)具有很強的魯棒性.

5結(jié)論

提出了采用分水嶺算法輔以形態(tài)學(xué)特征和區(qū)域連接強度為指導(dǎo)的區(qū)域合并方法.該方法能夠最大限度地發(fā)揮分水嶺算法對微弱邊緣的檢測能力，基于形態(tài)學(xué)特征的區(qū)域合并算法則能夠有效消除過分割.實驗結(jié)果證明，在成像質(zhì)量不高的場合，該算法較現(xiàn)有方法優(yōu)越.

參考文獻：

[1]徐蔚欽. 基于分水嶺分割的粘連顆粒圖像分析技術(shù)研究[D].長沙：中南大學(xué)信息科學(xué)與工程學(xué)院，2011.

[2]林小竹,王彥敏,杜天蒼,等.基于分水嶺變換的目標(biāo)圖像的分割與計數(shù)方法[J].計算機工程,2006,32(15): 181-183.

[3]BORGEFORS G. Distance transformations in digital images [J]. Computer Vision, Graphics, And Image Processing, 1986, 34(3):344-371.

[4]LIN H, TSAI M, CHEN C, et al. Cell segmentation and NC ratio analysis of third harmonic generation virtual biopsy images based on marker-controlled gradient watershed algorithm[C]//Seoul, Korea (South)：2012 IEEE International Symposium on Circuits and Systems (ISCAS). 2012: 101-104.

[5]VINCENT L, SOILLE P. Watersheds in digital spaces: an efficient algorithm based on immersion simulations [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1991, 13(6): 583-598.

[6]SMET P D, PIRES R L, VLEESCHAUWER D, et al. Implementation and analysis of an optimized rain falling watershed algorithm[C]// Proceedings of SPIE. San Diego, USA: The International Society for Optical Engineering. 2000: 759-766.

[7]張立東, 畢篤彥.一種基于洪水消退模型的快速分水嶺算法[J].模式識別與人工智能, 2006, 19(3): 349-360.

[8]HARIS H, EFSTRATIADIS S, MAGLAVERAS N, et al. Hybrid image segmentation using watersheds and fast region merging [J]. IEEE Transactions on Image Processing, 1998, 7(12): 1684-1699.

[9]XU Lu-lu, LU Hua-xiang. Automatic morphological measurement of the quantum dots based on marker-controlled watershed algorithm[J]. IEEE Transactions on Nanotechnology, 2013, 12(1): 51-56.

[10]BALLAED D, BROWN C, Computer vision[M]. Upper Saddle River, New Jersey，USA: Prentice Hall, 1982：247-158.

[11]趙琳，羅漢，劉京. 加權(quán)馬氏距離判別分析方法及其權(quán)值確定[J].經(jīng)濟數(shù)學(xué)，2007，24(2): 185-188.

Segmentation and Counting of Granular Objects Based on Morphological Features

WANG Zhong-yong1, JIA Meng1, HOU Zhong-xin2, CHEN Zai-xu2

(1.School of Information Engineering，Zhengzhou University，Zhengzhou 450001，China； 2.Henan Gao-Yuan road maintenance equipment Corp., Xinxiang 453000, China)

Abstract:For the task of segmentation and counting of granular objects with low quality images, a new method which based on watershed and morphology is presented. First, an oversegmented result is gotten by conducting the watershed algorithm to the smoothed image. Then the region merge stage, which is directed by using the morphological features defined in this article, as well as weighted Mahalanobis distance and region adjacency graph (RAG), takes place. Experiments with real images demonstrate the validity of the method, that it is invariant under translation of gray scale, contrast and noise. And it performs better than any present method, in terms of precision-recall criterion.

Key words:segmentation and counting; watershed algorithm; morphological features; RAG

中圖分類號：TP391

文獻標(biāo)志碼：A

doi:10.3969/j.issn.1671-6833.2015.02.018

文章編號:1671-6833(2015)02-0080-05

作者簡介：王忠勇(1965-)，男，河南鄭州人，鄭州大學(xué)教授，研究方向為無線通信技術(shù)，數(shù)字信號處理.

基金項目：國家自然科學(xué)基金資助項目(61172086)

收稿日期:2014-11-25；

修訂日期：2015-01-03