細(xì)胞膜機(jī)制螢火蟲算法優(yōu)化多閾值Otsu圖像分割

劉鑫晶，劉彥隆，徐鑫鑫

(太原理工大學(xué) 信息與計算機(jī)學(xué)院，山西 晉中 030600)

1 引 言

數(shù)字圖像處理是一種將圖像轉(zhuǎn)換為數(shù)字形式并對其執(zhí)行某些操作的方法.其中，圖像分割這一步驟是非常重要的，其精確性對整個圖像處理得到的結(jié)果有著至關(guān)重要的影響.基于閾值的分割算法具有分割穩(wěn)定性強(qiáng)，模型簡單及魯棒性強(qiáng)的優(yōu)點而得到了研究者們十分廣泛地研究和應(yīng)用.最大類間方差法(Otsu)是由日本學(xué)者Nobuyuki Otsu在1979年提出的具有代表性的閾值確定方法[1]，Otsu算法的基本原理是最小二乘法.為了使該算法能夠適應(yīng)具有多個目標(biāo)的復(fù)雜圖像，有研究者將Otsu算法從單閾值擴(kuò)展到多閾值[2].多閾值Otsu算法解決了復(fù)雜圖像的分割問題，但隨閾值個數(shù)逐漸增多，算法運算時間會指數(shù)增長，這也是所謂的NP(nonlinear problems)難題[3].

近幾年，為了減少多閾值Otsu算法的運算時間，如布谷鳥搜索(CS)，人工蜂群(ABC)，螢火蟲算法(GSO)，粒子群優(yōu)化(PSO)，差分進(jìn)化(DE)等啟發(fā)式算法用于圖像分割[4-8].GSO算法不但可以找到局部最優(yōu)解，還能夠得到全局最優(yōu)解，則在多閾值Otsu圖像分割中有普遍應(yīng)用.Assistant等人[9]提出了基于GSO算法優(yōu)化的最大期望算法應(yīng)用于多閾值分割，Simone A.Ludwig等人[3]提出了一種適用于多閾值圖像分割任務(wù)的改進(jìn)的GSO算法(EGSO).國內(nèi)，于超杰等人[10]提出了基于反向?qū)W習(xí)策略的改進(jìn)型螢火蟲算法(OBLFA)應(yīng)用于Otsu圖像分割，申玄京等人[11]提出了基于GSO算法的遞歸多閾值Otsu算法，陳凱等人[12]提出了基于反向GSO算法針對缺陷圖像的多閾值分割.

標(biāo)準(zhǔn)GSO算法尋求最優(yōu)解的過程中其精度低，速度慢和容易陷入局部最優(yōu)值，即“早熟”現(xiàn)象[13].雖然國內(nèi)外針對標(biāo)準(zhǔn)GSO算法存在的這些缺點提出了許多改進(jìn)型算法，并研究應(yīng)用于多閾值圖像分割中，但是結(jié)果都不盡理想.

為了更好地解決標(biāo)準(zhǔn)GSO算法存在的收斂速度慢和易出現(xiàn)“早熟”現(xiàn)象等問題以及多閾值圖像分割的時間復(fù)雜度高的NP難題，本文提出了一種基于細(xì)胞膜和GSO混合優(yōu)化算法的多閾值Otsu圖像分割.細(xì)胞膜優(yōu)化算法采用了保留最優(yōu)解策略，將其應(yīng)用于GSO算法中，收斂速度提高并能夠有效避免GSO算法陷入局部的最優(yōu)值.采用混合優(yōu)化算法將對閾值的搜索問題轉(zhuǎn)換成對目標(biāo)函數(shù)求它的最優(yōu)解問題，很好地降低了其時間復(fù)雜度.

2 螢火蟲混合優(yōu)化算法的產(chǎn)生

2.1 標(biāo)準(zhǔn)的螢火蟲算法

標(biāo)準(zhǔn)的GSO算法[14]主要包括5個步驟：初始化GSO中參數(shù)，更新熒光素值，選擇移動的方向，移動到更新的位置，更新決策域的半徑大小.

1)初始化GSO中參數(shù)

