基于點集與互信息的肺部CT圖像三維彈性配準算法

夏 威，高 欣，王 雷，周志勇，張 冉

(1.中國科學院長春光學精密機械與物理研究所，吉林長春130000;2.中國科學院蘇州生物醫(yī)學工程技術研究所醫(yī)學影像室，江蘇蘇州215163;3.中國科學院大學，北京100049)

肺部組織運動是復雜、不均勻的，為了獲得肺部組織及病灶運動軌跡，文中借助影像信息，通過圖像配準的方法，構筑肺部三維運動模型.圖像配準是指確定兩幅圖像中對應點的位置關系，待配準的圖像可以是不同時間點采集的圖像或者屬于不同模態(tài)的圖像.彈性配準是醫(yī)學圖像配準的重要研究領域，常用于追蹤病灶變化情況、運動校正、計算機輔助影像導航［1-3］.對一個周期內多個時間點上的肺部圖像進行配準是術前計算機手術路徑規(guī)劃、模擬手術和手術導航的必要步驟.通過配準能夠確定肺部組織的運動軌跡，提高手術器械在肺內部的定位.

目前主要有3種彈性配準算法:①基于特征的配準算法;②基于灰度的配準算法;③ 混合配準算法.基于特征的配準算法使用圖像中明顯的幾何結構或者解剖特征做圖像配準，例如特征點等，常用的特征點方法有SIFT點，SURF點等，醫(yī)學圖像中的特征點常由醫(yī)生手動提取或者使用半自動的特征點標定系統(tǒng).基于特征點的方法計算效率高，但是特征點的自動提取與精確匹配是一個難點［4］;另外，由于特征點分布稀疏，遠離特征點的區(qū)域的變形場精度不高［5］.基于灰度信息的配準算法利用圖像的灰度信息構成相似性測度來衡量配準的程度［3，6］，這種方法能夠實現(xiàn)全自動配準，并且精度較高，但是會耗費更多的時間與內存.對于復雜的肺部圖像配準，混合配準算法受到越來越多的關注，混合配準算法整合了基于灰度和基于特征的方法［7-8］，首先使用手動或者半自動的特征點提取方法提取氣管或血管樹的分支點，將分支點的位移信息作為灰度信息的約束，再使用基于灰度的配準算法進行配準，這種方法能夠顯著提高配準的精度，但特征點提取復雜費時不適用于臨床［8］.

表面組織是肺部的顯著解剖結構，表面分割的自動算法快速簡單［9］.肺表面的膨脹與收縮反映了肺部組織的形變趨勢，因此文中擬利用肺表面信息，提出一種全自動混合配準算法.該方法整合了基于特征與基于灰度的配準算法，先使用一致點漂移算法(coherent point drift，CPD)［10］得到點集位移向量，將點集位移向量與基于B樣條變形模型［11］的位移向量之間的歐幾里得距離最小化，得到基于B樣條變形模型的大尺度變換函數(shù)，再將此大尺度變換函數(shù)作為基于互信息的配準算法［6］的初始變換函數(shù)進行細化配準.變換函數(shù)的計算中引入節(jié)點位移約束［11］，從而保證圖像不會發(fā)生折疊現(xiàn)象，確保配準成功.

1 變形模型

給定一組肺部三維CT圖像F和R，設F為浮動圖像，R為參考圖像，定義x=(x，y，z)為圖像中體素的坐標.配準就是確定一個能使F匹配到R的變換函數(shù)T(x).

對于肺部復雜的彈性形變，需要先確立彈性變形模型，文中使用基于B樣條的變形模型［11］.定義u為間距s=(sx，sy，sz)的均勻網(wǎng)格;di，j，k為網(wǎng)格u中第i，j，k個節(jié)點的位移.坐標x的變換函數(shù)定義如下［11］:

式中:i=?x/sx」-1，j=?y/sy」-1，k=?z/sz」-1，u=x/s-i-1，v=y/s-i-1，w=z/s-i-1;B0，B1，B2，B3是4個三次B樣條基函數(shù)，即B0(t)=(-t3+3t2-3t+1)/6，B1(t)=(3t3-6t2+4)/6，B2(t)=(-3t3+3t2+3t+1)/6，B3(t)=t3/6，其中 0≤t＜1.

