基于最大互信息和量子粒子群優(yōu)化算法的醫(yī)學(xué)圖像配準(zhǔn)研究

(江南大學(xué) a.教育學(xué)院; b.信息學(xué)院， 江蘇 無(wú)錫 214422)



摘 要：

研究了基于最大互信息的圖像配準(zhǔn)算法，在圖像配準(zhǔn)中引入了新的相似性測(cè)度，在分析具有量子行為的粒子群優(yōu)化算法基礎(chǔ)上，將量子粒子群算法作為優(yōu)化策略用于圖像配準(zhǔn)并與Powell算法和PSO算法進(jìn)行了仿真比較，對(duì)仿真結(jié)果進(jìn)行了分析。

 關(guān)鍵詞：互信息； 圖像配準(zhǔn)； 量子粒子群算法

 中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)志碼： A

文章編號(hào)：10013695(2008)12368503



Research of medical image registration based on max mutual information and quantumbehaved particle swarm algorithm

WANG Xiaogena， XU Wenbob



(a.School of Education， b. School of Information Technology， Jiangnan University， Wuxi Jiangsu 214422， China)



Abstract:

This paper researched the image registration algorithm and introduced a new similar measure into image registration. On the basis of quantumbehaved particle swarm optimization， proposed QPSO algorithm to used into image registration and shown the simulation comparison with Powell and PSO algorithm. Finally gave analysis of simulation results.

 Key words：mutual information; image registration; quantumbehaved particle swarm optimization



隨著現(xiàn)代數(shù)字醫(yī)療技術(shù)的飛速發(fā)展，圖像處理被廣泛應(yīng)用于臨床診斷和治療中^[1]，已經(jīng)成為現(xiàn)代醫(yī)療不可或缺的一部分。由于成像原理和設(shè)備的不同，在醫(yī)療圖像處理中存在著多種成像模式，但不同的成像模式對(duì)人體同一解剖位置得到的圖像信息通常具有很強(qiáng)的互補(bǔ)性。臨床上希望能對(duì)這些不同的圖像信息進(jìn)行適當(dāng)程度的整合，以便提供更為全面的診斷信息^[2]。 整合的第一步就是使這些圖像信息在空間域中達(dá)到空間幾何位置的完全一致，稱之為配準(zhǔn)；第二步是已配準(zhǔn)的圖像信息融合成一個(gè)新的圖像模態(tài)并顯示出來(lái)，這就是所說(shuō)的圖像融合。 

進(jìn)行圖像配準(zhǔn)的過(guò)程實(shí)際上就是尋求兩幅待配圖像間的一一映射的過(guò)程。它需要將不同圖像中對(duì)應(yīng)于空間同一位置的點(diǎn)通過(guò)某種幾何變換（旋轉(zhuǎn)和平移）使其在空間中嚴(yán)格對(duì)齊，所以圖像配準(zhǔn)問(wèn)題本質(zhì)上是多參數(shù)優(yōu)化問(wèn)題^[3，4]，優(yōu)化算法的選擇至關(guān)重要。本文采用基于最大互信息的相似性測(cè)度和量子粒子群算法進(jìn)行了醫(yī)學(xué)圖像配準(zhǔn)研究，取得了較好的配準(zhǔn)效果。

1 具有量子行為的粒子群優(yōu)化算法

粒子群優(yōu)化算法是在1995年由Kennedy等人^[5]提出的。粒子群算法源于對(duì)生物群體的研究，在群體中，個(gè)體（粒子）通過(guò)搜索多維空間，在每一輪迭代中評(píng)價(jià)自身的目標(biāo)位置信息（適應(yīng)值），在整個(gè)過(guò)程中，粒子共享它們最優(yōu)位置的信息，然后使用它們的記憶調(diào)整它們自己的速度和位置，不斷地比較和追隨候選的空間解，最終發(fā)現(xiàn)最優(yōu)解或局部最優(yōu)解。

在基本PSO算法中，粒子位置和速度分別表示為 X i=(x0i，x1i，…，xni)， V i=(v0i，v1i，…，vni)，粒子通過(guò)下列方程更替：

