(1.廣州軍區廣州總醫院 放射科 廣州 510010; 2.華南理工大學 自動化科學與工程學院 廣州 510640)
摘 要:為了實現胸部醫學圖像的自動配準,提出了一種基于人工免疫及最大互信息的配準方法。首先定義了最優配準的目標函數,接著運用人工免疫算法,結合最大互信息熵函數產生最優的仿射變換系數,從而實現醫學圖像自動配準。該人工免疫算法中,抗原是指最大互信息熵目標函數,而抗體是指最優的仿射變換系數。實驗證明該方法配準效果較好。
關鍵詞:人工免疫算法; 多模圖像配準; 最大互信息; 最優化; 胸部
中圖分類號:TP391.4文獻標志碼:A
文章編號:1001-3695(2009)05-1951-04
Thorax multimodal image registration based on artificial immune algorithm and maximization of mutual information
LI Bin1 OU Shan-xing1 TIAN Lian-fang2 MAO Zong-yuan2
(1.Dept. of Radiology Pediatrics Guangzhou General Hospital of Guangzhou Command Guangzhou 510010 China; 2.School of Automation Science Engineering South China University of Technology Guangzhou 510640China)
Abstract:In order to realize the automatic registration of pulmonary medical images this paper proposed an image registration method based on artificial immune algorithm and maximization of mutual information. Firstly presented the objective function for the optimal registration. Then,presented an artificial immune approach to automatically generate the optimal parameters of affine transform thus realized the multimodal image registration. In this approach regarded the objective function as antigens and regarded the optimal parameters of affine transform as antibodies. Experiments demonstrate the good performance of the proposed method.
Key words: artificial immune algorithm; multimodal image registration; maximization of mutual information; optimization; thorax
在胸部臨床診斷中通常需要對同一個病人進行多種模式或同一種模式的多次成像,即同時從幾幅圖像獲得信息,進行信息融合、綜合分析。單一模式成像只使用一種成像設備,可用于觀察病灶生長,對比手術前后的治療效果等。當一種成像設備所提供的信息不能滿足需要時,可以采用多種模式成像。醫學圖像的配準是融合的基礎。胸部醫學圖像的配準方法可分為剛性配準[1]、仿射配準[2]和彈性配準[3]。心肺是運動的器官,因此在胸部臨床診斷的應用中,配準方法是計算時間和精確性、魯棒性要求的一個折中。提供一個與人工指導配準精確性相當而且快速的自動配準方法,仍然是一個大的挑戰。文獻[4]比較了隨機采集圖像的剛性、仿射和彈性配準,得出結論:仿射配準是速度、精確性和魯棒性間較好的折中。
在配準中,最大互信息法[5]以互信息作為相似性測度,具有較強的魯棒性,是目前應用較多的一種方法。配準問題的本質是一個多參數優化問題,優化算法的優劣直接決定了配準是否成功。局部極值問題一直是最大互信息配準算法研究的中心問題。局部極值使優化算法過早收斂,導致誤配準的嚴重后果。文獻[6]基于遺傳算法,實現醫學圖像的配準。然而,遺傳算法容易出現早熟現象。文獻[7]采用人工免疫算法對遙感圖像進行配準。人工免疫算法模擬生物免疫系統中抗體對抗原的學習、記憶、自適應調節等功能,從處理數據中提取有用信息,并找出數據間的相關特性。同時,它有效利用了待求問題中的一些有用信息,可以抑制優化過程中產生的退化現象,迅速把問題的最優解限制在一個較小的空間范圍內[8,9]。
在研究上述文獻基礎上,本文提出一種基于人工免疫算法及最大互信息的胸部多模醫學圖像自動配準方法。
1 圖像配準問題描述及目標函數的定義
1.1 配準原理
將在不同時間、不同成像條件下對同一景物獲取的兩幅或多幅圖像在空間上對準,這就是配準。假設有兩幅圖像,定義IR為參考圖像,IF為待配準圖像,也即浮動圖像,IR(i,j)和IF(i,j)分別代表兩幅圖像在(i,j)位置的灰度值。將這兩幅圖像的配準定義為兩幅圖像之間在坐標位置和灰度級上的雙重映射關系:
IR(i,j)=g(IF(T(i,j)))(1)
其中:T是一個二維空間坐標變換;g是一個一維的灰度變換。
在圖像配準中,尋找兩幅圖像之間的最佳空間或幾何變換是配準的關鍵所在,而灰度變換g可以根據實際情況進行調整或忽略,通常是將映射函數表達為兩個單值函數Ti、Tj:
IR(i,j)=IF(Ti(i,j),Tj(i,j))(2)
1.2 仿射配準
用于圖像配準的空間變換關系有剛體變換、仿射變換、透視變換、非線性變換等。考慮到醫學圖像配準速度、精確性和魯棒性間的折中,本文采用仿射變換。其表示式如下:
x′y′1=e11 e12 e13e21 e22 e23001 xy1(3)
其中:(x′,y′)為對浮動圖像IF配準后的點;(x,y)為對浮動圖像IF配準前的點。令e11=k cos θ,e12=k sin θ,e13=tx,e21=-k cos θ,e13=ty。其中:θ是旋轉角;k是縮放比例因子;tx是X方向上的平移;ty是Y方向上的平移。
1.3 以互信息為相似性測度
本文采用基于灰度信息的配準方法。其主要特點是實現簡單、配準精度高,但其在最優變換的搜索過程中需要巨大的運算量。一般來說,基于灰度的圖像配準方法隱含地認為兩幅圖像對應空間位置的灰度是獨立同分布的。根據這一假設,互信息可以作為相似性度量。最大互信息法[5]在多模圖像配準中應用廣泛,其配準精度一般高于基于分割的圖像配準方法。該方法不需要對圖像進行分割、特征提取等預處理,且幾乎可以用于任意模態圖像的配準,并具有較強的魯棒性。
設參考圖像為IR,浮動圖像為IF,圖像IR、IF的信息熵H(IR)、H(IF)表示如下:
H(IR)=-∑IRP(IR)log P(IR)(4)
H(IF)=-∑IFP(IF)log P(IF)(5)
其中:P(IR)、P(IF)分別表示參考圖像IR和浮動圖像IF的邊緣概率密度。令H(IR,IF)代表IR、IF圖像的聯合信息熵:
H(IR,IF)=-∑IR,IFlog P(IR,IF)(6)
其中:P(IR,IF)表示聯合概率分布。
若H(IR|IF)表示已知圖像IF是IR的條件熵,那么H(IR)與H(IR|IF)的差值代表在圖像IF中所包含的IR的信息量,即互信息。因此兩個圖像間的互信息I(IR|IF)可用式(7)描述:
I(IR,IF)=H(IR)+H(IF)-H(IR,IF)=
H(IR)-H(IR|IF)=H(IF)-H(IF|IR)(7)
則圖像間的互信息定義如下:
MI(IR,IF)=H(IR)+H(IF)-H(IR,IF)(8)
當兩幅圖像之間嚴格匹配時,互信息MI(IR,IF)達到峰值。在胸部多模醫學圖像配準中,雖然兩幅圖像一般都來自于不同的成像設備,但它們基于共同的人體解剖信息,所以當兩幅圖像的空間位置達到一致時,其中一幅圖像表達另一幅圖像的信息,即其互信息應為最大。
通常用聯合概率分布P(IR,IF)和完全獨立時的邊緣概率分布P(IR)和P(IF)間的廣義距離來估計互信息MI(IR,IF):
MI(IR,IF)=∑IR,IFP(IR.IF)log(P(IR,IF)/(P(IR)#8226;P(IF)))(9)
其中P(IR,IF)為
P(IR,IF)=h(IR,IF)/∑IR,IFh(IR,IF)(10)
其中:h(IR,IF)是IR、IF的聯合灰度直方圖。
本文采用基于互信息的配準方法。其基本思想是:首先對待配準圖像作幾何變換;然后根據灰度信息的統計特性——互信息熵定義一個目標函數,作為參考圖像與浮動圖像之間的相似性測度,使得配準參數在目標函數的極值處得到,并以此作為配準的判決準則和配準參數最優化的目標函數,從而將配準問題轉換為多元函數的極值問題;最后通過一定的最優化方法求得正確的幾何變換參數。
Studholme等人[5]發現互信息本身的大小與待配準兩圖像間的重疊度具有一定的關聯性,為了消除這種關聯關系,本文應用的標準化互信息方法如式(11)。實驗證明它比式(8)的互信息法更具有魯棒性。
MI(IR,IF)=(h(IR)+h(IF))/h(IR,IF)(11)
1.4 問題描述及目標函數定義
綜合1.2節與1.3節,可以將圖像自動配準問題描述如下:當將兩幅圖像進行自動配準時,以互信息為相似性測度,尋找一組最優的參數P={tx,ty,θ,k},使得衡量配準后的圖像I′F和參考圖像IR的相似測度——互信息J(P)=maxV(J(P))最大(V表示可能的參數集合),也就是對任何P≠P都滿足J(P)<J(P)。
為了獲得最優的參數,可以對每一個確定的tx、ty、θ、k(tx∈Tx={-256,-255,…,255,256},ty∈Ty={-256,-255,…,255,256},θ∈Θ={0,1,…,359},k∈K={0,1,2,3,4,5}),求J(tx,ty,θ,k)=maxV(J(P)),直到J(tx,ty,θ,k)達到一個極大值為止,但這樣運算量太大,工程實際上是不可行的。因此,可以運用人工免疫算法,將本配準問題轉換為四維空間Tx×Ty×Θ×K,求J(tx,ty,θ,k)=maxV(J(P))的問題。
2 基于人工免疫及最大互信息的圖像配準算法
2.1 配準算法所用人工免疫算法的機理
配準算法優化過程主要利用了免疫系統的兩大重要特性[8],即免疫記憶能力和自適應能力。免疫系統產生的抗體與入侵的抗原相結合,同時產生記憶細胞,將與這種抗原的接觸記憶下來,從而產生再次免疫應答;同時,免疫系統通過對抗體產生促進和抑制、加入新產生的B淋巴細胞和清除無用的B淋巴細胞,自適應地調節高親和力抗體的數量,保證抗體的多樣性。
2.2 相關的基本概念[10,11]
1)抗原 算法中的抗原是指目標函數(式(11))。
2)抗體 算法中的抗體是指目標函數的優化解,本算法的抗體即是最優結果的一組仿射變換參數,使得圖像得到配準。
3)信息熵 算法中為了表明群體中抗體的多樣性,引入信息熵的概念[11]。在免疫算法中有一個抗體就用一個二進制位串來表示。假設有N個抗體(即位串),每個抗體有M位,每位可供選擇的字母表中共有S個字母k1,k2,k3,…,ks,則這N個抗體的信息熵為
H(N)=1/M∑Mj=1Hj(N)(12)
其中:Hj(N)=∑Sj=1-Pij log Pij,它是第N個抗體第j位的信息熵;Pij為N個抗體中的第j位為字母ki的概率。本算法中,用30位的二進制位串來表示一個抗體,即M=30。 其中:用9位表示旋轉角θ;9位表示tx,tx 為X方向上的平移,首位為正負;9位表示ty,ty為Y方向上的平移,首位表示正負;3位表示縮放比例因子k。令e為編碼的值,平移的解碼值d=-255+e;旋轉角度的解碼值d=(e/511)×360。
4)抗體間的親和力 用于表明兩抗體之間的相似度。抗體V和抗體W的相似度為
SVW=1/(1+H(2))(13)
5)抗原與抗體的親和力 用于表明抗體對抗原的識別程度。本算法中將抗體V和抗原V的親和力定義為
AV=fitness(xV)+praise(xV) (14)
其中:fitness(xV)表示抗原(問題)和抗體(解)之間的適應值函數,fitness=1/(J+1),J為目標函數;praise(xV)表示對局部極值或極值附件的點的獎勵函數[11]。
2.3 人工免疫算法描述
a)分析問題 對問題及其解的特征進行分析和了解,并設計解的合適表達形式,本文中的解采用二進制位串的形式表達(30位串)。
b)產生初始解群體 假設解群體中解個體數為n(100),記憶庫中解個體數為m(局部最小值設定為20)。隨機產生n個解個體并從記憶庫中提取m個個體構成初始群體。
c)對上述群體中各抗體進行評價 在本算法中對個體的評價是以個體的期望繁殖率為標準。具體計算過程如下:
(a) 計算抗體在群體中所占的比率:
RV=1/n∑nw=1KVW(15)
其中:KVW=1 SVW≥Tacl0 otherwise ,Tacl為一個預先確定的閾值[11]。
為了加快基于信息熵的免疫算法運行速度,抗體V和W的相似度計算公式可寫成[12]:SVW=M/[M+(M-ns)log 2]。其中,設抗體V和W有ns個相同的位。
(b)計算抗體V的期望繁殖率:eV=AV/RV。
d)形成新一代解群體。將上述群體按eV的降序排列,并取前n個個體構成新一代群體;同時取前m個個體存入記憶庫中。
e)判斷是否滿足結束條件。是,則結束;否,則繼續下一步操作。
f)解群體的產生。基于步驟d)的計算結構對抗體群體進行選擇、交叉、變異操作得到解群體,再從記憶庫中取出記憶的個體,共同構成解群體。
g)轉步驟c)。
2.4 基于人工免疫及最大互信息的圖像配準算法描述
a)采用3×3的模板對圖像做維納自適應預濾波處理。由于不同模態醫學圖像的成像原理不同,分辨率、信噪比等也各有差異,在配準之前必須對圖像進行預處理。通常功能圖像(如PET、SPECT)在采集時引入的噪聲比較大,需要對其進行去噪和增強,突出感興趣區域。
b)插值,實現多模態圖像的一致大小。由于解剖圖像與功能圖像的分辨率通常不同,解剖圖像分辨率較高,而功能圖像的分辨率較低,要進行多模醫學圖像的配準及融合,首先要求待處理圖像具有相同的圖像分辨率,為了保留高分辨率圖像的信息,必須對低分辨率圖像進行插值。
本文采用雙三次樣條插值方法[13]對PET圖像進行插值,提高了圖像分辨率。c)求出初始時,參考圖像與浮動圖像的互信息熵值。對所求問題編碼并初始化,各個體串長30。串的位號意義為:用9位表示旋轉角θ;9位表示tx,tx 為X方向上的平移,首位為正負;9位表示ty,ty為Y方向上的平移,首位表示正負;3位表示縮放比例因子k。令e為編碼的值,平移的解碼值d=-255+e;旋轉角度的解碼值d=(e/511)×360。
d)用上述2.1~2.3節的人工免疫算法求出最佳的仿射變換參數。
e)用d)的結果對圖像進行配準。
3 實驗結果
基于人工免疫算法及最大互信息的醫學圖像配準方法采用MATLAB編程語言進行軟件編寫,實驗平臺的機器配置為Pentium4 3 GHz CPU 1 GB內存,顯卡為GeForce 6800 GT 。
按照第2章所述的配準方法,對三組CT和PET醫學圖像進行配準,配準的結果如圖1、2所示。
圖1給出了具有平移失配的兩幅CT和PET圖像。使用基于人工免疫算法及最大互信息的配準方法得到配準參數為Δx=1,Δy=-10,Δθ=0,k=1。圖1(c)給出了相應參數集的配準結果,配準時間為20.7 min。
圖2給出了具有旋轉平移失配的兩幅CT和PET圖像。使用基于人工免疫算法及最大互信息的配準方法得到配準參數為Δx=0,Δy=0,Δθ=-10,k=1。圖2(c)給出了相應參數集的配準結果,配準時間為25.3 min。
4 醫學圖像配準效果的評價
4.1 醫學圖像配準的評價方法
醫學圖像配準,特別是多模醫學圖像配準結果的評估一直是一個難點。由于待配準的多幅圖像是在不同時間或條件下獲取的,不存在絕對的配準,也不存在金標準,只有相對的最優配準。配準算法的優劣可以從多個方面進行評估,如配準速度、魯棒性、配準精度等,對于醫學圖像配準而言,配準效果一直是首要的評價標準。常用的評估方法主要有體模、準標和目測法[2]。
本文提出了一種基于圖像統計特性的醫學圖像配準效果定量評價方法,該方法利用圖像的統計特性,包括均方根誤差(root mean square error RMS error)、相關系數(correlation coefficient,CC ),給出了配準算法的定量評價指標。通過實驗發現,統計特性法是一種客觀且實用的評價方法。
設有兩幅圖像I1、I2,圖像大小為M×N,其均方根誤差RMS定義如下:
RMS=(∑M-1i=0∑N-1j=0(I1(i,j)-I2(i,j))2)/(M×N)(16)
RMS越小,說明兩幅圖像之間的差異越小,配準效果就越好;評價配準效果的另一個統計特性相關系數CC定義如下:
CC=[∑M-1i=0∑N-1j=0(I2(i,j)-I2)(I1(i,j)-I1)]/
(∑M-1i=0∑N-1j=0(I2(i,j)-I2)2
∑M-1i=0∑N-1j=0(I1(i,j)-I1)2)(17)
其中:I1、I2為兩幅圖像的灰度均值,I1=(∑M-1i=0∑N-1j=0I1(i,j))/(M×N),I2=(∑M-1i=0∑N-1j=0I2(i,j))/(M×N)。CC的取值在0~1,當兩幅圖像之間不存在任何相關性時,取值為0;反之,如果兩幅圖像達到完全匹配時,CC趨向于1。這是一種非常理想的情況,在實際情況中,特別是多模圖像配準時,CC取值往往會比較小。
4.2 配準效果評價
為了驗證本文配準方法的配準精度,采用了如下評估方法:先對一幅PET圖像進行不同角度的旋轉和平移,然后用本文使用的配準方法進行配準;再計算兩幅圖像在配準前和配準后的均方根誤差和相關系數,如表1所示。
由表1可以看出,配準之前參考圖像與浮動圖像之間的均方根誤差均大于配準后的均方根誤差,而配準之前參考圖像與浮動圖像之間的相關系數均小于配準后的相關系數。可見,本文配準算法對多模醫學圖像配準均有效。
5 結束語
本文提出一種基于人工免疫算法及最大互信息的醫學圖像配準方法。在本算法中,運用了人工免疫網絡強大的免疫記憶能力和自適應調節能力的特點進行圖像配準,可以實現自動圖像配準,并不需要預先確定特征點。經過基于統計特性的客觀評價,表明圖像配準效果較好。
參考文獻:
[1] BETKE M HONG H KO J P. Landmark detection in the chest and registration of lung surfaces with an application to nodule registration[J]. Medical Image Analysis 20037(3):265-281.
[2]MAINTZ J B A VIGERGEVER M A. A survey of medical image registration[J]. Medical Image Analysis 19982(1) :1-36.
[3]DOUGHERTY L ASMUTH J C GEFTER W B. Alignment of CT lung volumes with an optical flow method[J].Academic Radiology 200310(3): 249-254.
[4]BLAFFERT T WIEMKER R. Comparison of different follow-up lung registration methods with and without segmentation[J]. Journal of Proceeding SPIE 20045370:1701-1708.
[5]STUDHOLME C HILL D HAWKES D J. An overlap invariant entropy measure of 3D medical images alignment[J]. Pattern Recognition 199932(1): 71-86.
[6]葛培明 陳虬 李堯臣. 混合遺傳算法在醫學圖像配準中的應用[J]. 中國生物醫學工程學報 200726(3):326-329.
[7]葉發茂 蘇林 李樹楷. 基于人工免疫網絡算法的圖像配準方法[J]. 計算機工程 200733(13):197-199.
[8]TIMMIS J. Artificial immune systems: a novel data analysis technique inspired by the immune network theory[D].Wales: Department of Computer Science University of Wales 2000.
[9]李彬,田聯房,毛宗源. 基于人工免疫的灰度圖像多閾值自動分割[J]. 計算機工程與設計 2007 28(1):106-108.
[10]陸益民,毛宗源,張波. 基于免疫算法的混沌多模型微擾控制[J]. 控制理論與應用 200421(1):30-34.
[11]葛紅,毛宗源. 免疫算法幾個參數的研究[J]. 華南理工大學學報:自然科學版 200230(12):15-18.
[12]鄭日榮,毛宗源. 一種改進的人工免疫算法[J]. 計算機工程與應用 200339(33):55-57.
[13]THOMAS M L CLAUDIA G KLAUS S. Survey:interpolation methods in medical image processing[J]. IEEE Trans on Medical Imaging 199918(11):1049-1075.