結(jié)合各向異性擴(kuò)散和分水嶺算法的圖像分割?

陳世文

(嘉應(yīng)學(xué)院計(jì)算機(jī)學(xué)院，廣東 梅州 514015)

圖像分割是根據(jù)圖像的灰度信息在空間上將具有不同類別的同質(zhì)區(qū)域劃分開來，采用一些數(shù)學(xué)方法把感興趣的區(qū)域提取出來。作為圖像處理的基礎(chǔ)，圖像分割的質(zhì)量好壞對(duì)圖像分析、圖像理解以及圖像識(shí)別等高層圖像理解產(chǎn)生直接影響，其地位十分重要。因此在圖像處理中圖像分割一直是一個(gè)難點(diǎn)和熱點(diǎn)。近幾十年來，提出了基于區(qū)域、基于邊緣和基于閾值等圖像分割算法[1－3]。

文獻(xiàn)[4]提出一種改進(jìn)FAST 特征點(diǎn)檢測(cè)算法，該算法提取邊緣特征點(diǎn)與分水嶺算法結(jié)合，提高了分割正確率，但分割仍然不夠理想。文獻(xiàn)[5]提出基于區(qū)域的分裂與合并的分割方法，能有效地獲取圖像大量區(qū)域特征信息，但怎么樣確定合并及終止合并的條件是一個(gè)難點(diǎn)。文獻(xiàn)[6]提出在分水嶺算法中結(jié)合空間變換和距離變換，分割結(jié)果較為精確，但仍然有過分割的問題。很多學(xué)者通過改進(jìn)分水嶺算法應(yīng)用在不同場(chǎng)合的分割[7－14]，取得了不同程度的效果。在分割方法中分水嶺方法具有獨(dú)特的優(yōu)點(diǎn)，其優(yōu)點(diǎn)主要是定位的精確性、輪廓線的封閉性、較快的計(jì)算速度和對(duì)圖像的變化高度敏感。因?yàn)閳D像存在局部不規(guī)則結(jié)構(gòu)以及圖像采集時(shí)受噪聲的影響，如果直接使用分水嶺分割容易引發(fā)過分割，即圖像分割得十分細(xì)小，產(chǎn)生過多的區(qū)域。區(qū)域數(shù)太多通常會(huì)給下一步區(qū)塊合并產(chǎn)生一定麻煩，導(dǎo)致區(qū)域合并錯(cuò)誤概率升高、耗時(shí)過于長，由于分割區(qū)域太細(xì)碎，降低了一些基于區(qū)域統(tǒng)計(jì)特性的準(zhǔn)確度。通常抑制過分割的方法是先對(duì)圖像進(jìn)行平滑濾波，然后再對(duì)圖像進(jìn)行分水嶺分割。傳統(tǒng)的平滑濾波器能夠較好地抑制圖像噪聲，但圖像邊緣特征同時(shí)受到嚴(yán)重?fù)p害，引發(fā)圖像邊界移位和邊界模糊。為了解決這個(gè)問題，本文提出了一種基于各向異性擴(kuò)散濾波的邊緣保持分水嶺圖像分割算法，將各向異性擴(kuò)散濾波與分水嶺算法相結(jié)合，應(yīng)用梯度模的變化設(shè)定閾值并結(jié)合圖像結(jié)構(gòu)張量構(gòu)造一個(gè)擴(kuò)散函數(shù)，這樣邊緣處沿切線擴(kuò)散保持邊緣細(xì)節(jié)，平坦處具有各向同性易于平滑噪聲，保證了良好的分水嶺結(jié)構(gòu)。然后對(duì)圖像的梯度信息進(jìn)行分層重構(gòu)，為了使梯度信息得到補(bǔ)償，采用數(shù)學(xué)形態(tài)學(xué)的開閉運(yùn)算對(duì)圖像梯度信息進(jìn)行處理。然后運(yùn)用形態(tài)學(xué)極小值標(biāo)定方法標(biāo)記處理后的圖像局部極小值，最后用分水嶺算法對(duì)圖像進(jìn)行分割。分割結(jié)果表明，該算法能夠有效地抑制噪聲和抑制過分割，同時(shí)保持物體邊緣信息，獲得準(zhǔn)確的區(qū)域輪廓信息和精確的區(qū)域分塊信息。

