一種穩健的非剛性醫學影像配準方法

李 琦,姬紅兵,臧 博,劉 靳

(西安電子科技大學電子工程學院,陜西西安 710071)

一種穩健的非剛性醫學影像配準方法

李 琦,姬紅兵,臧 博,劉 靳

(西安電子科技大學電子工程學院,陜西西安 710071)

提出了一種基于局部線性嵌入(LLE)和改進的L-BFGS優化算法的非剛性配準新方法.該方法首先計算影像不同方向上的定序特征,用于補充傳統互信息測度中缺失的空間結構信息;然后,運用LLE方法及其逆映射對高維定序特征進行降維和融合;進而結合影像灰度信息構造了一種基于混合熵的配準測度,有效保證了配準測度函數的光滑性和收斂性;最后,采用改進的L-BFGS優化方法搜索最優配準參數.多組仿真數據的測試結果表明,在噪聲情況下,所提方法具有精度高、魯棒性強的特點,優于現有幾種方法.

非剛性配準;局部線性嵌入;定序特征

醫學影像配準是醫學影像處理領域的一個重要的和基本的研究課題,是醫學影像融合、醫學影像重建、影像與標準圖譜的匹配等研究的基礎[1-2].醫學影像配準就是尋求一種最優空間變換,使兩幅醫學影像的對應點在給定相似性測度下達到空間位置和解剖位置的一致.醫學影像的配準方法可以分為剛性配準和非剛性配準兩大類.經過多年的發展,剛性配準算法已經成熟,可以達到很高的配準精度.然而,醫學影像配準中剛性體的配準只是很小的一部分,許多重要的臨床應用需要非剛性變換來描述影像之間的空間關系.

另外,因為噪聲環境的復雜性,含噪醫學影像的配準也受到很大關注.在醫學影像的成像過程中,通常將成像必需的但對人體有害的強磁場、示蹤劑和放射線等控制在一定范圍內,加上成像模式本身的一些物理限制,往往導致影像模糊,并伴有噪聲甚至偽影的出現.因此,在醫學影像配準方法研究中,分析相似性測度對噪聲的穩健性也是研究配準性能的一個重要方面.

當前,非剛性配準方法主要分為基于特征的方法和基于灰度的方法[1-4].基于特征的配準方法一般算法效率較高,但算法精度通常取決于特征提取的精度[3].基于灰度的配準度量影像灰度間的信息分享程度,具有較高的精度和魯棒性,但計算量較大[4].在此類方法中,互信息由于不需要獲取待配準影像的先驗知識、無需預處理及人工干預,適用性強,得到了高度關注和廣泛應用.然而,由于Shannon熵假設相鄰像素灰度值不相關,互信息忽略了像素灰度分布的空間位置關系,導致了相似性測度存在過多局部極值、對噪聲敏感等問題,最終影響了配準精度.當前的研究熱點是建立將影像空間結構特征與互信息相結合的混合測度[5-9].表示空間結構信息的特征描述符很多,例如影像梯度、尺度不變特征變換(SIFT)、灰度共現矩陣(GLCM)以及定序特征等.Pluim等[6]將一個梯度相似性系數與歸一化互信息相乘作為最終的配準測度,但影像梯度對噪聲十分敏感,無法有效增強配準相似性測度的抗噪能力.SIFT特征對影像旋轉、尺度縮放、亮度變化等保持不變性,但在其算法的關鍵步驟,利用關鍵點鄰域像素的梯度方向分布特性,為每個關鍵點賦予方向;同時, SIFT丟棄低對比度的關鍵點以致對邊緣模糊的目標無法準確提取特征點.因此,SIFT不適用于低分辨率的醫學影像的配準.Rueckert等[7]提出了二階互信息的概念,將相鄰像素對的GLCM作為概率密度函數,擴展了邊緣熵及聯合熵的定義.然而,由于GLCM是像素距離和角度的矩陣函數,因此完整的GLCM計算,其參數的選取范圍很廣,計算量很大.文獻[9]將影像定序特征引入到醫學影像配準算法中,作為一種新的圖像特征,定序特征包含了豐富的紋理、邊緣等空間結構信息,具有很強的抗噪能力和通用性,同時不易受到影像亮度、對比度等因素的影響[10-11],能夠有效補充互信息測度的不足.

