孫 正

（華北電力大學電子與通信工程系，河北 保定 071003）

血管內超聲(intravascular ultrasound，IVUS)成像是臨床診斷冠心病的一種基于導管的腔內成像技術，具有獨特的在活體中觀察血管壁、管腔及斑塊形態(tài)性質，甚至管壁功能狀態(tài)的特點，已經成為 X射線冠狀動脈造影(X-ray coronary angiography)的彌補影像技術。

由于心臟運動以及血壓變化所造成的運動偽像是影響IVUS圖像的視覺效果，進而影響定量分析和血管三維重建精度的主要因素。目前主要采用門控的方法減少 IVUS序列中的運動偽像，包括3種方式：

一是心電(electrocardiogram，ECG)門控的圖像采集方法：即采用 ECG門控的馬達回撤超聲導管，根據(jù) ECG信號只在每個心動周期的相同時相（一般是R波）采集圖像，改善心臟搏動所致的偽像[1]。該方法的缺點在于需要專門的ECG門控采集裝置，而目前臨床采用的多數(shù)IVUS系統(tǒng)不包含此功能。并且由于每個心動周期只采集一幀，因此與連續(xù)回撤導管的非門控采集方式相比，圖像采集時間延長至少3倍[2]。

二是離線 ECG門控方法：即連續(xù)回撤超聲導管采集圖像序列，同時記錄 ECG信號。待介入檢查過程結束后，由醫(yī)生根據(jù) ECG信號從圖像序列中選取在相同心臟時相采集的圖像。這種方式雖然不會延長介入檢查的時間，但是其不足之處在于：ECG信號是心臟的生理電特性通過骨骼、肌肉和皮膚等反映在體表的信號，而 IVUS圖像是通過插入血管腔內的超聲探頭采集的，兩者不可能完全同步；門控結果的客觀性和可重復性差，受操作者的臨床經驗和主觀因素影響較大；雖然多采用 0%點即舒張末期作為采樣點，但無法確定采用該點是否可以獲得最穩(wěn)定的門控圖像序列；目前尚無公認、明確的標準，來確定采用心動周期中的哪一時相作為最佳采樣點，可獲得最大幀間穩(wěn)定度[3]。

三是基于圖像的離線門控方法：基于 ECG門控技術在抑制運動偽像方面的不足，近年來，運用數(shù)字圖像處理技術，從連續(xù)回撤超聲導管采集的IVUS圖像序列中提取出隱含的心動周期信息，并據(jù)此從每個心動周期中選取一幀圖像，組成同步采樣序列，即基于圖像的門控(image-based gating)技術[2-8]，逐漸成為該領域的研究熱點。這種方法雖然可以避免離線 ECG門控方法的缺點，但是由于每個心動周期只選擇一幀，需要拋棄大量幀，因而可能會丟失很多有診斷價值的信息。

目前，對非門控IVUS圖像序列直接進行運動估計的方法包括光流法(optical flow，OF)[9]和塊匹配法(block matching，BM)[10]。此類方法既不需要門控圖像采集裝置，無需記錄ECG信號，也不需要拋棄有用幀，可保證圖像數(shù)據(jù)集合的完整性。但是光流法的估計結果極易受到血流隨機運動的干擾，且對噪聲十分敏感；塊匹配法需要明顯的標志物（例如鈣化，血管分叉等），而此類標志并非在每幀中都存在。

本文在綜合分析IVUS圖像序列中由周期性心臟運動所致運動偽像的產生機制和表現(xiàn)形式的基礎上，對覆蓋多個心動周期的非門控 IVUS圖像序列直接進行剛性運動偽像的定量估計和補償。首先建立橫向視圖中的血管壁在心動周期中的剛性運動模型；然后定量估計相鄰幀間的血管橫斷面剛性運動參數(shù)（包括位移和旋轉角），并采用譜分析的方法分離運動分量與幾何分量；最后通過補償運動分量，達到抑制運動偽像的目的。無需利用 ECG信號，也不需要拋棄任何一幀圖像，在抑制運動偽像的同時保證圖像數(shù)據(jù)集合的完整性。

