四維CT肺圖像多層次時空一致性配準方法

楊 浩 胡冰鋒 孫瑞芳 楊 烜*

1(西安電力高等?？茖W校 陜西 西安 710032)2(深圳大學計算機與軟件學院 廣東 深圳 510086)

四維CT肺圖像多層次時空一致性配準方法

楊 浩1胡冰鋒2孫瑞芳2楊 烜2*

1(西安電力高等專科學校 陜西 西安 710032)2(深圳大學計算機與軟件學院 廣東 深圳 510086)

四維CT肺圖像提供了肺部動態信息，可廣泛用于肺部疾病診斷和精確放射治療。提出一種針對四維CT肺圖像的多層次時空一致性配準方法。該方法首先從最大吸入相位的圖像中提取控制點，利用結構張量跟蹤算法對控制點進行跟蹤；構建多層時空一致性優化模型，對控制點進行調整，從整圖配準角度提高控制點配準精度；每次調整后根據配準的一致性誤差，添加新的控制點，進一步降低配準的一致性誤差，并對新的控制點集重新進行優化調整。該模型在空間上對形變進行約束，保持形變扭曲的平滑性；在時間上對控制點軌跡進行約束，保持控制點的軌跡平滑性；通過一致性約束，降低了正向和反向形變函數的一致性誤差。在測試數據集上的實驗結果表明，該方法可以得到平滑、穩定、合理的形變結果，可以用于估計肺部運動模型。

圖像配準 CT圖像 一致性形變 時空約束

0 引 言

CT成像技術已廣泛用于醫學圖像研究領域。近些年，三維CT成像技術已經發展到四維成像，即在三維圖像的基礎之上加上時間維度，構成一系列離散的三維圖像集合。四維CT成像可以用于研究肺部的生理特征，基于4DCT的圖像配準技術可以得到肺部的運動模型。該模型可以用于圖像引導放療領域，為圖像引導放療提供更準確的肺部運動信息(包括腫瘤運動)，結合精確放療技術，減少放射量對正常組織的損傷，提高肺癌的治療效果[1]。

圖像配準就是估計浮動圖像與目標圖像之間的形變函數，使浮動圖像經過形變之后與目標圖像的差異盡可能小。由于四維圖像在時間軸的關聯性，四維配準并不是簡單的多個三維圖像配準，其求解的形變函數不僅僅要考慮諸如空間平滑約束，還需考慮時間平滑、正反形變的一致性約束等。所謂空間平滑指的是某個特定時間的空間形變具有較低的扭曲能量；時間平滑指器官隨時間變化形成穩定的運動軌跡，沒有大幅度的震蕩；一致性約束是指正向和反向形變之間滿足相互可逆的條件，這樣得到的形變更具生物學意義[2]。目前空間平滑和空間平滑約束已有較成熟的研究成果[3-4]，而一致性約束是一個比較困難且仍在進行研究的問題。

一致性圖像配準要求建立浮動圖像和目標圖像之間的相互可逆的正向和反向形變函數，即浮動圖像中的像素點A能映射到目標圖像唯一像素點B。反之，目標圖像中的點B也都能通過反向形變映射到浮動圖像中的點A。而目前大量的圖像配準技術并不滿足這個約束[5]，這種不一致造成的結果是配準得到的對應關系并不能真正反映圖像目標間的對應關系。

目前一致性配準算法主要分為兩類，一類如Zhang等[6]提出的基于圖像強度信息的一致性配準方法，如一致性demons配準算法，該算法在demons原算法基礎上加入一致性的約束；Alvarez等[7]針對稠密光流場引入了雙向一致性規約；通過迭代優化一致誤差，求解得到一致性形變函數；Ye等[8]在變分彈性模型下提出一種一致性圖像配準方法，通過引入殘差圖像的一致性變換約束提高配準精度；Leow等[9]在線性彈性模型上提出一種同時迭代正反變換的一致性方法，其單向形變函數通過對稱代價泛函進行求解，避免計算正反變換的逆函數；Modat等[10]在FFD(Free-Form Deformation)基礎上提出了一種一致性對稱自由形變算法。

