基于改進(jìn)狼群算法的三維OTSU閾值法

曹 爽，賀玉珍，安建成

(1.運(yùn)城學(xué)院 數(shù)學(xué)與信息技術(shù)學(xué)院，山西 運(yùn)城 044000；2.太原理工大學(xué) 信息與技術(shù)學(xué)院，山西 晉中 030600)

0 引 言

一維OTSU閾值法是由日本學(xué)者大津于1979年提出的[1]，其求解閾值簡(jiǎn)單、分割速度較快、提取效果良好，因此得到了廣泛應(yīng)用。但一維OTSU法只使用圖像的灰度信息，沒(méi)有涉及空間鄰域信息[2]，當(dāng)圖像有噪聲干擾時(shí)無(wú)法得到滿(mǎn)意的分割結(jié)果。為提高抗噪性和分割準(zhǔn)確率，研究者們將OTSU算法拓展到二維和三維[3]。其中三維OTSU閾值法引入了鄰域中值及均值信息，可得到理想的分割效果。

由于中值濾波易濾除椒鹽噪聲，均值濾波易濾除高斯噪聲[1]，因此，三維OTSU法既可對(duì)椒鹽噪聲或高斯噪聲有較高的抗噪性，又可對(duì)疊加混合噪聲(椒鹽噪聲加高斯噪聲)的圖像有良好的提取結(jié)果。然而OTSU法采用窮盡搜索策略，對(duì)于三維OTSU法的計(jì)算需要三重循環(huán)，導(dǎo)致計(jì)算量過(guò)多，求解復(fù)雜度呈指數(shù)型增長(zhǎng)，實(shí)時(shí)性較差。為此，景曉軍等[4]提出了三維OTSU快速遞推法，使計(jì)算效率提高，但分割準(zhǔn)確率一般，在此基礎(chǔ)上，范九倫等[5]對(duì)遞推公式進(jìn)行修正，給出新的遞推算法，使求解精度提高。葉志偉等[6]對(duì)螢火蟲(chóng)算法(firely algorithm，F(xiàn)A)進(jìn)行研究，將萊維飛行引入其中，用改進(jìn)的螢火蟲(chóng)算法來(lái)搜尋最優(yōu)分割閾值，使分割準(zhǔn)確率顯著上升。仇國(guó)慶等[7]提出了改進(jìn)的三維OTSU法，使用一維OTSU法縮小搜索空間以及迭代空間，并采用布谷鳥(niǎo)算法(cuckoo search，CS)求解閾值，既降低了計(jì)算量又具有較高的抗噪性。

新的群智能算法—狼群算法(wolf pack algorithm，WPA)是由吳虎勝等[8]于2013年提出的，其參數(shù)設(shè)置簡(jiǎn)單，已被大量應(yīng)用在多個(gè)領(lǐng)域。對(duì)狼群算法進(jìn)行研究，將花授粉算法(flower pollination algorithm，F(xiàn)PA)[9]求解的最優(yōu)閾值加入到狼群算法的初始閾值，再在圍攻過(guò)程中引入高斯變異，將改進(jìn)的狼群算法與三維OTSU算法結(jié)合，可取得較好的分割結(jié)果。

1 三維OTSU閾值法

三維OTSU算法不僅考慮了像素點(diǎn)灰度值，還考慮了像素點(diǎn)鄰域均值和中值。設(shè)計(jì)包含像素灰度及其鄰域均值和中值的三維直方圖，根據(jù)三維直方圖求取分割閾值，從而實(shí)現(xiàn)自動(dòng)閾值選取。

對(duì)于m×n的圖像，用h(x,y)表示點(diǎn)(x,y)的像素灰度值。對(duì)于某一個(gè)選定點(diǎn)k×k區(qū)域內(nèi)的平均灰度值用f(x,y)表示，計(jì)算公式如下：

(1)

該選定點(diǎn)k×k區(qū)域內(nèi)的灰度中值是g(x,y)，計(jì)算公式如下：

g(x,y)=med{h(x+m,y+n)}

(2)

