魏本征　甘　潔　尹義龍

(1.山東中醫(yī)藥大學(xué)理工學(xué)院，濟(jì)南，250355；2.山東中醫(yī)藥大學(xué)計(jì)算醫(yī)學(xué)實(shí)驗(yàn)室，濟(jì)南，250355；3.山東中醫(yī)藥大學(xué)第二附屬醫(yī)院放射科，濟(jì)南，250001；4.山東大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，濟(jì)南，250100)

引　　言

醫(yī)學(xué)圖像配準(zhǔn)是醫(yī)學(xué)圖像處理領(lǐng)域中的一項(xiàng)重要技術(shù)，在醫(yī)學(xué)基礎(chǔ)研究和臨床診斷治療中起著越來(lái)越重要的作用。醫(yī)學(xué)圖像配準(zhǔn)是使包含成像數(shù)據(jù)集中的信息最大化的一個(gè)重要步驟，其任務(wù)是對(duì)不同獲取時(shí)間、不同成像設(shè)備及不同獲取角度的數(shù)據(jù)進(jìn)行空間上的匹配，通過(guò)匹配可以得到多源醫(yī)學(xué)圖像信息的空間對(duì)應(yīng)關(guān)系。因此，醫(yī)學(xué)圖像配準(zhǔn)算法的主要任務(wù)就是確定一個(gè)映射，使圖像集在空間上關(guān)聯(lián)，從而使這些圖像能夠在一個(gè)統(tǒng)一的坐標(biāo)系上表示出來(lái)。

醫(yī)學(xué)圖像配準(zhǔn)技術(shù)的基本框架主要包括特征空間、搜索空間、相似性測(cè)度和搜索策略4個(gè)方面。其中，作為圖像配準(zhǔn)依據(jù)的圖像信息集合稱之為特征空間，其主要是指從兩幅待配準(zhǔn)圖像(參考圖像和浮動(dòng)圖像)中提取具有區(qū)別力的圖像信息特征，這對(duì)配準(zhǔn)算法的運(yùn)行速度和魯棒性等性能具有重要影響。配準(zhǔn)用圖像特征通常應(yīng)具有較高的普適性和穩(wěn)定性，較小的特征匹配計(jì)算量，以及特征提取簡(jiǎn)單、快捷等特性。當(dāng)前圖像配準(zhǔn)方法一般分為兩種，一種是基于圖像特征的配準(zhǔn)方法，其特征空間主要包括點(diǎn)、邊緣、曲線、曲面、形狀或是特征描述子，如輪廓矩、Zernike矩、小波描述符及方向鏈碼等。該類(lèi)方法主要是利用如點(diǎn)、線段和面等圖像形狀共同特征構(gòu)建配準(zhǔn)目標(biāo)函數(shù)，配準(zhǔn)速度快，但特征提取困難且不穩(wěn)定，誤配準(zhǔn)率較高。另一種是基于灰度的配準(zhǔn)方法[1-2]，圖像的灰度值特征是其最主要的特征空間。此類(lèi)配準(zhǔn)方法主要是對(duì)圖像的灰度進(jìn)行操作，算法思想簡(jiǎn)單，且無(wú)需預(yù)處理，然而該配準(zhǔn)過(guò)程計(jì)算量過(guò)大，易受到圖像噪聲和形狀變換的影響。在基于灰度的各種配準(zhǔn)測(cè)度中，互信息熵配準(zhǔn)方法具有魯棒性好、精度高等優(yōu)點(diǎn)，但其計(jì)算量大，而且當(dāng)圖像受低采樣、噪聲、變形等因素影響時(shí)，該配準(zhǔn)方法易出現(xiàn)誤匹配、魯棒性差等問(wèn)題[3]。

在諸多改進(jìn)圖像配準(zhǔn)算法的研究中，將幾何特征與像素相似性方法相結(jié)合，實(shí)現(xiàn)兩者優(yōu)勢(shì)互補(bǔ)，從而設(shè)計(jì)出更加穩(wěn)定、性能更好的相似性配準(zhǔn)測(cè)度，是一種有效和可行的研究思路。Rangarajan等人[4]提出了一種利用互信息匹配形狀特征點(diǎn)進(jìn)行配準(zhǔn)的策略。Pluim等人[1]通過(guò)提取圖像的梯度信息，并將其視為空間信息加入互信息中，獲得了較高的匹配精度?；鹪彽热薣5]提出一種基于輪廓特征點(diǎn)最大互信息的多源醫(yī)學(xué)圖像配準(zhǔn)方法。盧振泰等人[6]則采用先整體后局部的剛性配準(zhǔn)策略，提出一種基于等效子午面與互信息量的三維醫(yī)學(xué)圖像快速配準(zhǔn)算法。趙海峰等人[7]提出一種基于特征點(diǎn)Rényi互信息的醫(yī)學(xué)圖像配準(zhǔn)算法，該算法適于單模和多模醫(yī)學(xué)圖像配準(zhǔn)，圖像配準(zhǔn)速度較快，魯棒性較強(qiáng)。

