摘 要：為實現相同個體在不同呼吸狀態下產生較大形變的三維肺部醫學影像配準，提出一種基于改進Demons 算法的精確有效配準方案。首先，對待配準影像進行全局非剛性配準。通過尺度不變特征變算法對影像進行特征點提取與匹配，根據匹配結果計算變換參數，完成全局配準；其次，利用改進Demons 算法對全局配準后的影像進行局部非剛性配準。使用改進的方案實現了人體肺部影像配準，并且肺部整體輪廓以及內部組織的配準結果較理想。配準前，影像間的均方誤差值為25835.3，經配準后影像間均方誤差值降為3726.31，均方誤差值下降率為85.58%。提出的方案能夠有效配準三維肺部影像，為對肺部呼吸運動估計以及呼吸功能分析提供良好的基礎。

關鍵詞：肺部；非剛性配準；三維醫學影像；尺度不變特征變換；Demons 算法；均方誤差

引言：醫學影像配準技術是手術計劃制定、圖像引導放射治療、以及醫學圖像融合處理等方面的基礎和前提，具有重要的臨床應用價值。例如，在基于錐形束CT 的圖像引導放射治療中，醫生可以使用設備自動完成將錐形束CT 得到的病灶部位影像與計劃的CT 影像的配準，計算目標區域的中心位置，從而得到治療床所需調整的參數[1]。同時，人體大部分組織和器官具有非剛性特征，因而針對非剛性醫學影像配準的研究具有重要意義[2]。

肺部疾病主要通過計算機斷層成像技術來進行探查。多層螺旋CT 掃描可形成高分辨率組織斷層影像，能夠清晰展現肺部病灶區域的組織結構[3]。四維計算機斷層掃描技術可快速準確形成斷層影像，確定肺部腫瘤位置，捕捉組織與器官的運動軌跡，從而更準確勾畫放療靶區和計算放射劑量等[4]。

根據計算的特征不同，可以分為基于特征和基于灰度的配準算法。前者通過特征的提取和匹配，計算最優配準參數。由于特征信息相對較少，算法的計算效率較高。后者需要利用影像全部信息，避免在提取匹配特征時造成誤差，算法的配準精度高，受影像質量影響較小，有些方法用于多模態醫學影像配準。王安娜等人提出一種基于特征點匹配的非剛性配準方法，通過雙向匹配對應點建立匹配關系，采用粒子群優化算法尋找最優變換參數[5]。王婕妤等人改進傳統變分光流模型，解決傳統算法不能很好處理影像亮度不均勻等問題，并增加特征匹配項解決大形變影像配準問題[6]。通過上述考慮，提出的層次配準模型結合以上兩種算法的特點，通過全局配準加快影像局部配準的收斂速度，并且采用Diffeomorphic Demons 算法使影像局部形變更加精確，實現對產生較大形變的三維肺部影像進行配準。

1 方法與設計1.1 基于特征提取配準算法尺度不變特征變換算法具有對影像間的平移、旋轉、縮放、亮度、噪聲、光照變化、視覺變化保持不變性和穩定性，能有效匹配影像間特征點[7,8]。針對三維影像的處理，實現過程主要由三步組成：

(1) 特征檢測與生成 利用高斯卷積建立三維影像尺度空間，并通過差值處理生成高斯差分尺度空間，尺度空間生成示意圖如圖1 所示。在尺度空間中，影像共分為M 階，每階S+3 層，生成第一階影像的高斯核分別為σ,kσ,…,kS+2σ。下一階第一層影像為降采樣上一階第S 層影像的1/8，每階以此類推。差分尺度空間共分為M 階，每階S+2 層，在差分尺度空間中比較每個體素點與相鄰體素點、相鄰尺度層對應體素點及其相鄰體素點的灰度值大小，選取候選極值點。對候選極值點進行精確定位后，再篩選除去低對比度點和邊緣響應點，得到真正的特征點。

----第二階第一階高斯尺度空間高斯差分尺度空間σkσk2σkσk2σk3σk2σkσkσσ圖1 影像尺度空間生產示意圖每一個確定的特征點需要生成相應的特征向量，如圖2 所示。在三維空間中，以特征點為中心劃分16×16×16 的區域。

然后將區域再劃分為4×4×4 的子區域，計算體素的梯度方向直方圖。通過直方圖統計梯度方向，使每個子區域形成一個種子點，最終生成4096 維的特征向量。