1 方 法

1.1 運動偽像的產生機制和表現(xiàn)形式

IVUS圖像序列中由周期性心臟運動和搏動血流所致的運動偽像的主要表現(xiàn)形式為：

1） 幀間的錯位，即橫斷面圖像序列（圖1(a)）中相鄰幀之間血管橫截面的平移和旋轉（圖1(b)）；

2） 在沿管腔長軸方向的縱向視圖中，血管壁邊緣呈現(xiàn)鋸齒形（圖1(c)）。

圖2說明了運動偽像的產生機制。假設對一段均勻直徑的直血管段進行IVUS成像，在回撤超聲導管的過程中，導管始終位于各幀斷層圖像的中心。若無運動偽像，且假設導管回撤路徑與管腔軸線重合，則將各幀圖像按照采集順序依次排列得到的縱向視圖中，血管壁的上下輪廓應近似為直線（圖 2(a)）。若血管處于周期性運動狀態(tài)，則可能導致導管和管腔之間的相對運動。因而在回撤導管的過程中，采集到的各幀斷層圖像記錄的并非血管的真實狀態(tài)，在縱向視圖中將出現(xiàn)“鋸齒樣”的血管壁上下輪廓，如圖2(b)所示。這種偽像不僅會影響對血管腔和斑塊形態(tài)的觀察，同時也會給形態(tài)參數(shù)的準確測量和血管的三 維重建帶來很大困難。

圖1 IVUS圖像序列

圖2 IVUS圖像序列示意圖

1.2 建立血管壁剛性運動模型

IVUS斷層圖像序列中相鄰切片之間管腔橫截面重心的位移和空間方向的改變主要由運動和幾何兩方面的因素造成[3,11-13]：

1） 運動：即外部因素，指由周期性心臟運動和搏動的血流造成的超聲導管相對于管腔的運動和血管形態(tài)的變化；

2） 幾何：即內部因素，指血管腔本身不規(guī)則的幾何形狀。

從各幀斷層圖像中分割出血管腔邊界之后，相鄰幀之間管腔橫截面的平移和旋轉可分別用管腔邊界曲線重心的位移(,)xyΔΔ和邊界曲線之間的旋轉角αΔ來表示，它們分別由兩部分組成

其中，腳標d和g分別表示運動分量和幾何分量。

如圖3所示，設相鄰兩幀中的管腔邊界曲線分別為1γ和2γ，其重心分別為，在運動和幾何兩方面因素的共同作用下，1γ和2γ之間的位移為，1γ和2γ之間的旋轉角為Δα，旋轉中心為導管(x0,y0)，即

其中，Δx、Δy和Δα如式(1)所示分別由運動和幾何分量兩部分組成。對于圖像序列中的各相鄰幀，本文通過計算出位移(Δx,Δy)和旋轉角Δα，并從中分離出Δxd、Δyd和Δαd，完成對運動偽像的補償。

圖3 相鄰幀IVUS圖像運動參數(shù)示意圖

1.3 運動偽像的估計和補償

1.3.1 提取IVUS圖像中的血管壁邊界

對IVUS圖像進行分割，提取血管內腔邊界和血管壁中-外膜邊界，不僅是本文抑制運動偽像方法的重要步驟，而且也是量化分析[14]、三維重建[15]以及 IVUS與其它冠狀動脈圖像融合[16]的基礎和保證其精度的關鍵。近年來，提出了很多計算機輔助的自動分割方法[17-18]。考慮到處理效率等因素，本文采用一種三維分割方法[19]，與二維分割方法是對圖像序列進行逐幀串行處理相比，三維分割技術可實現(xiàn)對圖像序列的并行處理，提高處理效率。

1.3.2 計算運動參數(shù)

