梁志剛，顧軍華,2

(1.河北工業(yè)大學(xué) 電氣工程學(xué)院，天津300401； 2.河北工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與軟件學(xué)院，天津 300401)

0 引言

圖像配準(zhǔn)是醫(yī)學(xué)圖像分析領(lǐng)域的一個(gè)基本問題。臨床上醫(yī)學(xué)圖像輔助疾病診斷、確定治療方案、手術(shù)模擬和術(shù)中導(dǎo)航等[1-2]很多應(yīng)用，都需要圖像的精確配準(zhǔn)作為先決條件。醫(yī)學(xué)圖像配準(zhǔn)是指通過求解兩幅或多幅圖像間的變換關(guān)系，使得同一患者不同時(shí)間、不同成像方式得到的圖像或者圖像中具有診斷意義的特征點(diǎn)在空間上對(duì)應(yīng)起來[3]。圖像配準(zhǔn)的目的是綜合利用多幅圖像對(duì)同一器官的表達(dá)，在一幅圖像中同時(shí)反映人體器官解剖結(jié)構(gòu)、代謝功能、病理發(fā)展等多方面的信息。隨著醫(yī)學(xué)影像技術(shù)的廣泛應(yīng)用，圖像配準(zhǔn)在現(xiàn)代醫(yī)學(xué)中發(fā)揮著越來越重要的作用，已經(jīng)成為當(dāng)前醫(yī)學(xué)圖像分析領(lǐng)域的一個(gè)熱點(diǎn)問題。

醫(yī)學(xué)圖像配準(zhǔn)就是通過優(yōu)化算法搜索兩幅或多幅圖像之間變換參數(shù)的過程，本質(zhì)上是一個(gè)參數(shù)尋優(yōu)問題。圖像的互信息(Mutual Information, MI)值通常作為優(yōu)化問題的目標(biāo)函數(shù)，常用的優(yōu)化算法包括Powell算法[4]、單純形法[5]、遺傳算法(Genetic Algorithm, GA)[6]、蟻群優(yōu)化(Ant Colony Optimization, ACO)算法[7]、粒子群優(yōu)化(Particle Swarm Optimization, PSO)算法[8-9]等。由于互信息存在多個(gè)局部極值和計(jì)算負(fù)荷大的問題，因此近年來多選擇基于多分辨率分析的優(yōu)化策略，并采用全局和局部?jī)?yōu)化算法結(jié)合的混合算法以提高速度和精度。楊帆等[10]運(yùn)用蟻群算法和Powell 算法相結(jié)合的多分辨率搜索算法，對(duì)三維計(jì)算機(jī)斷層掃描成像(Computed Tomography, CT)和磁共振成像(Magnetic Resonance Imaging, MRI)圖像進(jìn)行配準(zhǔn)，得到了亞像素級(jí)的配準(zhǔn)精度。李超等[11]應(yīng)用Logistic混沌映射生成可控且可再現(xiàn)遺傳算法初始粒子，基于Powell算法與遺傳算法結(jié)合的混合算法在多分辨率策略下實(shí)現(xiàn)了配準(zhǔn)，優(yōu)化了配準(zhǔn)精度和性能。趙海峰等[12]提出了一種改進(jìn)的人工魚群算法和Powell算法結(jié)合的多分辨率醫(yī)學(xué)圖像配準(zhǔn)算法，提高了配準(zhǔn)的精確度和性能。呂曉琪等[13]利用PSO 算法在小波分解低分辨率層粗配準(zhǔn)后，將配準(zhǔn)參數(shù)作為Powell 算法在高分辨率層的配準(zhǔn)起始點(diǎn)，實(shí)現(xiàn)了快速精確配準(zhǔn)。由于蟻群算法和粒子群算法存在群聚效應(yīng)，容易出現(xiàn)早熟現(xiàn)象，遺傳算法和人工魚群算法計(jì)算復(fù)雜且精度較差，無法為局部搜索提供可靠的初始解，并且四種算法在極值點(diǎn)附近搜索易出現(xiàn)震蕩現(xiàn)象，存在收斂速度慢的問題。