(2) 特征向量匹配 通過特征向量匹配建立影像間特征點對應關系。比較一幅影像的特征向量u 與另一幅影像所有特征向量的歐氏距離，當另一幅影像的特征向量v 與u 的歐氏距離最近鄰時，最近鄰距離記為d(u,v)。同時，特征向量u 與v′之間距離為次近鄰時，次近鄰距離記為d(u,v′)。通過將特征點間的最近鄰距離與次近鄰距離的比值d(u,v)/d(u,v′)和設定的閾值T 作比較。當比值小于設定的閾值時，則特征向量v 是u 的最佳匹配。

(3) 薄板樣條變換 建立特征點的對應關系后，需要利用點集配準求解兩幅影像之間最優薄板樣條變換參數。通過薄板樣條插值計算得到變換的影像，完成對待配準影像的全局非剛性配準。

圖2 特征向量生成示意圖1.2 Diffeomorphic Demons 算法Demons 算法[9]是基于影像灰度的配準方法，根據光流場理論將待配準影像看作為連續運動影像序列中的兩幀，然后計算變形場大小，從而實現影像的配準[10,11]。為了使配準后的影像內部結構的連接性和鄰接關系與配準前影像保持一致，不會出現拓撲結構破壞的現象，配準過程中所求解的變形場應具有一定的拓撲保持性。而Diffeomorphic Demons 算法[12]的空間變換具有微分同胚性，通過李群復合原理替代Demons 算法中變形場的更新方式。李群是一種光滑流形也是群的數學空間，每個李群G 對應一個李代數g，李代數可以刻畫李群單位元處的局部形狀，并通過指數映射提升李代數性質至李群，從而能夠得到具有微分同胚性的形變向量，保證最終配準后影像具有拓撲保持性。

Diffeomorphic Demons 算法通過最小化能量函數以求得最終的形變向量。而在能量函數中以相似性測度均方誤差作為其一部分，所以能量函數表達式如式(1)所示。

22 22 21( , ) Dist( , ) Reg( ) ix tE c t F M c c t t     (1)其中，c 為算法當前迭代時未規則化變形場，t 為算法前一次迭代規則化變形場， Reg(t) 為限制 t 的規則項，Dist(c, t)  c  t ，表示距離。

待配準影像在三維空間p 處的灰度差為( ) p f t ，灰度差表達式如式(2)所示。

f p (t)  F( p) M t( p) (2)根據李群指數映射，exp( ) d 是從李代數g 中0 鄰域到李群G 中e 鄰域的微分同胚映射，所以當前迭代時未規則化變形場更新表達式如式(3)所示。

c  t exp(d) (3)根據公式(1)、(2)、(3)及Sherman-Morrison 公式并最小化能量函數可以得到具有微分同胚性的形變向量，如式(4)所示。

222( ) ( )( )( )TpipxF p M t pd p JpJ(4)算法中的參數有ε、x  、fluid  、diff  。exp( ) d 的計算精度由參數ε控制，ε過小時會增加計算量，文中ε取值為0.125。x  可以影響()dp 的大小，與最大迭代步長有關，避免影像不同分辨率對計算()dp 的影響。fluid  、diff  大小為影像分辨率值。Diffeomorphic Demons 算法實現過程如下：

1）設定初始變形場t 為零向量場，算法的迭代次數為N；2）根據公式(4)計算每一體素的形變向量d(p)，然后更新位移場d；3）對位移場d 高斯卷積；4）根據單參數子群性質和scaling and squaring 算法計算出exp( ) d 后，再根據公式(3)更新c；5）對c 高斯卷積，令變形場t=c。若算法迭代次數為N，則停止迭代并輸出t；反之，算法返回至步驟2)繼續迭代。

1.3 配準方案設計實驗所用數據為一個呼吸運動周期內十個時相的肺部CT影像，選擇其中的最大吸氣相和最大呼氣相。首先，通過肺部呼吸運動規律大致確定所需影像,初步選取第10_P 和20_P 時相作為最大吸氣相，第50_P 和60_P 時相作為最大呼氣相；其次，根據均方誤差最終確定實驗所需的時相影像。從表1 中可見，第10_P 和60_P 時相之間的均方誤差值最大，所以選取第10_P時相作為最大吸氣相，第60_P 時相作為最大呼氣相。各時相間均方誤差值如表1 所示。