另一類是以Johnson等[2]為代表的基于控制點對的一致性配準方法。這類方法在基于薄板樣條的形變函數下加入一致性形變約束，通過迭代逼近求解一致性形變。但是該方法存在控制點對應關系存在誤差的問題，同時計算量巨大。Yang等[12]提出了一種快速一致性形變算法，該算法省略了傳統算法[2]的求逆步驟，可節省大量的運算。Shen等[13]則提出了代價相對較低的一致性配準優化模型，該模型在代價函數中同時加入一致性約束，通過多次迭代優化使正向形變和反向形變滿足互逆條件。

本文提出了一種時空一致性約束優化模型。該模型首先提取控制點，對控制點進行跟蹤，在跟蹤的基礎上利用時空約束優化模型對控制點進行調整，同時加入了空間平滑約束、時間平滑約束和簡化的一致性約束。進一步地，在每一次優化之后，又根據形變的一致性誤差增加控制點，并對新增控制點進行優化調整，重復該過程，直至一致性誤差滿足要求。該方法可以提高四維圖像配準精度，得到平滑、穩定、滿足一致性約束的形變函數。

1 時空一致性配準模型

優化的思想就是構建代價函數，利用跟蹤準確的控制點對應關系調整跟蹤不準確的控制點。如圖1所示，跟蹤算法得到所有控制點在不同相位圖像中的預測位置，進而根據控制點對應關系，利用薄板樣條算法[3]求解浮動圖像到所有目標圖像的形變函數ft。由于大部分控制點跟蹤較為準確，ft可以較好地表達相位間的形變關系，為后期優化調整提供一個理想的初值。優化的思想就是建立目標函數，通過調整控制點的跟蹤位置使目標函數達到最優，進而根據調整后的控制點確定形變函數。

圖1 4D肺圖像配準示意圖ft是正向形變，ht反向形變

在優化過程中，代價函數容易陷入局部極值，而控制點的數量越多，優化過程陷入局部極值的可能性越大。針對該問題，本文提出多層次配準方法。該方法在初始控制點基礎上進行優化配準，每一次優化之后再根據一致性誤差增加新的控制點，并對新增的控制點進行優化調整，直至滿足一致性誤差的條件，如圖2所示。這樣就避免了一次性配準中控制點數量過多可能陷入局部極值的問題，同時逐層增加控制點也有助于減少一致性誤差。

圖2 多層次配準流程示意圖

1.1 控制點提取與跟蹤

1.2 時空約束優化

根據上述的肺運動特點，本文提出的優化模型考慮以下幾個因素：1) 形變圖像和目標圖像盡可能相似；2) 正向形變和反向形變之間盡可能滿足一致性約束；3) 空間形變的扭曲能量盡可能??；4) 控制點的時間軌跡要盡可能平滑。這里的一致性約束并沒有嚴格要求正反形變互逆，而是通過正向和反向形變的配準誤差表達一致性約束程度。這種方法是一種簡化的一致性誤差約束方法，可以避免求解正向形變函數的逆函數，從而減少優化的時間消耗。

本文的時空一致性約束配準模型由三部分組成：

CE(μ)=CESim(μ)+α·CEB(μ)+β·CES(μ)

(1)

(4)

優化模型中第一部分CESim是相似性度量，用來描述目標圖像It和形變圖像I0(ft)之間的相似度，以及浮動圖像I0和反向形變圖像It(ht)之間的相似度，體現了一種簡化的一致性約束。如果二者能同時達到最小，也能達到正向形變和反向形變趨于一致的效果。需要說明的是，本文比較的是控制點周圍的一個鄰域W之間的相似程度，而不是整個三維圖像的相似程度，這樣處理的目的，主要是為了提高計算效率，同時實驗分析也表明這種局部相似度是可行的。

第二部分用來計算正向形變函數ft和反向形變ht的扭曲能量[3]，形變扭曲能量越高表明形變越劇烈。

第三部分計算控制點軌跡的平滑度，肺圖像上的點隨時間變化的軌跡應該是比較平緩的，劇烈震蕩的軌跡往往是配準精度下降的表現。本文采用拉普拉斯算子度量軌跡的平滑度。

