基于多目標(biāo)交叉變異粒子群算法的多模醫(yī)學(xué)圖像非剛性配準(zhǔn)

許鴻奎 江銘炎 楊明強(qiáng)（山東大學(xué)信息科學(xué)與工程學(xué)院，濟(jì)南 5000）（山東建筑大學(xué)，濟(jì)南 500）

引言

不同模態(tài)的醫(yī)學(xué)圖像提供不同的信息，如CT（computer-assisted tomography）和MRI（magnetic resonance image）圖像提供清晰的解剖結(jié)構(gòu)信息，而PET（positron emission tomography）和 fMRI（functional MRI）則提供各種功能信息。圖像融合可以充分利用多模態(tài)醫(yī)學(xué)圖像的各種相互補(bǔ)充的信息，在醫(yī)療診斷、治療方案的制定等方面發(fā)揮重要的作用。

圖像融合的前提是圖像配準(zhǔn)，而以互信息為相似性測(cè)度的圖像配準(zhǔn)不需要對(duì)圖像進(jìn)行分割或其他預(yù)處理，特別適合于多模態(tài)醫(yī)學(xué)圖像的配準(zhǔn)?；诨バ畔⒌呐錅?zhǔn)最早由Maes提出［1］，基于互信息的非剛性配準(zhǔn)最早由Meyer［2］采用薄板樣條變換實(shí)現(xiàn)［2］；文獻(xiàn)［3-5］則基于互信息采用B樣條變換完成了非剛性配準(zhǔn)，其中文獻(xiàn)［5］相對(duì)［3］提高了配準(zhǔn)速度，但是配準(zhǔn)精度有所下降。雖然基于互信息的配準(zhǔn)已經(jīng)得到了很大發(fā)展和應(yīng)用，但在提高配準(zhǔn)精度上還面臨許多挑戰(zhàn)，主要問題在于傳統(tǒng)的優(yōu)化方法如梯度下降法、擬牛頓法（limited-memory Broyden-Fletcher-Goldfarb-Shanno，LBFGS）等很容易陷入互信息的局部極值，不能很好地搜索到最佳配準(zhǔn)參數(shù)。對(duì)圖像進(jìn)行預(yù)濾波處理，改善互信息曲面，使其更為平滑，一定程度上可以提高配準(zhǔn)精度。然而實(shí)驗(yàn)證明，使用不同的濾波器得到的配準(zhǔn)結(jié)果差別很大，而不同的圖像對(duì)濾波器的要求也不一樣，這就限制了基于互信息的配準(zhǔn)方法的應(yīng)用。為此，文獻(xiàn)［6-7］分別使用遺傳算法、粒子群算法去全局搜索最大互信息以獲得最佳剛性變換參數(shù)，取得了良好的結(jié)果；文獻(xiàn)［8］使用海明窗濾波對(duì)圖像進(jìn)行預(yù)處理，再利用粒子群算法進(jìn)行全局搜索。然而，這些方法都是針對(duì)剛性配準(zhǔn)的。

基于互信息的非剛性配準(zhǔn)（如基于B樣條的配準(zhǔn)），特別適合于具有局部差異的多模態(tài)醫(yī)學(xué)圖像之間的配準(zhǔn)。相對(duì)剛性配準(zhǔn)，更多的配準(zhǔn)參數(shù)及其相互作用，使得互信息曲面更加復(fù)雜；各種因素產(chǎn)生的互信息局部極值對(duì)搜索策略提出了更高的要求，若采用傳統(tǒng)的以及上述優(yōu)化方法，難以得到理想的配準(zhǔn)結(jié)果。針對(duì)這一問題，本研究提出一種新的配準(zhǔn)參數(shù)的全局搜索方法，用于基于互信息的多模態(tài)醫(yī)學(xué)圖像的B樣條配準(zhǔn)。它以LBFGS的優(yōu)化結(jié)果為參考來構(gòu)造初始種群，初始種群更優(yōu)良；采用粒子群算法進(jìn)行多目標(biāo)優(yōu)化，結(jié)合交叉變異策略來增加種群的多樣性，克服粒子群算法受初始值的選取等因素的影響易陷入局部最優(yōu)的缺點(diǎn)。實(shí)驗(yàn)證明，不必對(duì)圖像進(jìn)行預(yù)平滑處理也能得到理想的配準(zhǔn)精度，提高了互信息配準(zhǔn)的魯棒性。

