具有大形變特征的顱腦CT圖像的非剛性配準

徐 峰劉 偉李傳富，2馮煥清*

1(中國科學技術大學電子科技系，合肥 230027)2(安徽中醫(yī)學院第一附屬醫(yī)院影像中心，合肥 230031)

具有大形變特征的顱腦CT圖像的非剛性配準

徐 峰1劉 偉1李傳富1，2馮煥清1*

1(中國科學技術大學電子科技系，合肥 230027)2(安徽中醫(yī)學院第一附屬醫(yī)院影像中心，合肥 230031)

Demons算法是一種基于光流場模型的小形變非剛性配準算法，大形變情況下不具有拓撲保持性，將它用于顱腦CT圖像配準時效果不理想。為此，本研究對它進行了改進。首先建立Demons算法目標能量函數(shù)，將形變場求解轉化為目標函數(shù)優(yōu)化問題;然后通過增加sKL距離作為正則項來優(yōu)化目標函數(shù)，消除了形變場的不適定性，并使形變場更加光滑。對高分辨率顱腦CT圖像的實驗結果表明，改進算法不僅能夠處理大形變問題，還能在處理大形變時通過光滑的形變場得到更精確的配準結果。

醫(yī)學圖像配準;Demons算法;拓撲保持性;Jacobian行列式;sKL距離

Abstract:Demons is a non-rigid image registration algorithm which is derived by assuming small deformations.One of the limitations of the original Demons is that it can not produce topology preserving maps for the large deformations.Aiming to solve this problem，an improved Demons algorithm was proposed in this paper.First，the equation of force in the original Demons was regarded as the result of minimizing the energy function.Then，Demons algorithm was improved by adding a regularization term into the function.The symmetric Kullback-Leibler(sKL)distance in information theory was used as the regularization term.The experiment results with high resolution CT cerebral images demonstrated that the improved algorithm could not only handle large deformations，but also obtain more accurate registration results using smooth deformation fields.

Key words:medicalimage registration;Demonsalgorithm;topology preservation;Jacobian matrix;sKL distance

引言

隨著醫(yī)學影像輔助診斷在臨床上的廣泛應用，各種醫(yī)學圖像處理技術都得到了深入研究，其中圖像配準技術是研究熱點之一。圖像配準的本質(zhì)是尋找一種幾何變換使得待分析圖像和參考圖像達到空間位置上的對齊，已應用于圖像分割［1］、圖像融合［2］、病變檢出［3］等諸多方面。醫(yī)學圖像配準方法可分為剛性配準和非剛性配準兩大類，其中剛性配準屬于全局配準的方法，一般只是對圖像的旋轉、平移、縮放等幾個參數(shù)進行仿射變換，達到兩幅圖像在大體位置上的對齊，常用于相同個體不同次檢查或不同模式檢查數(shù)據(jù)之間的配準。而非剛性配準主要應用于異體配準 (inter-subject registration)，由于異體配準在統(tǒng)計圖譜的建立、病變檢出等方面具有重要的意義，因此非剛性配準方法是醫(yī)學圖像配準研究的重點。

非剛性配準算法可以分為兩類:模型驅動算法和灰度驅動算法。Demons算法是一種灰度驅動的算法［4］。Hellier等對常用的六種配準算法進行的回顧性對比研究表明，在多種配準精確度評價指標中Demons算法優(yōu)于其他五種算法［5］。但是 Demons算法只適用于小形變的情況，在出現(xiàn)大形變時Demons算法將不能保持圖像的拓撲結構，并產(chǎn)生物理上不合理的形變。然而，在異體配準過程中，由個體之間的生理差異以及占位性病變所導致的大形變問題是不可避免的。因此本研究針對此缺陷改進Demons算法，使其在大形變下具有拓撲保持性，以應用于具有大形變特征的顱腦CT圖像的配準。

1 Demons算法原理

Demons非剛性配準算法是 Jean-Philipe Thirion提出的［4］，基本思想源于19世紀 Maxwell為了解決熱動力學難題而提出的一種實驗假設，其假設在參考圖像S(x)的目標輪廓上分布著具有選擇性的驅動點，由它們驅動待配準圖像M(x)上的像素來實現(xiàn)兩幅圖像的配準。驅動力的形式取自光流場理論方程，而光流場理論的前提是假設圖像運動過程中灰度保持不變，即為計算驅動力，將式(1)微分化

進一步簡化式(2)得

其中，x=(x1，x2)，S(x)為S(x)的梯度向量，fD為驅動力。從而得到作用力

但是，在參考圖像比較均勻一致的區(qū)域ΔS(x)→0，式(4)非常不穩(wěn)定。為解決此問題，對式(4)進行調(diào)整，在其分母上增加一個分量，得

