改進的分水嶺圖像分割算法

孫惠杰，鄧廷權，李艷超

(1.哈爾濱工程大學計算機科學與技術學院，黑龍江哈爾濱150001;2.哈爾濱師范大學計算機科學與信息工程學院，黑龍江哈爾濱150025;3.哈爾濱工程大學 理學院，黑龍江 哈爾濱150001;4.哈爾濱師范大學 黑龍江省智能教育與信息工程重點實驗室，黑龍江哈爾濱150001)

圖像分割是把圖像分成各具特性的區域并從中 提取感興趣目標的技術和過程，它是圖像處理的重要環節。分割結果的好壞直接影響后續的圖像分析、理解及識別。圖像分割的方法很多，其中閾值分割、區域生長、分水嶺分割算法、聚類算法、神經網絡分割算法等都是一些經典方法。

分水嶺算法起源于Digabel等對簡單二值圖像的分析，是一種基于二值數學形態學的分割算法［1］。為了增加模型的適用范圍，Beucher等［2］建立了灰度圖像分割的分水嶺理論和算法。由于圖像中噪聲點和過多局部極小點的存在，使得分水嶺分割算法存在嚴重的過分割現象。因此降低過分割現象是當今研究熱點之一。目前，分水嶺算法的改進算法［3-4］很多，大體可以分為3類:1)對原始圖像在空間域進行預處理(區域平滑等)以消除過多的區域極值，進而減少連通分割類的個數［5］;2)對原始空間域圖像進行特定的有意義的變換，根據其性質在變換后的圖像上運用分水嶺算法進行圖像的連通分割，如在標記圖像上進行分水嶺分割［6］，將小波變換與分水嶺算法相結合對分水嶺算法進行改進［7］等;3)在原始圖像進行分水嶺變換并利用合適的區域合并方法對鄰接的同質連通區域進行融合以消除過分割現象。

本文旨在結合區域生長原理和圖像的視覺效果，基于優化方法對分水嶺算法進行改進，提出一種新型的圖像分割算法。

1 分水嶺算法與區域生長

1.1 分水嶺算法

分水嶺算法［8］是基于模擬浸水過程實現的。將圖像看成3D地形表示，即2D的地基(對應圖像空間)加上第3維的高度(對應圖像灰度值)。圖像中存在一些局部極小點。假設在每個極小點打一個孔，水將從這些孔中慢慢浸入表面，從最低的極小值點開始，水逐漸淹沒圖像的積水盆。當來自2個不同極小值點區域的水面不斷升高要匯集到一起時，在此筑起一道堤壩。在流溢過程的最后，每個極小值點被相應積水盆的堤壩所包圍，整個堤壩集合構成分水嶺。不同的積水盆代表圖像的不同分區，實現圖像分割。

實驗表明在梯度圖像上進行分水嶺分割比在原始圖像上做分割得到的結果更加準確，因此本文利用如下多尺度形態學梯度算子求取梯度圖像［9］進行分水嶺變換:

式中:GRAD為形態學梯度，Bi結構元素可以是任何形狀，只要滿足關系B0?B1?…?Bn即可。

1.2 區域生長

區域生長是一種經典的圖像分割方法，其基本思想是將具有相似性質的像素集合起來構成區域。具體來說就是先對每個需要分割的區域找一個種子像素作為生長的起始點(種子點)，將種子點鄰域中與種子像素有相同或相似性質的像素(根據某種事先確定的生長或相似準則來判定)合并到種子像素所在區域。將該區域中合并的像素當作新的種子點繼續進行上面的過程，直到沒有滿足條件的像素被包括進來為止，這樣一個區域就生長成了。這樣的過程一般用四叉樹算法可很快實現。

2 基于區域生長的分水嶺分割算法

經典分水嶺圖像分割算法會在每2個相鄰的極小點區域之間設定一個分水線，導致過分割現象的出現，影響圖像分割的質量和視覺效果。

2.1 基于區域生長的分水嶺算法

本文將區域生長算法和分水嶺算法融合實現圖像分割。借助于地形學解釋，當某2個相鄰的極小點區域的相鄰灰度值小于某一給定閾值T時，不在這2個積水盆之間筑堤壩，而把它們合并成一個區域。

設灰度圖像f的最小灰度值為hmin，最大灰度值為hmax，對于某一確定的閾值T，詳細的融合圖像分割算法可描述如下。

水流最先到達圖像的灰度最小點hmin，將這樣的點記為Mi(i=1，2，…，n)，n為像素值為hmin的點的個數。以Mi為種子點進行區域生長，用N8(Mi)表示點Mi的8鄰域，經過區域生長之后，點Mi所對應的區域為

