基于snake模型的IVUS圖像序列三維分割方法

孫 正， 楊 宇

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

近年來，血管內(nèi)超聲(Intravascular Ultrasound,IVUS)成像技術(shù)已經(jīng)越來越廣泛地被應(yīng)用到血管疾病，特別是冠心病的臨床診斷和介入治療中。它不僅能夠反映血管內(nèi)腔的變化，同時也能反映含斑塊在內(nèi)的管腔橫斷面結(jié)構(gòu)、血管壁厚度、形態(tài)以及斑塊成分等[1]。血管壁在IVUS圖像上呈現(xiàn)回聲不同的三層環(huán)形結(jié)構(gòu)：內(nèi)層為一強回聲亮環(huán)，在組織學(xué)上由內(nèi)膜及內(nèi)彈性膜組成；中層相當于中膜，呈現(xiàn)為低回聲暗區(qū)；外層為較亮的強回聲帶，由外彈性膜和外膜構(gòu)成。圖像序列中的前后幀之間非常相似，具有很強的相關(guān)性。

對IVUS圖像進行分割，提取出血管內(nèi)腔邊界和中-外膜邊界（包括可能存在的斑塊邊緣），是定量分析（包括測量血管腔徑、容積和斑塊體積、估計狹窄百分比、測量血流動力學(xué)參數(shù)等）[2]和血管三維重建[3]的重要步驟，分割質(zhì)量直接決定定量測量和三維重建的精度。目前臨床常采用手動分割，由富有經(jīng)驗的醫(yī)生用鼠標在屏幕上逐幀勾畫。由于IVUS圖像序列常包括上千幀，因而該方法不但耗時，工作量大，而且可重復(fù)性差，分割結(jié)果在很大程度上取決于操作者的經(jīng)驗和主觀因素。若采用數(shù)字圖像處理的方法進行自動提取，則由于IVUS圖像一般受噪聲污染嚴重（包括斑點噪聲、圖像偽影和部分血管壁鈣化陰影等的影響），常規(guī)的灰度圖像輪廓提取方法，如微分法等，難以獲得滿意結(jié)果。已經(jīng)提出的半自動提取方法主要包括圖搜索[4]、活動輪廓(snake)模型[5]及其相關(guān)改進算法[6]等。這些方法的分割效果較常規(guī)方法有很大提高，但仍需操作者一定程度的手動參與，從而延長處理時間或引入誤差。全自動的分割方法目前已成為該領(lǐng)域的研究熱點，例如，Brusseau等[7]的基于統(tǒng)計活動輪廓模型的方法、Filho等[8]的基于模糊聚類的方法、Giannoglou等[9]的自動獲得初始輪廓的 snake方法、Ruiz等[10]的概率統(tǒng)計分割內(nèi)腔方法、Papadogiorgaki 等[11]的基于紋理分析的方法、Bucher等[12]的形狀驅(qū)動方法、Bovenkamp等[13]的多智能體方法、Zhang等[14]的基于 contourlet變換的snake方法等。上述方法都是在二維空間中進行串行處理，即對IVUS序列中的各幀圖像進行逐幀分割。近年來，三維分割技術(shù)[15]成為該領(lǐng)域新的研究方向。三維分割與二維分割相比具有以下優(yōu)點：第一，可以利用整個序列的信息；第二，可同時完成對各幀圖像的分割，實現(xiàn)對整個序列的并行處理，從而大大提高處理效率，縮短處理時間。

本文提出一種基于snake模型的IVUS圖像序列三維分割方法。在對原始圖像進行濾波去噪和去除環(huán)暈偽像的預(yù)處理的前提下，首先獲取IVUS序列的四個縱向視圖，并從中提取出內(nèi)腔和中-外膜邊界；然后，將各縱向視圖中的邊界曲線映射到各IVUS幀中，得到各橫向視圖中的初始輪廓；將此初始輪廓作為snake模型的初始形狀，通過使預(yù)先設(shè)定的能量函數(shù)最小，模型在內(nèi)外力的共同作用下不斷變形，最終得到各IVUS幀中的內(nèi)腔邊界和中-外膜邊界，從而實現(xiàn)對各幀圖像的并行分割。

1 方 法