將M只螢火蟲隨機(jī)地分在一個D維的目標(biāo)空間內(nèi)，每只的初始熒光素大小l0和感知半徑r0一樣.并初始化GSO算法的最大迭代次數(shù)Tmax.

2)更新熒光素值

所有螢火蟲在下一次迭代前，都依據(jù)式(1)得到新的熒光素濃度.

li(t+1)=(1-ρ)li(t)+γJ[xi(t+1)]

(1)

3)選擇移動的方向

在每次迭代時，每只螢火蟲i都會在決策域半徑內(nèi)尋找滿足以下兩個前提的同伴j的集合Ni(t)：一是螢火蟲j在i的決策域范圍中；二是螢火蟲j的熒光素濃度值必須大于i.然后，螢火蟲i在集合Ni(t)中使用輪盤賭法選擇一只螢火蟲.Ni(t)內(nèi)每只螢火蟲被選的概率如式(2)所示.

(2)

4)移動到更新的位置

螢火蟲i向第3)步選好的螢火蟲j進(jìn)行第t+1次位置更新，如式(3)所示.

(3)

5)更新決策域半徑

每次迭代完后，決策域半徑進(jìn)行更新，如式(4)所示.

(4)

2.2 自適應(yīng)步長的螢火蟲算法

2018年李順等[15]提出一種自適應(yīng)步長GSO算法(adaptive-step GSO，A-GSO).A-GSO中引入熒光因子Hi，使得在搜索時動態(tài)改變螢火蟲的移動步長，使算法在一定程度上改善了“早熟”現(xiàn)象，提高求解速度.

熒光因子Hi的大小由式(5)計算.

(5)

式中，Xi,Xm分別表示螢火蟲i和當(dāng)代熒光素濃度最大的螢火蟲所在的位置；dmax表示當(dāng)代最優(yōu)的螢火蟲到其余螢火蟲距離的最大值.

移動步長s按照式(6)進(jìn)行調(diào)整.

si=smin+(smax-smin)×Hi

(6)

2.3 基于細(xì)胞膜和螢火蟲混合優(yōu)化算法

雖然A-GSO在一定程度上改善了標(biāo)準(zhǔn)GSO的“早熟”問題，但是該算法在運行時很容易出現(xiàn)一些螢火蟲的同伴集合Ni(t)為空集，這些螢火蟲在迭代過程中將不會移動到更新的位置，導(dǎo)致算法的收斂速率降低并且容易陷入局部最優(yōu).2011年林土勝、譚世恒等人[16]提出了一種新型的全局優(yōu)化算法，即細(xì)胞膜算法(CM)，并且通過數(shù)值實驗，已經(jīng)驗證該算法有很好的全局尋優(yōu)能力和快速收斂能力.但是該算法目前沒有得到廣泛的應(yīng)用.本文借鑒該算法中在迭代前將物質(zhì)分類和保留最優(yōu)解的思想將CM算法運用在A-GSO算法上，將兩個算法進(jìn)行結(jié)合，得到細(xì)胞膜和螢火蟲混合優(yōu)化算法(CMA-GSO).結(jié)合了細(xì)胞膜算法思想和自適應(yīng)步長思想的螢火蟲算法很好地改善了標(biāo)準(zhǔn)GSO算法的“早熟”問題，收斂速度慢以及求解精度不高的問題.CMA-GSO算法流程圖如圖1所示.

圖1 CMA-GSO算法流程圖Fig.1 Flow chart of CMA-GSO algorithm

CMA-GSO算法是對標(biāo)準(zhǔn)的GSO算法的“使用輪盤賭注法選擇移動方向”這一步驟使用了細(xì)胞膜算法中對物質(zhì)類型進(jìn)行劃分以及每種物質(zhì)有各自不一樣的運動方式來代替了.CMA-GSO算法的詳細(xì)步驟如下所示：

1)螢火蟲類型劃分