本文采用擬牛頓(quasi-Newton)優化算法的改進算法BFGS(BroydenFletcherGoldfarbShanno)求解代價函數CE(μ)。BFGS優化算法具有高精確度和快速收斂的優勢，在計算機內存不足、計算量大的情況下可以用L-BFGS(Limited-Memory Quasi-Newton Method)近似。在優化過程中，如果控制點的個數比較多，就選擇L-BFGS優化算法。

1.3 形變模型

本文方法中采用薄板樣條[3]對形變場進行估計，其中圖像中非控制點的形變要靠所有控制點的拉伸作用力去控制。如果控制點過少就會造成遠離控制點區域形變誤差較大，同時一致性誤差會變得更大。如果在一致性誤差較大的區域加入新的控制點，既可以保證一致性誤差有所改善，同時也可以使新控制點附近區域的形變更準確。

1.4 計算一致性誤差

假設當前正向形變函數為f，反向形變函數為h，點x處的正向一致性誤差ef和反向一致性誤差eh定義如下：

ef=f(h(x))-xeh=h(f(x))-x

(5)

1.5 增加控制點

多層次配準方法的思想就是在少量初始控制點基礎上，逐漸增加新的控制點進行更細致的配準，直到滿足結束條件為止。增加新的控制點的原則就是基于當前形變的一致性誤差。由于控制點上的一致性誤差為0，通過在一致性誤差較大的區域加上控制點，可以減少新增控制點鄰域內的一致性誤差，同時也可以提高形變函數的精度。我們把局部范圍內一致性誤差最大的點列為新增控制點候選，如果候選新增控制點處于非特征結構上(如非血管或非血管分叉等)，可能會導致跟蹤困難。所以本文提取一致性誤差較大且處于血管結構上的點作為本次迭代新增控制點。

在確定了新增控制點之后，我們采用SEST算法預測新增控制點在所有相位目標圖像中的位置，然后通過時空一致性優化模型對這些新增控制點進行調整。將調整后的新增控制點加入到當前層次控制點集中得到新的控制點集，再進行下一層的迭代，這樣逐層增加控制點，直至一致性誤差滿足結束條件。

1.6 配準流程

整個配準過程是通過迭代增加控制點，并調整控制點的位置，以完成對目標的優化，完整過程見2.1節。層次時空一致性配準步驟如下：

Step1在浮動圖像I0中提取一定數量的初始控制點集ψ0。設置迭代收斂閾值參數ε，令迭代數i=0。

本文算法是在MATLAB R2014a開發環境， Windows Server 2008操作系統環境下實現的。下面附上算法的部分核心代碼：

While (iterε) //iter為迭代次數，MaxIter為最

//大迭代次數，一致性誤差大于收斂誤差

iter=iter+1;

MCE_all(iter)=MCE;

//MCE是平均一致性誤差

BCEpointset_all=selectControlPoints(I_ref,CE,iter,init_grid,MCE);

//選擇新增加的控制點

var_added_all_5 = pointsCorrdinatesTrack(I_ref,I_tar,

BCEpointset_all);

//對新增控制點進行跟蹤

[BCEpointset,var_added_5]=deleteExceptinoalTrackedPoints

(init_grid,var_5,BCEpointset_all,var_added_all_5);

//刪除異常跟蹤點

var_added_other = getLMOtherPhaseCors(I_ref,BCEpointset,4,cases);

//獲取控制點在其他相位的坐標位置

…

var=fminlbfgs(@(x)costfunction(I_ref,x,init_grid,sizeV_new,

weight_cost,cases),var,optim);

//對控制點進行優化調整

[init_grid,var] = changeSamePoints(init_grid,var);

CE = getConsistentError(sizeInit,init_grid,round(var),

size_enlarge,C);

//計算一致性誤差

I_CE = double(I).*CE; MCE=sum(I_CE(:))/num_lung_tissue;

//計算平均一致性誤差

end

2 實驗結果及分析