如圖1所示，本文的分割方法包括以下步驟：① 原始圖像的預(yù)處理，包括濾波去噪和去除環(huán)暈偽像；② 獲取 IVUS圖像序列的四個縱向視圖；③ 提取縱向視圖中的內(nèi)腔邊界和中-外膜邊界；④ 將步驟③中獲得的邊界曲線映射到各幀IVUS切片圖像中，得到各 IVUS幀的初始內(nèi)腔邊界和中-外膜邊界；⑤ 在各IVUS幀中，初始輪廓演化變形，得到最終的內(nèi)腔和中-外膜輪廓。

圖1 IVUS圖像序列三維分割方法流程圖

1.1 圖像預(yù)處理

首先，采用中值濾波和高斯平滑兩種通用預(yù)處理方法，減少IVUS圖像中的椒鹽噪聲和隨機噪聲。

在IVUS成像過程中，超聲導(dǎo)管在圖像中心區(qū)域產(chǎn)生一個無回聲暗區(qū)，圍繞導(dǎo)管黑色影像的厚薄不一的白色環(huán)狀影像稱為環(huán)暈偽像(ring-down artifact)[16]，即圖2中箭頭所指。它與傳感器聲波振蕩引起的高振幅信號有關(guān)，不是導(dǎo)管本身的影像，應(yīng)該在預(yù)處理階段加以去除，否則會影響圖像分割的準確度。根據(jù)環(huán)暈偽像在IVUS極坐標視圖中的分布特性，本文采用極坐標變換法對其加以去除。

圖2 環(huán)暈偽像在兩種視圖上表示

首先，對IVUS圖像進行極坐標變換，得到其極坐標視圖，如圖2(b)所示?？梢姡h(huán)暈偽像固定的位于極坐標視圖的上部。采用閾值法，按照下式去除極坐標視圖中的環(huán)暈偽像

其中 I( r,θ)和I'(r,θ)分別為原極坐標視圖和去除環(huán)暈偽像后的極坐標視圖中像素點(r,θ)處的灰度值；Icatheter是原 IVUS圖像中導(dǎo)管區(qū)域的像素灰度值；r為像素點的極徑；t=ε×ImageHeight 2；ε為權(quán)重參數(shù)，經(jīng)實驗確定其取值范圍為[0.1,0.35]。在極坐標視圖中消除環(huán)暈偽像后，再經(jīng)過極坐標逆變換，即可得到直角坐標系下去除環(huán)暈偽像后的IVUS圖像。

1.2 獲得IVUS序列縱向視圖

采集IVUS圖像時，隨著超聲導(dǎo)管的緩慢回撤，獲得的是一系列血管橫截面切片圖像。IVUS縱向視圖是指IVUS圖像序列沿平行于導(dǎo)管方向的縱向截面視圖。如圖3所示，本文沿血管長軸（縱向）方向取四個縱向視圖，即垂直切面A、水平切面B、左對角線切面C和右對角線切面D。

圖3 獲取IVUS縱向視圖

1.3 縱向視圖中血管壁邊緣的提取

如圖4所示，IVUS縱向視圖中的每一行像素即為一幀切片圖像與該縱向視圖的交線。本文根據(jù)縱向視圖中的內(nèi)腔邊界和中-外膜邊界的灰度同其周圍像素灰度相差較大的特征，采用灰度閾值法進行提取。算法的基本原理是：遍歷縱向視圖中的各像素，判斷識別目標邊界點。由于縱向視圖中的目標邊界分布于圖像中部，所以從圖像中軸線開始，分別向左和向右進行逐行遍歷，如圖 5所示。對于當前像素(i, j)，其灰度值為I( i, j)，若 I ( i, j + 1 )- I ( i, j) ≥η ，η為閾值，則為目標邊界點，否則不是。將每行左右兩部分的像素中，第一個符合上述條件的像素記為內(nèi)腔邊界點，第二個記為中-外膜邊界點。經(jīng)實驗確定η的取值區(qū)間為[10, 20]。

圖4 縱向與橫向視圖的對應(yīng)關(guān)系示意圖

圖5 逐行遍歷縱向視圖像素

1.4 橫向視圖中初始輪廓的獲取

縱向視圖的輪廓線對應(yīng)到橫向視圖中為輪廓點，四條縱向視圖輪廓線映射到橫向視圖中為四個輪廓點，如圖6所示。在得到各幀橫向視圖中的內(nèi)腔和中-外膜初始邊界點以后，采用DDA(Digital Differential Analyzer)算法[17]依次連接各初始邊界點，如圖7所示，同時獲得各幀圖像中的初始內(nèi)腔邊界和中-外膜邊界。