對于圖像序列中的各相鄰幀，例如第k-1和k幀（k=2,3,…,M，M是IVUS圖像序列的總幀數(shù)），計算血管內腔邊界曲線的幾何中心作為對其重心的近似，并如§1.2中所述方法計算出重心之間的位移量和旋轉角和kαΔ包含運動分量和幾何分量)兩部分，需從中分離出運動分量,dkxΔ、,dkyΔ和,dkαΔ。

1.3.3 分離運動分量和幾何分量

本文對計算出的kxΔ、kyΔ和kαΔ（k=2,3,…,M）分別進行傅立葉變換，得到其幅度譜ΔXk、ΔYk和ΔAk。由于由血管本身的不規(guī)則幾何形狀所引起的相鄰幀之間管腔橫截面空間方向和重心位置的變化速度，遠小于由周期性心臟運動所致的變化速度，因此在ΔXk、ΔYk和ΔAk中，高頻分量對應于剛性運動分量，且其變化頻率應等于心率，而低頻分量對應于幾何分量 ( Δxk,g, Δyk,g, Δαk,g)。分別對ΔXk、ΔYk和ΔAk進行高通濾波，則濾波器的輸出即是剛性運動分量的幅度譜(ΔXk,d, ΔYk,d, Δ Ak,d)，對其進行逆傅立葉變換即可得到 ( Δxk,d, Δyk,d, Δαk,d)的估計值。

由于運動分量 ( Δxk,d,Δyk,d, Δαk,d)主要由周期性心臟運動引起，因此本文將高通濾波器的通帶截止頻率設定為病人的心率值（單位：次/秒，即Hz）。采用文獻[6]中的方法從IVUS圖像序列中估計平均心率的近似值，原理如下：

首先，對IVUS圖像序列進行逐幀比較，通過計算兩幀圖像之間灰度值的歸一化互相關(NCC)

計算各幀圖像之間的差異值

式中，i,j= 1 ,2,… ,M；Ii和Ij分別表示IVUS圖像序列中的第i幀和第j幀，其大小均為N1×N2像素，平均灰度值分別為μi和μj；M是圖像序列的總幀數(shù)。

然后，計算間隔為k幀的兩幀圖像Im和Im+k（m= 1,2,… ,M-k）的平均差異值

得到曲線～k，k=0,1,2,… ,M-1，且= 0 。該曲線具有近似周期的形狀，這是由心臟的周期性運動所造成的，該曲線的重復頻率，即其幅度譜曲線峰值所對應的頻率即為病人心率的近似值R（單位：次/分鐘，即Hz）。

1.3.4 補償運動偽像

2 實驗結果與分析

為了驗證方法的可行性并估計精度，本文分別利用計算機模擬的血管腔輪廓序列和臨床采集的IVUS圖像序列進行實驗，并對實驗結果進行分析。

2.1 模擬圖像實驗

圖4 7幀模擬管腔輪廓序列

為了評價§1.3中所述對血管腔輪廓序列剛性運動參數(shù)的估計精度，本文將用封閉曲線表示的管腔輪廓序列按照已知的位移和旋轉角進行平移和旋轉，模擬生成運動參數(shù)已知的管腔輪廓序列。圖4(a)和圖4 (b)分別為7幀模擬管腔輪廓的橫向視圖和完成重心計算的縱向排列圖，圖4(c) ～ (e)為各偏移量的柱狀圖。

圖5和圖6分別為對20幀和35幀模擬管腔輪廓曲線序列進行剛性運動補償前后的橫截面視圖和縱向視圖?？梢钥闯?，完成運動補償后的橫截面視圖中相鄰幀之間的扭動大大減小，在縱向視圖中體現(xiàn)為管腔邊界的鋸齒形尖峰與水平基線之間的距離大大縮小。

圖6 35幀模擬圖像實驗結果

2.2 臨床圖像實驗

臨床圖像序列是采用Jomed Endosonic超聲成像儀采集的，采用2.9F 30MHz單軌機械超聲探頭，由馬達驅動自動回撤探頭導管，回撤速度為0.5mm/s，當探頭以1800轉/分作360°旋轉時，可以30幀/秒的速率連續(xù)獲得血管橫軸實時切面圖像。