本文實驗環境采用雙處理器服務器，Intel Keon CPU E5620，主頻2.4 GHz，內存8 GB。通過dir-lab[16]網站提供的四維CT肺數據集驗證本文提出方法的有效性。該數據集共有10組不同病人的4D肺CT數據，每組數據提供了從T00到T50相位的6個時間點75個標志點的專家標注位置，這些標志點可用于評估本文方法的配準精度。表1列出了該數據集中圖像大小及分辨率。

表1 DIR-lab數據

為了避免胸腔數據的干擾，本文對原始圖像進行分割，得到肺實質數據。進一步地，采用文獻[17]的方法對肺圖像進行血管增強。實驗中所用優化參數見表2，其中代價函數空間平滑權重值大于軌跡平滑權重，是因為扭曲能量的計算結果比較大；平均一致性誤差閾值為迭代收斂條件，在實驗過程中發現，一致性誤差閾值過小會導致算法多次迭代，而一致性誤差提高并不明顯，經過多次實驗，得出該參數的合適經驗值為0.01。控制點局部鄰域窗口大小主要是用于計算目標函數中的局部區域度量。該區域不易過大，過大會帶來較大的計算量。該區域也不能過小，過小其統計結果不夠穩定；搜索區域大小主要是SEST跟蹤時的搜索范圍，過大的搜索范圍容易發生誤匹配，而過小的搜索范圍容易找不到對應位置。

表2 配準方法參數列表

2.1 配準誤差

本文利用專家標志點來評價多層次時空一致性配準HSCR(HierarchicalSpatio-temporal Consistent Registration,)方法的配準精度，以專家標志點與算法得到的標志點之間的距離(TRE)為評價標準。我們把HTSC的配準結果和文獻[18-19]的結果進行對比，表3列出75個控制點的TRE比較?？梢钥闯觯琀SCR的配準結果在平均誤差方面都大幅度優于文獻[18]的配準結果，而在標準差方面稍微高于[18]。HSCR在case1-case5的配準精度與文獻[19]相當，而在case6至case10實驗數據方面，HSCR的配準結果明顯好于文獻[19]的配準結果。綜上分析，從整體的配準誤差統計分析來看，HSCR取得了比較理想的配準效果。

表3 以75個標志點的配準誤差，括號外為平均誤差，括號內為標準差，“*”為沒有數據

圖3為10組數據中75個控制點的HSCR配準誤差分布圖，其中每組數據統計的是所有標志點在T10-T50 5個時刻的全部誤差。圖中可以看出這10組數據中75%的控制點配準誤差在2 mm以下，50%的標志點跟蹤誤差在1 mm以下。表4列出了75個控制點箱線圖統計中的離群點信息，num表示離群點個數，mean表示離群點的平均TRE，percent表示離群點所占的比例。從該表中可以看出case1的離群率較高，但是離群平均誤差較小；case6、case8和case9的離群平均誤差較大，但是離群率還是較小的。HSCR的配準誤差標準差較大，但是離群率還是較小的，所以總體的配準性能還是良好的。

圖3 case1-case10 10組數據上75個標志點配準平均誤差boxplot分布

Case12345678910num277785129161714mean2.953.493.506.316.659.626.737.543.327.94percent7.2%1.8%1.8%2.1%1.3%3.2%2.4%4.2%4.5%3.7%

為了進一步評價HSCR的配準效果，圖4展示了4個樣本的配準結果，這4個樣本選自呼吸運動較劇烈的case5。其中a為參考相位(T00)一個控制點的局部鄰域切片圖像；b所在行有兩幅圖，從左至右分別為配準前所有6個相位(T00-T50)沿著a中橫、縱虛線采樣的圖像強度疊加；c所在行有兩幅圖，從左至右分別為配準后所有6個相位沿著a中橫、縱虛線采樣的圖像強度疊加。從圖中可以看出，配準前同一目標位置有較大差異，導致在b中顯示出位移的效果，而配準后同一目標在c中基本是對齊的，這表明配準的結果比較理想。

圖4 四維肺CT圖像配準結果

2.2 一致性誤差和平滑性

