曹國剛，朱信玉，陳 穎，曹 聰，孔德卿

（上海應(yīng)用技術(shù)大學(xué)計算機科學(xué)與信息工程學(xué)院，上海，201418）

引 言

由于單一模態(tài)的醫(yī)學(xué)圖像難以提供充足的病灶信息，通常將多模態(tài)的醫(yī)學(xué)圖像進行融合，從而輔助醫(yī)生診斷。在進行圖像融合前，需要對患者的醫(yī)學(xué)圖像進行配準(zhǔn)，即在參考圖像與浮動圖像之間通過尋優(yōu)算法找到最優(yōu)的空間變換參數(shù)，使兩幅圖像在誤差最小的情況下對齊，其本質(zhì)是參數(shù)優(yōu)化問題[1-3]。

經(jīng)典醫(yī)學(xué)圖像配準(zhǔn)框架主要包括測度函數(shù)、優(yōu)化算法、空間變換、插值4 個方面。其中，相似性測度是圖像配準(zhǔn)結(jié)果的衡量指標(biāo)，測度函數(shù)值用來表示兩幅圖像的對齊程度[4]。常見的測度函數(shù)包括互信息(Mutual informational, MI)、條件方差和(Sum of conditional variance, SCV)[5]、歸一化互信息(Normalized MI, NMI)[6]、交叉累計剩余熵(Cross cumulative residual entropy, CCRE)[7-8]、歸一化互相關(guān)(Normalization cross-correlation, NCC)[9]等。尋優(yōu)算法在圖像配準(zhǔn)的過程中對配準(zhǔn)精度與速度起到?jīng)Q定性作用。當(dāng)前優(yōu)化算法主要分為局部優(yōu)化算法與全局優(yōu)化算法兩種，常用的局部優(yōu)化算法包括梯度下降法、Powell 算法、單純形搜索法，全局優(yōu)化算法包括遺傳算法、差分進化算法[10-11]、粒子群算法[12-13]等。測度函數(shù)普遍存在多極值問題，因而僅使用局部搜索算法難以逃離局部最優(yōu)，雖然全局搜索算法可以跳出局部最優(yōu)，但此類算法存在計算量大、收斂速度慢等問題。在醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域，采用多分辨率策略可以充分結(jié)合局部優(yōu)化算法和全局優(yōu)化算法的優(yōu)點，有效地縮短運算時間，同時保證較高的配準(zhǔn)精度。如文獻[14]將粒子群優(yōu)化算法（Particle swarm optimization, PSO）與單純形搜索法（Simplex）相結(jié)合使用多分辨率策略對醫(yī)學(xué)圖像配準(zhǔn)，其精度優(yōu)于Powell、Simplex、PSO 三種算法。文獻[15]將自適應(yīng)差分算法與Powell 相結(jié)合應(yīng)用于醫(yī)學(xué)圖像配準(zhǔn)，配準(zhǔn)精度達到了亞像素級。文獻[16]將頭腦風(fēng)暴優(yōu)化算法與Powell 相結(jié)合，與遺傳算法、粒子群算法、蟻群算法與Powell 結(jié)合算法進行比較，配準(zhǔn)結(jié)果均方根誤差與配準(zhǔn)時間均得到了一定程度的降低。

為了進一步提高配準(zhǔn)算法的精度與收斂速度，本文基于多分辨率策略開展了以下3 方面的工作：（1）改進了頭腦風(fēng)暴優(yōu)化算法，提出一種改進頭腦風(fēng)暴優(yōu)化（Improved brain storm optimization ，IBSO）算法；（2）提出一種基于多分辨率策略的醫(yī)學(xué)圖像配準(zhǔn)算法，該算法在低分辨率層使用IBSO 粗配準(zhǔn)，高分辨使用單純形搜索法精配準(zhǔn)；（3）對提出的配準(zhǔn)算法分別完成了MRI-MRI單模態(tài)與CT-MRI多模態(tài)實驗。

1 醫(yī)學(xué)圖像配準(zhǔn)方法

1.1 相似性測度

MI 源自于信息論，用于表示浮動圖像F包含參考圖像R的信息量。因互信息作為測度函數(shù)無需特征提取且配準(zhǔn)精度高，因而本文在單模態(tài)配準(zhǔn)實驗中采用其作為測度函數(shù)評價配準(zhǔn)效果。R與F的互信息定義為