1 基于互信息的圖像配準(zhǔn)

盡管多模醫(yī)學(xué)圖像的灰度并不相似，甚至差別很大，但對(duì)同一人體組織而言，對(duì)應(yīng)像素點(diǎn)之間的灰度在統(tǒng)計(jì)學(xué)上并非獨(dú)立而是相關(guān)的；而互信息用于描述兩個(gè)系統(tǒng)間的統(tǒng)計(jì)相關(guān)性，或者一個(gè)系統(tǒng)中所包含的另一個(gè)系統(tǒng)中信息的多少，因而以互信息為測(cè)度進(jìn)行配準(zhǔn)是可行的［2］。

給定圖像A和B，它們的互信息量定義為：

其中，H（A）和H（B）分別為圖像A和B的熵，H（A，B）為二者的聯(lián)合熵?；バ畔⒌挠?jì)算可以通過兩幅圖像的聯(lián)合直方圖統(tǒng)計(jì)來獲得：

式中，pAB（a，b）為兩圖像的灰度對(duì)（a，b）在兩圖像相同坐標(biāo)位置的出現(xiàn)概率，pA（a）是圖像A中像素取灰度值a的概率，pB（b）是圖像B中像素取灰度值b的概率。

基于互信息的圖像配準(zhǔn)可以描述為：

式中，A為參考圖像，B為浮動(dòng)圖像，T表示對(duì)浮動(dòng)圖像施加的某種變換，采用三次B樣條變換。配準(zhǔn)過程就是對(duì)互信息函數(shù)的最優(yōu)求解過程。

1.1 基于B樣條的圖像變換

基于B樣條自由變形的圖像配準(zhǔn)是通過操控覆蓋在圖像上的控制網(wǎng)格來實(shí)現(xiàn)的。B樣條函數(shù)的局部特性使得一個(gè)控制點(diǎn)的移動(dòng)只會(huì)影響它附近的像素點(diǎn)，而不是整個(gè)圖像所有的像素點(diǎn)都發(fā)生移動(dòng)，因而特別適合于具有局部變形的圖像之間非剛性配準(zhǔn)。

定義域Ω為包含圖像所有像素的最小平面矩形域，用Φ表示由nx×ny個(gè)控制點(diǎn)Φij（0≤i≤nx，0≤j≤ny）所組成的網(wǎng)格，網(wǎng)格步長(zhǎng)為δ。將網(wǎng)格Φ覆蓋在平面矩形域Ω上，Ω上的一個(gè)點(diǎn)P（x，y）經(jīng)過三次B樣條變換后的位置用下式表示：

某點(diǎn)的位移由其周圍的16個(gè)網(wǎng)格控制點(diǎn)及其位移決定。0≤u＜1，表示網(wǎng)格控制點(diǎn)對(duì)T貢獻(xiàn)的權(quán)重，該權(quán)重基于該點(diǎn)到各網(wǎng)格頂點(diǎn)的距離。

基于互信息的B樣條配準(zhǔn)，就是找到一組最合適的網(wǎng)格控制點(diǎn)Φij對(duì)浮動(dòng)圖像進(jìn)行變換，使得變換后的浮動(dòng)圖像與參考圖像之間的互信息最大，而最優(yōu)變換參數(shù)由搜索策略完成。

1.2 傳統(tǒng)優(yōu)化算法對(duì)于互信息配準(zhǔn)尚存在的問題