為了解決以上算法精度較差、易陷入局部極值和收斂速度慢的問題，本文結(jié)合多分辨率分析，提出改進(jìn)頭腦風(fēng)暴優(yōu)化(Modified Brain Storm Optimization, MBSO)算法和Powell算法結(jié)合的圖像配準(zhǔn)算法。頭腦風(fēng)暴優(yōu)化(Brain Storm Optimization, BSO)算法是Shi[14]于2011年提出的群體智能優(yōu)化算法，本文提出的MBSO算法對(duì)BSO算法全局和局部搜索的調(diào)節(jié)機(jī)制進(jìn)行了改進(jìn)，并采用可變步長(zhǎng)，能夠以較好的精度快速定位問題全局最優(yōu)解。本文算法在圖像小波分解的低分辨率層采用MBSO算法進(jìn)行全局搜索，快速確定全局配準(zhǔn)參數(shù)，然后將此參數(shù)作為Powell算法初始參數(shù)，在高分辨率層進(jìn)行搜索確定圖像最優(yōu)配準(zhǔn)參數(shù)。單模和多模圖像配準(zhǔn)實(shí)驗(yàn)表明，本文算法能夠快速、準(zhǔn)確確定圖像配準(zhǔn)參數(shù)。

1 互信息和多分辨率圖像配準(zhǔn)

圖像配準(zhǔn)過程中坐標(biāo)和角度保持不變的圖像稱為參考圖像，記為R，待配準(zhǔn)圖像稱為浮動(dòng)圖像，記為F。令I(lǐng)R(x,y)和IF(x,y)，(x,y)∈Ω?R2，分別表示參考圖像和浮動(dòng)圖像在(x,y)處的灰度值，Ω表示圖像空間。令Tp(x,y)表示作用在浮動(dòng)圖像F坐標(biāo)(x,y)處像素的剛體變換，p=(tx,ty,θ)表示此變換參數(shù)的向量形式，其中(tx,ty)表示x和y方向的平移，θ為旋轉(zhuǎn)角度，變換如式(1)所示：

(1)

對(duì)R和F進(jìn)行配準(zhǔn)，實(shí)質(zhì)上尋找最優(yōu)幾何變換T，使得圖像間相似性測(cè)度函數(shù)取得極值：

(2)

其中MI()表示相似性測(cè)度。

1.1 互信息原理

基于最大互信息的圖像配準(zhǔn)方法無需對(duì)圖像進(jìn)行預(yù)處理和特征點(diǎn)提取，具有魯棒性強(qiáng)、精度高的優(yōu)點(diǎn)，已經(jīng)被廣泛應(yīng)用于醫(yī)學(xué)圖像配準(zhǔn)研究中?；バ畔?MI)是信息論中基于熵的測(cè)度，描述兩個(gè)隨機(jī)變量相關(guān)性的度量函數(shù)。醫(yī)學(xué)圖像的互信息定義為兩幅間相似程度，即一幅圖像中包含另一幅圖像信息量的多少，互信息值越大，則兩幅圖像相關(guān)性越高。一幅圖像R的熵表示為：

(3)

其中:N表示圖像的最大灰度級(jí)別，每個(gè)灰度級(jí)別i=1,2,…,N出現(xiàn)的概率pi=hi/N，hi為灰度級(jí)別i的個(gè)數(shù)。圖像R和F的聯(lián)合熵由式(4)給出：

(4)

其中：pR,F(r,f)是將R和F兩幅圖像的灰度值作為隨機(jī)變量聯(lián)合概率分布。圖像R和F的互信息由式(5)表示：

MI(R,F)=H(R)+H(F)-H(R,F)

(5)

當(dāng)兩幅圖像達(dá)到最佳配準(zhǔn)位置時(shí)，聯(lián)合熵達(dá)到最小值，互信息取得最大值。MI并沒有解決重疊區(qū)域的問題，尤其是當(dāng)重疊區(qū)域灰度值較小時(shí)會(huì)影響MI的均衡性能，錯(cuò)誤的配準(zhǔn)同樣會(huì)使互信息值增大。為了抵消重疊區(qū)域的影響，Studholme引入歸一化互信息，如式(6)所示：

(6)

歸一化互信息抵消了重疊區(qū)域的影響，對(duì)于配準(zhǔn)結(jié)果的量化表示更加魯棒和精確，本文采用歸一化互信息作為圖像配準(zhǔn)的測(cè)度。

1.2 多分辨率分析

小波變換是一種時(shí)域和頻域綜合圖像分析方法，應(yīng)用多分辨率將圖像分解為不同尺寸、不同分辨率的圖像，在較低分辨率圖像上進(jìn)行配準(zhǔn)，可以大大降低配準(zhǔn)過程中計(jì)算聯(lián)合直方圖和插值的時(shí)間、空間等計(jì)算資源的消耗。對(duì)于分辨率為m×n二維圖像IR(x,y)，通過式(7)可以得到圖像的低頻部分：

Ik(x,y)=∑∑Ik-1(m,n)h(2m-x)l(2n-y)

(7)

其中:k表示分解層數(shù)，h和l分別表示小波分解高通濾波和低通濾波尺度方程系數(shù)。