其中，m=-k/2,…,k/2;n=-k/2,…,k/2。

根據(jù)上面的一系列推導(dǎo)可知，若L是圖像的灰度等級(jí)，則f(x,y)和g(x,y)的灰度等級(jí)也都為L(zhǎng)。將h(x,y)、f(x,y)以及g(x,y)組成的三元組(i,j,k)稱(chēng)作三維直方圖，如圖1(a)所示。對(duì)于直方圖上某一點(diǎn)代表的向量(i,j,k)，其出現(xiàn)的概率是pi,j,k，計(jì)算公式為：

(3)

對(duì)于一個(gè)可行的三維閾值向量(s,t,q)可將三維直方圖分成八個(gè)部分，如圖1(b)～(d)所示。其中，將區(qū)域1和0分別劃分成圖像的目標(biāo)和背景，區(qū)域2～7可看成噪聲和邊緣。由于目標(biāo)和背景區(qū)域的像素點(diǎn)多于邊緣和噪聲區(qū)域的像素點(diǎn)，從而將區(qū)域2～7的像素點(diǎn)忽略不計(jì)。

圖1 三維OTSU直方圖

令C1和C0分別代表目標(biāo)和背景這兩個(gè)區(qū)域，則目標(biāo)與背景發(fā)生的概率如下：

(4)

(5)

目標(biāo)與背景均值矢量分別為μ1和μ0：

(7)

直方圖總均值矢量μT為：

μT=(μTi,μTj,μTk)T=

(8)

由于區(qū)域2～7的像素點(diǎn)幾乎為0，則：

p0+p1=1

(9)

p0μ0+p1μ1=μT

(10)

那么目標(biāo)與背景區(qū)域之間的離散度矩SB為：

SB=p0[(μ0-μt)(μ0-μt)T]+

p1[(μ1-μt)(μ1-μt)T]

(11)

將離散度矩陣跡的值TrSB作為兩類(lèi)區(qū)域的類(lèi)間距度量函數(shù)：

TrSB=p0[(μ0i-μTi)2+(μ0j-μTj)2+(μ0k-μTk)2]+p1[(μ1i-μTi)2+(μ1j-μTj)2+(μ1k-μTk)2]

(12)

當(dāng)Tr(SB(s*,t*,q*))=max{Tr(SB(s,t,q))},0≤s,t,q≤L-1時(shí)，對(duì)應(yīng)的(s*,t*,q*)為圖像最佳分割閾值。

2 基于改進(jìn)狼群算法的三維OTSU法

2.1 狼群算法原理

狼群算法是通過(guò)個(gè)體狼與整個(gè)狼群進(jìn)行信息交互來(lái)搜索全局最優(yōu)解[10-11]。一個(gè)狼群體系包括頭狼、猛狼以及探狼這三種成員。頭狼負(fù)責(zé)指揮其他個(gè)體狼搜索獵物。探狼隨機(jī)游走去搜尋獵物，根據(jù)獵物氣味濃度判斷獵物源的遠(yuǎn)近，濃度越大(離散度矩陣跡的值越大)則狼離獵物越近。探狼發(fā)現(xiàn)獵物會(huì)馬上通知頭狼，頭狼向猛狼發(fā)起召喚，猛狼收到指令進(jìn)行圍攻追捕獵物。為保證種群的生存和發(fā)展，狼群會(huì)不斷更新群體，除去瘦弱的個(gè)體狼，保存種群實(shí)力。

2.2 改進(jìn)的狼群算法

2.2.1 花授粉算法優(yōu)化狼群算法初值

花授粉算法是由Yang Xinshe于2012年提出的。由于該算法參數(shù)少、結(jié)構(gòu)簡(jiǎn)單、全局搜索性能良好，將其引入到狼群算法中，提高狼群算法的全局尋優(yōu)能力。