借鑒上述研究方法，為進(jìn)一步提高互信息熵類(lèi)醫(yī)學(xué)圖像配準(zhǔn)方法的效率和精度，基于醫(yī)學(xué)圖像局部特征信息，本文提出了一種基于邊緣特征點(diǎn)互信息熵的醫(yī)學(xué)圖像配準(zhǔn)方法(Feature point based mutual information method，F(xiàn)PMI)。

1　邊緣特征點(diǎn)互信息熵配準(zhǔn)方法

為達(dá)到提高醫(yī)學(xué)圖像配準(zhǔn)質(zhì)量和配準(zhǔn)速度，降低圖像噪聲對(duì)圖像配準(zhǔn)影響的目的，基于圖像輪廓、邊緣特征點(diǎn)和互信息熵等圖像局部信息，所提出的FPMI算法主要包括：計(jì)算圖像熵及互信息量，提取待配準(zhǔn)圖像邊緣，構(gòu)建邊緣特征點(diǎn)互信息能量函數(shù)，以及選取配準(zhǔn)優(yōu)化策略等。

1.1　計(jì)算圖像熵及互信息量

熵表達(dá)的是一個(gè)系統(tǒng)的復(fù)雜性或者是不確定性。對(duì)于一幅圖像I來(lái)說(shuō)，它含有信息量的多少可以用熵表示為

H(I)=-∑pxlog2px

(1)

式中px表示醫(yī)學(xué)圖像中圖像灰度強(qiáng)度值為I(x)的像素點(diǎn)在所在圖像中的出現(xiàn)頻率，即相關(guān)統(tǒng)計(jì)概率。

圖像互信息量是指兩幅圖像中的共有信息量，當(dāng)兩幅圖像實(shí)現(xiàn)了正確匹配的時(shí)候，它們相關(guān)的互信息量最大。所以，在選擇醫(yī)學(xué)圖像配準(zhǔn)用的相似性測(cè)度函數(shù)時(shí)，可以選取用兩幅圖像的互信息熵作為它們圖像配準(zhǔn)的依據(jù)[8]。

在表示兩幅醫(yī)學(xué)圖像X和Y的互信息熵MI(X,Y)時(shí)，可選用聯(lián)合信息熵的形式,即

MI(X,Y)=H(X)+H(Y)-H(X,Y)

(2)

由Dobrushin公式[9]可推導(dǎo)得互信息熵的計(jì)算公式為

(3)

式中邊緣概率分布計(jì)算公式pX(x)和pY(y)描述了兩幅待配準(zhǔn)圖像X和Y的邊緣概率分布情況。另外，待配準(zhǔn)圖像X和Y的聯(lián)合概率分布是pXY(x,y)，直接計(jì)算比較困難，可借助于兩幅圖像的聯(lián)合直方圖h(x,y)進(jìn)行計(jì)算。h(x,y)表示兩幅圖像重疊部分圖像的灰度值為(x,y)的像素對(duì)總數(shù)，其中x，y分別表示兩圖中灰度等級(jí)，(x,y)表示兩圖中同一個(gè)位置的像素值，如果浮動(dòng)圖像溢出邊界，則兩幅圖像中不重合的部分不參加計(jì)算，或是用圖像背景灰度填充因浮動(dòng)圖像幾何變換而出現(xiàn)的空白區(qū)域。設(shè)兩幅圖像的重疊部分像素點(diǎn)為(i,j) (i= 0,1,2,…，M-1;j= 0,1,2,…，N-1)，其聯(lián)合直方圖可以表示為[10]

(4)

pY(y)可以近似為兩幅圖像重疊部分像素對(duì)的數(shù)量(圖像X中第x灰度級(jí)和圖像Y中第y灰度級(jí)構(gòu)成的像素對(duì))與總的像素對(duì)數(shù)(所有的灰度的像素對(duì)數(shù))之比。利用聯(lián)合直方圖，兩幅圖像的聯(lián)合概率分布為

(5)

邊緣概率分布pX(x)和pY(y)也可以通過(guò)聯(lián)合概率密度求得

(6)

(7)

通過(guò)上述公式的推導(dǎo)和計(jì)算，不難發(fā)現(xiàn)對(duì)于互信息熵的計(jì)算，可以轉(zhuǎn)化成為對(duì)兩幅圖像的聯(lián)合概率分布和邊緣概率分布的計(jì)算。

假設(shè)兩幅待配準(zhǔn)圖像X(i,j)和Y(i,j)的灰度級(jí)數(shù)都是256，尺寸大小為M×N，其中(i,j)是兩幅圖像重疊部分像素點(diǎn)的坐標(biāo)，兩幅圖像中未重疊的部分與計(jì)算無(wú)關(guān)。則計(jì)算兩幅圖像互信息熵MI(X,Y)的算法如下。

算法1圖像互信息熵的計(jì)算算法

(1) 確定參考圖像、浮動(dòng)圖像的坐標(biāo)；初始化聯(lián)合直方圖hXY(x,y)數(shù)據(jù)，并令hXY(x,y)=0；