2 Demons算法的優(yōu)化

2.1 配準能量函數(shù)構造

圖像配準過程就是尋找一個位移場u(x)，使待配準圖像形變后M(x-u)與參考圖像S(x)在某種含義上相似。最通用的相似性測度(能量函數(shù))是L2范數(shù):

只需最小化式(6)，利用梯度下降法求解相應的Euler-Lagrange 方程得到

其中，

式(8)即為驅動待配準圖像的力場。M|x-u為形變后圖像的梯度，它決定了力的方向;M(x-u)-S(x)為兩幅圖像的灰度差。

這種優(yōu)化的問題在于其不適定性，尤其是形變場u的不唯一，因此需要對u進行正則化。一般是通過在式(6)中增加正則項來解決［6］，則能量函數(shù)變?yōu)?

可以看出式(8)與式(5)很相似，因此本研究假設式(5)是最小化相似函數(shù) ED(M，S，u)得到的驅動力場，由上述知需要對解進行正則化。而Thirion原來的Demons只是通過簡單的高斯平滑來正則化，這正是 Demons算法局限性的原因所在。像式(9)那樣引入正則項來改進Demons算法，所采用的正則項是配準變換的Jacobian行列式分布與恒等變換之間的sKL距離。

2.2 正則項

令待配準圖像M(x)到參考圖像S(x)的變換為h:M→S，則位移場u(x)與h(x)滿足h(x)=x-u(x)。只要變換h滿足兩個條件就可保持拓撲結構［7］:①變換h的Jacobian行列式為正值;②變換h是一一對應即雙射變換。

設Dh為h的Jacobian矩陣相應的 Jacobian行列式為令，則 2D圖像變換的Jacobian行列式為

對于二維空間坐標變換，Jacobian行列式的物理含義為說明面積微元膨脹1則收縮面積不變?yōu)槠纥c。當面積微元ds(x)處 Jacobian行列式都等于1時，則變換前后微元ds(x)的面積沒有變化，即變換前后是恒等映射，是一一對應的。因此，若正向變換和逆向變換各點的Jacobian行列式分布足夠接近恒等映射，就可以滿足局部坐標系的拓撲保持性。采用信息論中 symmetric Kullback-Leibler(sKL)距離［8］來衡量這種接近程度。

設X與Y是兩概率密度函數(shù)，則KL距離為sKL距離定義為

令 h-1為h的逆變換，則滿足

對式(13)兩邊求導并設y=h(x)得所以

由式(10)與式(15)可以推導得

由此可以得到

因此，最小化式(18)就使得正向變換與逆向變換的Jacobian行列式分布最接近恒等映射。

這樣，由式(9)與式(18)就可得到新的能量函數(shù)

令 u=(u1，u2)及函數(shù)L(J)=(J-1)log(J)，則驅動力的公式更新為

3 算法實現(xiàn)

3.1 雙向多分辨率配準策略

為了克服Demons算法因其完全自由變形特點可能導致的空間結構不連續(xù)，并增加該算法的魯棒性，在實現(xiàn)改進算法時采用了雙向配準策略［9］。即以M(x)為待配準圖像，S(x)為參考圖像得到正向的形變矩陣Tf;再以S(x)為待配準圖像，M(x)為參考圖像得到逆向的形變矩陣 Tr。最后 M(x)向S(x)配準的形變矩陣T為Tr與Tf的均值，由于Tr與Tf的方向相反，得出T=(Tf-Tr)/2。

多分辨率配準策略就是從粗到精的分級配準。較低分辨率的圖像更能反映圖像的輪廓特征并且更加穩(wěn)定，高分辨率圖像配準是在低分辨率配準結果的基礎上進行的，因此采用多分辨率配準策略不僅可以顯著提高算法的運算速度，而且能夠達到更精確的配準效果。

3.2 形變場疊加

在配準過程中，每次迭代都是以上次迭代的結果為待配準圖像，即 Mi+1=Ti°Mi，其中 Mi為第 i次迭代的待配準圖像，Ti為第i次迭代的形變場，Mi+1為第i次迭代的配準結果，同時也是第i+1次的待配準圖像。所以要進行多次Ti°Mi插值操作，這會導致配準后圖像細節(jié)丟失和模糊。

本研究通過將形變場疊加來解決這一問題，即將多次形變場Ti合并為T，再進行一次T°M0操作得到配準結果，具體的疊加方法如圖1所示。

圖1 形變場疊加Fig.1 Deformation fields stacking