使用B樣條變形模型有兩個優(yōu)點［12］:① 與薄板樣條等相比，B樣條模型是局部可控的，一個節(jié)點只會影響鄰近4×4×4體素，因此在計算變換函數(shù)T(x;d)時效率較高;②B樣條基函數(shù)是C2連續(xù)的，能夠計算變換函數(shù)T(x;d)的解析梯度，在優(yōu)化過程中可使用基于梯度的優(yōu)化方法.傳統(tǒng)的B樣條模型有可能產生不可逆的變換函數(shù)，即當某些節(jié)點位移過大時，會導致圖像的折疊，從而導致配準失敗［11］.為了避免圖像折疊，需要對節(jié)點位移進行約束.根據(jù)三次 B 樣條可逆性條件［11］，第i，j，k個節(jié)點分別在x，y，z方向上的位移dx，dy，dz應滿足條件:

其中k是一常數(shù)，在三維空間值為2.479 472 335［11］.在這一約束條件下，圖像中體素坐標x的最大位移被限制為(sx/k，sy/k，sz/k)，當且僅當周圍4×4 ×4 節(jié)點的位移都是(sx/k，sy/k，sz/k)時才能達到.

網(wǎng)格節(jié)點的位移是變換函數(shù)的參數(shù)，因此變換函數(shù)的精度取決于網(wǎng)格的間距，小間距的網(wǎng)格能夠構建更加精細的變換函數(shù).然而肺部組織在呼吸過程中的位移較大，由于限制了節(jié)點位移的大小，必須使用大間距的網(wǎng)格以滿足肺部組織的大尺度位移，但又需要小間距的網(wǎng)格以保證變換函數(shù)的精度，為了解決這一矛盾，需要使用n層B樣條模型框架.n層B樣條模型框架依次使用由疏到密的n個網(wǎng)格求得n個變換函數(shù)，大間距的網(wǎng)格能夠得到大尺度變換函數(shù)，滿足大尺度形變的要求，小間距的網(wǎng)格能夠得到細化變換函數(shù)，同時網(wǎng)格節(jié)點位移的大小滿足上述可逆性約束條件，最終的變換函數(shù)是各層求得的變換函數(shù)的組合.

文中采用2層B樣條模型框架，依次使用兩個網(wǎng)格u1和u2，第1層網(wǎng)格u1間距取為60個體素，第2層網(wǎng)格u2間距取為15個體素.整個配準過程從第1層開始，由第1層得到大尺度變換函數(shù)Tp，再將Tp作為第2層的初始參數(shù)，得到細化變換函數(shù)Ti，每一層網(wǎng)格的節(jié)點位移都滿足可逆性條件，最終的變換函數(shù)T是Ti與Tp的組合:

2 基于點集的彈性配準

基于點集的彈性配準目的是生成大尺度變換函數(shù)Tp，基于點集的彈性配準分為兩步:① 點集提取與CPD點集彈性配準，得到點集坐標位移向量P(Xt，Xf);② 由點集坐標位移向量P(Xt，Xf)得到基于u1的大尺度變換函數(shù)Tp.

2.1 點集提取與點集彈性配準

肺表面需要被分割出來以生成表面點集.由于身體部分密度較高，而肺部組織充滿空氣密度較低，因此CT圖像上的身體部分與肺部組織灰度值差異明顯，文中采用了一種快速簡單的基于閾值和形態(tài)學的自動分割方法來分割得到肺的表面區(qū)域［9］.

將分割得到的表面區(qū)域使用Canny邊緣提取算法提取出輪廓面，由于輪廓面上的點數(shù)過多，只需要其中的一部分，依據(jù)點的坐標分布，均勻地取出輪廓上1/50的點生成肺表面點集.設浮動圖像F中取出的肺表面輪廓點數(shù)量為N，則浮動圖像F的表面點集SF=(xF1，xF2，…，xFi，…，xFN)，其中xFi=(xFix，xFiy，xFiz)，各分量分別表示點xFi在x，y，z方向上的坐標值;設參考圖像R中取出的肺表面輪廓點數(shù)量為M，則參考圖像R中的表面點集為SR=(xR1，xR2，…，xRi，…，xRM)，其中xRi=(xRix，xRiy，xRiz).