式中：H(R)、H(F)分別表示參考圖像R、浮動圖像F包含的信息量，如式（2）和（3）所示；H(R,F)為R和F的聯(lián)合熵，其定義如式（4）所示。

多模態(tài)圖像配準(zhǔn)時，采用互信息作為測度函數(shù)容易受到參考圖像R與浮動圖像F重疊程度的影響，導(dǎo)致測度函數(shù)曲線不光滑。Studholme 等[6]提出NMI 可以解決上述問題，因而在多模圖像配準(zhǔn)中使用廣泛，其定義為

NCC 是另一種常用的測度函數(shù)，其取值范圍為[0,1]，NCC 系數(shù)越趨近于1 代表配準(zhǔn)精度越高，其定義為

Wang 等[8]基于累計剩余熵理論提出一種CCRE 測度,其通過實驗證明了CCRE 具有較高的魯棒性，CCRE 定義為

式中：ε(R)為累計剩余熵函數(shù)，表達式為

式中：P(R＞u)代表圖像R的像素點比灰度值u高的概率。

1.2 多分辨率策略

局部搜索算法容易陷入局部最優(yōu)，而全局優(yōu)化算法運算時間長。為了平衡配準(zhǔn)的效率與精度，在醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域常常對待配準(zhǔn)圖像進行兩層小波分解[17]，由原圖像及小波分解后的圖像構(gòu)成圖像金字塔，其模型如圖1 所示。

小波變換將時域信號變換成頻域信號，小波分解后圖像的低頻成分保留原始圖像的基本信息，高分辨率圖像保留圖像的細節(jié)部分以及夾雜了高頻噪聲。因此在低分辨率的圖像上使用全局尋優(yōu)算法進行粗配準(zhǔn)，得到尋優(yōu)參數(shù)作為局部尋優(yōu)算法的初始點，再對高分辨率圖像進行精配準(zhǔn)是醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域一種常用的配準(zhǔn)策略。

圖1 圖像金字塔模型Fig.1 Model of image pyramid

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

基于人類頭腦風(fēng)暴過程建模，Shi 于2011 年提出了頭腦風(fēng)暴優(yōu)化算法(Brain storm optimization algorithm, BSO)[18]，并于2016 年又進行了綜述[19]。BSO 算法將種群中每個想法個體看作是d維問題的一個潛在解，種群中第i個想法個體可表示為

頭腦風(fēng)暴優(yōu)化算法首先根據(jù)問題的規(guī)模產(chǎn)生想法種群并初始化，然后在每次迭代中對想法種群進行聚類、變異、產(chǎn)生新個體和選擇操作。

（1）聚類。將想法種群聚成m個子群，在每個子群中根據(jù)個體的適應(yīng)度值進行排序，選擇最優(yōu)適應(yīng)度值的個體作為每個子群的聚類中心。

（2）變異。生成隨機數(shù)p1，其取值范圍為[0,1]，如果p1小于預(yù)設(shè)概率ppre1則生成一個新個體取代隨機選擇的聚類中心。

（3）產(chǎn)生新個體。生成隨機數(shù)p2，其取值范圍為[0,1]，將其與預(yù)設(shè)概率ppre2進行比較：

① 如果p2≤ppre2，隨機選擇一個子群同時生成[0,1]之間隨機數(shù)p3，與預(yù)設(shè)概率ppre3比較，（a）如果p3≤ppre3，選擇聚類中心作為線索用于產(chǎn)生新個體；（b）如果p3＞ppre3，隨機選擇該子群中的一個想法個體作為線索用于產(chǎn)生新個體。

② 如果p2＞ppre2，隨機選擇兩個子群同時生成[0,1]之間隨機數(shù)p4，與預(yù)設(shè)概率ppre4比較，（a）如果p4≤ppre4，選擇兩聚類中心作為線索用于產(chǎn)生新個體；（b）如果p4＞ppre4，隨機在兩個子群中選擇兩個個體作為線索用于產(chǎn)生新個體。

（4）選擇。將產(chǎn)生的新個體與被選擇個體進行適應(yīng)度值比較，選擇其中適應(yīng)度值高的個體進入新一輪迭代。

（5）若達到最大迭代次數(shù)或滿足最優(yōu)解條件則停止迭代，輸出結(jié)果，否則轉(zhuǎn)步驟1，開始新一輪迭代。