提出的配準方案采用混合層次配準模型實現對三維肺部醫學影像配準。首先，方案中分別以最大吸氣相作為參考影像和最大呼氣相作為浮動影像，使用尺度不變特征變換算法提取影像特征點后實現薄板樣條變換，從而完成全局非剛性配準，其中，M=4,S=2；其次，全局配準后采用Diffeomorphic Demons算法對影像局部區域優化配準。算法采用多分辨率策略，提高算法的運行速度。在粗糙分辨率層中迭代時間短，再根據迭代結果初始化高分辨率層，使高分辨率層也能在較短時間內完成迭代，從而快速實現肺部影像精確配準。實驗中設置多分辨率層數為3 層，最大步長為2，每層的迭代次數分別為20、15、10。配準方案流程圖如圖3 所示。

讀取體數據特征檢測與生成特征向量匹配點集配準計算變換參數薄板樣條變換插值生成配準后影像對全局配準后影像進行多分辨率處理全局配準輸出最終配準影像局部配準滿足優化條件或達到迭代次數是否達到最大分辨率層數是是初始化變形場計算形變向量并對向量高斯卷積對形變向量進行指數映射并更新變形場否否圖3 配準方案流程圖2 配準實驗結果2.1 全局非剛性配準通過實驗驗證提出的方案對人體肺部影像配準的有效性。

實驗配置環境為：32 位Windows7 操作系統，Microsoft VisualStudio 2008 開發環境。實驗中使用原始基于光流場模型的Demons 算法與提出的配準方案進行對比，分別采用均方誤差、相關系數、Mattes 互信息作為配準結果的相似性評價指標。相似性評價指標公式分別如式(5-7)所示。

211( ( ) )Ni iiMSE F T MN    (5)12 21 1( )( ( ) ( ))( ) ( ( ) ( ))Ni iiN Ni ii iF F T M T MCCF F T M T M    (6)( , )( , ) - ( , ) log( ) ( ) x Yp x yS x y P x yP x P y(7)其中，N 為體素總數， i F 為參考影像， ( )i T M 為配準后的浮動影像, F 和 ( ) T M 分別為參考影像和配準后浮動影像中各體素點的平均灰度值，p(x, y)為影像的聯合概率密度，p(x)和 p(y) 分別為邊緣概率密度。

表1 三維肺部時序影像各時相間均方誤差值10_P 20_P 30_P 40_P 50_P 60_P 70_P 80_P 90_P00_P 1685.9 3769.45 8381.81 15424.4 21016.1 23224.2 18614.4 10837.5 3187.4610_P 1821.14 8361.07 16855.8 23170.2 25835.3 21786 14529.9 6795.8120_P 3978.05 12099 18777.5 21813.4 18518 12590.8 6820.8930_P 3921.45 10244.6 13960.4 12476.6 9448.11 7558.8540_P 3105.83 7152.06 8252.98 9231.61 11529.8配準前肺部二維剖面如圖4-圖5 所示，其中圖a 為參考影像，圖b 為浮動影像。全局非剛性配準后肺部二維剖面如圖6-圖7 所示,其中圖a 為參考影像，圖b 為全局配準后浮動影像。

分別使用原始光流場模型和提出的方案對肺部進行配準，其中全局非剛性配準后重采樣影像與參考影像相似性評價參數見表2。

a b a b圖4 配準前三維肺部影像第56 層橫切面 圖5 配準前三維肺部影像第119 層冠狀面a b a b圖6 全局配準后三維肺部影像第56 層橫切面 圖7 全局配準后三維肺部影像第119 層冠狀面表2 全局非剛性配準后結果評價參數對比均方誤差 相關系數 Mattes 互信息配準前 25835.3 0.954453 -0.986138全局配準后 11790.9 0.976352 -1.07627由圖3~6 以及表2 可知，三維肺部影像經全局非剛性配準后均方誤差下降為11790.9，下降率為54.36%，相關系數由0.954453 增大為0.976352，Mattes 互信息絕對值由0.986138 上升為1.07627。說明參考影像與配準后影像相似度提高，但肺部結構復雜，影像間還有一定差距，在肺的內部細節不是很好，需要進一步處理。

2.2 局部非剛性配準全局配準后并不能完全滿足其對精確性的要求，所以應繼續對影像進行局部配準。全局配準的肺部浮動影像經局部非剛性配準后,需重采樣得到所需的肺部體數據。肺部體數據的二維剖面如圖8-圖9 所示，其中圖a 為參考影像剖面圖，圖b 為局部非剛性配準后影像剖面圖。局部非剛性配準后與參考影像相似性評價參數對比如表3。