醫學影像多方向的定序特征是高維的,包含了很多不重要的冗余信息,文獻[9]中采用主成分分析(PCA)方法進行降維處理.但是由于高維定序特征本質上是非線性的,PCA方法很難發掘出其中非線性分布的幾何結構和相關性.局部線性嵌入(LLE)是一種無監督的非線性學習方法,以保存數據局部鄰域間相互關系的方式把高維數據映射到一個低維全局坐標系下,能夠揭示數據的全局非線性結構[12-13].因此,筆者采用LLE及其逆映射對多方向的定序特征進行降維和融合,并在此基礎上結合影像灰度信息構造了一種基于混合熵的配準測度[9],利用定序特征對噪聲的強魯棒性有效保證了配準測度函數的光滑性和收斂性,大大提高了配準精度.

1 空間結構特征

使用一組和差定序濾波器[11](相位分別為θ=0,π/12,π/6,π/4,π/3,5π/12)對一幅T2-MR影像進行6個方向濾波,如圖1所示.圖1(a)為待濾波影像,圖1(b)～(g)分別為6個方向的濾波結果.可以看到,一幅醫學影像多個方向上的定序特征所包含的信息不盡相同,但也存在極強的相關性.為了最小化不同方向的定序特征中所包含的冗余信息,同時降低配準算法的計算代價,文獻[9]采用PCA方法進行降維處理.但是由于高維定序特征是影像多方向的濾波結果,因此本質上是非線性的,而PCA方法作為一種線性的數據降維方法并不能準確揭示出非線性數據的內在幾何結構和相關性.若采用PCA方法處理非線性數據,很有可能會丟失一些重要的信息,從而影響最終的配準精度.因此,文中采用更為適合的LLE方法將高維觀測空間中的多方向定序特征映射為低維嵌入空間中的對應點集.接下來,如何通過這些低維嵌入空間中的對應點集獲得高維觀測空間中醫學影像的空間結構特征是一個關鍵問題.眾所周知,樣本的聚類中心通常最具代表性,因此與低維嵌入空間中的中心點對應的高維觀測空間中的影像也應該最具代表性,即能夠描述醫學影像整體空間結構特征的多方向定序特征的融合.文獻[13]在LLE方法的基礎上提出一種從低維嵌入空間向高維空間映射的方法,并在多姿態人臉圖像的重構實驗中得到有效驗證.通過這一逆向映射方法,與低維嵌入空間中的中心點對應的高維觀測空間中的影像得以重構生成.具體步驟如下:

(1)降維.設對一幅大小為M×N的醫學影像I(x,y)進行n個方向的和差定序濾波,得到的n個方向的影像定序特征Y(x,y,θi)(i=1,2,…,n),大小也為M×N.將這些影像定序特征作為n個樣本,分別按列連接成M×N維的n個影像向量Yi∈RM×N(i=1,2,…,n).參數k為最近鄰數,d為低維嵌入空間的維數.

圖1 影像定序特征的提取及融合

定義代價函數為

固定權值W,由式(3)計算d維嵌入矢量{Xi}i=1,2,…,n∈Rd,即

(2)計算低維嵌入空間中Xi(i=1,2,…,n)的中心為

(3)融合.用其k鄰域的加權和表示Xcenter,優化代價函數

得到權值矩陣

在LLE方法中,低維空間中的局部鄰域關系在高維空間中保持不變.因此,與低維嵌入空間中的Xcenter對應的高維空間中的影像YF,可由式(6)得到的低維空間中的權值W重構,即

將YF按列還原為M×N的影像,該影像即為I(x,y)融合后的定序特征,如圖1(i)所示.為了便于比較,文獻[9]通過PCA方法得到的定序特征在圖1(h)給出.可以明顯看出,通過LLE方法得到的定序特征比PCA方法得到的定序特征包含更清晰、更豐富的空間信息.采用均值、標準差等圖像融合的指標客觀評價兩種方法,比較結果如表1所示.除了熵的值幾乎相等,在其他指標上均是LLE方法明顯優于PCA方法的.因此, LLE方法比PCA方法更能揭示影像多方向定序特征的非線性分布,由文中方法得到的定序特征中,冗余信息被最小化,互補信息得以保留并加強,同時不易受到噪聲的影響,魯棒性強.

將待配準的兩幅醫學影像中分辨率較高的一幅影像設為參考影像R,其表示影像整體空間結構的定序特征設為X,另外一幅影像設為浮動影像F.這里采用一種基于混合熵的配準相似性測度[9],定義為