肺表面點集反映的是肺的解剖形態(tài)，在配準的過程中點集不應改變拓撲結構，一致點漂移算法(CPD)［10］加入了一致性約束，可以保證點集的拓撲結構，非常適用于肺表面點集配準.CPD是一種魯棒性很強的點集配準算法，能夠得到高精度的點集彈性配準結果.它基于混合高斯模型(GMM)，將點集匹配看做概率密度估計問題，給定一對點集a與b，將a看做高斯混合模型的質點，通過最大似然概率完成a與b的匹配，a中的點在配準過程中分別移動到合適的位置以生成配準后的點集c，點集c再與b對齊.

通過CPD算法對肺表面點集SF和SR進行點集彈性配準，將SF作為高斯混合模型的質點，計算SF中各個點的位移，生成配準后的點集ST=(xT1，xT2，…，xTi，…，xTM)，再與SR的對齊.考慮到點集中點的數(shù)量較大，使用CPD算法的快速實現(xiàn)［10］:使用低秩矩陣估計與快速高斯變換對算法進行加速.CPD算法中含有3個自由參數(shù):ω，λ和β.其中ω是對噪聲數(shù)量的假設值，λ和β是平滑性約束參數(shù)，文中參數(shù)取值為ω=0.7，λ=3和β=4.

2.2 大尺度變換函數(shù)Tp

定義點集坐標位移向量P(xT，xF)為

式中:xF=(xFx，xFy，xFz)是點集同SF中的點;xT=(xTx，xTy，xTz)是點集ST中對應于點xF的點.

定義大尺度變換函數(shù)Tp為式(1)中的變換函數(shù)，大尺度變換函數(shù)Tp反映的是包括肺表面組織與內部組織在內的整個肺部組織的位移，點集坐標位移向量P(xT，xF)反映的是肺表面組織的位移，在肺表面組織處由Tp確定的位移向量應該與P(xT，xF)相等，因此使用Tp擬合P(xT，xF):設ΔxF是xF點處由Tp確定的位移向量，則ΔxF與點集坐標位移向量P(xT，xF)相等，即ΔxF到P(xT，xF)的歐幾里得距離等于0.

定義點集位移向量P(xT，xF)與ΔxF的歐幾里得距離的平方和為

ΔxF是節(jié)點位移d的函數(shù)，則有

以E(SF，ST;ΔxF)為目標函數(shù)，將式(7)代入式(6)可得

其中下標x，y，z表示變量在x，y，z方向上的分量.式(9)的優(yōu)化問題是求d以便使目標函數(shù)E(Sf，St;d)達到極小值.d為所有節(jié)點位移向量，假設有m個節(jié)點，則d=(d1，d2，…，di，…，dm)，其中di=(dix，diy，diz)是節(jié)點i的位移向量.目標函數(shù)E(Sf，St;d)的梯度為

如果Tp的分量的方向與dix的方向不一致，則對應的導數(shù)值為0，即為

采用有約束的有限容量擬牛頓優(yōu)化方法LBFGS-B［13］來最小化目標函數(shù)E(SF，ST;d).這一優(yōu)化方法適用于高維度的參數(shù)空間，并且能夠在優(yōu)化的過程中對位移向量d的分量di的大小進行約束，我們通過這種方法對節(jié)點位移施加前述的可逆性約束條件(2)，以保證配準成功.

3 基于互信息的配準

基于點集的配準得到了大尺度變換函數(shù)Tp，這一變換函數(shù)反映的是不夠精確的大尺度形變，還需對此變換結果進一步細化.很多試驗表明，基于灰度的方法對局部形變比較敏感，在圖像之間差異比較小的時候能夠取得精確的結果［6-7］.由于呼吸運動，肺部組織的膨脹與收縮會導致組織密度的改變，浮動與參考圖像中相同的肺部組織對應的灰度值并不一定相等，基于互信息的方法并不直接依賴于灰度值的對應相等關系，因此文中選取基于互信息(mutual information)的Mattes方法［6］對大尺度變換函數(shù)Tp進行細化.Mattes方法中互信息目標函數(shù)定義為