首先把所有的螢火蟲的熒光素濃度從小到大排序，排在前P比例的螢火蟲劃分為當(dāng)代優(yōu)質(zhì)螢火蟲，排在后面的螢火蟲劃分為當(dāng)代劣質(zhì)螢火蟲.然后對劣質(zhì)螢火蟲再劃分成兩類：高濃度和低濃度.對于螢火蟲i，它所處的濃度定義為其決策域內(nèi)同伴螢火蟲數(shù)量占總螢火蟲數(shù)的百分比，如式(7).

(7)

式中，n表示螢火蟲i決策域內(nèi)同伴螢火蟲數(shù)量；m表示目標(biāo)空間總螢火蟲數(shù)量.所有濃度的平均值為MeanCon，如式(8).

(8)

若某劣質(zhì)螢火蟲所處的濃度大于MeanCon，那么該螢火蟲劃分為高濃度劣質(zhì)螢火蟲，否則劃分為低濃度劣質(zhì)螢火蟲.

2)當(dāng)代優(yōu)質(zhì)螢火蟲自由擴(kuò)散

優(yōu)質(zhì)螢火蟲即代表細(xì)胞膜算法中的脂溶性物質(zhì).自由擴(kuò)散是指脂溶性物質(zhì)由膜的高濃度側(cè)向低濃度側(cè)的擴(kuò)散過程，整個過程不需要載體和能量.對于每一個高濃度優(yōu)質(zhì)螢火蟲sXi將向低濃度螢火蟲群LX方向擴(kuò)散locn次，選取最好的擴(kuò)散位置處的螢火蟲作為最終的結(jié)果.其詳細(xì)計算方法如式(9)所示.

(9)

經(jīng)過此過程后，由原來的優(yōu)質(zhì)螢火蟲群sX轉(zhuǎn)變?yōu)樾碌膬?yōu)質(zhì)螢火蟲群newsX.

3)當(dāng)代高濃度劣質(zhì)螢火蟲協(xié)助擴(kuò)散

高濃度劣質(zhì)螢火蟲即代表細(xì)胞膜算法中的高濃度非脂溶性物質(zhì)，在細(xì)胞膜中，高濃度非脂溶物質(zhì)做協(xié)助擴(kuò)散由細(xì)胞膜的高濃度側(cè)轉(zhuǎn)移到低濃度側(cè).此過程只需要載體，不需要能量，因此引入載體因子cf.高濃度劣質(zhì)螢火蟲hiXi的載體因子cfi如式(10)所示計算.

(10)

高濃度劣質(zhì)螢火蟲存在兩種運動方式.當(dāng)存在載體時，螢火蟲進(jìn)行協(xié)助擴(kuò)散，向低濃度方向運動，同上面自由擴(kuò)散的計算方法類似.當(dāng)不存在載體時，其向當(dāng)前最優(yōu)螢火蟲Xbest方向擴(kuò)散locn次，并選取最好的擴(kuò)散位置處的螢火蟲作為最終的結(jié)果.其擴(kuò)散方法如式(11)所示.

(11)

經(jīng)過此過程后，由原來的高濃度劣質(zhì)螢火蟲群hiX轉(zhuǎn)變?yōu)樾碌母邼舛攘淤|(zhì)螢火蟲群newhiX.

4)當(dāng)代低濃度劣質(zhì)螢火蟲主動運輸

低濃度劣質(zhì)螢火蟲即代表細(xì)胞膜算法中的低濃度非脂溶性物質(zhì).在細(xì)胞膜中，低濃度物質(zhì)主動運輸，從細(xì)胞膜的低濃度側(cè)轉(zhuǎn)移到高濃度側(cè)，須要載體和能量.需要引入載體因子cf和能量因子ef.載體因子cfi的計算方法與上一步類似.能量因子efi的計算方法為：按照低濃度劣質(zhì)螢火蟲的熒光素值從小到大進(jìn)行排序，熒光素值最大的螢火蟲其能量因子為efmin，熒光素值最小的螢火蟲其能量因子為efmax，其它螢火蟲的能量因子ef介于efmin和efmax之間，并按照排序的順序線性計算.其中，efmin和efmax為[0，1]內(nèi)的常數(shù)，在這里取efmin和efmax分別為0和1.