圖7 62幀臨床圖像實驗結果

圖7和圖8分別為62幀和80幀臨床圖像序列的運動偽像補償前后的橫向視圖（血管腔邊界曲線）及縱向視圖?？梢钥闯?，完成運動偽像抑制后的橫截面視圖中相鄰幀管腔輪廓之間的扭動大大減小，在縱向視圖中體現(xiàn)為外膜輪廓的鋸齒形邊緣變的相對平滑，因而很大程度上抑制了圖像序列中存在的運動偽像。

圖8 80幀臨床圖像實驗結果

由于目前還沒有一種客觀、定量評價抑制IVUS運動偽像方法精度的體系[6]，因此，我們將采用本文方法進行運動偽像抑制的結果和采用基于圖像的離線門控[8]結果進行比較，評價本文方法的精度。圖9是分別對120幀和200幀的臨床圖像序列進行實驗的結果，可見雖然基于圖像的離線門控法可實現(xiàn)對運動偽像的抑制，但是門控序列的長度相對于原序列而言明顯變短，即丟棄了大量的幀數(shù)據(jù)。而本文方法在完成運動偽像抑制的同時，不丟棄任何一幀數(shù)據(jù)，保證了原始圖像數(shù)據(jù)集合的完整性。

2.3 誤差源分析

本文的運動偽像抑制方法包括圖像的采集和預處理、血管腔輪廓提取、管腔輪廓重心及偏移量的計算、運動分量的濾出、對運動分量的補償?shù)炔襟E。除去計算誤差之外，誤差主要來自兩方面：

一是圖像的分割：由于本文方法是直接對分割出的管腔輪廓進行運動估計的，因而IVUS圖像的分割精度直接影響對運動偽像的補償效果。因此，進一步提高二維輪廓提取的準確性有助于提高運動偽像的補償精度。

二是濾波器的設計：由于本文采用高通濾波的方法分離由心臟的周期性運動引起的偏移量，因而高通濾波器的設計（包括原型的選擇和參數(shù)的設定）直接影響對運動偽像的抑制效果?？紤]到運算復雜度等因素，本文選擇了 Butterworth濾波器，并將通帶截止頻率設定為從IVUS圖像序列中估計出來的平均心率值。對于濾波器原型的選擇，文獻[3]首次探討了濾波器的通帶形狀對采樣精度的影響，并提出采用統(tǒng)計的方法確定哪種濾波器更適合于心臟時相信息的提取。在后續(xù)工作中，我們將著重分析濾波器的性能對于運動偽像抑制效果的影響。對于平均心率值的估算，若圖像采集過程中同時記錄了 ECG信號，則可直接利用ECG記錄的心率值；若無ECG信號可供利用，則采用文獻[6]的方法，在病人心率保持恒定的情況下，可得到精確的估算結果。但是若病人的心率發(fā)生波動，則其估計結果存在一定誤差，詳細討論參見文獻[6]。

圖9 臨床IVUS圖像序列運動偽像抑制實驗結果

3 結 論

本文提出了一種抑制冠狀動脈內超聲圖像序列中由心臟運動所致運動偽像的方法，與現(xiàn)有的心電門控方法（包括聯(lián)機和脫機兩種方式）相比較，既不需要專門的心電門控圖像采集裝置，也不需要同步記錄心電信號，而是運用圖像分析技術，對連續(xù)回撤超聲導管采集到的IVUS圖像序列數(shù)據(jù)進行處理，因而適用于沒有記錄 ECG信號的IVUS記錄序列。與基于圖像的門控方法相比，本文方法不需要拋棄圖像序列中的有用幀，可在抑制運動偽像的同時，保證圖像數(shù)據(jù)集合的完整性。

[1]胡廣祿, 陳興新. 心電門控技術及其應用[J]. 醫(yī)療裝備, 2000, 13(1): 17.

[2]de Winter S A, Hamers R, Degertekin M, et al. A novel retrospective gating method for intracoronary ultrasound images based on image properties [C]//Proceedings of IEEE International Conference on Computers in Cardiology, Kassandra Chalkidiki,Greece: IEEE Press, 2003, 30: 13-16.