其中,H(F)為F的邊緣熵;H(R,X)為R和X的二階聯合熵;H(R,F,X)表示R、F和X的三階聯合熵.該測度將R和X的聯合(R,X)作為一個整體與F進行配準,配準過程視為使得(R,X)和F之間互信息最大的參數求解過程.在基于混合熵的配準測度中,影像定序特征的引入不僅增加了互信息中缺失的空間信息,使得配準測度能夠同時衡量待配準影像對應灰度信息和空間信息的一致性,而且能夠有效抑制噪聲的影響,使配準測度函數更加平滑,有效避免了在配準搜索過程中陷入局部極值,提高了配準的精度.

表1 PCA方法和LLE方法的性能比較

2 配準方法

任何非剛性配準方法都可以描述為3個部分:聯系浮動影像與參考影像的空間變換,測量浮動影像與參考影像相似程度的相似性測度,確定相似性測度中最優變換參數的優化方法.文中配準方法的相似性測度在上節中已給出相應的數學描述,以下對文中方法的空間變換與優化方法進行介紹.

2.1 空間變換模型

在醫學影像非剛性配準算法中,為適合各種復雜的組織形變建立合理的空間變換模型是非常活躍的研究領域[1-2].一直以來,研究人員提出了各種方法.基于基函數的方法使用基函數的線性組合來描述變形域,然而它補償解剖形變的能力非常有限,通常用于全局粗配準.基于薄板樣條的方法是彈性形變中彈性薄片應力最小的一種插值模型,但是因為每一個控制點對變換都具有全局影響,所以很難模擬局部變形,不適合具有局部幾何差異的圖像配準.同薄板樣條相比,B樣條可以控制局部變形,改變控制點只影響它附近局部鄰域的形狀改變,其中基于B樣條的自由形變模型(FFD)是目前使用較多的一種樣條配準方法[14-15],其基本思想是:變形操作是作用于物體所嵌入的變形空間,而不是直接作用于物體;如果變形空間被改變了,則嵌入其中的物體自然隨之改變.基于B樣條的FFD模型是通過操縱分布在圖像區域上的控制點網格來完成非線性變形的,以下是對該模型的數學描述.

對于二維情形,定義平面上的矩形域為

令Φ表示覆蓋在矩形域Ω上由nx×ny個控制點φi,j所組成的網格,每個控制點之間的間距為δ,則該位移場定義為

其中,0≤u≤1.

B樣條具有局部特性,即控制點位置的改變僅影響該控制點周圍鄰域網格內像素點的位移.因此,更新某個控制點的位置時,僅需更新該控制點鄰域網格內的點;進行相似度計算時,也只需計算兩幅影像對應控制點鄰域的相似度值,因此算法的計算量大大減小.在文中配準方法中,首先對待配準影像進行剛性全局粗配準,再對發生局部形變的感興趣區域(ROI)用基于B樣條的FFD模型進行非線性空間變換.

2.2 改進的L-BFGS優化方法

在確定了空間變換類型和相似性測度函數之后,配準過程實質上就是一個多參數的優化求解過程,目的是尋找最佳空間變換參數,使得配準相似性測度達到最大值.所以優化方法的選擇十分重要,快速有效的優化算法可以大大節省運行時間,為實時圖像處理提供可能性.

L-BFGS優化方法是在BFGS算法的基礎上發展起來的一種算法,非常適合求解高維度和大規模的無約束優化問題,它并不指定Hessian矩陣的形式和存儲,直接計算搜索方向,大大降低了對內存的需求,提高了計算速度[16].在L-BFGS算法中,位置迭代公式為

其中,αk為步長因子,dk為搜索方向.

但在實際的非剛性醫學影像配準過程中,L-BFGS優化方法存在以下缺陷:

(1)當待配準的兩幅醫學影像中存在較大非線性形變、局部形變不連續、灰度突變或者存在噪聲等情況時,配準相似性測度會出現異常波動、起伏劇烈、局部極值增多、光滑性和收斂性惡化,從而導致優化算法在開始時,位置迭代過程中αk數值很大,最終很有可能會錯過最優空間變換參數.

(2)在優化過程中,L-BFGS優化方法無法擺脫局部極值的吸引,因此有必要在保留當前最優位置的同時引入變異算子.但是,若將dk全部重新初始化,將會完全破壞當前搜索,使收斂速度大大減緩.

針對上述情況,將L-BFGS優化方法進行如下改進:

(1)在位置迭代更新前,對αk設定閾值進行限制,防止誤配情況的發生.

(2)若當前位置xk變化極小時,按隨機幾率將dk中少部分維數重新初始化,以此增強全局搜索能力,使算法不易陷入局部極值,提高搜索精度.

最終,滿足改進的L-BFGS優化方法迭代終止條件的空間變換參數輸出為最優配準參數,完成配準.

3 實驗結果與分析