(2) 對(duì)于兩幅圖像的重疊區(qū)域上的每一點(diǎn)對(duì)(i,j)，(i= 0, 1,2,…,M-1;j= 0,1,2,…，N-1),計(jì)算hXY(X(i,j),Y(i,j))=hXY(A(i,j),B(i,j))+1；

(5) 利用式(3)，可求得兩幅圖像的互信息熵。

1.2　圖像邊緣特征提取

為得到一個(gè)更加準(zhǔn)確的醫(yī)學(xué)圖像輪廓，有必要消除圖像噪聲的干擾，從而得到盡可能完整的圖像邊緣輪廓。相比傳統(tǒng)的邊緣檢測(cè)方法，Matherom和Serra根據(jù)幾何代數(shù)及拓?fù)湔撎岢龅臄?shù)學(xué)形態(tài)學(xué)具有更好的抑制噪聲效果[11]。形態(tài)學(xué)梯度算子可表示為

Grad(f)=(f⊕g)-(f?g)

(8)

式中：f(x,y)為原始圖像，g(x,y)為結(jié)構(gòu)元素，⊕表示膨脹運(yùn)算，?表示腐蝕運(yùn)算。顯而易見(jiàn)，形態(tài)學(xué)梯度算子能加劇輸入圖像的灰度級(jí)階躍變。

為使得提取的醫(yī)學(xué)圖像邊緣曲線更平滑，同時(shí)保證邊緣特征的完整性和細(xì)節(jié)豐富性，借鑒廣義形態(tài)濾波器原理，設(shè)計(jì)邊緣檢測(cè)算子[12]和形態(tài)學(xué)梯度濾波算子[13]，如式(9,10)所示。

IGrad(f)=(f°g)⊕g-(f·g)?g

(9)

FGrad(f)=c1IGrad1(f)+c2IGrad2(f)

(10)

式中：c1和c2的權(quán)系數(shù)大小采用最小均方自適應(yīng)方法確定；IGrad1()和IGrad2()為分別與3×3的十字形和交叉形結(jié)構(gòu)元素相對(duì)應(yīng)的邊緣檢測(cè)算子。

1.3　邊緣特征點(diǎn)互信息能量函數(shù)

由于待配準(zhǔn)的醫(yī)學(xué)圖像幾乎都基于共同解剖結(jié)構(gòu)信息，因此在盡量不減少配準(zhǔn)圖像特征信息的前提下，為減少互信息熵的計(jì)算量，提高配準(zhǔn)效率，基于互信息的配準(zhǔn)策略可采用相同特征點(diǎn)集的互信息熵來(lái)取代整幅醫(yī)學(xué)圖像的互信息熵作為配準(zhǔn)目標(biāo)函數(shù)[5]。

假設(shè)從兩幅待配準(zhǔn)醫(yī)學(xué)圖像中分別提取出的形狀特征點(diǎn)的集合為X={Xi,i=1,2,…,N1}和Y={Yi,i=1,2,…,N2},其中Xi和Yi是點(diǎn)在二維平面中的坐標(biāo)位置。則基于變換參數(shù)T，X和Y的距離測(cè)度可表示為

(11)

式中:N=min{N1,N2}，參照K-means 聚類(lèi)算法[14]，選取參數(shù)q=3.5。根據(jù)互信息的定義，可得點(diǎn)集的互信息熵

(12)

式中：Pij為兩幅圖像邊緣特征點(diǎn)Xi和Yj的聯(lián)合概率，即同時(shí)從X中選取Xi和從Y中選取Yj的聯(lián)合概率，Pij=P{I=i,i∈(1,2,…,N1),J=j,j∈(1,2,…,N2)}。由式(12)可知，當(dāng)兩幅圖像配準(zhǔn)時(shí)，對(duì)應(yīng)匹配點(diǎn)的空間位置一致，如果有N個(gè)Xi和Yj是對(duì)應(yīng)匹配點(diǎn)，N的取值應(yīng)為最大，此時(shí)圖像間的相關(guān)性最大，互信息量也取得最大值。因此實(shí)際計(jì)算時(shí)，可采用如下經(jīng)驗(yàn)計(jì)算公式[4]為

Pij(T)=exp(-αDij(T)-λ)

(13)

式中：α和λ是兩個(gè)拉格朗日常數(shù)，作為輔助參數(shù)分別用于約束點(diǎn)匹配相似測(cè)度和概率總值。特征點(diǎn)互信息能量函數(shù)可表示為

(14)

(15)

上述能量函數(shù)設(shè)計(jì)過(guò)程中，K-means聚類(lèi)算法被用于提取醫(yī)學(xué)圖像特征點(diǎn)集，基于邊緣點(diǎn)特征點(diǎn)集X和Y，選取聚類(lèi)中心K=200，即選取200個(gè)點(diǎn)。之后，需找到一個(gè)最優(yōu)的配準(zhǔn)參數(shù)，使得EMI最小。

1.4　配準(zhǔn)優(yōu)化策略