傳統(tǒng)的優(yōu)化算法如牛頓法、梯度下降法具有理論完備、算法成熟、效率較高和結(jié)果穩(wěn)定性較好等諸多優(yōu)點(diǎn)，因而在實(shí)際應(yīng)用過程中獲得了廣泛應(yīng)用。但是它們對(duì)函數(shù)的連續(xù)性、導(dǎo)數(shù)的存在性都有相對(duì)較高的要求，采用梯度下降法、LBFGS法尋找上述最優(yōu)配準(zhǔn)參數(shù)，需要預(yù)先對(duì)圖像進(jìn)行高斯平滑處理，以減少互信息的局部極值點(diǎn)，否則無法收斂到最優(yōu)，配準(zhǔn)效果不理想。

互信息局部極值的產(chǎn)生原因是多方面的［9］，比如圖像本身固有的特性，圖像變換時(shí)插值的影響，圖像重疊區(qū)域大小的變化，噪聲的影響等等。抑制局部極值有多種方法，比如可以采用更高階的插值方法，或?qū)D像先進(jìn)行圖像分割等預(yù)處理，再進(jìn)行配準(zhǔn)。高階的插值意味著計(jì)算量的增大，采用分割預(yù)處理與圖像本身的特性有較大關(guān)系，不具有普遍適用性。對(duì)于基于B樣條的非剛性配準(zhǔn)，互信息局部極值的產(chǎn)生原因更加復(fù)雜，更多的配準(zhǔn)參數(shù)及其相互作用使得互信息函數(shù)更加復(fù)雜。互信息函數(shù)的非凸性、不規(guī)則性使得傳統(tǒng)優(yōu)化算法在尋優(yōu)過程中很容易陷入局部極值。實(shí)驗(yàn)證明，若不進(jìn)行高斯平滑預(yù)處理，梯度下降法和LBFGS找不到最優(yōu)配準(zhǔn)參數(shù)；而不同的圖像，不同階次的B樣條，甚至不同的B樣條網(wǎng)格尺寸對(duì)高斯平滑濾波器的參數(shù)要求也都各不相同。這種對(duì)圖像預(yù)處理的要求，影響了互信息配準(zhǔn)的魯棒性。為此，擬采用全局搜索的智能優(yōu)化方法尋找最優(yōu)配準(zhǔn)參數(shù)。

2 粒子群搜索策略

從20世紀(jì)60年代開始，以遺傳算法（genetic algorithm，GA）和粒子群優(yōu)化方法（particle swarm optimization，PSO）為代表的智能算法得到了很大發(fā)展，其主要特點(diǎn)是全局隨機(jī)搜索策略，不用考慮目標(biāo)函數(shù)本身是否連續(xù)或者可微。而粒子群優(yōu)化算法相對(duì)遺傳算法又具有本質(zhì)的并行性，具有更高的搜索效率。PSO算法起源于對(duì)一個(gè)簡(jiǎn)化社會(huì)模型的仿真，它和鳥類或魚類的群集現(xiàn)象有明顯的聯(lián)系。作為種群進(jìn)化方法，每個(gè)優(yōu)化問題的解被認(rèn)為是搜索空間的一個(gè)飛行的粒子，所有的粒子都有一個(gè)由目標(biāo)函數(shù)決定的適應(yīng)度值，每個(gè)粒子由一個(gè)速度決定其飛行方向和距離。PSO算法初始化為一群隨機(jī)粒子，粒子們根據(jù)對(duì)個(gè)體和群體的飛行經(jīng)驗(yàn)的綜合分析來動(dòng)態(tài)調(diào)整各自的飛行速度，在解空間中通過迭代尋找最優(yōu)解。每一次迭代中，粒子通過跟蹤兩個(gè)“極值”來更新自己，一個(gè)是個(gè)體極值pbest，即粒子目前自身所找到的最優(yōu)解；另一個(gè)是全局極值gbest，即整個(gè)種群目前找到的最優(yōu)解。在基本粒子群算法中，速度和位置的更新按下面公式進(jìn)行［10］：

式中，下標(biāo)i代表第i個(gè)粒子，下標(biāo)j代表速度或位置的第j維，上標(biāo)k代表迭代次數(shù)；c1和c2是學(xué)習(xí)因子，通常c1、c2∈［0，4］；r1和r2是介于［0，1］之間的隨機(jī)數(shù)；是粒子Pi在第j維的個(gè)體極值的坐標(biāo)；是群體在第j維的全局極值的坐標(biāo)。