為了對比分析文中方法的有效性,選取相同算法框架的互信息測度[4]、互信息與影像梯度結合測度[6]、二階互信息測度[8]和文中式(8)定義的測度(以下簡稱為測度1～測度4)對醫學影像進行配準實驗.實驗環境為Pentium(R)Dual-Core CUP 2.50 GHz和RAM 2 GB的PC機,操作系統為Windows XP,實驗用Matlab7.6編程實現.

3.1 測試不同優化方法的搜索性能

將同一病人腦部的不同切片分別作為浮動影像和參考影像,大小均為256×256,如圖2所示.對測度2分別采用L-BFGS方法和文中改進的L-BFGS方法進行最優配準參數搜索,對比結果如圖3所示.可以看出,由于測度2光滑性較差,L-BFGS優化方法的位置迭代更新步長過大,并且最終陷于局部極值,導致配準結果失真;而改進的L-BFGS優化方法,由于能及時判斷是否已收斂于局部最優點,并通過變異操作迅速擺脫它的束縛,提高了搜索成功率.實驗表明,該改進算法是非常有效的.

圖2 配準數據1

3.2 測試不同配準測度的精度和魯棒性

為驗證4種配準測度對噪聲的魯棒性,在圖2兩幅影像中均加入均值為0、方差為0.01的高斯噪聲,采用改進的L-BFGS優化方法進行最優配準參數搜索,配準結果如圖4所示.

因受噪聲影響,直觀上并不能明顯看出4種測度配準結果的優劣.為客觀評價4種測度的配準精度,采用灰度差平方和(SSD)和相關系數(CC)作為評價標準[14],SSD數值越大,表明配準后影像與參考影像差值越大,配準精度低;CC值越大,表明配準后影像與參考影像相似度越高,配準精度高.4種配準測度的精度比較結果如表2所示.可以看出,在相同噪聲條件下,測度2配準效果最差,測度1和測度3相差不大,文中測度始終保持較小的誤差,配準精度高于其他測度.

表2 N(0,0.01)高斯噪聲條件下4種配準測度的精度對比

圖3 兩種優化方法的實驗結果比較

圖4 N(0,0.01)高斯噪聲條件下4種配準測度的實驗結果比較

3.3 測試不同配準測度的運行時間

實驗采用一對CT/MR影像,如圖5所示,大小均為256×256.以CT影像作為參考影像,MR影像作為浮動影像.4種配準測度單次運行時間比較結果見表3.可以看到,測度1用時最少,其他改進測度的運算時間均有所增加,其中測度1、測度2定義中只包含二維聯合概率分布,節約了較多的內存,運行時間較短;測度3定義中包含四維聯合概率分布,運行時間最長;文中測度定義中雖然只包含三維聯合概率分布,但LLE算法的正向與反向映射還是消耗了一些運算時間,故文中測度算法效率僅優于測度3,仍需進一步提高.

圖5 配準數據2

表3 4種配準測度的運行時間對比

4 結束語

提出了一種基于LLE和改進的L-BFGS優化的非剛性醫學影像配準方法,運用基于B樣條的FFD模型進行空間變換;針對高維定序特征,采用LLE及其逆映射進行降維和融合處理;得到的醫學影像整體定序特征與影像灰度構成特征空間,基于混合熵的配準測度通過衡量兩幅影像對應灰度信息和空間信息的一致性,使配準測度函數更加平滑,減少了配準過程中可能陷入的局部極值;在最優配準參數的優化搜索過程中,改進的L-BFGS優化方法采用變異操作增強了跳出局部最優解的能力.實驗表明,文中方法能有效抑制噪聲影響,配準精度較高.后續作者將致力于有效降低算法執行時間,使其能夠滿足臨床的需要.

[1]Sotiras A,Daratzikos C,Paragios N.Deformable Medical Image Registration:a Survey[J].IEEE Transactions on Medical Imaging,2013,32(7):1153-1190.

[2]Liao S,Chung A C S.Nonrigid Brain MR Image Registration Using Uniform Spherical Region Descriptor[J].IEEE Transactions on Image Processing,2012,21(1):157-169.

[3]Wang Shanhu,You Hongjian,Fu Kun.BFSIFT:a Novel Method to Find Feature Matches for SAR Image Registration [J].IEEE Geoscience and Remote Sensing Letters,2012,9(4):649-653.

[4]Pluim J P W,Maintz J B A,Viergever M A.Mutual-information-based Registration of Medical Images:a Survey[J]. IEEE Transactions on Medical Imaging,2003,22(8):986-1004.