在具體應(yīng)用過(guò)程中，因?yàn)閰⒖紙D像或浮動(dòng)圖像的灰度特征信息以及圖像含有噪聲的干擾，互信息測(cè)度函數(shù)在配準(zhǔn)過(guò)程中存在大量局部極值點(diǎn)。文獻(xiàn)[10]已經(jīng)對(duì)配準(zhǔn)過(guò)程中局部極值的產(chǎn)生原因和可能采取的克服辦法進(jìn)行了專門(mén)的研究，對(duì)此本文不再進(jìn)行討論。在多維空間，因?yàn)榇罅繕O值點(diǎn)的存在使得參數(shù)優(yōu)化成為一個(gè)難題。從計(jì)算的角度看，在確定了圖像配準(zhǔn)目標(biāo)函數(shù)和特征點(diǎn)集之后，圖像配準(zhǔn)就變?yōu)榱艘粋€(gè)配準(zhǔn)參數(shù)尋優(yōu)的過(guò)程。相應(yīng)算法即是基于某種搜索策略和圖像變換方法，在搜索空間中找最優(yōu)解，求解目標(biāo)函數(shù)EMI最小值的過(guò)程。最優(yōu)解的求解過(guò)程就是求參考圖像和浮動(dòng)圖像之間變換模型參數(shù)的過(guò)程[15]。一般來(lái)說(shuō)，相似性測(cè)度的計(jì)算比較復(fù)雜，并且會(huì)產(chǎn)生大量數(shù)據(jù)，所以常采用最優(yōu)化求解法來(lái)減少運(yùn)算量，加快配準(zhǔn)速度。在實(shí)際求解過(guò)程中，幾乎每一種配準(zhǔn)過(guò)程都需要最優(yōu)化算法，用來(lái)搜索使目標(biāo)函數(shù)取得最小值(或使用相似性測(cè)度取得最大值)的最優(yōu)變換。對(duì)一個(gè)特定的配準(zhǔn)算法，配準(zhǔn)過(guò)程可表示為

(16)

式中：T為配準(zhǔn)變換，XR為參考圖像，XF為浮動(dòng)圖像，f為需要找到最小值的目標(biāo)函數(shù)。

圖像配準(zhǔn)參數(shù)的優(yōu)化主要有兩個(gè)要求：全局尋優(yōu)和快的優(yōu)化速度。目前常用的優(yōu)化方法主要有：梯度下降法、Powell算法和下降單純形法，以及模擬退火算法、遺傳算法和粒子群優(yōu)化算法等。其中前3種算法的優(yōu)點(diǎn)是收斂速度快、計(jì)算量小，缺點(diǎn)是易收斂到局部最優(yōu)；后3種算法的優(yōu)點(diǎn)是全局尋優(yōu)能力強(qiáng)，缺點(diǎn)是收斂速度慢、計(jì)算量大[16]。

本文選用梯度下降法作為配準(zhǔn)測(cè)度函數(shù)的優(yōu)化策略，其優(yōu)化算法求解過(guò)程和執(zhí)行步驟如下。

基于梯度的最優(yōu)化方法通常用來(lái)確定搜索方向，在該方向目標(biāo)函數(shù)的值局部遞減。設(shè)x為n維向量x=[x1,x2,…,xn]T，目標(biāo)函數(shù)f(x)的梯度向量能夠表示為

(17)

其二階偏導(dǎo)數(shù)能夠用一個(gè)Hessian矩陣表示

(18)

函數(shù)f(x)以xk的泰勒展開(kāi)近似表示為