1.5 基于snake的初始輪廓演化

snake模型是Kass等[5]在1989年提出的一個能量最小化樣條模型，它通過將由圖像數(shù)據(jù)直接獲得的約束力和關(guān)于目標的先驗知識結(jié)合起來，完成對圖像中目標的分割、匹配和跟蹤[18]。本文將1.4節(jié)中獲得的各幀IVUS圖像中血管壁的初始輪廓作為snake模型的初始形狀，通過使預(yù)先設(shè)定的能量函數(shù)最小，snake不斷變形，當能量函數(shù)取得全局最小值時，snake模型即停留在目標輪廓處，從而完成對各幀圖像中內(nèi)腔和中-外膜邊界的并行提取。

圖6 縱向與橫向視圖輪廓對應(yīng)關(guān)系示意圖

圖7 由初始點獲取初始輪廓

將初始輪廓離散成由 N個點組成的有序點集，則snake能量函數(shù)的離散表達式為

其中 內(nèi)部能量Eint使 snake在變形過程中保持連續(xù)和平滑，其歸一化表達式為

其中 ci( xi, yi)(i=1,2,…, N-2)是第 i個 snake點；和maxd分別是相鄰snake點之間的平均距離和最大距離，在每次迭代結(jié)束時，它們的值都被更新；α和β是權(quán)重參數(shù)。式(3)中的第一項保證變形過程中snake的連續(xù)性，使snake點均勻分布，不致產(chǎn)生收縮的現(xiàn)象；第二項是曲線二階導(dǎo)數(shù)的離散形式，保證變形過程中snake的光滑性。由于進行了歸一化，故式(3)中兩項的取值范圍都在[0, 1]區(qū)間內(nèi)，權(quán)重α和β的取值區(qū)間也是[0, 1]。

外部能量Eext是保證snake收斂的外部力，它決定snake的移動方向。本文將其定義為

其中 I( xi, yi)和 ? I ( xi, yi)分別是像素(xi, yi)的灰度和灰度梯度值。由于8位灰度圖像的灰度和灰度梯度的取值范圍分別為[0, 255]和別進行了歸一化，使其取值范圍為[0, 1]區(qū)間。γ, λ ∈ [ 0,1]是權(quán)重參數(shù)。

通過使能量函數(shù)E最小，snake在內(nèi)外力的共同作用下從初始形狀開始不斷變形，最終停留在能量函數(shù)取得全局最小值的最優(yōu)位置，即為目標輪廓。能量函數(shù)的全局最優(yōu)采用Williams快速算法[19]完成，不僅計算速度快，所需存貯空間小，而且非常適合于編程實現(xiàn)。

2 實驗結(jié)果及分析

2.1 實驗結(jié)果

本文采用臨床采集的冠狀動脈IVUS圖像序列對所提出的分割方法進行了實驗驗證，在這一節(jié)中給出部分實驗結(jié)果。實驗序列是采用Jomed Endosonic超聲成像儀在臨床獲取的，超聲探頭為2.9F Jomed單軌機械探頭，頻率為30MHz。在回撤導(dǎo)管的過程中，探頭導(dǎo)管以1800轉(zhuǎn)/分作360°旋轉(zhuǎn)時可連續(xù)獲得30幀/秒的血管橫軸實時切面圖像。采用Visual C++6.0編程實現(xiàn)算法，所用計算機配置為2.6GHz Intel Pentium III Xeon Tanner processor、2 GB內(nèi)存。

圖8和圖9是對IVUS圖像去除環(huán)暈偽像的實驗結(jié)果。圖10是從一個IVUS序列中選出的7幀，均已完成了對環(huán)暈偽像的去除。圖 11是從該序列獲取的四個縱向視圖，其中提取出了血管壁的內(nèi)外膜邊緣；將其映射到各幀切片圖像上，即可得到血管壁內(nèi)外膜的初始輪廓點，如圖12(a)～(d)所示，然后采用 DDA算法用直線段依次連接各輪廓點即可得到初始輪廓，如圖12(e)～(h)所示。之后，通過使 snake能量函數(shù)最小，初始輪廓不斷變形，最終停留在使能量函數(shù)取得全局最小值的地方，即內(nèi)腔和中-外膜邊界，如圖13和圖14所示，其中式(3)和式(4)中的權(quán)重參數(shù)設(shè)定為 α=β=γ=λ=1，即各能量項的權(quán)重相等。