1 各向異性擴(kuò)散模型函數(shù)及改進(jìn)方法

各向異性擴(kuò)散濾波方法是一種基于邊緣增強(qiáng)的方法，有效去除噪聲的同時(shí)能夠保持圖像邊緣結(jié)構(gòu)信息，有利于圖像邊緣輪廓提取。Perona 和Malik于1990 年提出了一種邊緣保護(hù)的方法，即非線性各向異性擴(kuò)散方程[15](P－M 擴(kuò)散模型)，方程如式(1)所示:

式中:｜?u｜為梯度模，u0(x，y)為原始圖像；｜?u｜作用是區(qū)分該區(qū)域是邊緣區(qū)域還是平坦區(qū)域，處在平坦區(qū)時(shí)｜?u｜的值相對(duì)比較小，而處在邊緣區(qū)時(shí)｜?u｜值相對(duì)比較大；f(｜?u｜)為擴(kuò)散系數(shù)函數(shù)，其值用來表示擴(kuò)散的強(qiáng)度，其條件必須滿足f(0)＝1，＝0，且f(｜?u｜)∈(0，1)。

考慮以上各種條件，兩個(gè)擴(kuò)散系數(shù)函數(shù)被提出了，如式(2)和式(3)所示:

式中:滿足｜?u｜≥0；k為邊緣閾值，作用是用來區(qū)分邊緣區(qū)域還是平坦區(qū)域，k的取值有兩種方法，一是人工估算一個(gè)值，一是通過Canny 直方圖估計(jì)方法求得。

由式(2)和式(3)可推出:當(dāng)｜?u｜?k時(shí)，該區(qū)域?qū)儆趫D像均勻區(qū)域，擴(kuò)散會(huì)變得很慢；當(dāng)｜?u｜?k時(shí)，該區(qū)域?qū)儆趫D像邊緣處，擴(kuò)散也會(huì)變得很慢；而在邊緣閾值k位置處，該區(qū)域?qū)儆趫D像噪聲區(qū)域，因此擴(kuò)散的速度會(huì)變得比較快。由上推理可知，P－M擴(kuò)散模型的適定性問題仍有待解決。

雖消除噪聲和保留邊緣能被P－M 擴(kuò)散模型兼顧，但仍有一些問題:①噪聲點(diǎn)具有較大的梯度和較小的擴(kuò)散系數(shù)，噪聲消除不完全，容易破壞帶噪聲的邊緣；②模型在數(shù)學(xué)上屬于病態(tài)的模型，解的唯一性和存在性得不到保證。以上問題容易導(dǎo)致處理圖像后出現(xiàn)“階梯”效應(yīng)，因此提出改進(jìn)P－M 模型的擴(kuò)散系數(shù)[9]，經(jīng)過坐標(biāo)變換，可將P－M 模型變換成如式(4)所示:

式中:uηη和uζζ為u在η和ζ方向上的二階導(dǎo)數(shù)，由式(5)和式(6)求得:

式中:ux和uy分別為u在x和y方向上的偏導(dǎo)數(shù)。

為應(yīng)對(duì)P－M 模型解的不適定性，將P－M 模型轉(zhuǎn)換為式(7):

對(duì)于擴(kuò)散函數(shù)f(x)，在平坦區(qū)要求f(x)逼近0.5，梯度方向和切線方向達(dá)到同等擴(kuò)散；而在邊緣區(qū)域時(shí)f(x)逼近0，只沿著切線方向擴(kuò)散，為此本文采用式(8)所示的擴(kuò)散函數(shù)[10]:

圖像的局部信息不能很好地被圖像的梯度獲得，而更豐富的圖像局部信息能被結(jié)構(gòu)張量獲得，因此擴(kuò)散函數(shù)的自變量通過結(jié)構(gòu)張量引入[11]，結(jié)構(gòu)張量可通過海森矩陣所得到的，如式(9)所示:

這里J的兩個(gè)特征值對(duì)應(yīng)為式(10)和式(11):

因此，將獲取的結(jié)構(gòu)張量｜I1－I2｜代入到引用的擴(kuò)散函數(shù)中。

因在擴(kuò)散函數(shù)中k是屬于給定的確定閾值，在圖像中較大的噪聲幅度會(huì)引起較大的梯度值，邊緣和噪聲將難于區(qū)分，這樣會(huì)降低濾波性能，因此將閾值k設(shè)為隨著梯度｜?u｜變換而改變的函數(shù)，如式(12)所示:

式中:km邊緣閾值，經(jīng)過迭代后趨于穩(wěn)定，此處km＝3；?為常量。

基于以上幾點(diǎn)對(duì)相關(guān)參數(shù)進(jìn)行改進(jìn)，獲得最終P－M模型的改進(jìn)擴(kuò)散函數(shù)，如式(13)所示，式中e為自然常數(shù)。

2 控制標(biāo)記符的建立

一種抑制過分割的方法是對(duì)圖像中物體進(jìn)行標(biāo)記，關(guān)鍵是找尋一群與圖像局部信息相關(guān)的種子點(diǎn)。把獲取的局部區(qū)域最大值進(jìn)行標(biāo)記，把標(biāo)記結(jié)果轉(zhuǎn)化為二值圖像，即為目標(biāo)標(biāo)記圖像。對(duì)梯度圖像進(jìn)行閾值分割，經(jīng)進(jìn)一步分割可獲得背景的標(biāo)記圖像。

為實(shí)現(xiàn)對(duì)目標(biāo)的細(xì)化，通過歐氏距離變換對(duì)背景標(biāo)記圖像進(jìn)行操作，可以用一個(gè)二維數(shù)組表示背景標(biāo)記圖像，其大小與原圖像同樣，記為Am×n＝[aij]，其中aij＝1 的像素表示為目標(biāo)像素，aij＝0 的像素表示為背景像素。設(shè)B＝{(x，y)｜axy＝0} 為背景點(diǎn)集合，對(duì)A中所有像素點(diǎn)(i，j)求距離變換，如式(14)所示:

式中:D[(i，j)，(x，y)]＝獲取歐氏距離變換圖像DM×N。

為了求極值，引入數(shù)學(xué)形態(tài)學(xué)運(yùn)算。式(15)是腐蝕運(yùn)算，表示圖像f被結(jié)構(gòu)元素b腐蝕；式(16)是膨脹運(yùn)算，表示圖像f被結(jié)構(gòu)元素b膨脹。接著對(duì)腐蝕、膨脹后的圖像分別進(jìn)行開重構(gòu)運(yùn)算和閉重構(gòu)運(yùn)算，即對(duì)圖像f的n個(gè)灰度級(jí)逐點(diǎn)比較最大值，記為fc，即為局部最大值，進(jìn)行標(biāo)記，然后把標(biāo)記圖像轉(zhuǎn)換成二值圖像，記為fgm。

采用數(shù)學(xué)形態(tài)學(xué)中求極小值的標(biāo)定方法，融合距離變換生成的圖像和目標(biāo)提取的標(biāo)記圖像整合生成梯度圖像的局部極小值，舍棄原來圖像中的全部局部極小值。重新生成的梯度圖像中，強(qiáng)制對(duì)不為零的僅且對(duì)應(yīng)有二值標(biāo)記的圖像信息有局部極小值，重新生成的梯度圖像如式(17)所示:

式中:IMMIN(?)是形態(tài)學(xué)極小值標(biāo)定操作。

3 分水嶺算法