對于某低濃度劣質(zhì)螢火蟲liXi存在三種可能的運動方式，首先判斷是否滿足能量條件，其次判斷是否滿足載體條件.因此低濃度劣質(zhì)螢火蟲liXi的三種可能的運動方式分別為：一是能量條件和載體條件同時滿足時，螢火蟲liXi進(jìn)行主動運輸，向高濃度螢火蟲群HX方向運動locn次，詳細(xì)運動公式與優(yōu)質(zhì)螢火蟲自由擴(kuò)散的方法類似，這里不再贅述；二是滿足能量條件而不滿足載體條件時，螢火蟲liXi向當(dāng)前最優(yōu)螢火蟲Xbest方向擴(kuò)散locn次；三是不滿足能量條件時，螢火蟲liXi在其決策域半徑內(nèi)隨機(jī)運動locn次，運動方式如式(12)所示.最后需要選取最好的擴(kuò)散位置的螢火蟲作為最終的結(jié)果.

newliXi=ld+rand()×(ud-ld)

(12)

經(jīng)過此過程后，由原來的低濃度劣質(zhì)螢火蟲群liX轉(zhuǎn)變?yōu)樾碌牡蜐舛攘淤|(zhì)螢火蟲群newliX.

5)當(dāng)前最優(yōu)螢火蟲循環(huán)更新移動位置

在優(yōu)質(zhì)螢火蟲附近往往會存在更優(yōu)螢火蟲，充分發(fā)掘當(dāng)前最優(yōu)螢火蟲決策域空間，有利于提高尋優(yōu)精度.以當(dāng)前最優(yōu)螢火蟲Xbest為中心進(jìn)行搜索，向更亮的螢火蟲運動，不斷地縮小決策域半徑，進(jìn)行反復(fù)運動.最后得到新的最優(yōu)螢火蟲newXbest即為當(dāng)前全局最優(yōu)螢火蟲.

3 混合優(yōu)化算法應(yīng)用于多閾值Otsu圖像分割

3.1 多閾值Otsu圖像分割

Otsu指出最大類間原則和最小類內(nèi)原則基本等價，本文采用最大類間原則.

設(shè)圖像中設(shè)有n個待區(qū)分的類，有n-1個閾值T1,T2,…,Tn-1將圖像分為n個類.這些類分別表示為:

C0={0,1,…,T1}，…，Cm={Tm+1,…,Tm+1}，…，Cn-1

={Tn-1+1,…,L-1}

(13)

(14)

(15)

式中，k=0,1,…,n-1，T0=0，Tn=L.

整幅類間方差為:

(16)

3.2 混合優(yōu)化的多閾值Otsu圖像分割算法

多閾值Otsu計算過程的偽代碼如下所示：

Fori0 toL

Forji+1 toL

Forkj+1 toL

?

begin

maxVal=max(maxVal,cal(4))

t=canSwap(i,j,k,…)

end

可見，如果將一幅圖像使用n-1個閾值分割，則這個算法的解空間是n-1維的.在解空間中任意一個n-1維向量都需要計算式(16)和式(14)，算法的時間復(fù)雜度為O(Ln).

因此，本文考慮在對多閾值Otsu圖像分割的閾值搜索過程中使用CMA-GSO優(yōu)化算法將對最佳閾值的尋找過程轉(zhuǎn)變?yōu)閷δ繕?biāo)函數(shù)的尋優(yōu)問題，將多閾值Otsu圖像分割的最大類間方差準(zhǔn)則下的函數(shù)式(16)作為目標(biāo)函數(shù)帶入到CMA-GSO優(yōu)化算法中，得到最優(yōu)閾值.從而降低算法的時間復(fù)雜度.

4 實驗結(jié)果與分析

為了驗證本文提出的基于CMA-GSO優(yōu)化的多閾值Otsu圖像分割算法在復(fù)雜多目標(biāo)圖像分割上的分割效果和在算法運行速率上的優(yōu)越性，本文選取Boats圖、Camera圖、MRI圖和Saturn圖四幅圖像作為測試對象分別進(jìn)行雙閾值、三閾值、四閾值、五閾值圖像分割.并將本文提出的方法與經(jīng)典多閾值Otsu法、GSO優(yōu)化多閾值Otsu法進(jìn)行對比.MRI圖是腦部核磁共振圖[17]，對MRI圖進(jìn)行多閾值分割有助于對腦部病情的診斷，具有現(xiàn)實意義.實驗在2.80GHz CPU和8.00GB內(nèi)存的PC機(jī)、window10.0操作系統(tǒng)和Matlab R2015b環(huán)境中進(jìn)行的.