為了保證新產(chǎn)生的想法盡可能地利用現(xiàn)有的想法，通常采用在線索個體中添加噪聲的方式來產(chǎn)生新想法個體。通過已存在的想法xold產(chǎn)生一個新想法xnew，可以表示為

式中：T為最大迭代次數(shù)；t為當(dāng)前迭代數(shù)；k用于改變 logsig(?)函數(shù)斜率，文獻[14]證明參數(shù)k為 25 是個不錯的選擇。

1.4 單純形搜索法

單純形搜索法是一種高效的局部搜索算法，常用于求解無約束最優(yōu)化問題。單純形指的是n維空間Rn中具有n+1 個頂點的凸多面體。單純形搜索法的基本步驟如下：給定Rn中的一個單純形后，分別求出n+1 個頂點處的函數(shù)值，確定出其中函數(shù)最大值點及函數(shù)最小值點，然后通過反射、擴展、壓縮等幾種方式求出一個較好的點來替代目前函數(shù)最大值點組成新的單純形，或者向函數(shù)最小點方向收縮構(gòu)成新的單純形。通過多次的迭代，逐漸收縮單純形范圍，逼近函數(shù)最小值點。

單純形搜索法在局部搜索時具有搜索效率高、速度快等優(yōu)勢。但初始點對其十分重要，當(dāng)初始點距離局部最優(yōu)點較近時，算法難以逃離局部最優(yōu)。為了避免單純形搜索法的這一缺陷，本文算法將單純形搜索法初始點設(shè)置為全局尋優(yōu)算法求得的全局最優(yōu)的粗略位置。

2 IBSO 算法及其在醫(yī)學(xué)圖像配準(zhǔn)中的應(yīng)用

2.1 IBSO 算法

BSO 在每一次迭代中通過聚類、變異操作使種群個體向最優(yōu)解收斂。在聚類過程中，每個子群中適應(yīng)度值最優(yōu)的個體被選擇作為聚類中心，其在算法的迭代過程中具有優(yōu)勢地位。BSO 通過向線索個體添加噪聲的方式生成新個體，而聚類中心相對普通個體具有較高的概率被選擇作為線索個體。

觀察算法的迭代過程可以發(fā)現(xiàn)：迭代后期階段每個子群的聚類中心基本穩(wěn)定，最優(yōu)的聚類中心個體是當(dāng)前迭代的最優(yōu)解，而最差聚類中心在搜索的后期階段基本處于停滯更新狀態(tài)。

為了保證所有個體均參與當(dāng)前的優(yōu)化搜索過程，加快收斂速度，本文提出IBSO 算法。IBSO 在每次迭代聚類過程結(jié)束后，對當(dāng)前最差的聚類中心進行改進。選擇當(dāng)前最優(yōu)聚類中心Cen(best)作為線索個體，根據(jù)式（10）生成新個體，計算新生成個體的適應(yīng)度值，與當(dāng)前最差聚類中心個體的適應(yīng)度值進行比較，選擇適應(yīng)度值高的個體進入下一次迭代過程，選擇條件如式（14）所示，IBSO 算法流程如圖2 所示。

式中：Cen(worst)為適應(yīng)度值最差的聚類中心；Cen(gen)為將當(dāng)前最優(yōu)聚類中心作為線索生成的新個體；Fit()為適應(yīng)度函數(shù)，F(xiàn)it(Cen(worst))為原最差聚類中心適應(yīng)度值，F(xiàn)it(Cen(gen))為新產(chǎn)生個體的適應(yīng)度值。

使用最優(yōu)聚類中心作為線索生成的新個體是極具潛力的個體，將其替換最差聚類中心個體，一方面，可以使種群所有個體在算法的后期搜索階段均處于活躍狀態(tài)。另一方面，IBSO 算法每次迭代均把最差聚類中心替換掉，可以加快種群個體向最優(yōu)解的逼近速度。

圖2 IBSO 算法流程圖Fig.2 Flowchart of IBSO algorithm

2.2 IBSO 算法在醫(yī)學(xué)圖像配準(zhǔn)中的應(yīng)用