CBT(Mi)={p|p∈N8(Mi)，|f(p)-f(Mi)|≤T}式中:CB為積水盆;將CBT(Mi)中每個點p又作為種子點進行區域生長，生長準則和上述方法一樣，直到沒有滿足生長條件的點為止。假設共對N個種子點進行了區域生長，hmin對應的積水盆為

從而得到hmin對應的積水盆。

將當前灰度值增加1，記為h，則h=hmin+1，重復以上步驟，直到h=hmax。得到所有積水盆，將積水盆的邊緣作為分水線。

在該融合算法中，區域生長閾值T的選取決定了融合算法的效果。一般情況下，結合圖像的具體特點和人的視覺特性，通過經驗和試驗獲得。這種處理方法帶有很大的主觀性和盲目性。

為了獲得區域生長算法中的最優閾值T，本文結合圖像分割常用評價指標確定一個目標函數，通過尋求該目標函數的最優值確定其閾值。

2.2 目標函數的確定

目前評價圖像分割結果好壞的準則很多［10-12］，其中區域間對比度，區域內部均勻性等是經典且常用的圖像分割質量評判指標。

設Ri、Rj表示運用閾值T進行區域生長分水嶺分割之后的2個相鄰連通部分，Ai表示區域Ri的大小(面積)，μ 表示區域Ri的平均灰度，ki表示與Ri相鄰的區域個數，定義歸一化系數:

它們分別表示分割圖像連通部分的最大方差和相鄰連通部分間的灰度平均值的最大差。記

式中:Ti表示Ri與其相鄰區域的平均對比度，而Mi表示區域Ri內部的均勻性。

為說明閾值對分割質量的影響，令區域生長閾值在［1，50］逐漸變化，對圖1(a)進行區域生長分水嶺算法分割，利用式(3)、(4)分別計算分割結果圖像的區域內部均勻性和區域間對比度，以及它們各自的平均值，所得結果如圖1(b)、(c)所示。從圖中可以看出，區域間的對比度均值隨著分割閾值的增大呈現增大的趨勢，而區域內部均勻性則隨著分割閾值的增大呈現減小的趨勢。

經過實驗和分析發現，當生長閾值增大時，分割出的連通部分面積將越來越大，其區域內部的均勻性將越來越小，而區域間的對比度越來越大。

好的圖像分割效果意味著連通區域內部的一致性較好(均勻性較大)，相鄰區域間的對比度也較大。如何限制合適的閾值進行區域增長以實現上述目的變得十分關鍵。最優閾值T的選擇正是基于一致性和對比度這對矛盾評價指標的平衡和折中。本文運用熵來刻畫這些評價指標的平衡，將分割圖像區域內灰度值的一致性和區域間灰度對比度的香農熵作為評判圖像分割質量優劣的目標函數，以此確定區域生長的最佳閾值。

熵是信息論中一個非?；静⒂兄匾獞玫母拍睿坍嬕粋€對象所蘊含的平均信息量。當分割區域內灰度值的一致性和區域間灰度對比度大時，每一個分割區域都是有意義的，蘊含的信息最多，從而可以用香農熵來量化分割結果所蘊含的信息量，進而評價圖像分割結果的好壞。對于圖像分割而言，其蘊含的信息量越多，分割效果越好。香農熵定義如下:

當pi為0或1時，熵值取最小值0;pi為1/e時熵值取得最大值，此時研究對象含有最多信息。香農熵在［1，1/e］之間為單增區間。因此為了方便且不引起歧義，將Ti和Mi規范化到得到［1，1/e］之間，得到Mi’=Mi/e，Ti’=Ti/e。設K表示分割圖像連通區域的個數，則分割圖像連通區域內部均勻性的熵和對比度的熵為

為了更準確地度量分割效果，確保分割圖像連通區域不致過大，也不致過小，更加符合人的視覺特性，將分割區域的內部均勻性和區域間對比度的極端度量(最大值的熵和最小值的熵)也作為評判標準的一部分，形成如下目標函數H(T):

使H(T)達到最大的T所對應的分割結果蘊含的信息最多，分割結果是最有意義的。T即為區域生長的最優閾值。

2.3 基于粒子群的分水嶺圖像分割算法

基于熵最大的目標函數(6)不具有解析解，通常需要優化算法而非窮舉法在有效的時間范圍內可求得最優值。本文采用粒子群優化算法確定式(6)最優解，從而得到圖像分割算法的最優閾值參數T。