在花授粉算法初始時(shí)，隨機(jī)產(chǎn)生M個(gè)種群個(gè)體p(t)={xit}，其中xit=[xti,1,xti,2,…,xti,j,…,xti,Q]，i=1,2,…,M;Q為優(yōu)化問(wèn)題的維度，j=1,2,…,Q;t是迭代次數(shù)?；ㄊ诜鬯惴ɡ棉D(zhuǎn)換概率p∈[0,1]來(lái)把控自花授粉和交叉授粉之間的轉(zhuǎn)換。

自花授粉為局部授粉，其過(guò)程為：

x(t+1)i=x(t)i+ε*(x(t)j-x(t)i)

(13)

其中，i,j為[1,M]區(qū)間內(nèi)的隨機(jī)整數(shù)，ε為變異因子，是[0,1]上服從均勻分布的隨機(jī)數(shù)。

交叉授粉為傳粉者采用萊維飛行的全局授粉，其過(guò)程為：

x(t+1)i=x(t)i+γL(x(t)i-g*)

(14)

其中，γ為步長(zhǎng)控制因子，g*為當(dāng)前群體最佳值，參數(shù)L是步長(zhǎng)，服從萊維分布。

轉(zhuǎn)換概率p的計(jì)算公式為：

p=0.8+0.2×rand

(15)

其中，rand為[-1,1]之間的隨機(jī)數(shù)。

花授粉算法主要步驟如下：

(1)初始化參數(shù)。設(shè)置群體數(shù)量M和種群個(gè)體p(t)，優(yōu)化問(wèn)題維度Q，迭代次數(shù)t，變異因子ε，步長(zhǎng)控制因子γ。

(2)計(jì)算個(gè)體適應(yīng)度值，確定全局最佳適應(yīng)度值和對(duì)應(yīng)的最佳位置g*。

(3)根據(jù)式(15)計(jì)算轉(zhuǎn)換概率p。

(4)在[0,1]區(qū)間生成一個(gè)隨機(jī)數(shù)rand1，若p≤rand1，則根據(jù)式(13)進(jìn)行局部授粉。

(5)若p> rand1，則根據(jù)式(14)進(jìn)行全局授粉。

(6)判斷算法是否滿(mǎn)足循環(huán)結(jié)束的條件，若滿(mǎn)足則輸出結(jié)果；否則，返回步驟(2)。

狼群算法初始時(shí)位置是隨機(jī)的，若位置分布得均勻，則容易尋找全局最佳值；若位置分布得不均勻，則容易陷入局部極值，降低尋優(yōu)準(zhǔn)確率。為了提高尋優(yōu)精度，加快全局搜索速度，將花授粉算法求解的較優(yōu)值作為狼群算法的初始值。

具體實(shí)現(xiàn)過(guò)程如下：令個(gè)體狼的總數(shù)為N，待求解的閾值變量數(shù)為D，那么狼群被抽象在N×D的歐氏空間中。設(shè)待尋優(yōu)的初始閾值向量為Xi=(xi1,xi2,…,xiD)，則Yi=f(Xi)(Y是目標(biāo)函數(shù)，即離散度矩陣跡的值)為每條人工狼感知的獵物氣味濃度。將花授粉算法求得的最佳解對(duì)應(yīng)的位置加入到狼群算法的初始位置。花授粉算法中，若狼群算法的群體數(shù)量N大于M，可將花授粉算法求解的最佳位置全部隨機(jī)加入狼群算法的初始位置；否則，只將花授粉算法求解的個(gè)別最佳位置隨機(jī)加入狼群算法的初始位置。

待初始位置確定后，選出頭狼(具有最優(yōu)目標(biāo)函數(shù)的人工狼)記作Ylead。頭狼不參加游走行為、召喚行為和圍攻行為，它直接進(jìn)入下一輪迭代，直到其被更兇猛的人工狼取代。

2.2.2 游走過(guò)程

先選取L_Num只人工狼作為探狼，L_Num取[n/(α+1),n/α]中的整數(shù)(α是探狼比因子)。探狼群體從各自的當(dāng)前位置隨機(jī)向m個(gè)方向游走stepa(游走步長(zhǎng))，同時(shí)記錄游走后的獵物氣味濃度(目標(biāo)函數(shù)值)并回到初始位置。向第h(h=1,2,…,m)個(gè)方向游走后，探狼i在第d維求解空間中的位置公式如下:

xhid=xid+sin(2π×h/m)*stepha

(16)

Yih(h=1,2,…,m)為探狼i感知的獵物氣味濃度。探狼個(gè)體在h個(gè)方向中選擇氣味濃度最強(qiáng)的方向游走，對(duì)探狼的位置Xi和感知的氣味濃度Yi=max{Yih}進(jìn)行更新。探狼群體重復(fù)上述游走行為，當(dāng)某只個(gè)體狼i感知的氣味濃度Yi>Ylead時(shí)，將這只人工狼作為頭狼發(fā)起召喚，否則重復(fù)游走行為直到T(游走次數(shù))達(dá)到上限。

2.2.3 召喚過(guò)程

將余下的(N-1-L_Num)只人工狼作為猛狼。猛狼聽(tīng)到頭狼的召喚后迅速追捕stepb(奔襲步長(zhǎng))逼近獵物源，那么對(duì)于第k+1次迭代第i只猛狼在第d維求解空間中的位置公式如下:

xk+1id=xkid+stepdb*(gkd-xkid)/|gkd-xkid|

(17)

其中，gkd為第k代頭狼在第d維求解空間的位置。追捕過(guò)程中，如果猛狼嗅到的氣味濃度Yi>Ylead，那么Ylead=Yi，將該猛狼作為頭狼進(jìn)行召喚，反之繼續(xù)追捕直至其與頭狼Q之間的距離diq小于dnr時(shí)進(jìn)入圍攻行為。dnr的計(jì)算公式如下:

(18)

其中，xmaxd與xmind分別是第d維求解空間的上下限，w是對(duì)距離判定的因子。

2.2.4 改進(jìn)的圍攻過(guò)程

令Gkd為第k次迭代，獵物在第d維求解空間里的位置，那么圍攻過(guò)程的公式如下:

xk+1id=xkid+λ*stepdc*|Gkd-xkid|

(19)

其中，stepc為圍攻步長(zhǎng)，λ為隨機(jī)步長(zhǎng)因子([-1,1]中均勻分布的隨機(jī)數(shù))。如果探狼和猛狼嗅到的獵物氣味濃度超過(guò)它們?cè)恢眯岬降臍馕稘舛?，則對(duì)人工狼位置進(jìn)行更新，反之位置不變。stepa、stepb以及stepc之間關(guān)系如下:

stepda=stepdb/2=2*stepdc

(20)

stepda=|xmaxd-xmind|/s

(21)

原始圍攻過(guò)程中，狼群搜索速度逐漸趨于0，使得算法出現(xiàn)“假死”狀態(tài)，降低算法搜索效率。為此，將高斯變異引入圍攻行為，使得算法后期能夠從局部最佳值附近深度搜索全局最佳值，加快尋優(yōu)速度，并減少狼群陷入局部最佳的次數(shù)。

根據(jù)高斯分布概率密度函數(shù)構(gòu)造高斯變異函數(shù)：

(22)

其中，σ為高斯分布的標(biāo)準(zhǔn)差，x*為全局最佳值。

通過(guò)GMFi(t)的擾動(dòng)，可使狼群在圍攻行為中跳出局部極值。根據(jù)高斯變異移動(dòng)的位置公式如下：

x(t)i=GMF(t)i*x(t)i

(23)

設(shè)置變異函數(shù)是為了讓所有個(gè)體狼都有機(jī)會(huì)發(fā)生變異移動(dòng)。變異函數(shù)的公式如下：

(24)

其中，N為種群數(shù)量，n為函數(shù)維度。當(dāng)0

2.2.5 群體更新

整個(gè)狼群根據(jù)“強(qiáng)多弱少”的原則分配獵物。將算法中目標(biāo)函數(shù)值最差的L匹人工狼除去，并隨機(jī)生成L匹狼，L取[n/(2×ξ),n/ξ](ζ為狼群更新因子)內(nèi)的隨機(jī)數(shù)。