本文將IBSO 算法的全局搜索能力與單純形搜索法的局部搜索能力進行優(yōu)勢互補，使用兩種算法協(xié)作完成醫(yī)學(xué)圖像配準(zhǔn)任務(wù)。首先，將原始圖像進行小波分解，分解后的圖像構(gòu)成圖像金字塔；然后，在低分辨率層使用IBSO 算法進行全局尋優(yōu)；最后，將IBSO 算法的尋優(yōu)值進行倍率縮放后作為單純形搜索法的搜索起點，在金字塔的第2 層及第1 層使用單純形搜索法進行局部搜索。具體算法步驟如下：

步驟1將待配準(zhǔn)圖像R和F進行兩次小波分解，源圖像作為圖像金字塔的第1 層圖像，第一次小波分解圖像作為金字塔第2 層，對第一次小波分解的圖像進行第二次小波分解作為圖像金字塔的第3層。若高一層的尋優(yōu)參數(shù)為(x,y,θ)，則(2x,2y,θ)作為下一層初始參數(shù)，平移參數(shù)翻倍，旋轉(zhuǎn)參數(shù)值保持不變。

步驟2使用IBSO 算法作為尋優(yōu)函數(shù)對圖像金字塔頂層圖像配準(zhǔn)，在單模圖像配準(zhǔn)時使用MI 作為配準(zhǔn)的測度函數(shù)，多模配準(zhǔn)將MI、NMI、NCC、CCRE 分別作為配準(zhǔn)的測度函數(shù)進行實驗。

步驟3將IBSO 尋優(yōu)算法得到的參數(shù)做相應(yīng)的倍率縮放后作為單純形搜索的搜索起點，搜索圖像金字塔的第2 層圖像。

步驟4將步驟3 中尋優(yōu)結(jié)果進行相應(yīng)的倍率縮放后，在圖像金字塔第1 層即原始圖像上進行尋優(yōu)。尋優(yōu)結(jié)束，得到配準(zhǔn)參數(shù)(x,y,θ)，使用上述配準(zhǔn)參數(shù)對浮動圖像F進行空間變換，融合變換后的浮動圖像與參考圖像（圖3）。

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

為了驗證本文所提算法的穩(wěn)定性與有效性，分別對MRI-MRI 單模態(tài)與CT-MRI 多模態(tài)醫(yī)學(xué)圖像進行配準(zhǔn)。實驗環(huán)境為：Windows10 操作系統(tǒng)，Matlab R2018a 實驗平臺。硬件平臺為：Intel(R)Core(TM)i7-9750H CPU@2.60 GHz、內(nèi)存16 GB。本文方法與3 種主流配準(zhǔn)算法比較算法誤差、耗時、配準(zhǔn)精度等參數(shù)，3 種算法分別為 PSO 與 Simplex 的結(jié)合[14]、差分進化（DE）算法與 Powell 算法的結(jié)合[15]、BSO 算法與 Powell 算法的結(jié)合[16]。以上 3 種算法分別記為 PSO+Simplex、DE+Powell、BSO+Powell。PSO、DE、BSO 與 IBSO 種群規(guī)模均為 50，PSO、DE、BSO、IBSO、Simplex、Powell 算法最大迭代次數(shù)均為200，其他參數(shù)與文獻[14-16]保持一致。

3.1 單模態(tài)配準(zhǔn)實驗與結(jié)果分析

單模態(tài)醫(yī)學(xué)圖像配準(zhǔn)采用BrainWeb 醫(yī)學(xué)圖像數(shù)據(jù)集。參考圖像R選取MRI-T1 加權(quán)圖像的第90層切片，切片厚度為1 mm，噪聲水平0%，如圖4(a)所示。浮動圖像F是在參考圖像的x軸方向平移7 像素、y軸方向平移3 像素，繞中心旋轉(zhuǎn)5°得到。對每種優(yōu)化算法重復(fù)進行100 次實驗，統(tǒng)計其平均誤差、最大誤差和平均耗時。

圖4 單模態(tài)配準(zhǔn)融合圖像Fig.4 Mono-modality registration fusion images

單模態(tài)配準(zhǔn)結(jié)果如表1 所示，其中Δx、Δy分別表示x、y方向的平移量誤差，Δθ表示旋轉(zhuǎn)的角度誤差。maxΔX、maxΔY分別表示 100 次實驗中x、y方向的最大平移量誤差，maxΔθ表示 100 次實驗中旋轉(zhuǎn)的角度最大誤差。當(dāng)平移量誤差小于1 像素，旋轉(zhuǎn)角誤差小于1°稱本次配準(zhǔn)達到了亞像素級[13]。