實驗中對CMA-GSO優(yōu)化算法的參數(shù)設(shè)置如下：最大迭代次數(shù)Tmax=200；螢火蟲個數(shù)M=100；初始熒光素l0=5；初始感知半徑r0=10；劃分比例P=50%；空間維數(shù)D=10；揮發(fā)因子ρ=0.4；介質(zhì)吸收系數(shù)γ=0.6；擴(kuò)散次數(shù)locn=10.圖2給出了4幅圖像的原始灰度圖及采用CMA-GSO優(yōu)化多閾值Otsu算法得到的二、三、四和五閾值分割結(jié)果.可見，分割閾值越大，分割得到的圖像的灰度信息更加豐富，更能充分體現(xiàn)出原圖中的目標(biāo)對象.

表1體現(xiàn)了對于復(fù)雜多目標(biāo)圖像的分割，在分割閾值較小時，本文方法與其余兩種方法的分割精度較為一致，而隨閾值的增多，本文方法分割更加精確.

圖2 從左向右依次是原始灰度圖，二閾值、三閾值、四閾值、五閾值分割圖；從上向下依次是Boats圖、Camera圖、MRI圖和Saturn圖Fig.2 From left to right arethe picture of original grayscale map,two thresholds,three thresholds,four thresholds,and five thresholds segmentation;from top to bottom are the picture of Boats,Camera,MRI and Saturn

為了體現(xiàn)本文算法在計算時間和圖像分割效果上的優(yōu)越性，采用圖像的峰值信噪比PSNR、計算時間來進(jìn)行對比，結(jié)果如表2所示.可見本文方法在二、三、四、五閾值分割的速率分別比經(jīng)典多閾值Otsu算法快5.52倍，237.01倍，1184.24倍和9887.46倍.雖然與GSO優(yōu)化多閾值Otsu算法對比，本文算法的計算速率基本一樣，但是本文算法的圖像峰值信噪比(PSNR)值要優(yōu)于GSO優(yōu)化多閾值Otsu算法，說明本文算法在連續(xù)運行時分割質(zhì)量較高.

表1 3種算法得出的圖像分割閾值與Otsu函數(shù)值
Table 1 Segmentation thresholds and Otsu function from three methods

圖像閾值個數(shù)閾 值經(jīng)典多閾值Otsu法GSO優(yōu)化多閾值Otsu法CMA-GSO優(yōu)化多閾值Otsu法Otsu函數(shù)值/103經(jīng)典GSO優(yōu)化CMA-GSO優(yōu)化Boats260,11460,11260,1122.31682.31682.3168354,120,17554,118,17654,118,1742.44352.40532.4053463,119,148,17063,120,151,17060,120,150,1702.51032.54842.5488561,85,115,158,19761,85,115,160,20261,85,115,160,2002.75732.84292.8745Camera2109,147109,147109,1473.64563.64563.6456351,152,20051,156,20251,150,2023.74843.74883.7489454,91,144,20054,94,147,20054,97,147,2003.84693.85903.8706537,51,100,154,19837,51,100,154,20031,51,100,154,1973.96873.97073.9726MRI2159,182159,184159,1841.85241.85241.85243132,159,181136,163,183136,161,1831.92631.93551.93554142,162,183,205142,160,190,208142,166,190,2152.15322.18462.1846566,99,146,161,20066,99,145,161,20266,99,148,161,2002.25412.29402.2948Saturn251,11751,11351,1132.87562.87562.8756339,81,13839,86,13439,86,1382.96782.96822.9682422,54,95,14210,54,92,1428,54,97,1423.15933.16123.1612537,54,82,113,15337,54,86,113,14939,54,86,113,1543.25693.26963.2752