2.3 改進(jìn)狼群算法優(yōu)化的三維OTSU法

2.3.1 改進(jìn)算法原理

利用改進(jìn)的狼群算法優(yōu)化三維OTSU算法的本質(zhì)即人工狼個(gè)體表示一個(gè)可行的三維閾值向量[12]。在像素點(diǎn)灰度以及該像素點(diǎn)鄰域灰度均值和中值構(gòu)成的三維向量空間中[13-14]，狼群通過(guò)游走、召喚、圍攻、種群更新這幾種行為的不斷迭代以及狼群間的信息交互來(lái)求取一組最優(yōu)解(s*,t*,q*)，使得式(12)的值最大即以目標(biāo)函數(shù)最大化為基礎(chǔ)來(lái)尋找最優(yōu)閾值。

2.3.2 改進(jìn)算法步驟

(1)利用花授粉算法求取最佳解和對(duì)應(yīng)的最佳位置。

(2)設(shè)置狼群算法參數(shù)。確定人工狼位置Xi以及人工狼的數(shù)目N，迭代次數(shù)上限kmax和游走次數(shù)上限Tmax，探狼比因子α，步長(zhǎng)因子s，距離判定因子w以及種群更新因子ζ。

(3)種群初始化。將花授粉算法得到的最佳位置加入狼群算法的初始位置。

(4)人工狼的位置為初始化的二維閾值向量，通過(guò)式(12)計(jì)算人工狼嗅到的獵物氣味濃度。根據(jù)氣味濃度大小選出頭狼、探狼以及猛狼。

(5)若某只探狼嗅到的氣味濃度Yi>Ylead，將這只探狼作為頭狼或達(dá)到游走次數(shù)上限Tmax，轉(zhuǎn)步驟(6)。

(6)猛狼聽(tīng)到頭狼召喚追捕獵物。如果其嗅到的氣味濃度Yi大于頭狼Ylead，則Ylead=Yi，該猛狼作為頭狼召喚狼群。如果Yi≤Ylead，那么猛狼繼續(xù)追捕直至dis≤dnear，轉(zhuǎn)步驟(7)。

(7)根據(jù)公式(24)計(jì)算變異函數(shù)，若0

(8)計(jì)算人工狼個(gè)體對(duì)應(yīng)的目標(biāo)函數(shù)值Yi，同時(shí)對(duì)狼群群體進(jìn)行更新。

(9)判斷算法是否達(dá)到尋優(yōu)精度或最大迭代次數(shù)，如果達(dá)到則將頭狼位置作為圖像最優(yōu)分割閾值，同時(shí)二值化分割圖像，否則返回步驟(3)。

3 實(shí)驗(yàn)仿真與分析

實(shí)驗(yàn)運(yùn)行環(huán)境為Windows 10系統(tǒng)下的MATLAB R2016a。首先，選取Benchmarks函數(shù)集中的Sphere函數(shù)和Griewank函數(shù)對(duì)改進(jìn)的狼群算法進(jìn)行尋優(yōu)性能測(cè)試。測(cè)試函數(shù)如表1所示。

表1 測(cè)試函數(shù)

由于算法的隨機(jī)性會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生較大的影響，因此，對(duì)于標(biāo)準(zhǔn)狼群算法(WPA)和改進(jìn)狼群算法分別獨(dú)立運(yùn)行20次。參數(shù)設(shè)置：狼群群體個(gè)數(shù)N為30，向量空間維數(shù)D為10；游走次數(shù)上限Tmax=10；探狼比因子α=4，更新因子ζ=6，步長(zhǎng)因子s=1 000，距離判定因子w=600。表2是兩種算法在不同測(cè)試函數(shù)下的最優(yōu)值和平均值。圖2、3分別是兩種算法在Sphere函數(shù)和Griewank函數(shù)上尋優(yōu)的收斂性能對(duì)比。

表2 函數(shù)測(cè)試結(jié)果