圖像經(jīng)過小波分解后低頻近似分量保存了圖像的主要信息，在低分辨率層運(yùn)用全局搜索算法確定配準(zhǔn)參數(shù)粗略位置后，上升到高分辨率層運(yùn)用局部搜索算法進(jìn)一步細(xì)化配準(zhǔn)結(jié)果。低分辨率層運(yùn)用全局算法搜索，較少的數(shù)據(jù)量提高了搜索效率，并且使算法跳出局部極值的干擾，為高分辨率層局部搜索提供準(zhǔn)確初始點(diǎn)信息，增強(qiáng)整個(gè)搜索算法的魯棒性、提高精度。圖像經(jīng)過小波分解后，配準(zhǔn)參數(shù)關(guān)系如下：高分辨率層平移參數(shù)是下一級(jí)低分辨率層變換參數(shù)的2倍，縮放和旋轉(zhuǎn)參數(shù)在各分辨率層相等?；诨バ畔⒌亩嗄D像配準(zhǔn)中，對(duì)于圖像只需分解2～3層即可達(dá)到配準(zhǔn)精度和時(shí)間的要求。

2 優(yōu)化算法

2.1 Powell算法

Powell算法可用于求解一般無約束優(yōu)化問題，無約束最小化優(yōu)化問題可以描述為：

minσ=f(x);x∈Rn

(8)

其中:σ=f(x)是優(yōu)化問題的目標(biāo)函數(shù)，x是n維優(yōu)化自變量，Rn是優(yōu)化問題的可行域空間。

Powell算法是一種十分有效的局部直接搜索算法，算法將多維優(yōu)化問題轉(zhuǎn)化為一維迭代尋優(yōu)問題，利用共軛方向可以加快收斂速度的性質(zhì)，在每一維中使用Brent算法進(jìn)行循環(huán)搜索。在尋優(yōu)過程中，Powell算法不需要對(duì)目標(biāo)函數(shù)進(jìn)行求導(dǎo)，當(dāng)目標(biāo)函數(shù)的導(dǎo)數(shù)不連續(xù)時(shí)也能應(yīng)用。

對(duì)于(8)的最小化優(yōu)化問題，給定的初始解，首先沿給定的n個(gè)共軛方向進(jìn)行搜索，接著沿初始點(diǎn)與搜索終點(diǎn)連線方向進(jìn)行n+1輪搜索，最后用第n+1次方向替換前n個(gè)方向中其中一個(gè)，開始新的迭代過程。Powell算法在局部搜索中具有很高的精度，并且搜索速度非?？?，但是作為一種能力極強(qiáng)的局部搜索算法，依賴于初始點(diǎn)的選擇，當(dāng)初始點(diǎn)選擇在局部最優(yōu)附近時(shí)，算法無法跳出局部最優(yōu)的干擾，很容易陷入局部最優(yōu)。為了克服Powell算法的這一缺點(diǎn)，首先運(yùn)用頭腦風(fēng)暴優(yōu)化算法進(jìn)行全局搜索，確定全局極值的粗略位置，作為Powell算法的初始點(diǎn)進(jìn)行局部細(xì)致搜索，兩種算法的結(jié)合可以實(shí)現(xiàn)快速、精確、魯棒的圖像配準(zhǔn)。

2.2 頭腦風(fēng)暴優(yōu)化算法

BSO算法是對(duì)人類社會(huì)行為進(jìn)行抽象產(chǎn)生的群體智能算法，算法的靈感來自于人類創(chuàng)造性解決問題的方法——頭腦風(fēng)暴(Brain storming)法。頭腦風(fēng)暴法指具有不同專業(yè)背景的人或小組圍繞一個(gè)特定的興趣或領(lǐng)域，進(jìn)行發(fā)散思維、思想碰撞、觀點(diǎn)融合，以找到問題的最優(yōu)解決方案，是一種激發(fā)集體智慧產(chǎn)生和提出創(chuàng)新設(shè)想的思維方法。BSO算法將頭腦風(fēng)暴法的不同階段抽象為聚類、變異、生成新個(gè)體、選擇四種操作的不斷迭代過程，能夠跳出局部極值，快速定位優(yōu)化問題的最優(yōu)解。BSO已經(jīng)在許多優(yōu)化問題中取得了成功，算法流程如下。

在優(yōu)化問題的定義域內(nèi)隨機(jī)產(chǎn)生n個(gè)個(gè)體組成的種群，每個(gè)個(gè)體都是d維空間的一個(gè)向量，第i個(gè)個(gè)體表示為：

Xi=(xi1,xi2,…,xid);i=1,2,…,n

(9)