f(xk+x)≈f(xk)+(

(19)

算法2基于梯度的最優(yōu)化算法

(1) 初始化。設(shè)置迭代器k=0，并初始化向量x，計(jì)算f(xk)。

(2) 驗(yàn)證收斂標(biāo)準(zhǔn)。如果滿足了收斂標(biāo)準(zhǔn)，則停止最優(yōu)化過(guò)程，得到解xk。

(3) 計(jì)算搜索方向。計(jì)算搜索方向向量Pk。

(4) 計(jì)算步長(zhǎng)lk。確定滿足f(xk+lkPk)

(5) 更新變量。令xk+1=xk+lkPk，計(jì)算f(xk+1)，令k=k+1，返回步驟(2)。

搜索方向Pk的計(jì)算和步長(zhǎng)lk的尋找是基于梯度優(yōu)化算法的兩個(gè)主要問(wèn)題，用不同的方法計(jì)算搜索方向就產(chǎn)生了各種不同的基于梯度的方法。

2　圖像配準(zhǔn)過(guò)程

Rangarajan算法[4]在定義了基于最大互信息熵后，又以其為基礎(chǔ)，通過(guò)引入輔助變量作代數(shù)變換定義了一個(gè)能量函數(shù)EMI。兩幅圖像配準(zhǔn)時(shí)分別作為參考圖像和浮動(dòng)圖像。通過(guò)對(duì)配準(zhǔn)圖像的空間變換，使得兩幅圖像的像素點(diǎn)盡可能對(duì)應(yīng)匹配，此時(shí)兩幅圖像的互信息達(dá)到最大值，且EMI取得最小值。在計(jì)算過(guò)程中，一般選取邊緣輪廓清晰、圖像特征提取穩(wěn)定的圖像作為參考圖像，另一幅作為浮動(dòng)圖像。

結(jié)合兩幅圖像邊緣特征點(diǎn)互信息熵的計(jì)算方法，采用梯度下降法作為圖像配準(zhǔn)優(yōu)化算法。算法流程如下。

算法3FPMI圖像配準(zhǔn)算法

(1) 用形態(tài)學(xué)梯度濾波算法獲取兩幅圖像輪廓，并采用K-means聚類(lèi)算法提取圖像邊緣特征點(diǎn),得到形狀特征點(diǎn)集X和Y，X為參考圖像點(diǎn)集,Y為浮動(dòng)圖像點(diǎn)集；

(2) 初始化，Ρij=(N1N2)-1,α=αmin，λ=λ0，κ=κmin，T=0；

(3) 確定算法的搜索空間，ΔT=(Δx,Δy,Δθ)，Δα，Δκ；

(4) 初始化由ΔT，Δα和Δκ構(gòu)成梯度優(yōu)化算法的5維向量x;

(5) 計(jì)算待配準(zhǔn)圖像的梯度損失函數(shù)值(兩幅圖像的EMI)；

(6) 根據(jù)式(15,19)更新x；

(7) 迭代結(jié)束條件通常為全局最優(yōu)解小于最小允許誤差,或達(dá)到一個(gè)預(yù)設(shè)最大代數(shù)Gmax，如未達(dá)到該結(jié)束條件則返回步驟(5)；

(8) 達(dá)到迭代結(jié)束條件，配準(zhǔn)完成。

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

本文將FPMI算法與基于最大互信息熵等算法在配準(zhǔn)精度、運(yùn)行時(shí)間等算法性能方面進(jìn)行比較。用Matlab 7.0在PC機(jī)上(P4 2.8 GHz CPU,1 GB內(nèi)存)分別對(duì)多源醫(yī)學(xué)圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)。利用直方圖估計(jì)概率密度分布時(shí)采用64位灰度級(jí)。配準(zhǔn)成功率設(shè)定為配準(zhǔn)成功的次數(shù)與實(shí)驗(yàn)總數(shù)的比值，且每套數(shù)據(jù)獨(dú)立進(jìn)行100次實(shí)驗(yàn)。利用形態(tài)學(xué)梯度濾波算法提取邊緣時(shí)采用的結(jié)構(gòu)元素為{0 1 0; 1 1 1;0 1 0}T和{1 0 1;0 -1 0; 1 0 1}T。優(yōu)化算法均采用梯度下降法，最大迭代步數(shù)200次，算法初始化的搜索空間為：x方向的搜索范圍x∈[-20, 20]，y方向的搜索空間y∈[-20, 20]，旋轉(zhuǎn)的搜索空間角度θ∈[-15°,15°]。

3.1　模擬配準(zhǔn)測(cè)試實(shí)驗(yàn)

3.1.1人工合成圖像配準(zhǔn)實(shí)驗(yàn)

為檢驗(yàn)配準(zhǔn)算法的正確性，本文采用人工合成圖像及其變換圖像作為配準(zhǔn)實(shí)驗(yàn)的參考和浮動(dòng)圖像，對(duì)所提FPMI算法做正確性驗(yàn)證配準(zhǔn)實(shí)驗(yàn)。待配準(zhǔn)圖像如圖1所示,圖1(a)為參考圖像，圖1(b)為浮動(dòng)圖像，兩幅圖像中白色矩形區(qū)域有一定的旋轉(zhuǎn)角度和平移量偏差。其變換參數(shù)(Δx,Δy,Δθ)為(40,20,20°)。

實(shí)際計(jì)算過(guò)程中，配準(zhǔn)成功的標(biāo)準(zhǔn)設(shè)定為：如果計(jì)算出的變換參數(shù)(Δx,Δy,Δθ)和平移量與真實(shí)值相差1個(gè)像素，且旋轉(zhuǎn)角度偏差值在1°范圍內(nèi)，即認(rèn)為配準(zhǔn)成功。100次算法驗(yàn)證實(shí)驗(yàn)，成功率為100%。配準(zhǔn)結(jié)果如表1所示，F(xiàn)PMI算法平移配準(zhǔn)精度在0.2像素以內(nèi)，旋轉(zhuǎn)角度精度在0.1°以內(nèi)。實(shí)驗(yàn)結(jié)果表明了FPMI算法的正確性。

圖1　人造合成配準(zhǔn)測(cè)試實(shí)驗(yàn)圖像Fig.1　Synthetic images for registration test

指標(biāo)Δx/像素Δy/像素Δθ/(°)設(shè)定變換值402020FPMI配準(zhǔn)值39．8220．1519．88配準(zhǔn)誤差均值0．180．150．12配準(zhǔn)參數(shù)標(biāo)準(zhǔn)差0．070．060．06

3.1.2CT-MRI配準(zhǔn)實(shí)驗(yàn)