[3]Hernàndez-Sabaté A, Gil D, Garcia-Barnés J, et al.Image-based cardiac phase retrieval in intravascular ultrasound sequences [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control,2011, 58(1): 60-72.

[4]Zhu H, Oakeson K D, Friedman M H. Retrieval of cardiac phase from IVUS sequences [C]//Proceedings of SPIE Conference on Medical Imaging 2003:Ultrasonic Imaging and Signal Processing, Bellingham WA: SPIE Press, 2003, 5035: 135-146.

[5]Nadkarni S K, Boughner D R, Fenster A. Image-based cardiac gating for three-dimensional intravascular ultrasound [J]. Ultrasound in Medicine and Biology,2005, 31(1): 53- 63.

[6]O’Malley S M, Granada J F, Carlier S, et al.Image-based gating of intravascular ultrasound pullback sequences [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3):299-306.

[7]Carlo G, Oriol P, Oriol R L, et al. Robust image-based IVUS pullbacks gating [C]//Proceedings of 11th International Conference on Medical Image Computing and Computer-Assisted Intervention Part II(MICCAI 2008), New York City, USA: MICCAI Society Press, 2008, 5242: 518-525.

[8]孫 正. 抑制冠狀動脈內超聲圖像序列運動偽影的離線門控方法[J]. 光電子·激光. 2010, 21(4):632-638.

[9]Danilouchkine M G, Mastik F, Steen A F W. Accuracy in prediction of catheter rotation in IVUS with feature-based optical flow——a phantom study [J].IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 356-365.

[10]Danilouchkine M G, Mastik F, Steen A F W. Two methods for catheter motion correction in IVUS palpography [C]// Proceedings of 2007 IEEE International Ultrasonic Symposium (IUS), New York NY: IEEE Press, 2007: 1724-1727.

[11]Nadkarni S K, Mills G, Boughner D R, et al. A pulsating coronary vessel phantom for two-and threedimensional intravascular ultrasound studies [C]//Proceedings of the 4th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2001), Utrecht, Netherlands:MICCAI Society Press, 2001, 2208: 1164-1165.

[12]Holzapfel G, Gasser T, Stadler M. A structural model for the viscoelastic behavior of arterial walls:continuum formulation and finite element analysis [J].European Journal of Mechanics A-Solids, 2002, 21:441-463.

[13]Rosales M, Radeva P, Rodriguez-Leor O, et al.Modeling of image-catheter motion for 3-D IVUS [J].Medical Image Analysis, 2009, 13(1): 91-104.

[14]Baldewsing R A, Danilouchkine M G, Mastik F, et al.An inverse method for imaging the local elasticity of atherosclerotic coronary plaques [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 277-289.

[15]Paul S. 3-D Intravascular ultrasound (IVUS) and IVUS-Palpography: insights into the mechanical behavior of the coronary vessel wall [J]. The International Journal of Cardiovascular Imaging,2006, 22: 153-155.

[16]Sun Zheng. Three-dimensional reconstruction of vessels from intravascular ultrasound sequence and X-ray angiograms [C]//Proceedings of 2nd International Conference on BioMedical Engineering and Informatics (BMEI 2009), Tianjin,China: IEEE Press, 2009, 1: 444-448.

[17]Bovenkamp E G P, Dijkstra J, Bosch J G,et al.User-agent cooperation in multiagent IVUS image segmentation [J]. IEEE Transactions on Medical Imaging, 2009, 28(1): 94-105.

[18]Zhang Qi, Wang Yuanyuan, Wang Weiqi, et al.Automatic segmentation of calcifications in intravascular ultrasound images using snakes and the contourlet transform [J]. Ultrasound in Medicine &Biology, 2010, 36(1): 111-129.

[19]孫 正, 楊 宇. 一種基于snake模型的IVUS圖像序列三維分割方法[J]. 工程圖學學報, 2011, 22(6):25-32.