粒子群算法這種全局搜索策略，不必考慮目標(biāo)函數(shù)本身是否連續(xù)或者可微，因而用于基于互信息的多模醫(yī)學(xué)圖像非剛性配準(zhǔn)時(shí)，就沒有對(duì)圖像進(jìn)行平滑預(yù)處理的特別要求。實(shí)驗(yàn)表明，高斯平滑濾波對(duì)其收斂性影響很小，因而采用粒子群優(yōu)化算法來搜索B樣條配準(zhǔn)參數(shù)。但是實(shí)驗(yàn)中也發(fā)現(xiàn)，即使采用帶慣性權(quán)重因子的粒子群算法［11］，受初始值隨機(jī)選取的影響，結(jié)果也不穩(wěn)定，出現(xiàn)了早熟現(xiàn)象，配準(zhǔn)結(jié)果不甚理想，為此本研究提出多目標(biāo)交叉變異粒子群優(yōu)化算法加以改進(jìn)。

3 多目標(biāo)交叉變異粒子群優(yōu)化算法

隨著基本粒子群算法的不斷應(yīng)用，出現(xiàn)了很多改進(jìn)方法。除了上述引入慣性權(quán)重和壓縮因子外，典型的還有自適應(yīng)模糊粒子群算法、雜交微粒群優(yōu)化算法、免疫粒子群算法、小生境粒子群算法［12］等等，在不同的應(yīng)用領(lǐng)域取得了良好的效果。

針對(duì)基于互信息的多模態(tài)醫(yī)學(xué)圖像的非剛性配準(zhǔn)，采用粒子群算法進(jìn)行多目標(biāo)優(yōu)化［13］，結(jié)合交叉變異策略［14］來跳出局部極值，克服粒子群算法早熟的缺點(diǎn)；同時(shí)對(duì)種群的初始化加以改進(jìn)，以LBFGS的優(yōu)化結(jié)果作為參考值構(gòu)造初始種群，這樣構(gòu)造的初始種群離目標(biāo)函數(shù)的全局最優(yōu)更近，陷入局部最優(yōu)的機(jī)會(huì)大大減少。雖然這樣會(huì)消耗一定的時(shí)間，但可以進(jìn)一步提高配準(zhǔn)精度。

3.1 多目標(biāo)粒子群優(yōu)化算法

多目標(biāo)優(yōu)化最早由意大利經(jīng)濟(jì)學(xué)家V.F.Pareto提出，通常一個(gè)多目標(biāo)優(yōu)化問題可以表示如下［13］：

式中，S?Rn稱為可行解區(qū)域，x為決策變量，E=｛F（x）|x∈Rn｝為目標(biāo)向量，fi（x）i=1，2，…n是第i個(gè)目標(biāo)函數(shù)。在大多數(shù)情況下，各個(gè)目標(biāo)函數(shù)可能是有沖突的，這就使得多目標(biāo)優(yōu)化問題不存在惟一的全局最優(yōu)解，使所有目標(biāo)函數(shù)同時(shí)最優(yōu)。

由于誤差的平方和（sum of square differences，SSD）也是配準(zhǔn)圖像的常用的相似性測(cè)度，采用互信息為第一目標(biāo)函數(shù)f1（x），而SSD為第二目標(biāo)函數(shù)f2（x）。采用下述交叉變異的方法，同時(shí)優(yōu)化兩個(gè)目標(biāo)函數(shù)，以期增加種群的多樣性。

3.2 子群間的交叉變異

在兩個(gè)子群中按一定的雜交概率選取指定數(shù)量的粒子放入繁殖池中隨機(jī)地兩兩進(jìn)行雜交，產(chǎn)生相同數(shù)目的子代粒子。子代粒子的位置按式（9）和式（10）根據(jù)父代粒子的位置計(jì)算得到，子代粒子的速度按式（11）和（12）根據(jù)父代粒子的速度計(jì)算得到［15］。雜交后的子代粒子回歸到各自的子群并代替父代粒子，從而保持種群的規(guī)模不變。