V i(t+1)= V i(t)+c1×rand1( P i- X i(t))+

c2×（rand2( P g- X i(t))(1)

X i(t+1)= X i(t)+ V i(t+1)(2)

其中： V i(t)是粒子i在t次迭代中的速度信息； X i(t)是粒子i在t次迭代中的位置信息;c1和c2是學(xué)習(xí)因子，分別調(diào)節(jié)向全局最優(yōu)粒子和個(gè)體最優(yōu)粒子方向飛行的最大步長(zhǎng);rand1和rand2是（0，1）中的隨機(jī)數(shù)。為了防止粒子飛離解空間，粒子的速度被限定在[-Vmax，+Vmax]中。矢量 P i=(Pi1，Pi2，…，PiD)是粒子i的最好的早先位置，稱為pbest（這個(gè)位置給出了最優(yōu)的適應(yīng)值）；矢量 P g=(Pg1，Pg2，…，PgD)是所有粒子中最優(yōu)粒子的位置，稱為gbest。

自從1995年P(guān)SO算法提出以后，在此基礎(chǔ)上提出了許多改進(jìn)算法。其中，慣性權(quán)重（inertia weight）法^[6]和壓縮因子法^[7]的改進(jìn)效果比較好。Jun Sun等人^[8]認(rèn)為，在上述基本PSO及其改進(jìn)算法中，由于pbest與隨機(jī)數(shù)一樣都是常量，都是線性的方法，生物群體中個(gè)體的思維模式要比方程（1）（2）描述的復(fù)雜得多。實(shí)際上，個(gè)體的思維是不確定的，在很大程度上就像粒子有量子行為，因此將量子理論應(yīng)用于粒子群算法，提出了量子行為粒子群算法（QPSO）。實(shí)驗(yàn)結(jié)果表明QPSO算法在幾個(gè)基準(zhǔn)函數(shù)的實(shí)驗(yàn)中要優(yōu)于基本PSO算法。

在QPSO算法中，粒子的速度和位置信息均歸結(jié)為一個(gè)參數(shù)，算法方程如下：

p=(φ1×pi+φ2×pg)/(φ1+φ2)（3）