個(gè)體初始化后開始1)～5)的迭代過程，每輪迭代過程中生成n個(gè)新個(gè)體。

1)種群聚類：種群的個(gè)體采用k-means聚類為m個(gè)子群，根據(jù)適應(yīng)度函數(shù)對(duì)個(gè)體進(jìn)行排序，每個(gè)子群中適應(yīng)度值最好的定義為中心個(gè)體，其余為一般個(gè)體；

2)個(gè)體變異：隨機(jī)產(chǎn)生一個(gè)[0,1]的隨機(jī)數(shù)p11，若大于預(yù)先設(shè)定的概率參數(shù)p1，隨機(jī)選中一個(gè)子群的中心個(gè)體，替換為隨機(jī)生成的新個(gè)體；

3)生成新個(gè)體：預(yù)先設(shè)定概率p2，隨機(jī)生成[0,1]之間的隨機(jī)數(shù)p2i，與p2進(jìn)行比較，并通過以下兩種方式產(chǎn)生候選個(gè)體：

①如果p2i≥p2,隨機(jī)選中一個(gè)子群，選擇子群中心個(gè)體或隨機(jī)選擇一個(gè)一般個(gè)體，作為候選個(gè)體Xs；

②如果p2i

Xs=ω×Xs1+(1-ω)×Xs2

(10)

候選個(gè)體Xs加上高斯隨機(jī)擾動(dòng)生成新個(gè)體Y，如式(11)和(12)所示：

Y=Xs+ξ×N(μ,σ)

(11)

ξ=logsig((0.5×Tmax-t)/S)×rand(0,1)

(12)

其中，t為當(dāng)前迭代次數(shù)，Tmax是最大迭代次數(shù)，N(μ,σ)是中心在μ方差為σ的高斯隨機(jī)函數(shù)，ξ為對(duì)數(shù)S變換函數(shù)；

4)個(gè)體選擇：新個(gè)體Y的適應(yīng)度與當(dāng)前需要更新的個(gè)體進(jìn)行比較，適應(yīng)度好的個(gè)體進(jìn)入下一代，n個(gè)個(gè)體全部更新完后轉(zhuǎn)5)，否則轉(zhuǎn)3)進(jìn)行個(gè)體更新；

5)若個(gè)體達(dá)到最優(yōu)解條件或達(dá)到預(yù)先設(shè)定的迭代次數(shù)時(shí)，算法結(jié)束，否則轉(zhuǎn)1)開始新的迭代過程。

2.3 改進(jìn)頭腦風(fēng)暴算法

由BSO算法迭代過程可以看出，算法新個(gè)體生成方式和步長(zhǎng)是算法收斂速度關(guān)鍵因素。MBSO從兩個(gè)方面進(jìn)行改進(jìn)：一是調(diào)節(jié)新個(gè)體生成方式，二是生成新個(gè)體過程中的步長(zhǎng)調(diào)節(jié)方法。

在BSO算法迭代過程3)生成新個(gè)體時(shí)：①中選擇一個(gè)子群的中心或一般個(gè)體作為候選個(gè)體Xs，稱由此候選個(gè)體加隨機(jī)擾動(dòng)生成新個(gè)體的方式為a方式；②中選擇兩個(gè)子群中心或兩個(gè)一般個(gè)體融合產(chǎn)生新個(gè)體作為候選個(gè)體Xs，稱由此候選個(gè)體加隨機(jī)擾動(dòng)生成新個(gè)體的方式為b方式。分析可知：b方式生成的個(gè)體大概率出現(xiàn)在兩個(gè)子群的中間，適合于搜索初期進(jìn)行大范圍的全局搜索，使算法快速收斂至全局極值區(qū)域；a方式生成的新個(gè)體大概率出現(xiàn)在候選個(gè)體的附近，適合于后期局部細(xì)致搜索定位極值點(diǎn)。因此，通過調(diào)節(jié)每輪迭代過程中b和a方式生成新個(gè)體的比例來控制算法全局和局部搜索的能力，前期以b方式生成新個(gè)體為主加強(qiáng)全局搜索，后期以a方式生成新個(gè)體為主加強(qiáng)局部細(xì)致搜索。通過調(diào)節(jié)概率參數(shù)p2可以達(dá)到此目的，改進(jìn)后的概率參數(shù)p2變化形式如下：

p2=k×p20

(13)

k=exp(-α(t/Tmax)β)

(14)

其中：p20表示初始概率參數(shù)，k是更新因子，t為當(dāng)前迭代次數(shù)，Tmax是最大迭代次數(shù)，α和β是正整數(shù)。