待配準圖像Mi經(jīng)兩次形變得到配準結果Mi+2=Ti+1°(Ti°Mi)。所以只需要知道 Ti(p)和Ti+1(p)，就可以得到從Mi中p位置到Mi+2中p位置的形變T(p)=Ti(p)+Ti+1(p)。而Ti只是已知整數(shù)坐標點的形變，因此需要通過插值來得到非整數(shù)坐標點的形變，采用雙線性插值方法來計算圖1中的形變 Ti+1(p)。

3.3 實現(xiàn)步驟

改進算法的具體步驟如下:

初始化:u(x，0)=0;

步驟1:按式(20)計算出待配準圖像M(x)每一個 Demons點的力 f(x，u(x，t));

步驟2:根據(jù)力 f(x，u(x，t))和 Demons算法的參數(shù)項力衰減系數(shù)α，計算出形變的位移場u(x，t)= α·f(x，u(x，t));

步驟3:將位移場 u(x，t)與高斯核 Gσ卷積u(x，t)=Gσ*u(x，t)，將位移場 u(x，t)應用于待配準圖像M(x)使其形變，得到M(x-u);

步驟4:若式(19)能量函數(shù)的值與前一次迭代的值之間的差小于某閾值則停止，否則返回步驟1。

4 實驗結果及討論

采用512像素×512像素的高分辨率顱腦 CT圖像來比較 Demons算法與本改進算法的配準性能。Demons算法是文獻［9］優(yōu)化后的 Demons算法，并采用該文獻提供的參數(shù)。在實現(xiàn)改進算法時采用了雙向配準和多分辨率策略。圖2中(a)和(b)分別為待配準圖像與參考圖像，兩幅圖像來源于不同個體的同一層片，且已去除了顱骨和支架。圖2中(c)和(d)分別為Demons算法與改進算法的配準結果，均進行了3次迭代，其中改進算法式(20)中的參數(shù)λ取λ=1.2。

圖2 CT圖像配準。(a)待配準圖像;(b)參考圖像;(c)Demons配準結果;(d)改進算法配準結果Fig.2 CT image registration.(a)template image;(b)reference image;(c)registration result of Demons;(d)registration result of proposed algorithm

圖3 配準結果與參考圖像的差值圖像。(a)Demons算法;(b)改進算法Fig.3 Differencebetween registration resultand reference image.(a)Demons;(b)proposed algorithm

由圖2中(c)和(d)可以看到兩種方法都得到了精確的配準結果，從圖3的差值圖像也可以看出兩種方法的結果都比較精確。圖4是將兩種方法的形變場施加于網(wǎng)格得到的結果，通過對比圖4的(a)和(b)中圓圈標記處的網(wǎng)格，可明顯看出改進算法的網(wǎng)格更加規(guī)則，這說明改進算法有更加平滑的形變;而Demons算法的網(wǎng)格在某些區(qū)域出現(xiàn)了嚴重的彎曲，甚至幾乎交叉折疊。這是因為在這些區(qū)域有Jacobian行列式為負的情況，表明拓撲結構發(fā)生了變化?？梢?，由于獲得了更加平滑的形變場，改進算法得到的結果更精確。也就是說，在大形變的情況下，該算法能更好地保持圖像的拓撲結構。

圖4 形變場網(wǎng)格對比圖。(a)Demons算法;(b)改進算法Fig.4 Deformed grids Comparison.(a)Demons;(b)proposed algorithm

圖5為兩種方法在3次迭代中顱腦區(qū)域內(nèi)Jacobian行列式的統(tǒng)計直方圖，其中圖5(a)、(c)、(e)為Demons算法的結果，圖 5(b)、(d)、(f)為改進算法的結果。很明顯，本算法的Jacobian行列式統(tǒng)計分布更加尖銳，Jacobian行列式更集中于1附近，這也說明改進算法具有更加平滑的形變。從圖6(圖5(c)的放大圖)中的標記處可以看到，Demons算法有少量的Jacobian行列式為負值，這正是圖4(a)中Demons網(wǎng)格出現(xiàn)嚴重彎曲的原因。這說明了Demons算法在大形變的情況下改變了拓撲結構，而改進算法得到了很好的結果。圖7是兩種方法第一次迭代后形變的拓撲改變程度圖(其他迭代的結果類似)，可以看到本算法的結果明顯地優(yōu)于原Demons算法的結果。

圖5 顱腦區(qū)域內(nèi)Jacobian行列式統(tǒng)計直方圖。(a)和(b)Demons和改進算法第1次迭代;(c)和(d)Demons和改進算法第2次迭代;(e)和(f)Demons和改進算法第3次迭代Fig.5 Histograms of Jacobian values of deformations inside brain.(a)and(b)the first iteration of Demons and the proposed algorithm;(c)and(d)the second iteration of Demons and the proposed algorithm;(e)and(f)the third iteration of Demons and the proposed algorithm

