曹國(guó)剛，朱信玉，陳 穎，曹 聰，孔德卿

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

引 言

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

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

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

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

1.1 相似性測(cè)度

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

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

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

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

Wang 等[8]基于累計(jì)剩余熵理論提出一種CCRE 測(cè)度,其通過(guò)實(shí)驗(yàn)證明了CCRE 具有較高的魯棒性，CCRE 定義為

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

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

1.2 多分辨率策略

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

小波變換將時(shí)域信號(hào)變換成頻域信號(hào)，小波分解后圖像的低頻成分保留原始圖像的基本信息，高分辨率圖像保留圖像的細(xì)節(jié)部分以及夾雜了高頻噪聲。因此在低分辨率的圖像上使用全局尋優(yōu)算法進(jìn)行粗配準(zhǔn)，得到尋優(yōu)參數(shù)作為局部尋優(yōu)算法的初始點(diǎn)，再對(duì)高分辨率圖像進(jìn)行精配準(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)暴過(guò)程建模，Shi 于2011 年提出了頭腦風(fēng)暴優(yōu)化算法(Brain storm optimization algorithm, BSO)[18]，并于2016 年又進(jìn)行了綜述[19]。BSO 算法將種群中每個(gè)想法個(gè)體看作是d維問(wèn)題的一個(gè)潛在解，種群中第i個(gè)想法個(gè)體可表示為

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

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

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

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

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

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

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

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

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

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

1.4 單純形搜索法

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

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

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

2.1 IBSO 算法

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

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

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

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

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

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

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

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

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

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

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

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

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

為了驗(yàn)證本文所提算法的穩(wěn)定性與有效性，分別對(duì)MRI-MRI 單模態(tài)與CT-MRI 多模態(tài)醫(yī)學(xué)圖像進(jìn)行配準(zhǔn)。實(shí)驗(yàn)環(huán)境為：Windows10 操作系統(tǒng)，Matlab R2018a 實(shí)驗(yàn)平臺(tái)。硬件平臺(tái)為：Intel(R)Core(TM)i7-9750H CPU@2.60 GHz、內(nèi)存16 GB。本文方法與3 種主流配準(zhǔn)算法比較算法誤差、耗時(shí)、配準(zhǔn)精度等參數(shù)，3 種算法分別為 PSO 與 Simplex 的結(jié)合[14]、差分進(jìn)化（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ù)與文獻(xiàn)[14-16]保持一致。

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

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

圖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 次實(shí)驗(yàn)中x、y方向的最大平移量誤差，maxΔθ表示 100 次實(shí)驗(yàn)中旋轉(zhuǎn)的角度最大誤差。當(dāng)平移量誤差小于1 像素，旋轉(zhuǎn)角誤差小于1°稱本次配準(zhǔn)達(dá)到了亞像素級(jí)[13]。

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

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

3.2 多模態(tài)配準(zhǔn)實(shí)驗(yà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 種相似性測(cè)度函數(shù)作為指標(biāo)對(duì)配準(zhǔn)結(jié)果進(jìn)行評(píng)價(jià)，相似性指數(shù)的值越高代表參考圖像與浮動(dòng)圖像配準(zhǔn)效果越好。表2 為4 種配準(zhǔn)算法對(duì)參考圖像與浮動(dòng)圖像分別使用MI、NMI、NCC、CCRE 作為測(cè)度函數(shù)進(jìn)行100 次配準(zhǔn)實(shí)驗(yàn)結(jié)果的平均值。配準(zhǔn)結(jié)果顯示本文算法MI、NMI、CCRE 與 NCC 均 優(yōu) 于 其 他 3 種 配 準(zhǔn) 算 法 。圖 5（a）為 CT 參考圖像，圖 5（b）為 MRI-T1 浮動(dòng)圖像，圖 5（c）—（f）分別為 PSO+Simplex、DE+Powell、BSO+Powell 和本文算法對(duì)參考圖像與浮動(dòng)圖像配準(zhǔn)融合后的圖像。

表2 多模態(tài)配準(zhǔn)實(shí)驗(yàn)結(jié)果對(duì)比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 中每次迭代過(guò)程均使用當(dāng)前最優(yōu)解作為線索生成新個(gè)體替代本次迭代最差聚類中心，使算法收斂速度大大提高。雖然，在每次迭代中替換最差聚類中心個(gè)體的操作可以使得算法收斂速度顯著提高，但這個(gè)操作同時(shí)降低了原始BSO 的種群多樣性，增大了算法早熟收斂的風(fēng)險(xiǎn)。從配準(zhǔn)結(jié)果來(lái)看，本文算法各項(xiàng)指標(biāo)均有不同幅度的提升，較上述3 種配準(zhǔn)算法性能優(yōu)越，具有更高的臨床使用價(jià)值。

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

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