式中，x是粒子的位置矢量，v是粒子的速度矢量；Ci（x）和Pi（x）i=1，2是子代粒子和父代粒子的位置；Ci（v）和Pi（v）i=1，2是子代粒子和父代粒子的速度；p是雜交概率。在迭代過程中，后代不斷地繼承來自兩個(gè)子群的父代的特征，然后又尋找使自己的適應(yīng)度函數(shù)最優(yōu)的位置。這種交叉變異機(jī)制，保證了種群的多樣化，使得算法在解空間搜索的遍歷性得到改善，可以有效地避免求解陷入局部最優(yōu)。

3.3 算法流程

步驟1：子群劃分：根據(jù)目標(biāo)函數(shù)f1（x）和f2（x）把微粒群平均分成2個(gè)子群ChiledS1和ChiledS2。

步驟2：設(shè)置子群參數(shù)：設(shè)置粒子的雜交概率p，設(shè)置其它參數(shù)學(xué)習(xí)因子c1、c2，慣性參量ω、速度限制vmax等。

步驟3：以LBFGS的優(yōu)化結(jié)果為參考，初始化兩個(gè)子群ChiledS1和ChiledS2中各粒子的初始位置。

步驟4：按照慣性權(quán)重粒子群算法更新子群中粒子的速度和位置。

步驟5：子群間的交叉變異：從兩個(gè)子群中按參數(shù)選擇一定的粒子進(jìn)行雜交，根據(jù)式（9）和式（10）以及式（11）和式（12）分別產(chǎn)生子代的速度和位置。

步驟6：滿足迭代中止的條件，則停止迭代，否則返回步驟4繼續(xù)迭代。

由兩個(gè)目標(biāo)函數(shù)的搜索結(jié)果取均值作為最終配準(zhǔn)參數(shù)，然后對(duì)浮動(dòng)圖像進(jìn)行B樣條變換，完成和參考圖像的配準(zhǔn)。

4 實(shí)驗(yàn)結(jié)果與分析

腦部實(shí)驗(yàn)圖像來自BrainWeb［16］。圖1（a）是217像素×181像素的MR-T2圖像，作為配準(zhǔn)的參考圖像；圖1（b）是相對(duì)應(yīng)的MR-PD圖像，在左上方和右下方經(jīng)過了局部變形，作為待配準(zhǔn)的浮動(dòng)圖像。此變形是對(duì)兩個(gè)網(wǎng)格控制點(diǎn)進(jìn)行一定的移位產(chǎn)生的。

首先以互信息為相似性測(cè)度，采用LBFGS優(yōu)化方法進(jìn)行三次B樣條變換實(shí)現(xiàn)配準(zhǔn)，B樣條變換的網(wǎng)格大小為32像素×32像素。正如前面1.2節(jié)分析，由于互信息曲面存在局部極值，配準(zhǔn)前需要預(yù)濾波平滑處理。圖1（c）是使用模板大小為20像素×20像素和尺度為5的高斯濾波器得到配準(zhǔn)圖像，表1列出了使用不同參數(shù)的高斯平滑濾波器得到的配準(zhǔn)結(jié)果?？梢?，使用不同的濾波器配準(zhǔn)的結(jié)果差別很大。LBFGS優(yōu)化方法對(duì)圖像預(yù)處理的依賴，影響了基于互信息配準(zhǔn)方法的魯棒性。

粒子群算法是一種全局優(yōu)化方法，不需要圖像預(yù)處理，也能得到較高的配準(zhǔn)精度。圖2是采用慣性權(quán)重粒子群優(yōu)化方法的配準(zhǔn)結(jié)果，配準(zhǔn)前沒有進(jìn)行高斯濾波；其中，圖2（a）是圖1（a）和圖1（b）配準(zhǔn)前的差，圖2（b）是圖1（a）和圖1（b）配準(zhǔn)后的差。對(duì)20次配準(zhǔn)的結(jié)果進(jìn)行平均計(jì)算，互信息和誤差的平均值分別是2.049和0.045，這是不錯(cuò)的配準(zhǔn)結(jié)果。然而從圖2（c）所示控制網(wǎng)格的變化可以看出仍有早熟的跡象，所得配準(zhǔn)參數(shù)距離真正的網(wǎng)格點(diǎn)變化值仍有一定的差距。