圖8 一幀IVUS圖像的環(huán)暈偽像去除結(jié)果(ε = 0.21)

圖9 四幀IVUS圖像的環(huán)暈偽像去除結(jié)果(ε=0.21)

圖10 去除環(huán)暈偽像的7幀IVUS圖像

圖11 完成了管壁輪廓提取的IVUS序列縱向視圖

圖12 IVUS橫向視圖中血管壁的初始輪廓 (a)～(d)初始輪廓點；(e)～(h)初始輪廓

圖13 IVUS橫向視圖中初始輪廓的變形結(jié)果1

圖14 IVUS橫向視圖中初始輪廓的變形結(jié)果2

2.2 分割精度和速度

目前還沒有一種公認的客觀、定量評價IVUS圖像分割方法精度的體系或標準，一般采用將自動分割結(jié)果與醫(yī)學(xué)專家手動分割結(jié)果進行比較的方法來衡量自動方法的分割質(zhì)量[20]。作者從臨床獲取的冠狀動脈內(nèi)超聲圖像序列中隨機選取了一個序列，由200幀組成，由醫(yī)生手動描繪出各幀圖像中的內(nèi)腔邊界和中-外膜邊界，并將其作為“金標準”。然后，分別采用本文方法和二維方法對圖像進行自動分割，并記錄處理時間。對各幀圖像，計算自動分割出的輪廓與“金標準”之間的歐氏距離作為誤差[21]，定量估計分割算法的精度。其中二維方法采用 Giannoglou等[9]提出的自動獲取初始輪廓的 snake模型法。采用兩種方法進行分割之前，均對原始圖像進行1.1節(jié)中所述的預(yù)處理，即濾波去噪和去除環(huán)暈偽像（參數(shù)設(shè)置為 ε=0.21）。處理時間均是在2.6GHz Intel Pentium III Xeon Tanner processor、2GB內(nèi)存的計算機參數(shù)下測得的，且均不包含預(yù)處理（去噪和抑制偽影）的時間。

對200幀圖像的分割誤差以及處理時間的統(tǒng)計結(jié)果如表1所列。由表可知，無論對內(nèi)腔邊界還是中-外膜邊界，本文方法的平均誤差與二維方法均十分接近。而且，對內(nèi)腔邊界的分割結(jié)果均明顯優(yōu)于中-外膜邊界。這是由于IVUS圖像中血管內(nèi)腔與內(nèi)膜的灰度值差異通常比較明顯，因而在去除環(huán)暈偽像的情況下，便可較準確地提取出內(nèi)腔邊界。而中-外膜邊界的精確提取依賴于圖像中中膜的形態(tài)及初始輪廓的位置，若中膜在圖像中表現(xiàn)為較厚的低灰度帶，且初始輪廓較接近真實輪廓，則最終分割結(jié)果較好，否則易被偽邊緣吸引[22]。對于處理時間，二維方法是進行逐幀處理的串行分割，對每幀圖像的處理時間約為0.26s，則處理200幀共需約52s。而采用本文方法，由于對各幀的分割是同時進行的，因而對200幀圖像的處理時間僅需 16.8s，從而大大提高了處理效率。大大提高處理效率，縮短處理時間。

表1 對200幀圖像的分割誤差及處理時間統(tǒng)計結(jié)果

3 結(jié) 論

對IVUS圖像進行分割，提取出管腔橫截面圖像中的內(nèi)腔邊界和中-外膜邊界，是對臨床重要的血管形態(tài)參數(shù)進行定量測量和血管三維重建等的重要步驟，對血管疾病的臨床診治具有重要意義。本文提出了一種基于snake模型的IVUS圖像序列的三維分割方法，不僅可充分利用整個圖像序列的信息，而且可實現(xiàn)對圖像序列的自動、并行分割。該方法在對原始圖像進行去噪和抑制環(huán)暈偽像等預(yù)處理的基礎(chǔ)上，將從四個縱向視圖中提取出的血管壁邊界映射到各幀IVUS橫向視圖中，獲得各幀圖像中的初始輪廓；然后通過使預(yù)先設(shè)定的snake能量函數(shù)最小，初始輪廓不斷演化得到最終輪廓。實驗結(jié)果表明，在獲得與現(xiàn)有逐幀處理的二維串行分割方法相類似的分割精度的情況下，采用本文的三維分割方法可