mbest=∑ m i=1 pi/M=(∑ M i=1 pi1/M，∑ M i=1 pi2/M，…，∑ M i=1 piD/M(4)

x(t+1)=p±β×|mbest-x(t)|×ln(1/u)(5)

其中：φ1、φ2是在（0，1）中產(chǎn)生的隨機(jī)數(shù)；pi 是粒子i的pbest，pg是gbest；mbest稱為中值最優(yōu)位置;β是系數(shù)創(chuàng)造力；x(t)為粒子的相關(guān)位置信息;M是群體中所含粒子的數(shù)目;u是（0，1）中的隨機(jī)數(shù)。在迭代過(guò)程中，±是由（0，1）中隨機(jī)產(chǎn)生的隨機(jī)數(shù)的大小決定的，當(dāng)產(chǎn)生的隨機(jī)數(shù)大于05時(shí)，取“-”號(hào)；其他情況取“+”號(hào)。

量子粒子群算法在函數(shù)測(cè)試^[8]、濾波器設(shè)計(jì)^[9]、多階段金融規(guī)劃^[10]、神經(jīng)網(wǎng)絡(luò)優(yōu)化^[11]和H∞控制^[12]等應(yīng)用中顯示了其優(yōu)越性。

2 最大互信息法

互信息是香農(nóng)信息論中的一個(gè)基本概念，由信息論中的另外一個(gè)基礎(chǔ)理論“信息熵”引申而來(lái)。它是對(duì)兩個(gè)隨機(jī)事件（變量）統(tǒng)計(jì)相關(guān)性的測(cè)度，通常用于描述兩個(gè)系統(tǒng)間的統(tǒng)計(jì)相關(guān)性，或者是在一個(gè)系統(tǒng)中包含的另一個(gè)系統(tǒng)的信息的多少，可以用熵來(lái)描述。

基于互信息的醫(yī)學(xué)圖像配準(zhǔn)方法最早由Viola和Wells等人提出。由于該方法無(wú)須對(duì)不同成像模式中圖像灰度間的關(guān)系作任何假設(shè)，也無(wú)須對(duì)圖像作特征提取或其他預(yù)處理，可以避免特征提取造成的精度損失，被廣泛地用于CT/MR、PET/MR等多種配準(zhǔn)工作，認(rèn)為是目前最準(zhǔn)確和魯棒性最強(qiáng)的回溯性圖像配準(zhǔn)的度量之一^[13]。最大互信息法幾乎可以用在所有不同模式圖像的配準(zhǔn)，特別是當(dāng)其中一個(gè)圖像的數(shù)據(jù)部分缺損時(shí)也能得到很好的配準(zhǔn)效果。

在多模醫(yī)學(xué)圖像配準(zhǔn)中，雖然待配準(zhǔn)的兩幅圖像分別來(lái)自于不同的成像設(shè)備，但都基于共同的人體解剖信息，所以當(dāng)兩幅圖像的空間位置完全一致時(shí)，其中一幅圖像的信息表達(dá)的關(guān)于另一幅圖像的信息，也就是對(duì)應(yīng)像素灰度的互信息應(yīng)為最大^[14]。那么兩個(gè)系統(tǒng)A和B的互信息可以定義為



I(A，B)=∑ a，b pAB(a，b)log[pAB(a，b)/（pA(a)×pB(b)）](6) 

其中：pA(a)、pB(b)為邊緣概率分布，pAB(a，b)為聯(lián)合概率分布。由于聯(lián)合直方圖無(wú)須計(jì)算互信息的導(dǎo)數(shù)且易于實(shí)現(xiàn)，本文用它來(lái)統(tǒng)計(jì)聯(lián)合概率分布，然后再求得邊緣概率分布。

pAB(i， j)=h(i， j)/∑ i， j h(i， j)

pA(i)=∑ j pAB(i， j)， pB(j)=∑ i pAB(i， j)

將上式分別代入式(6)即可得到互信息的計(jì)算公式。上面的公式顯示：互信息就是聯(lián)合概率分布pAB(i， j)與各自的邊緣概率分布乘積pA(i)×pB(j)之間的相對(duì)熵。

最大互信息法算法流程如下：

a)采用灰度直方圖計(jì)算圖像A和B的邊緣概率分布PDF[i， j]；

b)計(jì)算圖像A和B的聯(lián)合熵H(A，B)=-∑ j， k PDF[j，k]log PDF[j，k]；

c)計(jì)算圖像A和B的邊緣熵H(A)=∑ j (∑ k PDF[j，k]log∑ l PDF[j，l]；

d) 計(jì)算圖像A和B的互信息I(A，B)=H(A)+H(B)-H(A|B)。 

最大互信息法具有如下優(yōu)點(diǎn)：

a）人工干預(yù)少，只依賴于圖像本身的信息，無(wú)須任何假設(shè)或先驗(yàn)醫(yī)學(xué)知識(shí)，也無(wú)須對(duì)圖像進(jìn)行特征點(diǎn)提取、組織分類等預(yù)處理，是一種自動(dòng)而有效的配準(zhǔn)算法。

b）精度高，可達(dá)到亞像素級(jí)。

c）可靠性高，對(duì)圖像中的幾何失真不敏感。

d）不依賴于成像設(shè)備，可應(yīng)用于多模醫(yī)療圖像配準(zhǔn)。

3 基于最大互信息和QPSO優(yōu)化策略的圖像配準(zhǔn)

為了比較QPSO、PSO以及Powell算法作為優(yōu)化策略時(shí)的配準(zhǔn)效果，將之分別應(yīng)用到相同兩幅圖片中。圖1(a)為MR圖，作為參考圖像；圖1(b)為CT圖，作為浮動(dòng)圖像。共設(shè)計(jì)六組實(shí)驗(yàn)，將CT事先移離原配準(zhǔn)位置10、15、20個(gè)像素點(diǎn)，以及旋轉(zhuǎn)10°、20°、25°的位置，用三種算法分別求配準(zhǔn)。在每種算法中，粒子群中粒子數(shù)均設(shè)為20，算法的主循環(huán)迭代次數(shù)均為100。每種算法運(yùn)行50次，記錄其平均最優(yōu)適應(yīng)值以作比較。在PSO算法中，本實(shí)驗(yàn)的參數(shù)設(shè)置如下：φ1＝φ2＝2.0，ω(t)=0.8(maxiter-t)/maxiter+0.6;在QPSO算法中，只需設(shè)定惟一的一個(gè)參數(shù)β=0.5(maxiter-t)/maxiter+0.5。