式(12)中，ξ的兩個(gè)分量logsig函數(shù)和rand函數(shù)取值都在(0,1)，乘以均值為0、方差為1的高斯隨機(jī)函數(shù)后，取值范圍為[-4,4]。logsig函數(shù)在Tmax/2處存在一個(gè)階躍，且高斯函數(shù)特點(diǎn)決定取值均值附近概率較大。這些因素導(dǎo)致步長(zhǎng)與算法所處階段不匹配，造成算法前期收斂較慢，而后期在極值附近反復(fù)震蕩，無法快速收斂。因此，可以調(diào)節(jié)高斯隨機(jī)函數(shù)的方差以調(diào)節(jié)算法的搜索步長(zhǎng)，改進(jìn)后高斯函數(shù)方差變化形式為：

σ=kσ0

(15)

其中:σ0為高斯函數(shù)初始方差，k取值與式(14)一致。MBSO算法流程如圖1所示。

圖1 MBSO算法流程

3 本文算法

醫(yī)學(xué)圖像配準(zhǔn)的研究目標(biāo)是要尋找一種精度高、魯棒性強(qiáng)、速度快、計(jì)算代價(jià)較小的優(yōu)化算法，通過算法尋找配準(zhǔn)參數(shù)將兩幅圖像在空間上嚴(yán)格對(duì)應(yīng)起來。MBSO算法對(duì)起始點(diǎn)沒有要求，具有極強(qiáng)的魯棒性，能夠快速定位全局最優(yōu)的粗略位置，但是算法復(fù)雜度、高收斂緩慢。Powell算法具有極強(qiáng)的局部搜索能力且收斂速度快，但是依賴于初始點(diǎn)的選擇，易陷入局部極值。本文在對(duì)圖像進(jìn)行小波分解的基礎(chǔ)上，將MBSO的全局尋優(yōu)能力和Powell算法的局部搜索能力結(jié)合起來，采用二者的混合優(yōu)化算法實(shí)現(xiàn)圖像配準(zhǔn)任務(wù)。對(duì)原始圖像進(jìn)行三層小波分解后，在低分辨率層采用MBSO算法進(jìn)行全局搜索，此時(shí)由于圖像較小可以快速確定全局最優(yōu)值的大概位置，接下來將MBSO算法的優(yōu)化結(jié)果作為初始點(diǎn)，在第2層和原始圖像層采用Powell算法進(jìn)行局部細(xì)致搜索，Powell算法的局部搜索能力可以快速定位全局最優(yōu)的位置。多分辨率的加入可以有效加快MBSO算法的搜索速度，MBSO算法的搜索結(jié)果作為Powell算法的初始值，可以有效防止Powell算法陷入局部極值中。算法的具體步驟如下。

1)對(duì)于參考圖像R和浮動(dòng)圖像F進(jìn)行三層小波分解，原始圖像作為第1層圖像，小波分解第2、3層近似系數(shù)作為第2、3層圖像。原始圖像的配準(zhǔn)可以轉(zhuǎn)化為第2、3層圖像的配準(zhǔn)問題，相鄰兩層平移、旋轉(zhuǎn)參數(shù)關(guān)系為：低分辨率層圖像平移參數(shù)為(tx,ty)時(shí)，前一層圖像平移參數(shù)為(2tx,2ty)，而旋轉(zhuǎn)角度參數(shù)θ保持不變。

2)利用MBSO算法對(duì)第3層圖像進(jìn)行配準(zhǔn)，將兩幅圖像的歸一化互信息(Normalized Mutual Information, NMI)作為適應(yīng)度函數(shù)，圖像灰度直方圖灰度區(qū)間設(shè)置為128，采用最近鄰法作為插值方法，在低分辨率層采用這樣設(shè)置有助于提高搜索速度而不影響搜索的精度。MBSO算法搜索結(jié)束后，輸出最優(yōu)配準(zhǔn)參數(shù)。

3)將步驟2)中MBSO算法得到的配準(zhǔn)參數(shù)作為Powell算法搜索的起始點(diǎn)在第2層繼續(xù)搜索。Powell算法的迭代次數(shù)設(shè)為200，初始步長(zhǎng)設(shè)為1，允許誤差設(shè)為0.001，圖像灰度直方圖灰度區(qū)間設(shè)置為256，采用部分體積(Partial Volume, PV)插值方法，算法結(jié)束時(shí)輸出本層最優(yōu)配準(zhǔn)參數(shù)。

4)將步驟3)得到的最優(yōu)配準(zhǔn)參數(shù)作為Powell算法的初始點(diǎn)，重復(fù)步驟3在原始圖像層進(jìn)行搜索，得到最優(yōu)圖像配準(zhǔn)參數(shù)(tx,ty,θ)。