粒子群優化算法［13-14］是一種基于群智能的優化算法，其把優化問題的解抽象成粒子運動的最優狀態。1999年Clerc［15］提出帶有收縮因子的粒子群優化算法，使得粒子群算法具有更好的收斂性。本文正是采用這種優化算法實現式(6)的最優解。該算法需要設定初始種群規模M，初始種群的速度:

初始位置:

式中:Vi(0)和Xi(0)分別表示第i個粒子的初始速度和初始位置。粒子的速度向量的更新公式為

結合圖像分割特點，其粒子的位置向量的每個分量都應該是整數，其更新公式采用:

式中:［*］表示最接近* 的整數，rand表示［0，1］之間的隨機值，Vi(t)為粒子i在第t次迭代的更新速度，要求在［-Vmax，Vmax］，Vmax為事先給定的最大更新速度。參數k為收縮因子，C1、C2為學習因子，C1控制粒子向自身學習的能力，C2控制粒子向群體學習的能力。依實驗結果，取k=0.729，C1=C2=2.05。粒子群的初始狀態X(0)和初始速度V(0)可隨機選取，但要求X(0)在解空間內，V(0)在［-Vmax，Vmax］內。Pi(t)表示粒子i在第t次迭代時經過的最佳位置，由

確定，G(t)表示所有粒子在尋找最優解的過程中第t次迭代時所經過的最佳位置，其中:

基于粒子群分水嶺圖像分割算法的步驟如下:

1)粒子群初始化。確定閾值T:A＜T＜B，選取初始種群規模M，設定最大迭代次數tmax。設定粒子運動的最大速度為Vmax，在解空間中選擇M個數作為粒子群的初始位置X(0)，并在［-Vmax，Vmax］內隨機選擇初始速度V(0)。

2)對于初始種群X(0)中的每個分量，依據2.1節所描述的方法求得其分割圖像，利用式(1)、(2)計算Xi(0)對應的CT和DT。取Cmax=max{CT}，Dmax=max{DT}。以C，D為歸一化系數，根據式(6)分別計算Xi(0)對應的目標函數H(Xi(0))。

3)根據式(9)、(11)，得到每個粒子的初始最優位置Pi(0)和所有粒子的最優位置G(0)。

4)令t=1。

5)依據式(7)、(8)更新速度和位置V(t)，X(t)。

6)依據2.1節所描述的方法求得Vi(t)對應的分割圖像，利用式(1)、(2)計算Vi(t)對應的CT和DT。取Cmax=max{CT，C}，Dmax=max{DT，D}。以C、D為歸一化系數，根據式(6)分別計算Xi(t)對應的目標函數H(Xi(t))與Pi(t-1)的目標函數H(Pi(t-1))。

7)根據式(10)，得到每個粒子在第t次迭代時所經過的最佳位置Pi(t)。根據式(12)，得到所有粒子在第t次迭代時所經過的最佳位置G(t)。

8)若t≤tmax，t=t+1，返回 7);若t＞tmax，則將G(t)作為區域生長的最優閾值，對應的圖像為最優分割結果，每個區域的邊界為分水線。算法結束。

經該算法后，圖像的過分割現象會得到很大程度的緩解，但是分割區域仍然可能較多，主要原因是圖像存在噪聲，算法會將噪聲點及其鄰域像素聚類形成一個面積很小的區域。因此，依據實際情況，結合人的視覺感知特性對分割結果進行后處理。

設定一個臨界值L，對于面積小于臨界值L的區域Ri，計算Ri與其相鄰區域間的灰度差異度:

式中:dif為灰度差異度，將Ri合并到差異度最小區域中，得最終分割結果。

設圖像共有N個像素點，在基于區域生長分水嶺算法中，因每個像素點都要和其周圍的8個像素點比較大小，令k=8，則這一部分的時間復雜度為O(k×N)。當粒子群的規模為M，迭代次數為tmax時，本文算法總的時間復雜度為:O(M×tmax×k×N)。

3 實驗結果與分析

為驗證本文圖像分割算法的有效性，在MATLAB 7.1平臺上實現了該算法，并對實驗結果進行分析。

3.1 參數的設置對實驗結果的影響

本文算法中涉及的參數包括:種群規模M、粒子的最大運行速度Vmax、生長參數閾值T的取值范圍和迭代次數tmax以及臨界值L。

