基于LevelSet的超聲乳腺腫瘤圖像的輪廓提取

黃韞梔 HUANG Yunzhi

劉 奇 LIU Qi

四川大學電氣信息學院醫(yī)學信息工程系四川成都 610065

乳腺癌發(fā)病率居女性癌癥首位，近年呈不斷上升趨勢，因此開展乳腺癌的診斷和防治研究具有重要意義。超聲診斷以其無損傷、非侵入性、靈敏度高、適用于鑒別軟組織等特點，越來越多地應用于乳腺腫瘤的輔助檢測，而對乳腺腫瘤超聲圖像進行分割以獲取腫瘤的邊界是判斷乳腺腫瘤良惡性的關(guān)鍵。由于超聲圖像存在大量的斑點噪聲，圖像質(zhì)量明顯下降，造成醫(yī)學圖像的復雜性和多樣性，從而大大增加了乳腺腫瘤分割的難度。本研究旨在采用曲率各向異性擴散的方法，在保留主要邊界信息的前提下對超聲乳腺圖像加以平滑，并利用水平集的方法自動提取乳腺腫瘤輪廓，為后續(xù)研究打下基礎(chǔ)。

1 資料與方法

1.1 腫瘤粗定位 鎖定目標，提取腫瘤所在的方形區(qū)域（圖1），以減小后期進行自動輪廓提取的迭代時間，并增加迭代的成功率。

1.2 保留邊緣平滑濾波 傳統(tǒng)的圖像平滑算法如均值濾波、中值濾波和高斯濾波等，屬于各向同性擴散，這些平滑算子在去除噪聲的同時會模糊甚至破壞圖像的邊緣[1]。為了在去除噪聲的同時保持圖像的邊緣，Perona等[2]提出在不同的方向設(shè)置不同分布系數(shù)的各向異性擴散濾波（anisotropic diffusion），建立了各向異性擴散方程：

圖1 腫瘤粗定位

其中原始圖像f (x, y)作為其初始值，在不同時刻可以得到在此時刻下擴散濾波后的圖像。若設(shè)置：

為擴散程度的主控因子，其中X＝（x,y）代表圖像中像素點的位置，g（X）表示每個像素點的梯度強度[3]。顯然，C（X）是梯度幅值的減函數(shù)。采用該方法對乳腺超聲圖像的濾波效果見圖2。

圖2 A.包含腫瘤的區(qū)域；B.基于梯度各向異性擴散濾波（步長=0.25，迭代次數(shù)=5，傳導系數(shù)=5）

本次實驗中，以公式（3）作為主控因子表達式的曲率各向異性擴散（MCDE）[4]對圖像進行濾波。MCDE方程見公式（4）。

對比圖2B和圖3B可知，與經(jīng)典的Perona-Malik擴散方式相比，MCDE具有更低的敏感度，在去除噪聲和保留超聲圖像的邊界信息方面優(yōu)于梯度各向異性擴散，并能夠更好地保留圖像中的結(jié)構(gòu)細節(jié)。此外，對比圖2B和圖3C可知，MCDE使用更少的迭代即可以得到一個可接受的結(jié)果。因此，采用MCDE對超聲醫(yī)學圖像進行預處理，不僅能夠減少圖像中的斑點和噪聲，同時能夠保留乳腺腫瘤的對比度和形狀。1.3 邊緣的潛在圖像 圖像的梯度強度廣泛應用于圖像分析，幫助檢測物體的輪廓以及分離均勻區(qū)域。為了能夠得到更加準確的輪廓，先將經(jīng)過MCDE濾波后的圖像再次進行濾波——高斯濾波（選擇合適的標準差，控制圖像邊緣的影響范圍），然后計算圖像的梯度。濾波后的梯度圖像見圖4，其中σ＝1.0。

圖3 A.包含腫瘤的區(qū)域；B、C.基于曲率各向異性擴散濾波（B：步長=0.25，迭代次數(shù)=5，傳導系數(shù)=5；C：步長=0.25，迭代次數(shù)=3，傳導系數(shù)=5）