5)用步驟4)得到的配準(zhǔn)參數(shù)(tx,ty,θ)對(duì)浮動(dòng)圖像F進(jìn)行空間變換，使之與參考圖像完全對(duì)應(yīng)起來，完成圖像配準(zhǔn)過程。

步驟2)中在第3層采用MBSO搜索的流程如下：

①算法初始化：在解空間隨機(jī)生成均勻分布的50個(gè)個(gè)體組成初始種群，個(gè)體形為(tx,ty,θ)，算法最大迭代次數(shù)設(shè)為200，算法終止適應(yīng)度閾值設(shè)為0.95，預(yù)設(shè)概率p1為0.5，p20初始設(shè)置為0.85，高斯函數(shù)參數(shù)均值μ設(shè)置為0，初始方差σ0設(shè)置為3，更新因子k的參數(shù)α設(shè)置為3，β設(shè)置為2，當(dāng)前迭代代數(shù)t為1，用個(gè)體對(duì)浮動(dòng)圖像進(jìn)行變換，計(jì)算參考和浮動(dòng)圖像的NMI作為個(gè)體適應(yīng)度值。

②采用k-means將種群聚類為m個(gè)子群，對(duì)個(gè)體按適應(yīng)度值排序并確定子群中心個(gè)體。

③依概率選中任一子群，用隨機(jī)生成的個(gè)體替換該子群中心個(gè)體。根據(jù)當(dāng)前迭代次數(shù)計(jì)算概率參數(shù)p2和高斯函數(shù)方差σ。

④隨機(jī)生成[0,1]的隨機(jī)數(shù)p2i，如果p2i>p2，則按照a方式選擇一個(gè)子群中個(gè)體作為候選個(gè)體Xs，否則按照b方式選擇兩個(gè)子群中個(gè)體融合作為候選個(gè)體Xs，候選個(gè)體加上隨機(jī)擾動(dòng)生成新個(gè)體Y。

⑤用新個(gè)體Y對(duì)浮動(dòng)圖像進(jìn)行變換后，計(jì)算參考和浮動(dòng)圖像的NMI作為Y適應(yīng)度值，若該值小于待更新個(gè)體X的適應(yīng)度值，將Y舍棄，否則用Y替換X。如果當(dāng)前最大適應(yīng)度值大于閾值0.95，則轉(zhuǎn)⑦輸出結(jié)果，如果50個(gè)個(gè)體全部更新完畢，轉(zhuǎn)⑥，否則轉(zhuǎn)④生成新個(gè)體。

⑥如果已經(jīng)達(dá)到最大迭代次數(shù)，則轉(zhuǎn)⑦，否則轉(zhuǎn)②開始新一輪迭代。

⑦終止MBSO算法搜索過程，將當(dāng)前最優(yōu)適應(yīng)度函數(shù)值的個(gè)體作為結(jié)果輸出。

4 實(shí)驗(yàn)與分析

為了驗(yàn)證本文配準(zhǔn)算法的性能，采用加拿大McGill大學(xué)BrainWeb數(shù)據(jù)庫(kù)[15]的模擬腦MRI圖像，以及哈佛大學(xué)醫(yī)學(xué)院的腦部數(shù)據(jù)[16]分別進(jìn)行單模和多模圖像配準(zhǔn)實(shí)驗(yàn)。實(shí)驗(yàn)平臺(tái)選用Windows操作系統(tǒng)下的Matlab R2012B科學(xué)計(jì)算軟件，硬件配置為Intel Core i3-4170 CPU 3.70 GHz、內(nèi)存4.00 GB。配準(zhǔn)效果可以從定性和定量分析[17-18]兩個(gè)方面進(jìn)行評(píng)價(jià)，本文采用目測(cè)效果進(jìn)行定性評(píng)價(jià)，采用配準(zhǔn)成功次數(shù)、精準(zhǔn)度和配準(zhǔn)時(shí)間三個(gè)指標(biāo)對(duì)算法進(jìn)行定量評(píng)價(jià)。采用均方根誤差(Root Mean Square Error, RMSE)作為精度的度量，RMSE值越小則圖像配準(zhǔn)效果越好。

選擇文獻(xiàn)[10-11,13]中配準(zhǔn)算法作為對(duì)比算法，與本文算法進(jìn)行比較。三種算法分別是粒子群優(yōu)化(PSO)算法、蟻群優(yōu)化(ACO)算法和遺傳算法(GA)與Powell算法結(jié)合的混合算法，以下記三種算法為PSO+Powell、ACO+Powell和GA+Powell。ACO算法、GA、PSO算法種群個(gè)體數(shù)量設(shè)置為50，最大迭代次數(shù)設(shè)置為200，算法其他參數(shù)與文獻(xiàn)[10-11,13]中保持一致。Powell算法的迭代次數(shù)設(shè)為200，初始迭代步長(zhǎng)設(shè)為1，允許誤差設(shè)為0.001。