式中:d為節(jié)點位移向量;p，pT和pR分別為浮動圖像和參考圖像的聯(lián)合概率密度分布函數(shù)、浮動圖像邊緣概率密度分布函數(shù)、參考圖像邊緣概率密度分布函數(shù)［6］.Mattes方法使用Parzen窗方法進行連續(xù)聯(lián)合直方圖估計以計算互信息目標函數(shù)的解析梯度，使用單層B樣條模型作為變形模型，使用L-BFGS-B優(yōu)化方法對互信息目標函數(shù)進行優(yōu)化.

混合配準算法將Tp作為Mattes方法的初始變換函數(shù);設F'是對浮動圖像F應用Tp進行變換后生成的圖像，將F'作為Mattes方法的浮動圖像，將R作為Mattes方法的參考圖像;使用2層B樣條變形模型中的u2作為Mattes方法的變形模型;使用LBFGS-B優(yōu)化方法對互信息目標函數(shù)進行優(yōu)化，同時對節(jié)點位移施加可逆性約束條件(2)，最終得到細化變換函數(shù)Ti(x;d).

4 試驗

文中采用4套3D肺部CT圖像對該算法進行評估，每一套圖像都是肺部運動周期中的最大吸氣時刻與最大呼氣時刻圖像.CT圖像的大小為512×512 ×320，空間分辨率是0.97 mm ×0.97 mm ×2.00 mm.試驗平臺硬件環(huán)境是 CPU為 Intel Core i3 2350m@2.3 GHz和RAM為8 GB的計算機.

提取肺表面點集SR與SF，使用CPD對點集SR與SF進行彈性配準.圖1顯示了肺表面點集SR與SF的提取效果，圖2顯示了配準前后的點集.

圖1 肺表面點集SR與SF

圖2 點集彈性配準前后的肺表面點集

算法性能的量化指標由互信息值、標志點距離誤差以及計算時間組成.

互信息值反映了一對圖像的相似程度，一般互信息值越大則圖像越相似.初始的互信息值由浮動圖像與參考圖像的互信息計算得到，配準后的互信息值由經過變形后的浮動圖像與參考圖像的互信息計算得到.

標志點距離誤差反映了配準算法的精度［14］.使用半自動的標志點標定軟件系統(tǒng)(Murphy，2008)［13］對肺部區(qū)域內進行標志點標定，使用文獻［13］中描述的標定步驟，每一對圖像中都標定了250對均勻分布的標志點，每一對標志點表示相同部位在浮動與參考圖像上的對應位置.圖3以三維視角顯示了標志點的分布情況.圖4以綠色十字顯示標志點標定系統(tǒng)中的標志點位置，由于將所有標志點投影到一層圖像上，部分綠色十字顯示在該層肺部區(qū)域之外.初始的標志點距離誤差是浮動和參考圖像上對應的每對標志點的歐幾里得距離的均值.配準后的標志點距離誤差是浮動圖像上的標志點坐標與配準結果所預測的標志點坐標的歐幾里得距離的均值.

圖3 標志點分布示意圖(藍色小球為所選取的標志點)

圖4 標志點標定系統(tǒng)中標志點的位置(綠色十字形為標志點)

將混合配準算法與僅基于互信息的Mattes方法［6］對相同的圖像進行配準得到的結果加以比較，對比結果列于表1中，其中MI表示互信息值，LD表示標志點距離誤差，SD表示LD的標準差.Initial是初始值，Mattes代表僅基于互信息的Mattes方法的結果，PR代表混合配準算法中基于點集的配準算法的結果，HR代表混合配準算法的結果.

表1 配準算法的對比結果

由試驗結果可以看出，混合配準算法HR具有最小的標志點誤差和最大的互信息值，基于點集的配準算法PR的標志點誤差最大、互信息值最小，僅基于互信息的Mattes方法的標志點誤差與互信息值介于兩者之間;僅基于互信息的Mattes方法耗時最多，基于點集的配準算法PR耗時最少，而混合配準算法HR耗時介于兩者之間.相對于僅基于互信息的Mattes方法，混合配準算法HR的標志點誤差平均減少了5%，互信息值平均提升了1%，耗時平均減少了70%;基于點集的配準算法PR耗時比混合配準算法HR平均減少74%，但混合配準算法HR的標志點誤差平均減少了28%，互信息值平均提升了16%.