[1]楊青苗, 呂樹錚. 血管內(nèi)超聲在冠心病診斷中的臨床應(yīng)用進展[J]. 心肺血管病雜志, 2010, (01):73-76.

[2]Baldewsing RA, Danilouchkine MG, 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.

[3]孫 正. 應(yīng)用血管內(nèi)超聲與X射線造影圖像融合的血管三維重建[J]. 工程圖學(xué)學(xué)報, 2010, 31(1):116-123.

[4]Olszewski ME, Wahle A, Vigmostad SC, et al.Multidimensional segmentation of coronary intravascular ultrasound images using knowledgebased methods [C]//Proceedings of SPIE Conference on Medical Imaging and Image Processing, 2005,5747: 496-504.

[5]Kass M, Witkin A, Terzopoulos D. Snakes: active contour models [J]. International Journal of Computer Vision, 1987, 1(4): 321-331.

[6]Klingensmith J D, Nair A, Kuban B D, et al.Segmentation of three-dimensional intravascular ultrasound images using spectral analysis and a dual active surface model [C]//Proceedings of 2004 IEEE Ultrasonics Symposium, 2004, 3: 1765-1768.

[7]Brusseau E, de Korte C L, Mastik F, et al. Fully automatic luminal contour segmentation in intracoronary ultrasound imaging-a statistical approach [J]. IEEE Transactions on Medical Imaging,2004, 23: 554-566.

[8]Filho ES, Yoshizawa M, Tanaka A, et al. Detection of luminal contour using fuzzy clustering and mathematical morphology in intravascular ultrasound images [C]//Proceedings of 27th IEEE-EMBS Annual International Conference on Engineering in Medicine and Biology, 2005: 3471-3474.

[9]Giannoglou G D, Chatzizisis YS, Koutkias V, et al. A novel active contour model for fully automated segmentation of intravascular ultrasound images: in vivo validation in human coronary arteries [J].Computers in Biology and Medicine, 2007, 37(9):1292-1302.

[10]Mendizabal-Ruiz G, Rivera M, Kakadiaris IA. A probabilistic segmentation method for the identification of luminal borders in intravascular ultrasound images [C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 2008: 1-8.

[11]Papadogiorgaki M, Mezaris V, Chatzizisis Y S, et al.Image analysis techniques for automated IVUS contour detection [J]. Ultrasound in Medicine and Biology, 2008, 34(9): 1482-1498.

[12]Bucher S, Carlier S, Slabaugh G, et al. Shape-driven segmentation of the arterial wall in intravascular ultrasound images [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(5):335-347.

[13]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.

[14]Zhang Q, Wang Y, Wang W, et al. Automatic segmentation of calcifications in intravascular ultrasound images using snakes and the contourlet transform [J]. Ultrasound in Medicine and Biology,2010, 36(1): 111-129.

[15]Cardinal M H R, Meunier J, Soulez G, et al.Intravascular ultrasound image segmentation: a three-dimensional fast-marching method based on gray level distributions [J]. IEEE Transactions on Medical Imaging, 2006, 25(5): 590-602.

[16][美]Paul Schoenhagen, et al. 輕松掌握血管內(nèi)超聲[M].劉茜蒨, 劉 健, 陳 蕓譯. 北京: 人民軍醫(yī)出版社, 2009: 20-30.

[17]范幸義. 計算機圖形學(xué)[M]. 重慶: 重慶大學(xué)出版社, 2008: 110-120.

[18]孫 正. 變形模型技術(shù)在冠狀動脈造影圖像序列后處理中的應(yīng)用[J]. 工程圖學(xué)學(xué)報, 2007, 28(4):90-95.

[19]Williams D J, Shah M. A fast algorithm for active contours and curvature estimation [J]. CVGIP: Image Understanding, 1992, 55(1): 14-26.

[20]Papadogiorgaki M, Mezaris V, Chatzizisis Y S, et al.Image analysis techniques for automated IVUS contour detection [J]. Ultrasound in Medicine and Biology Journal, 2008, 34(9): 1482-1498.

[21]Gil D, Hernandez A, Rodriguez O, et al. Statistical strategy for anisotropic adventitia modelling in IVUS [J].IEEE Transactions on Medical Imaging, 2006, 25(6):768-778.

[22]張 麒, 汪源源, 王威琪, 等. 活動輪廓模型和Contourlet多分辨率分析分割血管內(nèi)超聲圖像[J].光學(xué)精密工程, 2008, 16 (11): 2303-2311.