通過(guò)對(duì)比表2中兩個(gè)測(cè)試函數(shù)的最優(yōu)值及平均值可知，利用改進(jìn)的狼群算法求解的全局最優(yōu)值更加接近兩個(gè)測(cè)試函數(shù)的實(shí)際最小值，改進(jìn)算法的求解精度要比標(biāo)準(zhǔn)狼群算法好，尋優(yōu)能力也更強(qiáng)。

從圖2、3中可發(fā)現(xiàn)，改進(jìn)狼群算法計(jì)算全局最優(yōu)解時(shí)所需的迭代次數(shù)比標(biāo)準(zhǔn)狼群算法要小，收斂性能更好，能夠快速求解最佳值。其中適應(yīng)度值的指數(shù)為求取兩個(gè)函數(shù)實(shí)際最小值的指數(shù)。因此，可以使用改進(jìn)狼群算法對(duì)三維OTSU法進(jìn)行優(yōu)化。

圖2 Sphere函數(shù)的收斂曲線(xiàn)

圖3 Griewank函數(shù)的收斂曲線(xiàn)

為驗(yàn)證提出算法的分割效果，將文中算法與傳統(tǒng)三維OTSU算法以及文獻(xiàn)[10]提出的CWPA算法優(yōu)化三維OTSU法進(jìn)行圖像分割對(duì)比。實(shí)驗(yàn)選擇三幅圖片Bubbles、Baboon和Couple，均為512×512像素。狼群群體規(guī)模N設(shè)置為50，空間維數(shù)D為3，最大迭代次數(shù)K=100，CWPA算法的信仰空間接受更新的概率門(mén)限為0.7。分割效果如圖4所示，運(yùn)行結(jié)果如表3～表5所示。

圖4 三種算法的分割結(jié)果

表3 分割閾值對(duì)比

從表中可觀察到，對(duì)于Bubbles圖片，文中算法的分割閾值與三維OTSU算法相差不大。而對(duì)于其他兩張圖片，文中算法的分割閾值均小于另外兩種分割算法。由于算法將高斯變異引入圍攻行為,不易陷入局部最優(yōu)，能精準(zhǔn)找到最佳閾值。

表4 分割用時(shí)對(duì)比 s

從表中可看出，在時(shí)間消耗方面，三維OTSU算法用時(shí)最長(zhǎng)，文中算法和CWPA算法用時(shí)較短。這是因?yàn)槲闹兴惴ú捎没ㄊ诜鄯▋?yōu)化狼群算法后，加快了全局搜索速度。

表5 區(qū)域一致性對(duì)比

區(qū)域一致性反映了算法對(duì)圖像的分割效果。在同一類(lèi)別區(qū)域內(nèi)，像素相似度越大，表明一致性越好，分割效果越佳。區(qū)域一致性可定義為：

(25)

其中，Ri是分割完的第i個(gè)區(qū)域，Ii是區(qū)域i的面積，C是圖像的像素點(diǎn)總個(gè)數(shù)，f(x,y)是圖像(x,y)處的灰度值。區(qū)域內(nèi)部的均勻性越好則V越大(V的范圍在[0,1]內(nèi))。從表5觀察到，對(duì)于Bubbles圖片，三維OTSU算法的區(qū)域一致性較好，但分割時(shí)間過(guò)長(zhǎng)。而對(duì)于另外兩張圖片，文中算法的區(qū)域一致性均優(yōu)于CWPA算法和三維OTSU法。因此，文中算法不僅減少了搜索閾值的運(yùn)行時(shí)間，而且提高了分割圖像的精度。

4 結(jié)束語(yǔ)

三維OTSU法引入了鄰域中值及均值信息，較好地解決了三維OTSU法對(duì)噪聲敏感的問(wèn)題，但該算法計(jì)算量龐大，計(jì)算復(fù)雜度較高。為此，提出應(yīng)用改進(jìn)的狼群算法來(lái)搜索最優(yōu)三維閾值向量。實(shí)驗(yàn)結(jié)果表明，該算法使計(jì)算量大大減少，比CWPA算法的求解能力更強(qiáng)，求解精度顯著提高。