M越大代表種群包含粒子越多，尋優能力越強，但計算量越大［16］。一般選取M為一個較小的整數值即可達到尋優目的。由于區域生長閾值參數T(灰度值，整數)的取值范圍有限，T很小時會使得區域生長的效果很不明顯，較大時有可能會使圖像產生欠分割現象，因此T不能太小，也不能太大。根據圖像分割的特點，如果圖像灰度變化比較劇烈，區域生長閾值參數T的取值范圍［A，B］應稍大一些，而對于灰度值變化比較平緩的圖像，T的取值范圍應稍小一些。Vmax控制粒子的運動速度，Vmax太大，粒子可能會飛過最優區域，太小可能導致粒子無法跳出局部最優。Vmax的值也與T的值有很大關系，由于T的值有限，而Vmax控制粒子從一個值向另一個值轉變的速度，因此Vmax也取的比較小。迭代次數tmax取的越大，計算時間越長，效果越好。初始種群V(0)是進行尋優的初始值，在解空間確定的情況下，標準粒子群算法規定可以隨機確定初始種群。如果選擇得當，會使算法具有更好的收斂性。在種群規模M和T的取值范圍［A，B］確定的情況下，為了盡量避免粒子群算法陷入局部最優，使初始粒子的值均勻分布在T的取值范圍內，依據:

確定位置向量中第j個粒子的初始值。

大量實驗表明，對于很多圖像，當T=50時，已經產生欠分割現象，如圖2(e)～(h)所示。另一方面，當T＜10時，大多數圖像存在過分割現象，如圖2(i)～(l)所示(其為圖2(a)～(d)在生長閾值為10時的分割結果)。

圖2 閾值分別為50和10時，lena、pepper、house、brain圖的分割結果Fig.2 Segmented results for images at different thresholds

因此一般可取10 ＜T＜50，在［11，49］求取最優值，當然也可以視不同的情況做一些調整。

當10 ＜T＜50 時，在式(7)中取k=0.729，C1=C2=2.05。令M=3，Vmax=4，tmax=40，利用式(13)計算初始值，基于本文算法，得到圖2中的4幅圖像的收斂情況如圖3所示。從圖中可以看到，對于house和brain圖像，利用本文算法可以很快的收斂到最優值，但是pepper圖像的收斂速度比較慢，lena圖并沒有收斂到最優值，其主要原因為house與brain圖像的目標相對單一，且圖像中不存在紋理區域，而pepper圖像目標相對比較復雜，lena圖像中存在紋理區域。

在M=8，Vmax=2，tmax=40的情況下，圖2中的4幅圖像的收斂情況如圖4所示。從圖3、4可以看出，對于比較復雜的圖像(比如說pepper圖像，lena圖像)，應該將種群規模取的大一些，以免算法陷入局部最優。當種群規模比較大時，應相應的將最大速度取的小一些，以免粒子跳過全局最優值，而對于比較簡單的圖像(如house、brain圖像)，種群規模取的較大，反而會增加運行時間，因此應將這一類圖像的種群規模取的較小一些。一般M在［3，8］取值即可。同時根據M調整Vmax的大小。迭代次數tmax也不用過大，一般迭代5～10次即可收斂到全局最優值。

基于如上所述分割算法及確定的參數值，得到圖5中圖像(a)～(d)的最優分割結果，對應的最優生長閾值分別為35、27、29、26。對于臨界值L的選擇，主要根據分割區域的分布及大小并兼顧實際需要設定L的值，將無意義的小區域合并到與之相鄰的區域中。圖5(e)～(h)為合并后的圖像分割結果，對應的臨界值分別為 200、200、50、100。

圖3 優化變量敏感性分析Fig.3 The sensitivity analyses of the optimization variables

圖 4 lena、pepper、house、brain 圖的收斂結果Fig.4 Convergence for lena，pepper，house and brain images

圖5 合并之后的最終分割結果Fig.5 Final results of segmengtation

3.2 與現有算法的比較

利用改進分水嶺算法對圖2(a)進行分割。圖6(a)為利用文獻［17］所提出的改進分水嶺算法對lena圖像進行分割得到的結果，該算法將lena圖像分成了64個區域，較傳統分水嶺和梯度分水嶺算法有了很大程度的改進。圖6(b)為利用本文算法對lena圖像進行分割的分割結果，將lena圖像分成了23個區域，其分割結果較傳統分水嶺算法、梯度分水嶺算法以及文獻［17］中提到的算法都有很大程度的改進，而且分割結果與人的直觀視覺相一致。圖6(d)為利用本文算法對6(c)的分割結果，(e)為在原圖上疊加分割線的結果，分割效果不錯，6(f)為利用文獻［18］中提到的算法進行分割的結果，可以看出，本文算法的分割結果更加符合人的直觀視覺。