4.1 單模圖像配準(zhǔn)

選取BrianWeb數(shù)據(jù)庫(kù)中MRI-T1腦部橫斷面圖像進(jìn)行單模圖像配準(zhǔn)實(shí)驗(yàn)。如圖2所示，圖2(a)是MRI圖像，作為參考圖像R，圖2(b)將R分別沿水平、豎直方向平移10個(gè)和8個(gè)像素，繞中心旋轉(zhuǎn)10°，得到浮動(dòng)圖像F，圖2(c)～(f)分別是PSO+Powell、ACO+Powell、GA+Powell和本文算法的配準(zhǔn)結(jié)果。

圖2 單模腦部圖像配準(zhǔn)結(jié)果

單模圖像配準(zhǔn)實(shí)驗(yàn)的定量分析中，對(duì)圖2中參考圖像R隨機(jī)施加±20個(gè)像素的平移和±20°的旋轉(zhuǎn)，重復(fù)進(jìn)行100次的配準(zhǔn)實(shí)驗(yàn)，從配準(zhǔn)成功次數(shù)、配準(zhǔn)精度和平均配準(zhǔn)時(shí)間對(duì)算法效果進(jìn)行評(píng)價(jià)。每次實(shí)驗(yàn)中計(jì)算配準(zhǔn)變換參數(shù)與真實(shí)變換參數(shù)的絕對(duì)差值(Δx,Δy,Δθ)，當(dāng)平移誤差小于1個(gè)像素、旋轉(zhuǎn)誤差小于1°時(shí)，稱此次配準(zhǔn)達(dá)到亞像素級(jí)，視為配準(zhǔn)成功。

表1列出了100次實(shí)驗(yàn)中本文算法和對(duì)比算法配準(zhǔn)成功的次數(shù)、平均RMSE和平均配準(zhǔn)時(shí)間。綜合來看，PSO+Powell算法精度在四種算法排在第二，平均運(yùn)行時(shí)間第二短。GA+Powell算法平均運(yùn)行時(shí)間最長(zhǎng)，配準(zhǔn)成功率排在第二。ACO+Powell算法配準(zhǔn)成功率最低，精度也最差，平均運(yùn)行時(shí)間排在第三。本文算法效果最好，相比PSO+Powell、ACO+Powell和GA+Powell算法，平均RMSE分別下降了28.81%、46.15%和38.24%，平均配準(zhǔn)時(shí)間分別縮短了18.08%、24.38%和27.47%。

表1 單模圖像配準(zhǔn)實(shí)驗(yàn)定量分析

4.2 多模醫(yī)學(xué)圖像配準(zhǔn)

選取哈佛大學(xué)醫(yī)學(xué)院的計(jì)算機(jī)斷層掃描/磁共振成像(Computed Tomography/Magnetic Resonance Imaging, CT/MRI)和正電子放射斷層成像/磁共振成像(Positron Emission Tomography/Magnetic Resonance Imaging, PET/MRI)圖像分別進(jìn)行多模配準(zhǔn)實(shí)驗(yàn)。第一組待配準(zhǔn)圖像選擇急性腦卒中患者CT/MRI圖像，如圖3所示，圖3(a)是CT圖像作為參考圖像，圖3(b)是MRI圖像作為浮動(dòng)圖像，圖3(c)～(f)分別顯示了PSO+Powell、ACO+Powell、GA+Powell和本文算法配準(zhǔn)結(jié)果。第二組圖像選擇腦部神經(jīng)膠質(zhì)瘤患者PET/MRI圖像，圖4(a)是PET圖像作為參考圖像，圖4(b)是MRI圖像作為浮動(dòng)圖像，圖4(c)～(f)分別顯示了PSO+Powell、ACO+Powell、GA+Powell和本文算法配準(zhǔn)結(jié)果。

在CT/MRI和PET/MRI配準(zhǔn)成功的基礎(chǔ)上，對(duì)本文算法多模圖像配準(zhǔn)效果進(jìn)行定量分析。實(shí)驗(yàn)中以MRI圖像配準(zhǔn)結(jié)果為基準(zhǔn)，對(duì)兩組圖像中MRI圖像分別隨機(jī)施加±20個(gè)像素的平移和±20°的旋轉(zhuǎn)，重復(fù)進(jìn)行100次的配準(zhǔn)實(shí)驗(yàn)，從配準(zhǔn)成功次數(shù)、配準(zhǔn)精度和平均配準(zhǔn)時(shí)間對(duì)算法效果進(jìn)行評(píng)價(jià)。表2列出了CT/MRI、PET/MRI實(shí)驗(yàn)中四種算法配準(zhǔn)成功的次數(shù)、RMSE平均值和平均配準(zhǔn)時(shí)間。