圖4 梯度圖像

在曲面演化的過程中，需要設(shè)置演化終止時的條件，本實驗中終止的條件就是曲面到達了圖像邊緣。為了產(chǎn)生邊緣的潛在圖像，本次實驗中選用Sigmoid函數(shù)[5,6]：

其中g(shù)是輸入像素的梯度，g'是輸出像素的梯度，gmax、gmin是經(jīng)過MCDE濾波以后所得圖像的最大梯度值和最小梯度值，α定義了輸入梯度值范圍的寬度，β定義了圍繞在范圍中心的梯度。所得潛在邊界見圖5。

圖5 潛在的腫瘤邊界

Osher等[7,8]于1982年提出了基于幾何變形模型的水平集方法。利用水平集實現(xiàn)圖像分割的典型用法為：由用戶先定義出待分割對象的初始輪廓，再利用水平集的進化功能，使水平集一直進化到符合待分割對象的解剖結(jié)構(gòu)為止。

在尋找腫瘤邊界的過程中，不能保證速度函數(shù)F恒正[9-12]，因而不能利用FastMarching[13]提取乳腺腫瘤輪廓，需采用動態(tài)水平集：

其中初始條件為距離函數(shù)：

2.1 簡單的初始輪廓 按照以下方式生成一個初始水平集曲面：①待分割區(qū)域內(nèi)用戶提供一個用于輪廓擴張的種子點（a，b），一個好的種子點集將增大分割一個復雜對象而不丟失數(shù)據(jù)的可能性。在腫瘤輪廓提取實驗中，由于腫瘤的特殊形狀，一般設(shè)置軸線上的中心位置為種子點。②以種子點為圓心，建立一個半徑為r的圓，該圓作為初始的輪廓曲線，初始水平集曲面的方程表達式為：

水平集的初始輪廓見圖6。

圖6 水平集的初始輪廓

2.2 乳腺輪廓的提取 曲面隨著時間而變化，嵌入其中的零水平集曲線的拓撲結(jié)構(gòu)也會隨之發(fā)生改變。在曲面演化過程中，可利用圖像信息，如灰度均值、梯度和邊緣等來控制曲面演化的過程。通用的曲面運動偏微分方程（PDE）為：

其中A（X）是水平對流系數(shù)，將水平線集引向?qū)ο蟮倪吔?；P（X）是傳播（膨脹）系數(shù)，用于初始輪廓的擴散；其中Z（X）是一個曲率均值的空間調(diào)節(jié)器系數(shù)：高曲率區(qū)域的平滑系數(shù)，進一步將噪聲平滑，曲率越大，分割結(jié)果越平滑，然而曲率縮放比例太大會把輪廓移出形狀邊界[15]。常數(shù)α、β、γ是每個系數(shù)在界面運作上相關(guān)影響的權(quán)值。因而曲面演化的方程表達式即為：

2.2.1 測地活動輪廓（GAC） Caselles等[16]闡述了公式（10）的具體實現(xiàn)，即GAC的求解。演化過程需要定義達到收斂的標準或達到迭代器的最大數(shù)量，和邊界的潛在圖像一起構(gòu)成水平集的停止條件。本實驗的收斂標準是以水平集函數(shù)的均方根變換的形式來定義的。如果均方根變換在用戶指定的一個門限之下，運動將會聚合到一點。具體參數(shù)設(shè)置：收斂標準：RMS=0.02，迭代器的最大數(shù)量為1200。

圖7顯示了利用GAC采用表1中的參數(shù)所得的圖像輪廓。實驗過程中測地線方法中的參數(shù)和傳播參數(shù)很難控制，并且從測試的結(jié)果也發(fā)現(xiàn)最后得到的速度圖像收斂到的邊緣處并不清晰，與原圖結(jié)構(gòu)有很大差異。圖7B中孔洞即是分割的泄露。

圖7 測地活動輪廓水平集分割腫瘤