圖1 待配準(zhǔn)圖像及采用LBFGS優(yōu)化的配準(zhǔn)結(jié)果。（a）MR-T2參考圖像；（b）MR-PD浮動(dòng)圖像；（c）配準(zhǔn)后的浮動(dòng)圖像Fig.1 Images to be made registration and the result using LBFGS optimization.（a）the MR-T2 reference image；（b）the MR-PD float image；（c）the float image after registration

表1 不同參數(shù)的高斯平滑濾波器的配準(zhǔn)結(jié)果Tab.1 Registration results using different filters

圖2 采用慣性權(quán)重粒子群優(yōu)化算法的配準(zhǔn)結(jié)果。（a）配準(zhǔn)前的差；（b）配準(zhǔn)后的差；（c）控制網(wǎng)格的變化Fig.2 Registration result using inertia weight PSO.（a）the difference before registration；（b）the difference after registration；（c）the changes of control grid

圖3 采用多目標(biāo)交叉變異粒子群算法的配準(zhǔn)結(jié)果。（a）配準(zhǔn)后的浮動(dòng)圖像；（b）配準(zhǔn)后的差；（c）控制網(wǎng)格的變化Fig.3 Registration result using proposed optimization in the paper.（a）the float image after registration；（b）the difference after registration；（c）the changes of control grid

圖3是采用所提出的多目標(biāo)交叉變異粒子群優(yōu)化算法得到的配準(zhǔn)結(jié)果；其中，圖3（a）是配準(zhǔn)后的輸出圖像，圖3（b）是兩圖像配準(zhǔn)后的差。對(duì)比圖3（b）和圖2（b），配準(zhǔn)效果比慣性權(quán)重粒子群算法得到了明顯改善。對(duì)20次配準(zhǔn)后的結(jié)果進(jìn)行平均計(jì)算，互信息和誤差的平均值分別是2.075和0.019，相比慣性權(quán)重粒子群算法的2.049和0.045，提高了配準(zhǔn)精度。圖3（c）是控制網(wǎng)格的變化，對(duì)比圖2（c）可以看出，所得配準(zhǔn)參數(shù)更接近真正的網(wǎng)格點(diǎn)變化值。

將分別采用梯度下降法［3，5］、慣性權(quán)重粒子群優(yōu)化方法以及本研究?jī)?yōu)化方法得到的配準(zhǔn)結(jié)果數(shù)據(jù)列于表2加以比較，本方法的配準(zhǔn)精度優(yōu)于另外兩種方法。

除了上述腦部MR-T2圖像和MR-PD圖像，從BrainWeb中又取了50個(gè)MR-T2圖像及序列號(hào)對(duì)應(yīng)的50個(gè)MR-PD圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)。按前述相同的方法，人為地對(duì)MR-PD圖像進(jìn)行局部變形。假設(shè)將配準(zhǔn)誤差大于和等于0.02定義為誤配準(zhǔn)，依此統(tǒng)計(jì)各種方法的配準(zhǔn)率，如表3所示?？梢?，本方法取得了較高的配準(zhǔn)率。

為了進(jìn)一步驗(yàn)證本算法的有效性，還采用其它模態(tài)的醫(yī)學(xué)圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)。圖4是CT和CBCT（Cone beam computer tomography）圖像的配準(zhǔn)實(shí)驗(yàn)結(jié)果，其中圖4（a）參考圖像，圖4（b）是浮動(dòng)圖像，圖4（c）是采用改進(jìn)的多目標(biāo)交叉變異粒子群優(yōu)化算法的配準(zhǔn)輸出，圖4（g）是配準(zhǔn)輸出與參考圖像的差；作為比較，給出其它方法配準(zhǔn)后兩圖像的差，圖4（d）是剛性配準(zhǔn)后浮動(dòng)圖像與參考圖像的差，圖4（e）是采用LBFGS優(yōu)化方法進(jìn)行非剛性配準(zhǔn)后的浮動(dòng)圖像與參考圖像的差，圖4（f）是采用慣性權(quán)重粒子群算法進(jìn)行非剛性配準(zhǔn)后的浮動(dòng)圖像與參考圖像的差。很明顯，采用本方法配準(zhǔn)后的浮動(dòng)圖像與參考圖像的差最小，配準(zhǔn)效果最好。配準(zhǔn)前，兩圖像的互信息及距離（像素值差的平法和，以下相同）分別是0.527和11 583，采用本方法配準(zhǔn)后互信息及距離分別是1.082和7 155，配準(zhǔn)效果顯著。

