黃韞梔 HUANG Yunzhi

劉 奇 LIU Qi

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

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

1 資料與方法

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

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

圖1 腫瘤粗定位

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

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

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

本次實(shí)驗(yàn)中，以公式（3）作為主控因子表達(dá)式的曲率各向異性擴(kuò)散（MCDE）[4]對(duì)圖像進(jìn)行濾波。MCDE方程見(jiàn)公式（4）。

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

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

圖4 梯度圖像

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

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

圖5 潛在的腫瘤邊界

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

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

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

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

水平集的初始輪廓見(jiàn)圖6。

圖6 水平集的初始輪廓

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

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

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

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

圖7 測(cè)地活動(dòng)輪廓水平集分割腫瘤

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

[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ī)學(xué)圖像編程技術(shù). 北京:電子工業(yè)出版社, 2010: 9-11.

[6] 周振環(huán), 王安明,王京陽(yáng). 醫(yī)學(xué)圖像分割與配準(zhǔn)Ⅰ. 成都:電子科技大學(xué)出版社, 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模型的腦血管圖像骨架提取研究. 蘇州:蘇州大學(xué)碩士學(xué)位論文, 2009.

[10] 陳波,賴(lài)劍煌.用于圖像分割的活動(dòng)輪廓模型綜述.中國(guó)圖象圖形學(xué)報(bào), 2007, 12(1): 11-20.

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

[12] 王崢, 楊新, 李俊, 等. 基于水平集人機(jī)交互模型的醫(yī)學(xué)圖像分割. 模式識(shí)別與人工智能, 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] 鄭國(guó)賢. LevelSet方程的求解與應(yīng)用. 杭州: 浙江大學(xué)碩士學(xué)位論文, 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ù)字圖像處理(修訂版). 南京:東南大學(xué)出版社,1999: 222-223.

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