使用 QPSO算法進(jìn)行圖像配準(zhǔn)，算法步驟如下：

a)初始化，在由X、Y、Z軸平移度量和繞X、Y、Z軸的旋轉(zhuǎn)角度構(gòu)成的解空間中隨機(jī)選取一個(gè)點(diǎn)T，并在以T為中心的空間內(nèi)隨機(jī)初始化一個(gè)粒子群；

b)使用QPSO算法開(kāi)始迭代一定步數(shù)，得到當(dāng)前最好的一組解T′;

c) 利用互信息計(jì)算公式求得T、T′位置的互信息目標(biāo)函數(shù)值分別為I(xiàn)、I′，若I(xiàn)

d)將T′作為初始點(diǎn)繼續(xù)進(jìn)行迭代優(yōu)化，得到T附近的極值T″;

e)令T=T″作為新一輪搜索中的初始最優(yōu)點(diǎn)，重新初始化粒子群，轉(zhuǎn)b)；

f)輸出當(dāng)前最優(yōu)解T′及相應(yīng)的互信息函數(shù)值I′。

將CT圖移至離配準(zhǔn)位置10個(gè)像素點(diǎn)處后，分別用Powell、PSO和QPSO算法進(jìn)行配準(zhǔn)。每種算法均作50次，記錄100次迭代所求出的最優(yōu)值的平均值。由圖2可以看出，QPSO算法在迭代中所表現(xiàn)出來(lái)的收斂速度比另外兩種算法快30%左右；而且用QPSO在100次迭代內(nèi)就可以將兩張圖基本配準(zhǔn)，另外兩種算法就無(wú)法做到這一點(diǎn)。圖3顯示的是一次配準(zhǔn)后的結(jié)果，其I(xiàn)N=1161 7。為了表明配準(zhǔn)后的效果，一、三象限為原先的MR圖，二、四象限為配準(zhǔn)后的CT圖，各象限的接縫處過(guò)渡平滑，可以認(rèn)為基本配準(zhǔn)。每種算法運(yùn)行50次后得到的I(xiàn)N的最小值、最大值以及平均值列于表1中。由圖4可以清楚地看出，QPSO算法在100次循環(huán)迭代內(nèi)配準(zhǔn)達(dá)到的兩幅圖互信息值比另外兩種算法均有明顯提高，對(duì)應(yīng)的也就是配準(zhǔn)度有大幅提高。

由上述結(jié)果可以看出，在時(shí)間上Powell算法明顯快許多，但同時(shí)配準(zhǔn)結(jié)果誤差很大，由于最大互信息函數(shù)的多峰值使得該算法很快陷入局部極值而停止。用這三種優(yōu)化算法進(jìn)行的配準(zhǔn)過(guò)程，QPSO和PSO算法配準(zhǔn)所耗時(shí)間基本相同，前者稍快；PSO算法較之Powell算法配準(zhǔn)精度提高了不少，但是它搜索到的仍只是靠近全局最優(yōu)解附近的局部極值；QPSO比PSO算法配準(zhǔn)結(jié)果更優(yōu)，基本接近全局最優(yōu)解。究其原因，在標(biāo)準(zhǔn)PSO系統(tǒng)中，粒子的收斂是以軌道的形式實(shí)現(xiàn)的，且粒子的速度總是有限的，因此在搜索過(guò)程中粒子的搜索空間是一個(gè)有限的區(qū)域，不能覆蓋整個(gè)可行的空間。所以標(biāo)準(zhǔn)PSO算法不能保證以概率1收斂到全局最優(yōu)解，甚至不能收斂于局部最優(yōu)解。因此要找到最佳配準(zhǔn)位置的速度較慢，甚至找不到最佳配準(zhǔn)位置。