表2 不同優(yōu)化方法的配準(zhǔn)結(jié)果比較Tab.2 Comparison of registration results using different optimizations

表3 不同優(yōu)化方法的配準(zhǔn)率比較Tab.3 Comparison of registration rate using different optimizations

圖4 對(duì)CT和CBCT圖像采用不同配準(zhǔn)方法進(jìn)行配準(zhǔn)的結(jié)果比較。（a）CT參考圖像；（b）CBCT浮動(dòng)圖像；（c）使用本文方法配準(zhǔn)后的浮動(dòng)圖像；（d）剛性配準(zhǔn)后的差；（e）LBFGS B樣條配準(zhǔn)后的差；（f）慣性權(quán)重PSO B樣條配準(zhǔn)后的差；（g）本方法配準(zhǔn)后的差Fig.4 Comparison of CT and CBCT image registration results using different registration methods.（a）the CT reference image；（b）the CBCT float image；（c）the float image after registration using the method proposed in this paper；（d）the difference after rigid registration；（e）the difference after LBFGS B-spline registration；（f）the difference after B-spline registration with inertia weight PSO；（g）the difference after registration proposed in this paper

圖5是CT和PET圖像的配準(zhǔn)實(shí)驗(yàn)結(jié)果。其中，圖5（a）參考圖像，圖5（b）是浮動(dòng)圖像，圖5（c）是采用改進(jìn)的多目標(biāo)交叉變異粒子群優(yōu)化算法的配準(zhǔn)輸出。配準(zhǔn)前，兩圖像的互信息及距離分別是0.884和48 254，采用本文方法配準(zhǔn)后互信息及距離分別是1.004和47 527，也取得了不錯(cuò)的配準(zhǔn)效果。

圖5 CT和PET圖像的配準(zhǔn)。（a）CT參考圖像；（b）PET浮動(dòng)圖像；（c）使用本方法配準(zhǔn)后的浮動(dòng)圖像Fig.5 Registration results of CT and PET image using the method proposed in this paper.（a）the CT reference image；（b）the PET float image；（c）the float image after registration

然而，配準(zhǔn)精度提高的代價(jià)是配準(zhǔn)速度的下降。采用圖1中的圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)時(shí)，本方法相比慣性權(quán)重粒子群算法，配準(zhǔn)后兩圖像的互信息改善了0.026，而配準(zhǔn)時(shí)間增加了近10％，大約19s。這是因?yàn)槎嗄繕?biāo)操作及交叉變異都需要額外占用一定的時(shí)間。

5 結(jié)論

基于多目標(biāo)交叉變異粒子群算法，以互信息為相似性測(cè)度采用B樣條變換對(duì)多模態(tài)醫(yī)學(xué)圖像進(jìn)行非剛性配準(zhǔn)，克服了傳統(tǒng)優(yōu)化方法如梯度下降法、LBFGS依賴于圖像預(yù)平滑處理的缺點(diǎn)，提高了基于互信息的非剛性配準(zhǔn)的魯棒性；該方法以LBFGS的優(yōu)化結(jié)果為參考構(gòu)造的初始種群更優(yōu)良，充分繼承了傳統(tǒng)優(yōu)化方法和智能優(yōu)化方法各自的優(yōu)點(diǎn)，所構(gòu)造的初始種群距離全局最優(yōu)更近；針對(duì)基本粒子群算法易陷入局部最優(yōu)的缺點(diǎn)，采用多目標(biāo)優(yōu)化方法引入交叉變異機(jī)制對(duì)粒子群算法加以改進(jìn)，使得算法在解空間搜索的遍歷性得到改善，優(yōu)化結(jié)果更接近全局最優(yōu)，從而減少了配準(zhǔn)誤差，提高了配準(zhǔn)的精度。如何在減少配準(zhǔn)誤差的同時(shí)進(jìn)一步提高配準(zhǔn)的速度將是以后研究的重點(diǎn)。