圖8(a)是一幅具有占位性病變的CT圖像，圖8(b)是圖8(a)(待配準配準圖像)向圖2(b)(參考圖像)配準的結果，圖9(a)(b)分別是兩種方法第一次迭代中顱腦區(qū)域內(nèi)Jacobian行列式統(tǒng)計直方圖(另外兩次迭代結果與之相似)，從圖9(a)中可看出原 Demons算法仍存在Jacobian行列式為負值的情況，而改進算法仍然得到很好的結果，這說明改進算法適用于由占位性病變引起的大形變問題。

圖6 圖5(c)的放大圖Fig.6 Larger image of Fig.5(c)

圖7 拓撲改變程度圖。(a)Demons算法;(b)改進算法Fig.7 Topology change images.(a)Demons;(b)the proposed algorithm

圖8 占位性病變圖配準結果。(a)待配準圖像;(b)配準結果Fig.8 CT image with occupying lesionsand result of registration.(a) template image;(b) result of registration

圖9 顱腦區(qū)域內(nèi)Jacobian行列式統(tǒng)計直方圖。(a)Demons算法;(b)改進算法Fig.9 Histograms of Jacobian values of deformations inside brain.(a)Demons;(b)proposed algorithm

5 結論

Demons算法是在假設小形變的前提下由光流場原理推導出的方程，因此在處理不適用于Demons算法的大形變問題時，拓撲結構不能得以保持，故本研究對 Demons算法進行了改進。先將 Demons算法中的驅動力方程看作是極值化能量函數(shù)的結果，再通過引入Jacobian行列式統(tǒng)計分布與恒等變換之間的sKL距離來更新能量函數(shù)，得到新的驅動力方程，從而實現(xiàn)了改進。改進算法實現(xiàn)時采用了雙向配準和多分辨率策略。用高分辨率顱腦CT圖像數(shù)據(jù)的實驗結果表明，改進算法在處理大形變時具有良好的拓撲保持性，并得到了比Demons更精確的配準結果，將被應用到基于HRCT圖像的腦部病變計算機輔助診斷系統(tǒng)設計中。

［1］Pitiot A，Delingette H，Thompson PM，et al.Expert knowledgeguided segmentation system for brain MRI［J］.NeuroImage，2004，23(1):S85-S96.

［2］Pietrzyk U，Herholz K，Schuster A，et al.Clinical applications of registration and fusion of multimodality brain images from PET，SPECT，CT，andMRI［J］.EuropeanJournalof Radiology，1996，21(3):174-182.

［3］Lee JM，Kim SH，Jang DP，et al.Deformable model with surface registration for hippocampal shape deformity analysis in schizophrenia［J］.NeuroImage，2004，22(2):831-840.

［4］Thirion JP.Non-rigid matching using demons［A］.In:Computer Vision andPatternRecognition(CVPR '96)［C］.San Francisco:IEEE，1996.245-251.

［5］Hellier P，Barillot C，Corouge I，et al.Retrospective evaluation of intersubject brain registration.Medical Imaging［J］.IEEE Transactions on Medical Imaging，2003，22(9):1120 –1130.

［6］Vercauteren T，Pennec X，Perchant A，et al.Non-parameteric Diffeomorphic Image Registration with the Demons Algorithm［J］.MICCAI2007，Part II.LNCS，2007，4792:319 – 326.

［7］Holden M.A review of geometric transformations for nonrigid body registration［J］.IEEE Transactions on Medical Imaging，2008，27(1):111–128.

［8］Yanovsky I，Thompson PM，Osher S，et al.Topology preserving log-unbiased nonlinear image registration: theory and implementation［A］.In:Computer Vision and Pattern Recognition［C］.Minneapolis:IEEE，2007.1-8.

［9］李傳富.顱腦CT病變的計算機自動化檢出方法研究［D］.合肥:中國科學技術大學，2005.

A Non-rigid Registration of the Cerebral CT Images with Large Deformations

XU Feng1LIU Wei1LI Chuan-Fu1，2FENG Huan-Qing1
*1(Department of Electronic Science and Technology，University of Science and Technology of China，Hefei 230027，China)2(Medical Imaging Center，F(xiàn)irst Affiliated Hospital of Anhui Traditional Chinese Medicine College，Hefei 230031，China)

R319

A

0258-8021(2010)02-0172-06

10.3969/j.issn.0258-8021.2010.02.003

2009-06-19，

2009-11-11

國家自然科學基金資助項目(60771007)

*通訊作者。 E-mail:hqfeng@ustc.edu.cn