[5]So R W K,Tang T W H,Chung A C S.Non-rigid Image Registration of Brain Magnetic Resonance Images Using Graph-cuts[J].Pattern Recognition,2011,44(10-11):2450-2467.

[6]Pluim J P W,Maintz J B A,Viergever M A.Image Registration by Maximization of Combined Mutual Information and Gradient Information[J].IEEE Transactions on Medical Imaging,2000,19(8):809-814.

[7]Rueckert D,Clarkson M J,Hill D L G,et al.Non-rigid Registration Using Higher-order Mutual Information[C]// Proceedings of SPIE—the International Society for Optical Engineering.Bellingham:SPIE,2000:438-447.

[8]Chen B J,Li J L,Chen G.Study of Medical Image Registration Based on Second-Order Mutual Information.[C]// Proceedings of the 1st International Conference on Bioinformatics and Biomedical Engineering.Piscataway:IEEE,2007: 956-959.

[9]李琦,姬紅兵,同鳴.一種多模態醫學影像魯棒配準方法[J].西安電子科技大學學報,2011,38(6):30-36.

Li Qi,Ji Hongbing,Tong Ming.Robust Registration of Multimodality Medical Images Based on the Principal Ordinal Feature and Hybrid Entropy[J].Journal of Xidian University,2011,38(6):30-36.

[10]Tan T N,Sun Z N.Ordinal Representations for Biometrics Recognition[C]//Proceedings of 15th European Signal Processing Conference.Poland:EURASIP,2007:35-39.

[11]Han Y F,Tan T N,Sun Z N,et al.Embedded Palmprint Recognition System on Mobile Devices.[C]//Lecture Notes in Computer Science:4642.Heidelberg:Springer Verlag,2007:1184-1193.

[12]Roweis S T,Saul L K.Nonlinear Dimensionality Reduction by Locally Linear Embedding[J].Science,2000,290 (5500):2323-2326.

[13]Zhang C S,Wang J,Zhao N Y,et al.Reconstruction and Analysis of Multi-pose Face Images Based on Nonlinear Dimensionality Reduction[J].Pattern Recognition,2004,37(2):325-336.

[14]Rueckert D,Sonoda L I,Hayes C,et al.Non-rigid Registration Using Free-form Deformations:Application to Breast MR Images[J].IEEE Transactions on Medical Imaging,1999,18(8):712-721.

[15]Díez Y,Oliver A,LladóX,et al.Revisiting Intensity-based Image Registration Applied to Mammography[J].IEEE Transactions on Information Technology in Biomedicine,2011,15(5):716-725.

[16]Morales J L,Nocedal J.Remark on“Algorithm 778:L-BFGS-B:Fortran Subroutines for Large-scale Bound Constrained Optimization"[J].ACM Transactions on Mathematical Software,2011,38(1):7.

(編輯:齊淑娟)

Non-rigid registration of medical images based on local linear embedding and improved L-BFGS optimization

LI Qi,JI Hongbing,ZANG Bo,LIU Jin
(School of Electronic Engineering,Xidian Univ.,Xi’an 710071,China)

Non-rigid registration of medical images has become a challenging task in medical image processing and applications.In this paper,we propose a local linear embedding(LLE)and improved LBFGS(limited-memory Broyden Fletcher Goldfarb Shanno)optimization based registration method.With abundant spatial information and good stability in noisy environment,the ordinal features are computed on different orientations to represent spatial information in medical images.For high dimensional ordinal features,the LLE algorithm is used for dimensionality reduction and the inverse mapping of LLE is used to fuse complementary information together.Then a hybrid entropy based similarity measure which integrates image intensity with ordinal feature is chosen as the registration function.Finally an improved L-BFGS algorithm is used to search for the optimal registration parameters.We evaluate the effectiveness of the proposed approach by applying it to the simulated brain image data.Experimental results show that the proposed registration algorithm is less sensitive to noise in images.Compared with some traditional methods,the proposed algorithm is of higher precision and better robustness.

non-rigid registration;local linear embedding;ordinal feature

TP391.4

A

1001-2400(2014)05-0054-07

2013-05-13< class="emphasis_bold">網絡出版時間:

時間:2014-01-12

國家自然科學基金資助項目(61101246);中央高校基本科研業務費專項資金資助項目(JB140209,72125748)

李 琦(1979-),女,講師,西安電子科技大學博士研究生,E-mail:qili@xidian.edu.cn.

http://www.cnki.net/kcms/doi/10.3969/j.issn.1001-2400.2014.05.010.html

10.3969/j.issn.1001-2400.2014.05.010