5 結論

1)相對于僅基于互信息的配準算法，所提方法減少了配準時間，提高了配準精度.

2)點集配準可以生成較精確的大尺度變換函數(shù)，為混合配準算法中基于互信息的方法提供了誤差較小的初始變換函數(shù)，因此經過細化后可以得到誤差更小的細化變換函數(shù).

3)混合配準算法能夠避免陷入局部極值并擁有更好的魯棒性，能夠獲得更小的誤差并且加快收斂速度.

4)所提算法對于確定肺部組織的運動軌跡、建立肺部呼吸運動模型、實現(xiàn)計算機輔助術前手術規(guī)劃與手術導航有很大的幫助.

References)

［1］Boldea V，Sharp G C，Jiang S B，et al.4D-CT lung motion estimation with deformable registration:quantification of motion nonlinearity and hysteresis［J］.Med Phys，2008，35(3):1008-1018.

［2］Dinkel J，Hintze C，Rochet N，et al.Computed tomography of the lungs［J］.Radiologe，2009，49:698-704.

［3］Rueckert D，Aljabar P.Nonrigid registration of medical images:theory，methods，and applications［J］.IEEE Signal Processing Magazine，2010，27(27):113-119.

［4］張俊雄，吳科斌，宋 鵬，等.基于BP神經網(wǎng)絡的玉米單倍體種子圖像分割［J］.江蘇大學學報:自然科學版，2011，32(6):621-625.Zhang Junxiong，Wu Kebin，Song Peng，et al.Image segmentation of maize haploid seeds based on BP neural network［J］.Journal of Jiangsu University:Natural Sci-ence Edition，2011，32(6):621-625.(in Chinese)

［5］Muenzing S E A，van Ginneken B，Murphy K，et al.Supervised quality assessment of medical image registration:application to intra-patient CT lung registration［J］.Medical Image Analysis，2012，16:1521-1531.

［6］Mattes D，Haynor D R，Vesselle H，et al.PET-CT image registration in the chest using free-form deformations［J］.IEEE Trans Med Imaging，2003，22(1):120-128.

［7］Gorbunova V，Durrleman S，Lo P，et al.Lung CT registration combining intensity，curves and surfaces［C］∥Proceedings of ISBI2010 7th IEEE International Symposium on Biomedical Imaging:From Nano to Macro.Rotterdam， Netherlands:IEEE ComputerSociety，2010:340-343.

［8］Yin Y B，Hoffman E A，Ding K，et al.A cubic B-spline-based hybrid registration of lung CT images for a dynamic airway geometric model with large deformation［J］.Physics in Medicine and Biology，2011，56:203-218.

［9］?zkan H，Osman O，?ahin S，et al.Lung segmentation algorithm for CAD system in CTA images［J］.World A-cademy of Science，Engineering and Technology，2011，77:306-309.

［10］Myronenko Andriy，Song Xubo.Point set registration:coherent point drift［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2010，32(12):2262-2275.

［11］Choi Y，Lee S.Injectivity conditions of 2D and 3D uniform cubic B-spline functions［J］.Graphical Models，2000，62:411-427.

［12］Yin Y B，Hoffman E A，Lin C L.Mass preserving nonrigid registration of CT lung images using cubic B-spline［J］.Med Phys，2009，36(9):4213-4222.

［13］Murphy K，van Ginneken B，Pluim J P W，et al.Semiautomatic reference standard construction for quantitative evaluation of lung CT registration［C］∥Proceedings of2008 11th Medical Image Computing and Computer-Assisted Intervention.New York:Springer Verlag，2008:1006-1013.

［14］Murphy K，van Ginneken B，Reinhardt J M，et al.E-valuation of registration methods on thoracic CT:the EMPIRE10 challenge［J］.IEEE Transactions on Medical Imaging，2011，30(11):1901-1920.