為檢驗(yàn)FPMI算法的有效性，本文對(duì)CT-MRI圖像進(jìn)行模擬配準(zhǔn)對(duì)比實(shí)驗(yàn)，結(jié)果如圖2所示，圖2(a)為CT圖像，作為參考圖像，圖2(b)為MRI圖像，作為浮動(dòng)圖像。每幅圖像的大小為512像素×512像素，圖像取自BIS醫(yī)學(xué)圖像庫(kù) (http://nova.nlm.nih.gov/Mayo)。在實(shí)驗(yàn)中，圖像配準(zhǔn)分別采用最大互信息熵算法(MI算法)、具有代表性的特征互信息量配準(zhǔn)算法(Rangarajan算法[4])以及本文的FPMI算法。各算法均獨(dú)立運(yùn)行100次，實(shí)驗(yàn)結(jié)果如圖2(c,d)及表2，3所示。其中，圖2(c)為FPMI算法求解的MRI變換后圖像，圖2(d)為配準(zhǔn)后的融合圖像；表2是算法模擬配準(zhǔn)變換參數(shù)(Δx,Δy,Δθ)的均值統(tǒng)計(jì)表；表3是3種算法的配準(zhǔn)指標(biāo)統(tǒng)計(jì)表，包括4個(gè)指標(biāo)：平均互信息量、互信息量標(biāo)準(zhǔn)差δ、單次平均計(jì)算時(shí)間和配準(zhǔn)成功率。

從實(shí)驗(yàn)結(jié)果可以看出，對(duì)模擬CT、MRI圖像進(jìn)行配準(zhǔn)時(shí)，F(xiàn)PMI的平均配準(zhǔn)偏差較小，配準(zhǔn)參數(shù)(Δx,Δy,Δθ)的精度高于另外兩種方法；FPMI算法的平均配準(zhǔn)互信息量略高于MI算法的平均互信息量，這在一定程度上保證了算法配準(zhǔn)的準(zhǔn)確度。MI算法成功率略低，為98%；Rangarajan算法和FPMI兩種算法成功率均可以達(dá)到100%，但FPMI算法運(yùn)算速度快，其運(yùn)算速度約為Rangarajan算法的1.5倍，MI算法的2.4倍。因此，F(xiàn)PMI算法的綜合性能優(yōu)于另外兩種算法。

圖2　腦部MRI-CT圖像模擬配準(zhǔn)示意圖Fig.2　Simulation registration of brain MRI and CT image

算　法Δx/像素Δy/像素Δθ/(°)MI算法11．2729．288．14FPMI算法11．2029．178．31Rangarajan算法11．1429．158．47

表3　CT-MRI模擬配準(zhǔn)各算法指標(biāo)統(tǒng)計(jì)表

3.2　測(cè)度函數(shù)曲線銳度測(cè)試實(shí)驗(yàn)

為對(duì)FPMI配準(zhǔn)測(cè)度函數(shù)進(jìn)行性能測(cè)試，本文選取MRI腦圖像進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)所用的兩幅MRI腦圖像取自BrainWeb(http://www.bic.mni.mcgill.ca/brainweb/.)，分別是大小為181×217×181的T1加權(quán)和PD加權(quán)的腦部MRI圖像，其體素大小為1.0 mm3。圖3為測(cè)試實(shí)驗(yàn)圖像及邊緣特征點(diǎn)集示意圖。如圖3(a)所示，將MRI-T1圖像沿水平方向上平移一定像素后與MRI-PD圖像(圖3(c))之間計(jì)算測(cè)度值，觀察峰值的變化，像素平移的范圍是[-10,10]，負(fù)號(hào)表示向左平移，正號(hào)表示向右平移；同樣將MRI-T1圖像繞中心旋轉(zhuǎn)一定的角度后與MRI-PD圖像計(jì)算測(cè)度值，再觀察峰值的變化，角度旋轉(zhuǎn)的范圍是[-100°,100°]，負(fù)號(hào)表示順時(shí)針旋轉(zhuǎn)，正號(hào)表示逆時(shí)針旋轉(zhuǎn)。圖3(b,d)分別為圖3(a,c)的圖像邊緣特征點(diǎn)集示意圖。

本文選用的配準(zhǔn)測(cè)度函數(shù)是互信息量和FPMI的邊緣特征點(diǎn)互信息量，結(jié)果如圖4所示。圖4(a,b)分別表示兩測(cè)度函數(shù)值隨圖像在水平方向平移像素和角度旋轉(zhuǎn)時(shí)的變化曲線，其中，橫軸分別表示水平方向像素的平移量和旋轉(zhuǎn)角度，縱軸均表示測(cè)度函數(shù)值。從圖4可得，兩圖像在水平平移為0像素及旋轉(zhuǎn)角度為0°時(shí)測(cè)量函數(shù)均有最大值，F(xiàn)PMI的配準(zhǔn)方法中測(cè)度函數(shù)的曲線峰值比MI更加尖銳，更有利于最優(yōu)變換參數(shù)的選取。

圖3　測(cè)試實(shí)驗(yàn)圖像及邊緣特征點(diǎn)集示意圖Fig.3　Test images and the edge feature point set

圖4　配準(zhǔn)測(cè)度函數(shù)變化示意圖Fig.4　The comparison chart of registration measure function

3.3　多模醫(yī)學(xué)圖像配準(zhǔn)實(shí)驗(yàn)

為了檢驗(yàn)FPMI新算法在多模圖像配準(zhǔn)上的配準(zhǔn)精度，本文選取兩組同一病人的腦部圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)，如圖5所示。每幅圖像的大小均為256×256，并且選取MRI圖像(圖5(a，c))作為參考圖像，CT(圖5(b))和SPECT(圖5(d))分別作為浮動(dòng)圖像，將本文算法與MI算法和Rangarajan算法進(jìn)行了對(duì)比，3種算法獨(dú)立測(cè)試50次。表4給出了MRI-CT和MRI-SPECT 的兩組圖像的配準(zhǔn)結(jié)果。