a b a b圖8 配準后三維肺部影像第56 層橫切面 圖9 配準后三維肺部影像第119 層冠狀面對比表3 中配準結果評價參數可知：提出的方案最終結果的精確度有顯著提高。配準結果的均方誤差值由配準前25835.3降為3726.31，下降率為85.58%，相關系數由0.954453 增大為0.992377，Mattes 互信息值由-0.986138 降為-1.38939，絕對值增大，說明參考影像與配準后影像間相似度增高，相關性增強，信息相關程度增高。傳統基于光流模型的配準算法均方誤差下降率為78.25%，相關系數變為0.98859，Mattes 互信息下降為-1.37972。相比傳統算法，提出的肺部配準方案能夠有效解決肺部因呼吸運動造成的大形變問題，配準后肺部浮動影像外部輪廓完整，復雜內部結構清晰，相似性評價參數更加精確。3 結束語由于相同個體在不同呼吸狀態下的肺部影像存在較大的形變，并且肺部內部組織結構復雜。為了能夠有效地對不同呼吸狀態下的肺部進行配準，提出了一種精確度更高、魯棒性較強的層次模型配準方案。該配準方案選擇肺部三維時序影像作為實驗數據，這樣能保證方案對其他時相配準也有效。首先對肺部三維影像進行全局配準，從整體上模擬肺部形態變化，再對其局部區域進行優化，從而實現肺部較大形變的配準。同時，算法中引進微分同胚變換，保證求解的變換具有光滑、可逆的性質，使得影像拓撲結構在配準前后保持一致。通過實驗對比發現，相比傳統算法該方案配準精度較高，能夠有效配準肺部邊緣處及內部細節，具有較高的臨床應用價值。

參考文獻

[1] 袁崢璽, 聶生東. 圖像配準技術在圖像引導放療中的應用[J]. 中國醫學物理學雜志, 2012, 29(5): 3628-3631, 3637.

[2] LORENZI M, AYACHE N, FRISONI G B, et al. LCC-Demons: a robustand accurate symmetric diffeomorphic registration algorithm[J].NeuroImage, 2013, 1(81): 470-483.

[3] 周穎玥, 馮煥清, 李傳富. 同體呼氣相和吸氣相肺部高分辨率CT 體積圖像的配準[J]. 航天醫學與醫學工程, 2009, 22(3): 198-204.

[4] 呂曉琪, 張傳亭, 任國印, 等. 基于圖形處理器的肺部及病灶區域四維可視化技術[J]. 中國醫學影像技術, 2013, 29(11): 1901-1905.

[5] 王安娜, 呂丹, 王哲, 等. 基于SIFT 特征提取的非剛性醫學圖像配準算法研究[J]. 生物醫學工程學雜志, 2010, 27(4): 763-768, 784.

[6] 王婕妤, 王加俊, 張靜亞. 基于改進光流場和尺度不變特征變換的非剛性醫學圖像配準[J]. 電子與信息學報, 2013, 35(5): 1222-1228.

[7] 汪道寅. 基于SIFT圖像配準算法的研究[D]. 合肥: 中國科學技術大學,2011.

[8] Zhang Rui, Zhou Wu, Li Yanjie, et al. Nonrigid Registration of Lung CTImages Based on Tissue Features[J]. Computational and MathematicalMethods In Medicine, 2013, 2013(2013): 1-7.

[9] Suh J W, Scheinost D, Dione D P, et al. A non-rigid registration method forserial lower extremity hybrid SPECT/CT imaging[J]. Medical ImageAnalysis, 2011, 15(1): 96-111.

[10] 周志勇. 醫學圖像非剛性配準方法研究[D]. 長春: 中國科學院研究生院(長春光學精密機械與物理研究所), 2013.

[11] Lu Xiaoqi, Zhao Yongjie, Zhang Baohua, et al. A non-rigid cardiac imageregistration method based on an optical flow model[J].Optik-International Journal for Light and Electron Optics, 2013,124(20): 4266-4273.

[12] Santos J, Chaudhari A J, Joshi A A, et al. Non-rigid registration of serialdedicated breast CT, longitudinal dedicated breast CT and PET/CT imagesusing the diffeomorphic demons method[J]. Physica Medica, 2014, 30(6):713-717.