Vincent 和Soille 將數(shù)學(xué)形態(tài)學(xué)中的分水嶺變換用于圖像處理中，其基本思想是把一幅圖像模擬成拓?fù)涞匦螆D，圖像中各個(gè)像素的灰度值對(duì)應(yīng)為該點(diǎn)的海拔高度，高灰度值看作為山峰，低灰度值看作為山谷，水總是從高處往低處流動(dòng)，分水嶺對(duì)應(yīng)著圖像中的脊或每一局部極小值。假設(shè)水從每一局部極小值涌出并且水位逐漸上升，水慢慢浸入整個(gè)集水盆地。當(dāng)處在不同的集水盆地的水位上升到相應(yīng)山脊時(shí)就將匯合，為了隔斷這些水匯合，就要建立堤壩。堤壩會(huì)因水位的升高而增加，直到整個(gè)圖像沉入水中，堤壩對(duì)應(yīng)分水線連通形成的路徑就成為分水嶺。

設(shè)給定待分割的一幅圖像為f(x，y)，其圖像梯度為g(x，y)，分水嶺是通過梯度空間進(jìn)行計(jì)算的。設(shè)M1，M2，…，MR表示梯度g(x，y)中各局部極小值的像素位置，與Mi對(duì)應(yīng)的區(qū)域中的像素坐標(biāo)集合記為C(Mi)。梯度閾值用n表示，(u，v)的像素集合記為T[n]，g(u，v)

與Mi相應(yīng)的區(qū)域，符合條件的坐標(biāo)集合Cn(Mi)定為一幅二值圖像，如式(19)所示:

圖像中全部符合梯度值小于閾值n的像素集合記為C[n]，則有式(20):

那么全部區(qū)域的并集是C[max＋1]，max 為灰度最大值。

分水嶺初始化算法如式(21)所示，接著算法進(jìn)行逐步迭代。

假設(shè)步驟n已創(chuàng)建了C[n－1]，令S代表T[n]中的連通組元集合，用3×3 全是1 的結(jié)構(gòu)模塊對(duì)S∩C[n－1]進(jìn)行膨脹，則可獲取像素為1 的邊界。由于噪聲影響，理想的分水嶺變換應(yīng)該為梯度圖像，因?yàn)樘荻葓D像能夠增強(qiáng)圖像中的噪聲，同時(shí)減弱了灰度平坦的區(qū)塊，這樣噪聲更容易被濾波器濾除掉。

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

本文圖像分割算法的實(shí)驗(yàn)步驟如下:①將原始圖像轉(zhuǎn)換為灰度圖像，如果是灰度圖像則保持；②利用基于梯度模和結(jié)構(gòu)張量融合的各向異性擴(kuò)散對(duì)灰度圖像進(jìn)行處理；③將經(jīng)各向異性擴(kuò)散處理后的圖像生成梯度圖像，再對(duì)梯度圖像首次進(jìn)行分水嶺算法；④對(duì)上一步生成的圖像進(jìn)行距離變換和數(shù)學(xué)形態(tài)學(xué)處理，得到局部極小值，標(biāo)記控制；⑤得到重構(gòu)后的梯度圖像，再次對(duì)其進(jìn)行分水嶺分割。

算法的實(shí)驗(yàn)采用仿真軟件為MATLAB 2016b。為了測(cè)試文中算法抑制噪聲的性能，引入峰值信噪比(PSNR)，如式(22)、式(23)所示:

式中:MN 是圖像像素個(gè)數(shù)，xi是無噪聲圖像像素值，ui是測(cè)試圖像像素值。加入噪聲，測(cè)試如表1 所示，一定程度上PSNR 反映了圖像數(shù)據(jù)被處理前后變化的統(tǒng)計(jì)平均，PSNR 越大，則數(shù)據(jù)失真越少。

表1 試驗(yàn)數(shù)據(jù)