在MRI-CT的配準(zhǔn)過(guò)程中，F(xiàn)PMI算法和Rangarajan算法均出現(xiàn)一次誤配，成功率達(dá)到了98%，而MI的成功率是94%；在運(yùn)算效率方面，F(xiàn)PMI的運(yùn)算效率約是MI算法的3.5倍，是Rangarajan算法的1.8倍。

圖5　多模醫(yī)學(xué)配準(zhǔn)圖像Fig.5　Multi-modal medical image registration

在MRI-SPECT的配準(zhǔn)實(shí)驗(yàn)中，3種算法均出現(xiàn)了不同程度的誤配準(zhǔn)現(xiàn)象，F(xiàn)PMI算法成功率達(dá)到了90%，而MI算法只取得了72%的成功率。這是由于SPECT圖像灰度分布較為模糊，且其形狀輪廓與MRI圖像差別較大，因此導(dǎo)致了分別基于灰度、輪廓及邊緣特征點(diǎn)的3種配準(zhǔn)算法的準(zhǔn)確率出現(xiàn)下降。在運(yùn)算效率方面，F(xiàn)PMI的運(yùn)算時(shí)間達(dá)到了1 027 s，其效率約為MI的2.3倍，Rangarajan算法的1.9倍。與CT-MRI配準(zhǔn)實(shí)驗(yàn)相比，在組織結(jié)構(gòu)更為復(fù)雜的醫(yī)學(xué)圖像配準(zhǔn)中MI和Rangarajan算法的運(yùn)算效率降低幅度更大，誤配率更高。

因此，對(duì)比實(shí)驗(yàn)結(jié)果可知，相對(duì)于另外兩種算法，F(xiàn)PMI算法在配準(zhǔn)精度和運(yùn)算效率上具有較為明顯的優(yōu)勢(shì)。且將其應(yīng)用于邊緣輪廓更為清晰的醫(yī)學(xué)圖像配準(zhǔn)時(shí)，其算法的性能優(yōu)勢(shì)會(huì)發(fā)揮得更好。

從上述實(shí)驗(yàn)可以看出，與MI算法及Rangarajan算法相比較，F(xiàn)PMI算法在配準(zhǔn)精度和運(yùn)算效率上都有一定優(yōu)勢(shì)。實(shí)驗(yàn)結(jié)果還表明， FPMI算法中采用MGZ算法提取圖像輪廓特征較為穩(wěn)定和準(zhǔn)確， 所以此方法對(duì)于圖像輪廓相對(duì)清晰的醫(yī)學(xué)圖像是一種較理想的配準(zhǔn)方法，但該方法對(duì)數(shù)據(jù)的缺失及圖像輪廓的形狀變異較為敏感。此外，相對(duì)于同類(lèi)算法，在配準(zhǔn)過(guò)程中本算法無(wú)需執(zhí)行粗配準(zhǔn)和精配準(zhǔn)兩個(gè)過(guò)程，因此算法的配準(zhǔn)過(guò)程簡(jiǎn)單、實(shí)用性更強(qiáng)。

表4　兩組多模醫(yī)學(xué)圖像配準(zhǔn)結(jié)果統(tǒng)計(jì)表

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

基于互信息熵配準(zhǔn)測(cè)度函數(shù)，本文設(shè)計(jì)了一種基于多特征信息融合的邊緣特征點(diǎn)互信息配準(zhǔn)測(cè)度，并利用梯度下降法優(yōu)化配準(zhǔn)過(guò)程。所提出的FPMI算法能反映醫(yī)學(xué)圖像的空間形狀特征，可大幅減少互信息的計(jì)算量，大大縮短配準(zhǔn)時(shí)間，并可有效克服局部極值問(wèn)題。此外，本文還研究用梯度下降法進(jìn)行醫(yī)學(xué)圖像配準(zhǔn)的參數(shù)尋優(yōu)，以進(jìn)一步提高配準(zhǔn)的效率和準(zhǔn)確度。研究結(jié)果表明，該方法能夠準(zhǔn)確、快速地處理剛性配準(zhǔn)問(wèn)題，在多源醫(yī)學(xué)圖像的配準(zhǔn)中具有一定的先進(jìn)性，是一種穩(wěn)健的自動(dòng)配準(zhǔn)方法，可用于臨床研究，具有較高的實(shí)用價(jià)值。但是，本文算法并不能適用于所有各類(lèi)醫(yī)學(xué)圖像，因此下一步工作將針對(duì)不同醫(yī)學(xué)圖像的配準(zhǔn)需求研究更具有針對(duì)性的配準(zhǔn)測(cè)度算法。