圖6 不同分割方法的比較Fig.6 Comparison of results of different segmentation methods

4 結束語

本文改進了現有的圖像分割算法，提出一種區域生長與分水嶺分割算法相結合的灰度圖像自動分割算法。1)將灰度值接近的極小區域合并成為一個區域，進而改善分水嶺算法的分割效果。2)引入生長閾值，使得分割區域內部具有較好的均勻性和對比度。3)結合香農熵提出一種圖像分割結果評判準則，依此尋找最優分割結果。

實驗證明，新算法有效地解決了過分割現象，是一種有效實用的圖像分割和目標提取方法。

［1］DIGABEL H，LANTUEJOUL C.Iterative algorithms［C］//Proc 2nd European Symp Quantitative Analysis of Microstructures in Material Science，Biology and Medicine.Sturrgart，West Germany:Riederer Verlag，1978:85-99.

［2］VINCENT L，SOILLE P.Watersheds in digital spaces:an efficient algorithm based on immersion simulation ［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1991，13(6):583-598.

［3］NG H P，ONG S H，FOONG K W C，et al.Masseter segmentation using an improved watershed algorithm with unsupervised classification［J］.Computers in Biology and Medicine，2008，38:171-184.

［4］RAO A R，SRINIVAS V V.Regionalization of watersheds by fuzzy cluster analysis［J］.Journal of Hydrology，2006，318:57-79.

［5］ZHANG M，ZHANG L，CHEN H D.A neutrosophic approach to image segmentation based on watershed method［J］.Signal Processing，2012，90:1510-1517.

［6］FLORES F C，LOTUFO R A.Watershed from propagated markers:An interactive method to morphological object segmentation in image sequences［J］.Image and Vision Computing，2010，28:1491-1514.

［7］JUNG C R.Combining wavelets and watersheds for robust multiscale image segmentation［J］.Image and Vision Computing，2007，25:24-33.

［8］王小鵬.形態學圖像分析原理與應用［M］.2版.北京:清華大學出版社，2008:58.WANG Xiaopeng.Morphological image analysis principles and applications［M］.2nd ed.Beijing:Tsinghua University Press，2008:58.

［9］WANG D.A multiscale gradient algorithm for image segmentation using watershed［J］.Pattern Recognition，1997，30(12):2043-2052.

［10］章毓晉.圖像分割評價技術分類和比較［J］.中國圖像圖形學報，1996，1(2):151-158.ZHANG Yujin.A classification and comparison of evaluation techniques for image segmentation［J］.China Journal of Image and Graphics，1996，1(2):151-158.

［11］POLAK M，ZHANG H，PI M H.An evaluation metric for image segmentation of multiple objects［J］.Image and Vision Computing，2009，27:1223-1227.

［12］CORDENESA R，LUIS-GARCIA R，BACH-CUADRAB M.A multidimensional segmentation evaluation for medical image data［J］.Computer Methods and Programs in Biomedicine，2009，96(2):108-124.

［13］王科俊，郭慶昌.基于粒子群優化算法和改進的Snake模型的圖像分割算法［J］.智能系統學報，2007，2(1):53-58.WANG Kejun，GUO Qingchang.Image segmentation algorithm based on the PSO and improved Snake model［J］.CAAI Transactions on Intelligent Systems，2007，2(1):53-58.

［14］EBERHART R C，KENNEDY J.A new optimizer using particle swarm theory［J］.Sixth International Symposium on Micro Machine and Human Science，1995(3):39-43.

［15］CLERC M.The swarm and queen:towards a deterministic and adaptive particle swarm optimization［J］.IEEE International Congress on Evolutionary Computation，1999(3):1951-1957.

［16］李艷超.基于分割區域的圖像壓縮方法研究［D］.哈爾濱:哈爾濱工程大學，2011.LI Yangchao.Research on image compression based on image segmentation［D］.Harbin:Harbin Engineering University，2011.

［17］LUIS P.Fuzzy relations applied to minimize over segmentation in watershed algorithms［J］.Pattern Recognition Letters，2005，26:819-828.

［18］李蘇祺，張廣軍.基于鄰接表的分水嶺變換快速區域合并算法［J］.北京航空航天大學學報，2008，34(11):1327-1348.LI Suqi，ZHANG Guangjun.Fast region merging algorithm for watershed transform based on adjacency list［J］.Journal of Beijing University of Aeronautics and Astronautics，2008，34(11):1327-1348.