為了測(cè)試邊緣信息保持性能，引入邊緣保持指數(shù)(FOM)作為評(píng)價(jià)標(biāo)準(zhǔn)，其定義如式(24)所示:

式中:NA表示含噪圖像邊緣信息中的非零像素?cái)?shù)，Nf表示去噪后邊緣信息中的非零像素?cái)?shù)，di是歐氏距離，表示含噪邊緣圖像第i個(gè)像素與濾波后的最近邊緣像素之間的幾何距離，α是懲罰因子，擬取α＝1/9。對(duì)原始米粒圖像添加高斯噪聲，變化方差進(jìn)行測(cè)試，測(cè)試數(shù)值如表2 所示，由數(shù)據(jù)可知本文算法能夠較好地保持邊緣信息。

表2 邊緣測(cè)試

經(jīng)去噪和邊緣性能檢測(cè)后，對(duì)米粒圖像進(jìn)行分割實(shí)驗(yàn)，分別使用傳統(tǒng)分水嶺算法、文獻(xiàn)[14]提出的分割算法和本文改進(jìn)的各向異性擴(kuò)散算法對(duì)圖像進(jìn)行分割，分割效果如圖1 所示?？梢园l(fā)現(xiàn)傳統(tǒng)算法對(duì)米粒的分割能力較弱，文獻(xiàn)[14]算法缺乏全局信息，邊緣信息提取較差。本文算法通過各向異性擴(kuò)散明顯保持了邊緣信息，去除了虛假極小值點(diǎn)，降低了過分割現(xiàn)象，同時(shí)亦實(shí)現(xiàn)了對(duì)目標(biāo)米粒的較完整分割。

圖1 米粒圖像分割實(shí)驗(yàn)

對(duì)原始硬幣圖像進(jìn)行分割實(shí)驗(yàn)，先給硬幣加入泊松噪聲，分別使用傳統(tǒng)分水嶺算法、文獻(xiàn)[14]提出的分割算法和本文改進(jìn)的各向異性擴(kuò)散算法對(duì)圖像進(jìn)行分割，分割效果如圖2 所示。傳統(tǒng)分水嶺算法分割的區(qū)域塊數(shù)較多，分割效果不理想，主要是因?yàn)樵肼曈绊?、邊緣過分割和虛假極小值點(diǎn)等因素影響。文獻(xiàn)[14]算法的分割結(jié)果得到了明顯改進(jìn)，效果較傳統(tǒng)方法好，但仍然存在虛假邊緣，仍有過分割等問題。本文所提算法，能夠抑制噪聲，同時(shí)使弱邊界保持較好，極小值點(diǎn)的標(biāo)記相對(duì)合適，因此，分割效果明顯優(yōu)于前兩種方法，具有較好的抗噪性能。

圖2 加入泊松噪聲硬幣圖像分割實(shí)驗(yàn)

5 結(jié)語

本文提出了一種基于各向異性擴(kuò)散的分水嶺分割算法。應(yīng)用梯度模的變化設(shè)定閾值取代P－M 算子人為設(shè)定邊界閾值的方法，并引入圖像結(jié)構(gòu)張量到擴(kuò)散函數(shù)里面，達(dá)到平滑噪聲的同時(shí)保持邊緣細(xì)節(jié)的目的，并使其具有良好的分水嶺結(jié)構(gòu)。然后，借助數(shù)學(xué)形態(tài)學(xué)的開閉運(yùn)算補(bǔ)償圖像梯度信息，實(shí)現(xiàn)梯度信息分層重構(gòu)，接著基于形態(tài)學(xué)最小值標(biāo)記方法對(duì)圖像局部極小值進(jìn)行標(biāo)記，最后用分水嶺算法對(duì)處理后的圖像進(jìn)行分割。實(shí)驗(yàn)對(duì)無噪聲圖像和加噪聲圖像進(jìn)行分割，結(jié)果表明該方法具有良好的分割效果，尤其對(duì)噪聲圖像有較好的魯棒性。