在保證配準精度達到一定程度的情況下，形變的一致性和平滑性越好，配準的效果也就更加合理。本文將從一致性誤差和形變平滑性兩方面評價HSCR配準方法。表5列出了一致性配準的統計結果，包括初始控制點數量(Initial Number)、最終收斂控制點數量(Final Number)、迭代次數(Iteration)、初始配準的平均一致性誤差(Initial Consistent Error)和最終收斂的配準平均一致性誤差(Final Consistent Error)。

表5 10組數據一致性誤差數據統計

從表5可以看出，HSCR收斂后的平均一致性誤差比初始控制點得到的平均一致性誤差明顯下降，而且控制點的數量也隨著迭代次數增加而增加，這表明層次增加控制點對于解決配準一致性的問題是有效的。

為了更加詳細地表明增加控制點的效果，選取了case4的配準結果進行分析。以case4的第73橫斷面切片一致性誤差為例，如圖5所示，可以看出新增的控制點都是位于一致性誤差較大的位置(圖像越亮，一致性誤差越大)“o”為初始控制點位置，“+”增加的控制點位置。(a)的平均一致性誤差為0.091 1，最大的一致性誤差為0.694 8；(b)的平均一致性誤差為0.028 3，最大的一致性誤差為0.326 4。上述分析表明HSCR增加的控制點可以使局部區域的一致性誤差減小。

圖5 增加控制點前后的一致性誤差對比

本文將從形變區域的雅可比行列式統計形變的平滑程度，原則上，形變的雅克比行列式越接近1，表明形變扭曲的情況越小，體保持特性更好。本文以case5為例，圖6為增加控制點前后形變結果的雅可比行列式直方圖統計，其中初始形變的雅可比行列式平均值為0.876 9，配準收斂后形變結果的雅可比行列式平均值為0.896 5。從圖6可以看出，經過HSCR優化后的形變結果的雅可比行列式的值更加靠近1，這表明增加控制點后得到的形變比初始形變更加平滑。且改善了形變的體保持特性。

(a) 初始控制點得到的形變的雅可比行列式統計直方圖 (b) 增加控制點后的雅可比行列式統計直方圖圖6 增加控制點前后配準結果得到的雅可比直方圖統計對比

3 結 語

為了得到連續、平滑且滿足一致性形變約束的肺形變模型，本文針對4DCT圖像提出了時空一致性配準方法。該方法提取初始點集，利用SEST算法對點集進行跟蹤，并從形變圖像和目標圖像盡可能相似，形變扭曲盡可能小，控制點軌跡盡可能平滑，以及正反形變之間滿足一致性約束等方面，提出了一個優化模型，對提取的控制點進行優化調整。采用層次時空一致性配準策略，逐層增加新的控制點集，并通過優化模型對控制點進行調整，以提高配準精度。降低正反向形變的一致性誤差。在DIR-lab數據庫中測試了算法性能，實驗結果表明，該配準方法可以達到良好的配準精度，同時可以保證形變的平滑度和一致性，并改善了形變的體保持特性。

[1] Tan C, Yi J B, Yang X. Transformation models for four-dimensional image registration: A survey [J]. Recent Patents and Topics on Imaging, 2015, 5: 88-96.

[2] Johnson H J, Christensen G E. Consistent landmark and intensity-based image registration [J]. IEEE Transactions on Medical Imaging, 2002, 21(5): 450-461.