對(duì)比仿真表明，從量子力學(xué)角度出發(fā)提出的QPSO算法在配準(zhǔn)過(guò)程中具有較大的優(yōu)勢(shì)。這是因?yàn)镼PSO所基于的量子系統(tǒng)具有非線性和不確定性，所以在此系統(tǒng)中，一個(gè)粒子能夠以某一確定的概率出現(xiàn)在整個(gè)可行的搜索空間中任意一個(gè)位置，而非限定在某一個(gè)固定區(qū)域內(nèi)，因而增加了找到更好適應(yīng)值的可能。

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

本文分析了醫(yī)學(xué)圖像匹配的一些方法，利用QPSO算法全局收斂性，提出使用QPSO算法作為搜索優(yōu)化策略算法來(lái)解決 圖像匹配過(guò)程中計(jì)算復(fù)雜的問(wèn)題。仿真表明，QPSO算法作為優(yōu)化策略，較好地解決了圖像配準(zhǔn)中連續(xù)變量全局優(yōu)化的問(wèn)題，對(duì)初始點(diǎn)的選取要求較低、不易落入局部極值、實(shí)現(xiàn)簡(jiǎn)單，在配準(zhǔn)結(jié)果上優(yōu)于其他算法，具有很好的精確性和魯棒性，且算法具有較好的實(shí)時(shí)性，具有很高的實(shí)際應(yīng)用價(jià)值。

參考文獻(xiàn)：

[1]MAINTZ J B A， VIERGEVER M A. A survey of medical image registration[J].Medical Image Analysis ， 1998， 2 (1):136.

[2] JENKINSON M， SMITH S. A global optimization method for robust affined registration of brain images[J].Medical Image Analysis ， 2001， 5 (2):143156.

[3] BROWN L G. A survey of image registration techniques[J].Computing Surveys ， 1992， 24 (4):325376.

[4] MOFRTDIYXKI J. Numerical methods for image registration[M]. New York: Oxford University Press， 2004.

[5] KENNEDY J， EBERHART R. Particle swarm optimization[C]//Proc of IEEE International Conference on Neural Network. Perth: IEEE Press， 1995:19421948.

[6] SHI Yuhui， EBERHART R. A modified particle swarm optimizer[C]//Proc of IEEE International Conference of Evolutionary Computation. Piscataway: IEEE Press， 1998:6973.

[7] CLERC M. The swarm and the queen: towards a deterministic and adaptive particle swarm optimization[C]//Proc of Congress on Evolutionary Computation. Piscataway: IEEE Press， 1999:19511957.

[8] SUN Jun， FENG Bin， XU Wenbo. Particle swarm optimization with particles having quantum behavior[C]//Proc of Congress on Evolutionary Computation. 2004:325331.

[9] FANG Wei， SUN Jun， XU Wenbo. Design IIR digital filters using quantumbehaved particle swarm optimization[C]//Proc of the 2nd International Conference on Natural Computation. Berlin: Springer， 2006:992998.

[10] SUN Jun， XU Wenbo， FANG Wei. Solving multiperiod financial planning problem via quantumbehaved particle swarm algorithm[C]//Proc of Internationd Conference on Intelligent Computing. Berlin: Springer， 2006:11581169.

[11] SUN Jun， XU Wenbo， LIU Jing. Training RBF neural network via quantumbehaved particle swarm optimization[C]//Proc of the 13th International Conference on Neural Information Processing. Berlin: Springer， 2006:11561163.

[12] XI Maolong， SUN Jun， XU Wenbo. Quantumbehaved particle swarm optimization for designing Hinfinity structured specified controllers[C]//Proc of DCABES Conference. 2006:4246.

[13] PLUIM J P W， MAINTZ J B A， VIERGEVER M A. Mutualinformationbased registration of medical images: a survey[J].IEEE Trans on Medical Imaging ， 2003， 22 (8):9861004.

[14] VIOLA P， WELLS W M. Alignment by maximization of mutual information[J]. Journal of Computer Vision ， 1997， 24 (2):137154.