表1 單模態(tài)圖像配準(zhǔn)實驗結(jié)果對比Table 1 Comparison of experimental results of mono-modality image registration

觀察表1 可以發(fā)現(xiàn)4 種配準(zhǔn)算法均可以達到亞像素級配準(zhǔn)。綜合來看，與PSO+Simplex 和BSO+Powell 算法相比，本文算法的平均誤差、最大誤差均得到了一定程度的降低，平均耗時較兩種配準(zhǔn)算法分別降低了32.89%和13.66%；與DE+Powell 算法相比，本文算法角度誤差略高，其余指標(biāo)均得到一定程度的優(yōu)化，且平均耗時降低了13.91%。

3.2 多模態(tài)配準(zhǔn)實驗與結(jié)果分析

多模態(tài)圖像配準(zhǔn)采用Vanderbilt 大學(xué)Retrospective Image Registration Evaluation Project,Version 2.0 數(shù)據(jù)集。配準(zhǔn)結(jié)果采用MI、NCC、NMI、CCRE 4 種相似性測度函數(shù)作為指標(biāo)對配準(zhǔn)結(jié)果進行評價，相似性指數(shù)的值越高代表參考圖像與浮動圖像配準(zhǔn)效果越好。表2 為4 種配準(zhǔn)算法對參考圖像與浮動圖像分別使用MI、NMI、NCC、CCRE 作為測度函數(shù)進行100 次配準(zhǔn)實驗結(jié)果的平均值。配準(zhǔn)結(jié)果顯示本文算法MI、NMI、CCRE 與 NCC 均 優(yōu) 于 其 他 3 種 配 準(zhǔn) 算 法 。圖 5（a）為 CT 參考圖像，圖 5（b）為 MRI-T1 浮動圖像，圖 5（c）—（f）分別為 PSO+Simplex、DE+Powell、BSO+Powell 和本文算法對參考圖像與浮動圖像配準(zhǔn)融合后的圖像。

表2 多模態(tài)配準(zhǔn)實驗結(jié)果對比Table 2 Comparison of experimental results of multimodality registration

圖5 多模態(tài)配準(zhǔn)融合圖像Fig.5 Multi-modality registration fusion images

單模態(tài)與多模態(tài)配準(zhǔn)結(jié)果表明，本文算法較PSO+Simplex、DE+Powell、BSO+Powell 配準(zhǔn)算法在速度方面得到了一定幅度的提升，這得益于在IBSO 中每次迭代過程均使用當(dāng)前最優(yōu)解作為線索生成新個體替代本次迭代最差聚類中心，使算法收斂速度大大提高。雖然，在每次迭代中替換最差聚類中心個體的操作可以使得算法收斂速度顯著提高，但這個操作同時降低了原始BSO 的種群多樣性，增大了算法早熟收斂的風(fēng)險。從配準(zhǔn)結(jié)果來看，本文算法各項指標(biāo)均有不同幅度的提升，較上述3 種配準(zhǔn)算法性能優(yōu)越，具有更高的臨床使用價值。

4 結(jié)束語

為了克服測度函數(shù)局部極值多、配準(zhǔn)消耗時間長等問題，本文采用多分辨率策略，將改進頭腦風(fēng)暴優(yōu)化算法與單純形搜索法相結(jié)合提出一種新的配準(zhǔn)方法，并將該方法與3 種主流配準(zhǔn)算法進行了單模態(tài)與多模態(tài)配準(zhǔn)的實驗對比。實驗結(jié)果表明，本文算法在有效提高醫(yī)學(xué)圖像配準(zhǔn)精度的同時縮短配準(zhǔn)所用的時間，具有更高的臨床使用價值。所提算法目前在顱腦MRI-MRI、顱腦CT-MRI 圖像配準(zhǔn)實驗中取得了較好的配準(zhǔn)效果。而臨床診斷與治療中CT-PET、US-MRI 等多模態(tài)圖像之間的配準(zhǔn)同樣具有很高的實用價值，這些多模態(tài)配準(zhǔn)圖像因其成像差異需要選定不同的測度函數(shù)評估其對齊程度。本文重點研究了配準(zhǔn)過程中的優(yōu)化算法，將其與更多相似測度結(jié)合應(yīng)用于其他多模態(tài)圖像配準(zhǔn)是下階段的主要工作。