表1 GAC水平集分割的參數(shù)設(shè)置（曲率縮放=1.0，水平縮放=1.0）

對于超聲乳腺腫瘤圖像的分割，需要一個相對較大的傳播縮放比例：①超聲圖像的邊界對比較低；②乳腺腫瘤生理結(jié)構(gòu)的邊界形狀比較復雜。但是這些縮放比例參數(shù)的最優(yōu)值只能通過不斷地實驗得到。

2.2.2 閾值水平集分割 由于速度函數(shù)的定義直接關(guān)系到水平集方法的分割質(zhì)量，針對GAC方法在連接成員方案中出現(xiàn)的“泄漏”，因此對GAC水平集的分割算法中的速度函數(shù)進行修改，并且定義一個灰度值范圍來對相關(guān)的組織類型繼續(xù)分類，然后求出對那個亮度范圍基于水平集等式上的傳播系數(shù)。任意點的速度函數(shù)F：

缺省了公式（10）中的水平對流系數(shù)。速度函數(shù)F中，需滿足圖8所示的要求傳播項：

其中g(shù)（X）為輸入的待分割圖像中各點的灰度值，U為閾值上限，L為閾值下限[17]。

對于公式（12）中閾值L和U的選取，可以參照以下2點：①乳腺腫瘤超聲中，目標處于低灰度值范圍，因此L值一般采用ROI區(qū)域中的最小灰度值。下門限設(shè)置成負數(shù)，以便確保被分割對象的內(nèi)部出現(xiàn)在二值區(qū)域內(nèi)。②U值采用大津閾值法（Ostu）[18]，即類間方差最大化或類中差異最小化進行估計。

圖8 傳播系數(shù)函數(shù)

實驗結(jié)果見圖9，其中各參數(shù)設(shè)置見表2。下門限設(shè)置成一個大的負數(shù)以便確保被分割對象的內(nèi)部將出現(xiàn)在二值區(qū)域內(nèi)。閾值水平集的方法改進了GAC和FastMarching分割方法在弱邊界處出現(xiàn)“泄漏”，并且提取到與圖10手工勾邊較為相似的輪廓，較好地包裹了腫瘤區(qū)域。

圖9 閾值水平集分割結(jié)果。圖10 醫(yī)師手動勾畫的輪廓

表2 閾值水平集分割的參數(shù)設(shè)置（曲率縮放=1.0）

2.2.3 閾值水平集分割定量分析 為了評價分割結(jié)果的精度，以醫(yī)師手工勾畫的結(jié)果作為“金標準”，圖11顯示了自動提取病灶區(qū)域與“金標準”之間的關(guān)系。本文采用Udupa[19]提出的3種不同的測量誤差對99幅乳腺腫瘤（61幅良性腫瘤，38幅惡性腫瘤）圖像的分割結(jié)果進行定量分析：假陽性（FP）、假陰性（FN）、真陽性（TP）：

圖12和表3顯示了99例腫瘤利用閾值水平集進行自動分割后得到的FP、FN、TP的測量結(jié)果。由于獲取的乳腺超聲圖像會出現(xiàn)邊界極不清晰的情況，使得自動分割未獲得滿意的結(jié)果（<5%，如圖12A中第42例和圖12B中第26例），此時不能完全依靠自動分割，還需要醫(yī)師進一步調(diào)整腫瘤輪廓以獲得比較滿意的結(jié)果。

圖11 自動分割區(qū)域和手動分割區(qū)域真陽性、假陽性、假陰性比較

圖12 閾值水平集獲取的超聲乳腺腫瘤輪廓（圖像均采自四川大學華西醫(yī)院超聲科）

表3 99例腫瘤分割的誤差測度（99例圖像均采自四川大學華西醫(yī)院超聲科）

3 結(jié)束語