參考文獻(xiàn):

[1]Pluim J P W,Antoine M J B,Viergever M A. Image registration by maximization of combined mutual information and gradient information[J]. IEEE Transactions on Medical Imaging, 2002,19(8): 809-814.

[2]陳慶芳, 吳小俊. 基于分塊互信息和量子粒子群算法的圖像配準(zhǔn)[J]. 數(shù)據(jù)采集與處理, 2011, 26(4):473-477.

Chen Qingfang, Wu Xiaojun. Image registration based on block mutual information and quantum-behaved particle swarm optimization[J]. Journal of Data Acquisition & Processing, 2011, 26(4)：474-477.

[3]何彥杰, 周焰, 劉超,等. 結(jié)合PCM和MIE的多模態(tài)圖像配準(zhǔn)方法[J]. 空軍預(yù)警學(xué)院學(xué)報(bào), 2010, 24(5):372-375.

He Yanjie, Zhou Yan, Liu Chao, et al．Method of multimode image registration based on PCM and MIE[J]. Journal of Air Force Radar Academy, 2010, 24(5)：372-375.

[4]Rangarajan A,Chui H,Duncan J S.Rigid point feature registation using mutual information[J].Med Image Anal,1999,3(4):425-440.

[5]火元蓮，齊永鋒，宋海聲.基于輪廓特征點(diǎn)最大互信息的多源醫(yī)學(xué)圖像配準(zhǔn)[J]. 激光與紅外, 2008,38(1):96-98.

Huo Yuanlian，Qi Yongfeng，Song Haisheng．Multi-modality medical image registration based on mutual information of feature points[J]. Laser & Infrared, 2008, 38(1)：96-98.

[6]盧振泰，馮衍秋，馮前進(jìn)，等. 基于等效子午面與互信息量的醫(yī)學(xué)圖像配準(zhǔn)[J].計(jì)算機(jī)學(xué)報(bào), 2009,32(8):1611-1617.

Lu Zhentai, Feng Yanqiu, Feng Qianjin, et al. Medical image registration using equivalent meridian plane and mutual information[J]. Journal of Computer, 2009, 32(8)：1611-1617.

[7]趙海峰，陸明，卜令斌，等. 基于特征點(diǎn)Rényi互信息的醫(yī)學(xué)圖像配準(zhǔn)[J].計(jì)算機(jī)學(xué)報(bào), 2015,38(6):1212-1220.

Zhao Haifeng, Lu Ming, Bu Lingbin, et al. Medical image registration based on feature points and Rényi mutual information[J]. Chinese Journal of Computers, 2015, 38(6)：1212-1221.

[8]Penney G P, Weese J, Little J A, et al. A comparison of similarity measures for use in 2D-3D medical image registration[J]. IEEE Trans Med Imag, 1999, 17(4):586-595.

[9]Wells W M,Viola P, Atsumi H, et al. Multi-modal volume registration by maximization of mutual information[J].Med Image Anal,1996,1(1)：35-51.

[10] Maes F, Collignon A, Vandermeulen D, et al. Multimodality image registration by maximization of mutual information[J]. IEEE Trans Med Imag, 1997, 16(2): 187-198.

[11] Serra J. Image analysis and mathematical morphology[M]. New York, NY, USA: Academic Press, 1982:1-610.

[12] 魏本征,趙志敏,宋一中.基于IPSO和綜合信息的醫(yī)學(xué)圖像配準(zhǔn)新方法[J].光電子·激光,2009, 20(9) : 1271-1274.

Wei Benzheng, Zhao Zhimin, Song Yizhong. A novel algorithm for multimodality medical image registration based on IPSO algorithm and hybrid information[J]. Journal of Optoelectronics Laser , 2009, 20(9): 1271-1274 .

[13] 魏本征，趙志敏，華晉.基于形態(tài)學(xué)梯度和Zernike矩的亞像素邊緣檢測(cè)方法[J].儀器儀表學(xué)報(bào)，2010，31(4):838-844.

Wei Benzheng, Zhao Zhimin, Hua Jin. Sub-pixel edge detection method based on improved morphological gradient and Zernike moment[J]. Chinese Journal of Scientific Instrument, 2010, 31(4): 838 - 844 .

[14] Hamerly G, Elkan C. Alternatives to the k-means algorithm that find better clusterings[C]∥Proc of the ACM Conference on Information and Knowledge Management, CIKM-2002．[S.l.]:CIKM, 2002:600-607.

[15] 王偉.醫(yī)學(xué)圖像非剛性配準(zhǔn)算法研究[D].大連：大連理工大學(xué)，2012.

Wang Wei. Research on non-rigid registration algorithm of medical images[D]. Dalian： Dalian University of Technology, 2012.

[16] 龔曉彥.基于互信息的醫(yī)學(xué)圖像配準(zhǔn)算法研究[D].秦皇島：燕山大學(xué)，2010.

Gong Xiaoyan. Research on medical image registration based on mutual information[D]. Qihuangdao: Yanshan University, 2010.