［1］Maes F.Multimodality image registration by maximization of mutual information［J］.IEEE Trans on Medical Imaging，1997，16（2）：187-198.

［2］Meyer CR，Leichtman GS，Brunberg JA，et al.Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin plane spline warped geometric deformations［J］.Medical Image Analysis，1997，1（3）：195-206.

［3］Rueckert D，Sonoda L，Hayes C，et al.Nonrigid registration using free-form deformations-application to breast MR images［J］.IEEE Trans on Medical Imaging，1999，18（8）：712-721.

［4］Kybic J.Fast parametric elastic image registration［J］.IEEE transaction on image processing，2003，12（11）：1427-1441.

［5］彭曉明，陳武凡，馬茜.基于B樣條的快速彈性圖像配準(zhǔn)方法［J］.計(jì)算機(jī)工程與應(yīng)用，2006，11：186-189.

［6］Ingole VT，Deshmukh CN，Anjali Joshi，et al.Medical image registration using genetic algorithm［C］∥Ajith A，Norway.2009 2nd International Conference on Emerging Trends in Engineering and Technology.USA：IEEE computer Society，2009：63-66.

［7］Chen Yenwei，Aya Mimori.Hybrid particle swarm optimization for medical image registration［C］∥Hua Zhang.2009 Fifth International Conference on Natural Computation.Tianjin：Institute of Electrical and Electronics Engineers（IEEE），2009：26-30.

［8］裴繼紅，田劍豪，楊烜.基于海明窗濾波及粒子群優(yōu)化搜索的醫(yī)學(xué)圖像配準(zhǔn)［J］.生物醫(yī)學(xué)工程學(xué)雜志，2007，24（02）：262-267.

［9］Ardizzone E，Gallea R.Effective and efficient interpolation for mutual information based multimodality elastic image registration［C］∥Takashi M.2009 IEEE 12th International Conference on Computer Vision Workshops.Washington：IEEE Computer Society，2009：376-381.

［10］Kennedy J，Eberhart R.Particle swarm optimization［C］∥Proceedings of IEEE International Conference on Neural Networks.New Jersey：IEEE Service Center，1995：1942-1948.

［11］Ememipour J，Nejad MS，Ebadzadeh MM，et.al.Introduce a new inertia weight for particle swarm optimization［C］∥Sungwon S.2009 Fourth International Conference on Computer Sciences and Convergence Information Technology.Seoul：IEEE Computer Society，2009：1650-1653.

［12］王萬良，唐宇.微粒群算法的研究現(xiàn)狀與展望［J］.浙江工業(yè)大學(xué)學(xué)報(bào)，2007，35（2）：136-141.

［13］Xiang Feng，F(xiàn)rancic CM.A parallel evolutionary approach to multi-objective optimization［C］∥Kay CT.2007 IEEE Congress on Evolutionary Computation.Singapore：IEEE Computational Intelligence Society，2007：1199-1206.

［14］Coelho Ls，Krohling Ra.Predictive controller tuning using modified particle swarm optimization based on cauchy and gaussian distributions［C］∥Mario K.8th On-Line World Conference on Soft Computing in Industrial Applications.Berlin：Springer，2005：287-298.

［15］張敏慧.改進(jìn)的粒子群計(jì)算智能算法及其多目標(biāo)優(yōu)化的應(yīng)用研究［D］.杭州：浙江大學(xué)，2005.

［16］EVANS AC.Brain Web：online simulated brain database［DB/OL］.http://www.bic.mni.mcgill.ca/brainweb，2006-06-08/2007-03-01.