本文首先定位到超聲圖像中的腫瘤區(qū)域，然后針對乳腺腫瘤超聲圖像自身的一些缺點如圖像對比度低、斑點噪聲大、部分邊緣缺失、腫瘤內(nèi)部微細結(jié)構(gòu)分布復雜（如血管、鈣化灶等），特別是惡性腫瘤還具有復雜形狀等，提出了針對性的MCDE的方法進行濾波處理，平滑了噪聲，但邊緣信息丟失甚少，對于后期圖像處理至關(guān)重要。乳腺腫瘤邊界的提取采用LevelSet方法，其中最關(guān)鍵的是速度函數(shù)的確定，不同速度函數(shù)導致不同的分割效果，不同圖像的目標分割需要的速度函數(shù)也不盡相同。本實驗首先設(shè)定腫瘤的初始輪廓，然后測試了GAC以及閾值水平集方法在腫瘤自動分割中的效果。然而GAC方法的實驗參數(shù)較難控制，容易出現(xiàn)分割泄露，分割效果不佳；相對而言，閾值水平集分割只需要通過Ostu方法找到閾值，通過小范圍的偏移就可以獲得與手動勾邊相似的輪廓，大幅度改善了GAC的分割效果。自動提取腫瘤邊界為后續(xù)腫瘤良惡性的判斷打下了重要基礎(chǔ)，也正朝著解放醫(yī)師讀片工作、自動診斷腫瘤的方向上前進。

[1] Gonzalez RC, Woods RE. Digital image processing. 2nd ed.Beijing: Publishing House of Electronics Industry, 2008: 183-215.

[2] Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis Machine Intelligence, 1990, 12(7): 629-639.

[3] Black MJ. Robust anisotropic diffusion. IEEE Trans Image Process, 1998, 7(3): 421-432.

[4] Huang YL, Jiang YR, Chen DR, et al. Level set contouring for breast tumor in sonography. J Digit Imaging, 2007, 20(3): 238-247.

[5] 周振環(huán), 伍云智, 趙明. 醫(yī)學圖像編程技術(shù). 北京:電子工業(yè)出版社, 2010: 9-11.

[6] 周振環(huán), 王安明,王京陽. 醫(yī)學圖像分割與配準Ⅰ. 成都:電子科技大學出版社, 2007: 111-113.

[7] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulation. J Computer Physics, 1988, 79(1): 12-49.

[8] Sethian JA. Level set methods and fast marching methods. 2nd ed. Cambridge: Cambridge University Press, 1999.

[9] 徐倩. 基于Level Set模型的腦血管圖像骨架提取研究. 蘇州:蘇州大學碩士學位論文, 2009.

[10] 陳波,賴劍煌.用于圖像分割的活動輪廓模型綜述.中國圖象圖形學報, 2007, 12(1): 11-20.

[11] 師為禮,苗語,顏雁. 基于ITK的水平集醫(yī)學圖像分割算法.長春工業(yè)大學學報(自然科學版), 2008, 29(5): 516-519.

[12] 王崢, 楊新, 李俊, 等. 基于水平集人機交互模型的醫(yī)學圖像分割. 模式識別與人工智能, 2002, 15(4): 392-396.

[13] Sethian JA. A fast marching level set method for monotonically advancing fronts. PNAS, 1996, 93(4): 1591-1595.

[14] 鄭國賢. LevelSet方程的求解與應用. 杭州: 浙江大學碩士學位論文, 2006.

[15] Sethian JA. Level set methods: evolving interfaces in geometry,fluid mechanics, computer vision and materials science.Cambridge: Cambridge University Press, 1996.

[16] Caselles V, Kimmel R, Sapiro G. Geodesic active contours.International Journal on Computer Vision, 1997, 22(1): 61-97.

[17] Ibanez L, Schroeder W, Ng L, et al. ITK software guide. 2nd ed. New York: Kitware Inc, 2005: 531-558.

[18] 夏良正. 數(shù)字圖像處理(修訂版). 南京:東南大學出版社,1999: 222-223.

[19] Udupa JK. A methodology for evaluating image segmentation algorithms. Proc SPIE, 2002, 4684(2): 266-277.