由表2可以看出，在多模圖像配準(zhǔn)中，本文算法在四種算法的結(jié)果具有最好的配準(zhǔn)表現(xiàn)：平均配準(zhǔn)時(shí)間在四種算法中是最少的，并且具有100%的配準(zhǔn)成功率，RMSE值在四種算法中也是最小的。由表2腦部CT/MRI圖像的各項(xiàng)指標(biāo)來看，在CT/MRI配準(zhǔn)中，相比PSO+Powell、ACO+Powell和GA+Powell算法，本文算法RMSE分別下降了15.80%、21.95%和8.89%，平均配準(zhǔn)時(shí)間分別縮短了19.01%、30.97%和25.88%。由表2腦部PET/MRI圖像各項(xiàng)指標(biāo)來看，在PET/MRI配準(zhǔn)中，本文算法RMSE分別下降了18.06%、23.28%和8.50%，平均配準(zhǔn)時(shí)間分別縮短了16.50%、25.79%和26.46%。

單模和多模醫(yī)學(xué)圖像配準(zhǔn)結(jié)果表明，本文算法相對(duì)于PSO+Powell、ACO+Powell和GA+Powell算法，平均RMSE分別下降了20.89%、30.46%和18.54%，平均配準(zhǔn)時(shí)間分別縮短了17.86%、27.05%、26.60%，并且達(dá)到了100%的成功率。實(shí)驗(yàn)結(jié)果說明，本文算法不僅能夠有效減少配準(zhǔn)的時(shí)間，跳出局部極值的干擾，還可以提高配準(zhǔn)的成功率和準(zhǔn)確率。MBSO算法特有的對(duì)于全局和局部搜索的調(diào)節(jié)機(jī)制和對(duì)局部細(xì)致搜索能力，不僅能在低分辨率層快速定位全局最優(yōu)解的范圍，同時(shí)能夠?qū)⒆顑?yōu)解范圍盡可能地縮小，為局部搜索算法Powell提供可信賴的初始點(diǎn)。本文方法的配準(zhǔn)結(jié)果優(yōu)于PSO算法、ACO算法和GA與Powell算法結(jié)合的優(yōu)化算法的配準(zhǔn)結(jié)果，能夠快速、準(zhǔn)確、自動(dòng)完成圖像配準(zhǔn)任務(wù)。

圖3 腦部CT/MRI圖像配準(zhǔn)結(jié)果

圖4 腦部PET/MRI多模圖像配準(zhǔn)結(jié)果

圖像算法成功次數(shù)RMSE平均時(shí)間/s腦部CT/MRI圖像腦部PET/MRI圖像PSO+Powell871.218349.6584ACO+Powell841.314258.2548GA+Powell911.125854.2624本文算法1001.025740.2157PSO+Powell851.258750.9825ACO+Powell831.344557.3657GA+Powell901.127456.3542本文算法1001.031542.5715

5 結(jié)語

針對(duì)現(xiàn)有圖像配準(zhǔn)算法精度較差、易陷入局部極值和收斂速度慢的問題，本文提出了基于多分辨率分析，MBSO算法與Powell算法結(jié)合的醫(yī)學(xué)圖像配準(zhǔn)算法。算法首先對(duì)圖像進(jìn)行小波分解，在低分辨率層采用速度快、精度高M(jìn)BSO算法進(jìn)行全局搜索，在高分辨率層采用適于局部搜索的Powell算法完成圖像配準(zhǔn)任務(wù)。以歸一化互信息為相似性測(cè)度，選擇粒子群優(yōu)化算法、蟻群優(yōu)化算法、遺傳算法與Powell算法結(jié)合的混合算法作為對(duì)比算法，運(yùn)用BrainWeb模擬腦數(shù)據(jù)和哈佛大學(xué)醫(yī)學(xué)院真實(shí)腦圖像進(jìn)行了對(duì)比實(shí)驗(yàn)，本文算法在配準(zhǔn)成功率、配準(zhǔn)時(shí)間和配準(zhǔn)精度上均優(yōu)于其他三種算法。實(shí)驗(yàn)結(jié)果表明，本文算法具有很強(qiáng)的魯棒性，能夠快速、準(zhǔn)確完成單模和多模醫(yī)學(xué)圖像配準(zhǔn)任務(wù)。