[3] Bookstein F L. Principal warps: Thin-plate splines and the decomposition of deformations [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(6): 567-585.

[4] Castillo E, Castillo R, Martinez J,et al. Four-dimensional deformable image registration using trajectory modeling[J]. Physics in Medicine and Biology, 2010, 55(1): 305-327.

[5] Sotiras A, Davatzikos C, Paragios N. Deformable medical image registration: A Survey [J]. IEEE Transactions on Medical Imaging, 2013, 32(7): 1153-1190.

[6] Zhang Z, Jiang Y, Tsui H. Consistent multi-modal non-rigid registration based on a variational approach [J]. Pattern Recognition Letter, 2006, 27(7): 715-725.

[7] Alvarez L, Deriche R., Papadopoulo T, et al. Symmetrical dense optical flow estimation with occlusions detection [J]. Lecture Notes in Computer Science, 2002, 2350: 721-735.

[8] Ye X, Chen Y. A new algorithm for inverse consistent image registration [J]. Advances in Visual Computing, 2009, 5875: 855-864.

[9] Leow A, Huang S, Geng A, et al. Inverse consistent mapping in 3D deformable image registration: Its construction and statistical properties [J]. Information Processing in Medical Imaging, 2005, 3565: 493-503.

[10] Modat M, Cardoso M J, Daga P, et al. Inverse-consistent symmetric free form deformation [J]. Information Processing in Medical Imaging, 2013, 4765: 323-329.

[11] Christensen G E, Rabbitt R D, Miller M I. Deformable templates using large deformation kinematics [J]. IEEE Transactions on Image Processing, 1996, 5(10): 1435-1447.

[12] Yang X, Zhang D, Yao S, et al. A fast algorithm to estimate inverse consistent image transformation based on corresponding landmarks[J]. Computerized Medical Imaging & Graphics the Official Journal of the Computerized Medical Imaging Society, 2015, 45:84-98.

[13] Shen D, Davatzikos C. HAMMER: Hierarchical attribute matching mechanism for elastic registration [J]. IEEE Transactions on Medical Imaging, 2002, 21(11): 1421-1439.

[14] Wu J H, Hu B F, Yang X. Deformable registration of 4D-CT image using landmark tracking [J]. Biomedical Research, 2016, 27(2): 797-807.

[15] Baboiu D M, Hamarneh G. Vascular bifurcation detection in scale-space[C]// Mathematical Methods in Biomedical Image Analysis (MMBIA), 2012 IEEE Workshop on. IEEE, 2012: 41-46.

[16] www.dir-lab.com.

[17] Yang X, Pei J. Lung deformation estimation using spatially mean shift for 4D-CT[C]// IEEE International Conference on Bioinformatics and Biomedicine. IEEE, 2013:237-242.

[18] Metz C T, Klein S, Schaap M, et al. Nonrigid registration of dynamic medical imaging data using nD+ t B-splines and a groupwise optimization approach [J]. Medical Image Analysis, 2011, 15(2): 238-249.

[19] Wu G, Wang Q, Lian J, et al. Reconstruction of 4D-CT from a Single Free-Breathing 3D-CT by Spatial-Temporal Image Registration[M]// Information Processing in Medical Imaging. Springer Berlin Heidelberg, 2011:686-698.

CONSISTENTREGISTRATIONMETHODFORTHEMULTILAYERSPATIO-TEMPORAL4DCTLUNGIMAGES

Yang Hao1Hu Bingfeng2Sun Ruifang2Yang Xuan2*1

(Xi’anElectricPowerCollege,Xi’an710032,Shaanxi,China)2(CollegeofComputerScienceandSoftware,ShenzhenUniversity,Shenzhen510086,Guangdong,China)

Four-dimensional CT lung images provide lung dynamic information, which can be widely used in the diagnosis of pulmonary diseases and precise radiotherapy. This paper proposes a hierarchic spatio-temporal consistent registration for 4D CT lung images. The proposed method extracts control points from the image at the maximum inhalation phase, and tracks these control points by using the structure tensor tracking algorithm; a hierarchical spatio-temporal consistent optimization model is proposed to adjust the control points finely, which has the advantages of improvement of the registration accuracy as a whole; new control points are selected from the region where consistent errors are large, and are added to the control point set after each optimization. These newly added control points are optimized further to reduce the consistent errors of registration. In the proposed method, the spatial smoothness constraint keeps the deformation smooth; the temporal smoothness constraint keeps the trajectories of points stable; the consistent constraint reduces the consistent errors introduced by the transformations in two directions. The experimental results of public dataset show that the proposed method can estimate smooth, stable, and reasonable deformation results.

Image registration CT images Consistent transformations Spatio-temporal constraints

TP391

A

10.3969/j.issn.1000-386x.2017.09.040

2016-09-16。國家自然科學基金聯合項目(U1301251)。楊浩，講師，主研領域：醫學圖像處理。胡冰鋒，碩士。孫瑞芳，碩士生。楊烜,教授。