表2 3種算法的PSNR與計算時間
Table 2 Comparison of PSNRs and computing time from three methods

圖像閾值個數(shù)經(jīng)典多閾值Otsu法PSNR/dB計算時間/sGSO優(yōu)化多閾值Otsu法PSNR/dB計算時間/sCMA-GSO優(yōu)化多閾值Otsu法PSNR/dB計算時間/sBoats213.521410.3213.85633.3313.85263.14315.2563912.4815.34874.2015.34873.85417.21686584.3517.56825.5217.56825.56518.357659423.6418.59846.3818.60356.01Camera214.315612.1414.33252.0214.33252.08317.2698845.6817.45313.5717.45313.45418.34924862.4918.58423.9918.54263.74519.154836514.8519.35495.3119.35495.38MRI215.604210.6915.25432.3115.41362.33318.3619863.2418.65453.9518.65243.64419.10344563.4719.59424.4819.59424.12520.630458642.1520.47235.0120.41295.12Saturn218.365211.9418.56142.1118.56142.16320.2540764.2920.65193.2120.65183.02421.23574982.3621.53164.3421.53163.99522.023363245.5922.23144.6322.23164.38

將使用本文所提出的算法和文獻(xiàn)[3]，文獻(xiàn)[11]所提出的算法分別對國際公共圖像測試庫中的Camera圖進(jìn)行迭代次數(shù)從0到1000的二閾值圖像分割實驗，并計算二閾值分割時間和與窮舉法得到的閾值相比閾值變化率，實驗結(jié)果如圖3所示.與文獻(xiàn)[3]和文獻(xiàn)[11]相比，本文所提出的算法隨迭代次數(shù)的增長分割時間增長最慢，即本文所提出的算法降低了時間復(fù)雜度.并且本文所提出的算法隨迭代次數(shù)的增長閾值變化率降低最快，即本文所提出的算法得到的閾值精確度高.

綜上所述，本文提出的基于CMA-GSO優(yōu)化算法的多閾值Otsu圖像分割方法，能夠有效快速并且精確地分割復(fù)雜多目標(biāo)圖像，并且在醫(yī)學(xué)，科學(xué)等領(lǐng)域具有現(xiàn)實意義.

圖3 對Camera圖進(jìn)行二閾值分割時，左圖表示計算時間隨迭代次數(shù)增長的變化曲線，右圖表示閾值變化率隨迭代次數(shù)增長的變化曲線Fig.3 For the two-threshold segmentation of the picture of Camera,the left graph shows the change curve of the computation time with the increase of iterations,and the right graph shows the change curve of the threshold change rate with the increase of iterations

5 結(jié) 論

本文主要針對復(fù)雜多目標(biāo)灰度圖像分割耗時多的缺點，提出了基于CMA-GSO優(yōu)化算法的多閾值Otsu圖像分割方法.首先，針對GSO算法易陷入局部最優(yōu)解的“早熟”問題，聯(lián)想到細(xì)胞膜算法的保留最優(yōu)解策略，提出了CMA-GSO優(yōu)化算法.其次，針對多閾值Otsu圖像分割算法在時間復(fù)雜度上的NP難題，提出了將混合優(yōu)化算法應(yīng)用于多閾值Otsu圖像分割，在最優(yōu)閾值的搜索過程中采用自然啟發(fā)式搜索而不是窮舉法搜索，降低閾值搜索的時間復(fù)雜度.最后，對多目標(biāo)灰度圖像進(jìn)行分割，并與經(jīng)典多閾值Otsu算法和GSO優(yōu)化多閾值Otsu算法進(jìn)行比較.實驗成果表明，在分割閾值數(shù)量較多時，CMA-GSO優(yōu)化多閾值Otsu算法的計算速度比經(jīng)典多閾值Otsu算法提高95%.同時，本文算法在分割的精確度上高于GSO優(yōu)化多閾值Otsu算法.該方法已經(jīng)應(yīng)用于腦部切片(MRI)圖，復(fù)雜風(fēng)景(Boats)圖，行星(Saturn)圖和高對比度人物(Camera)圖，并且取得了快速有效的分割結(jié)果.