改進頭腦風暴優化算法與Powell算法結合的醫學圖像配準

梁志剛，顧軍華,2

(1.河北工業大學 電氣工程學院，天津300401； 2.河北工業大學 計算機科學與軟件學院，天津 300401)

0 引言

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

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

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

1 互信息和多分辨率圖像配準

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

(1)

對R和F進行配準，實質上尋找最優幾何變換T，使得圖像間相似性測度函數取得極值：

(2)

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

1.1 互信息原理

基于最大互信息的圖像配準方法無需對圖像進行預處理和特征點提取，具有魯棒性強、精度高的優點，已經被廣泛應用于醫學圖像配準研究中。互信息(MI)是信息論中基于熵的測度，描述兩個隨機變量相關性的度量函數。醫學圖像的互信息定義為兩幅間相似程度，即一幅圖像中包含另一幅圖像信息量的多少，互信息值越大，則兩幅圖像相關性越高。一幅圖像R的熵表示為：

(3)

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

(4)

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

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

(5)

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

(6)

歸一化互信息抵消了重疊區域的影響，對于配準結果的量化表示更加魯棒和精確，本文采用歸一化互信息作為圖像配準的測度。

1.2 多分辨率分析

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

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

(7)

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

圖像經過小波分解后低頻近似分量保存了圖像的主要信息，在低分辨率層運用全局搜索算法確定配準參數粗略位置后，上升到高分辨率層運用局部搜索算法進一步細化配準結果。低分辨率層運用全局算法搜索，較少的數據量提高了搜索效率，并且使算法跳出局部極值的干擾，為高分辨率層局部搜索提供準確初始點信息，增強整個搜索算法的魯棒性、提高精度。圖像經過小波分解后，配準參數關系如下：高分辨率層平移參數是下一級低分辨率層變換參數的2倍，縮放和旋轉參數在各分辨率層相等。基于互信息的多模圖像配準中，對于圖像只需分解2～3層即可達到配準精度和時間的要求。

2 優化算法

2.1 Powell算法

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

minσ=f(x);x∈Rn

(8)

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

Powell算法是一種十分有效的局部直接搜索算法，算法將多維優化問題轉化為一維迭代尋優問題，利用共軛方向可以加快收斂速度的性質，在每一維中使用Brent算法進行循環搜索。在尋優過程中，Powell算法不需要對目標函數進行求導，當目標函數的導數不連續時也能應用。

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

2.2 頭腦風暴優化算法

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

在優化問題的定義域內隨機產生n個個體組成的種群，每個個體都是d維空間的一個向量，第i個個體表示為：

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

(9)

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

1)種群聚類：種群的個體采用k-means聚類為m個子群，根據適應度函數對個體進行排序，每個子群中適應度值最好的定義為中心個體，其余為一般個體；

2)個體變異：隨機產生一個[0,1]的隨機數p11，若大于預先設定的概率參數p1，隨機選中一個子群的中心個體，替換為隨機生成的新個體；

3)生成新個體：預先設定概率p2，隨機生成[0,1]之間的隨機數p2i，與p2進行比較，并通過以下兩種方式產生候選個體：

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

②如果p2i

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

(10)

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

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

(11)

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

(12)

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

4)個體選擇：新個體Y的適應度與當前需要更新的個體進行比較，適應度好的個體進入下一代，n個個體全部更新完后轉5)，否則轉3)進行個體更新；

5)若個體達到最優解條件或達到預先設定的迭代次數時，算法結束，否則轉1)開始新的迭代過程。

2.3 改進頭腦風暴算法

由BSO算法迭代過程可以看出，算法新個體生成方式和步長是算法收斂速度關鍵因素。MBSO從兩個方面進行改進：一是調節新個體生成方式，二是生成新個體過程中的步長調節方法。

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

p2=k×p20

(13)

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

(14)

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

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

σ=kσ0

(15)

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

圖1 MBSO算法流程

3 本文算法

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

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

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

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

4)將步驟3)得到的最優配準參數作為Powell算法的初始點，重復步驟3在原始圖像層進行搜索，得到最優圖像配準參數(tx,ty,θ)。

5)用步驟4)得到的配準參數(tx,ty,θ)對浮動圖像F進行空間變換，使之與參考圖像完全對應起來，完成圖像配準過程。

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

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

②采用k-means將種群聚類為m個子群，對個體按適應度值排序并確定子群中心個體。

③依概率選中任一子群，用隨機生成的個體替換該子群中心個體。根據當前迭代次數計算概率參數p2和高斯函數方差σ。

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

⑤用新個體Y對浮動圖像進行變換后，計算參考和浮動圖像的NMI作為Y適應度值，若該值小于待更新個體X的適應度值，將Y舍棄，否則用Y替換X。如果當前最大適應度值大于閾值0.95，則轉⑦輸出結果，如果50個個體全部更新完畢，轉⑥，否則轉④生成新個體。

⑥如果已經達到最大迭代次數，則轉⑦，否則轉②開始新一輪迭代。

⑦終止MBSO算法搜索過程，將當前最優適應度函數值的個體作為結果輸出。

4 實驗與分析

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

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

4.1 單模圖像配準

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

圖2 單模腦部圖像配準結果

單模圖像配準實驗的定量分析中，對圖2中參考圖像R隨機施加±20個像素的平移和±20°的旋轉，重復進行100次的配準實驗，從配準成功次數、配準精度和平均配準時間對算法效果進行評價。每次實驗中計算配準變換參數與真實變換參數的絕對差值(Δx,Δy,Δθ)，當平移誤差小于1個像素、旋轉誤差小于1°時，稱此次配準達到亞像素級，視為配準成功。

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

表1 單模圖像配準實驗定量分析

4.2 多模醫學圖像配準

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

在CT/MRI和PET/MRI配準成功的基礎上，對本文算法多模圖像配準效果進行定量分析。實驗中以MRI圖像配準結果為基準，對兩組圖像中MRI圖像分別隨機施加±20個像素的平移和±20°的旋轉，重復進行100次的配準實驗，從配準成功次數、配準精度和平均配準時間對算法效果進行評價。表2列出了CT/MRI、PET/MRI實驗中四種算法配準成功的次數、RMSE平均值和平均配準時間。

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

單模和多模醫學圖像配準結果表明，本文算法相對于PSO+Powell、ACO+Powell和GA+Powell算法，平均RMSE分別下降了20.89%、30.46%和18.54%，平均配準時間分別縮短了17.86%、27.05%、26.60%，并且達到了100%的成功率。實驗結果說明，本文算法不僅能夠有效減少配準的時間，跳出局部極值的干擾，還可以提高配準的成功率和準確率。MBSO算法特有的對于全局和局部搜索的調節機制和對局部細致搜索能力，不僅能在低分辨率層快速定位全局最優解的范圍，同時能夠將最優解范圍盡可能地縮小，為局部搜索算法Powell提供可信賴的初始點。本文方法的配準結果優于PSO算法、ACO算法和GA與Powell算法結合的優化算法的配準結果，能夠快速、準確、自動完成圖像配準任務。

圖3 腦部CT/MRI圖像配準結果

圖4 腦部PET/MRI多模圖像配準結果

圖像算法成功次數RMSE平均時間/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 結語

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