基于雙邊濾波-距離正則化水平集演化算法的甲狀腺超聲圖像分割*

冉冬梅，嚴(yán)加勇，崔崤峣，于振坤

（1.上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093；2.上海健康醫(yī)學(xué)院醫(yī)療器械學(xué)院，上海201318；3.中科院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，蘇州215163；4.南京同仁醫(yī)院，南京211102）

1 引 言

近年來(lái)，甲狀腺癌發(fā)病率在世界范圍內(nèi)逐年快速上升，早期發(fā)現(xiàn)、早期診斷是其治療關(guān)鍵。目前，臨床上廣泛使用B超來(lái)實(shí)現(xiàn)對(duì)甲狀腺相關(guān)疾病的篩查[1-2]。但是，傳統(tǒng)的二維超聲無(wú)法直接描述甲狀腺病變的精準(zhǔn)位置。因此，準(zhǔn)確分割甲狀腺具有重要臨床意義。

受超聲成像原理影響，甲狀腺圖像易產(chǎn)生斑點(diǎn)噪聲，灰度對(duì)比度低、邊緣模糊。傳統(tǒng)的圖像分割算法如邊緣檢測(cè)、閾值分割等[3-4]，都難以達(dá)到甲狀腺分割要求。臨床醫(yī)師手動(dòng)分割結(jié)果雖較為準(zhǔn)確，但費(fèi)時(shí)費(fèi)力，分割結(jié)果往往還帶有主觀因素。

為實(shí)現(xiàn)甲狀腺有效分割，國(guó)內(nèi)外學(xué)者做了很多關(guān)于分割甲狀腺超聲圖像的研究。Chang等[5]在2010年提出使用徑向基函數(shù)（radial basis function，RBF）神經(jīng)網(wǎng)絡(luò)自動(dòng)分割甲狀腺的方案；2013年，Kaur等[6]分析了三種用于甲狀腺分割的方法，即無(wú)邊緣的活動(dòng)輪廓、基于局部區(qū)域的活動(dòng)輪廓和距離正則化水平集；同年，翟建敏等[7]提出采用手動(dòng)繪圖進(jìn)行甲狀腺微小結(jié)節(jié)的定位；2016年，Poudel等[8]提出先使用活動(dòng)輪廓分割超聲圖像中的甲狀腺，然后借助三維重建工具，獲取甲狀腺的三維模型。為提高分割精度，2017年P(guān)oudel等[9]通過添加均方誤差比和直方圖之間的相關(guān)性這兩個(gè)區(qū)域相似性指標(biāo)，將分割精度提高到86.7%。這些算法雖然能分割出甲狀腺，但普遍計(jì)算復(fù)雜、耗時(shí)且分割效率低。

目前，圖像分割領(lǐng)域一大研究熱點(diǎn)是幾何活動(dòng)輪廓模型。根據(jù)處理圖像的不同特征，分為基于圖像區(qū)域[10-12]和邊緣信息的幾何活動(dòng)輪廓模型[13-15]。其中，Li等[16]提出的基于邊緣的無(wú)需初始化變分水平集方法受到廣泛關(guān)注。該方法無(wú)需重新初始化，極大提高了分割效率，被稱為距離正則化水平集演化，即DRLSE模型。

本研究針對(duì)甲狀腺超聲圖像的特點(diǎn)，通過雙邊濾波來(lái)增強(qiáng)圖像對(duì)比度、抑制噪聲并保留圖像邊緣信息，降低原始圖像復(fù)雜度。采用基于目標(biāo)邊界信息的變分水平集方法，即改進(jìn)邊緣指示函數(shù)的DRLSE模型，實(shí)現(xiàn)對(duì)甲狀腺超聲圖像的分割。通過與使用另外兩種邊緣指示函數(shù)[16-17]的DRLSE模型的對(duì)比，本研究算法明顯減少了迭代次數(shù)和運(yùn)行時(shí)間。

2 本研究算法

2.1 BF模型

BF是一種非線性濾波器，最早由 Tomasi[18]提出，該濾波器可以在濾除噪聲的同時(shí)，保留圖像的邊緣的信息，降低了圖像的復(fù)雜度。BF是采用基于高斯分布加權(quán)平均的方法，用周邊像素亮度值的加權(quán)平均代表某個(gè)像素的強(qiáng)度。相較于高斯濾波，BF模型增加了對(duì)圖像邊緣的保護(hù)。BF組合了空間域和值域核函數(shù)，使得輸出圖像與空間鄰近度和像素值相似度均有關(guān)，即輸出像素不僅與像素空間距離相關(guān)，還與像素點(diǎn)的像素差值相關(guān)。

BF中，輸出像素的值f BF(x，)y依賴于鄰域像素值的加權(quán)組合，定義為：

權(quán)重系數(shù)w(x，y，i，j)（即核函數(shù)）取決于空間域和值域，其中，空間域權(quán)重因子為：

值域權(quán)重因子為：

BF的權(quán)重w(x，y，i，j)等于空間域權(quán)重因子和值域權(quán)重因子的乘積：

式（2）中，σd是空間域方差。σd越小，圖像邊緣和細(xì)節(jié)越清楚，σd越大，圖像越模糊。d(x，y，i，j)等同于高斯濾波系數(shù)。

式（3）中，σr是值域方差。在圖像平滑區(qū)域，σr越大，平滑噪聲能力越強(qiáng)；而在邊緣，σr越小，能保留的邊界特征越多。r(x，y，i，j)與空間像素差值相關(guān)。

實(shí)際應(yīng)用中，在圖像非邊緣區(qū)域，像素差值較小，σr變大，此時(shí)，空間域方差σd起主要作用，等同于普通的高斯濾波，保邊性能下降；在圖像邊緣，像素差值較大，σr變小，w(x，y，i，j)減小，當(dāng)前像素受到影響就越小，從而保持了邊緣信息。

圖1是BF模型核函數(shù)產(chǎn)生過程，該圖是在大小為9×9（由0和1組成）的模板上構(gòu)造的二值圖像，圖1（a）、（b）、（c）分別是 BF空間域函數(shù)（設(shè) σd=2）、值域函數(shù)（σr=0.2）（在邊緣（5，5）處）及 BF核函數(shù)圖像（分別用 D、R、W 函數(shù)表示式（2）、（3）、（4）對(duì)應(yīng)函數(shù)）。其中，D函數(shù)表示區(qū)域的位置關(guān)系，R函數(shù)體現(xiàn)像素的灰度關(guān)系。當(dāng)輸入有噪聲和邊緣的圖像，通過核函數(shù)W 運(yùn)算后，就可以得到降低噪聲并增強(qiáng)邊緣信息的圖像。

圖1 BF模型核函數(shù)圖像Fig 1 BF model kernel function image

圖2是一幅甲狀腺超聲圖像濾波后的圖。其中，圖2（a）為原圖，圖2（b）、圖2（c）分別對(duì)圖2（a）進(jìn)行雙邊濾波、高斯濾波操作，為消除使用B超測(cè)量甲狀腺?gòu)骄€，圖2（d）在圖2（b）基礎(chǔ)上增加中值濾波操作。

圖2 甲狀腺超聲圖像雙邊濾波結(jié)果Fig 2 Bilateral filtering results of thyroid ultrasound images

2.2 DRLSE模型

為解決水平集演化重新初始化問題，Li[19]提出一種增加懲罰項(xiàng)的方法，極大提高了曲線演化速率。但引入的懲罰項(xiàng)會(huì)造成擴(kuò)散率趨于無(wú)窮大，導(dǎo)致曲線演化無(wú)法到達(dá)期望的邊界。因此，Li等[16]又在能量函數(shù)中添加正則化項(xiàng)，使擴(kuò)散率保持為一個(gè)有界常數(shù)，從而實(shí)現(xiàn)演化曲線無(wú)需重新初始化即能到達(dá)目標(biāo)區(qū)域邊界，該方法即DRLSE模型。此模型是基于能量泛函的活動(dòng)輪廓模型。令Ω為圖像區(qū)域，I（x，y）為Ω→R上的灰度圖像，能量函數(shù)定義為：

其中，φ（x，y）是定義在域Ω→R上的水平集函數(shù)，μ＞0，是正則化項(xiàng)權(quán)重參數(shù)，Rp(φ)是正則化項(xiàng)，定義為：

Ρ為勢(shì)函數(shù)，即：

引入Ρ是為了使φ（x，y）在零水平集附近保持符號(hào)距離特性，即同時(shí)，在遠(yuǎn)離零水平集的位置，保持εext(φ)為外部能量，定義為：

其中，λ＞0，α∈R。Lg(φ)是以g為權(quán)重的加權(quán)長(zhǎng)度項(xiàng)，Ag(φ ) 是以g為權(quán)重的加權(quán)面積項(xiàng)。Lg(φ)和Ag(φ )分別定義為：

其中，g為邊緣指示函數(shù)，其作用是使零水平輪廓能停止在目標(biāo)區(qū)域邊緣，定義為：

其中，Gσ為標(biāo)準(zhǔn)偏差為σ的Gaussian內(nèi)核函數(shù)。

式（9）、（10）中，δ和H分別是 Dirac函數(shù)和Heaviside函數(shù)。一般可使用近似的δε和Hε，即：

綜上，能量泛函為：

其中，λ＞0，α∈R。

能量泛函的最小化可由以下梯度流來(lái)求解：

式中第一項(xiàng)對(duì)應(yīng)正則化項(xiàng)。第二項(xiàng)對(duì)應(yīng)加權(quán)長(zhǎng)度項(xiàng)，用于調(diào)整φ的零水平輪廓，驅(qū)使曲線按照平均曲率的方式朝目標(biāo)邊界演化。當(dāng)曲線某一點(diǎn)曲率為正，向內(nèi)收縮，為負(fù)則向外擴(kuò)張。第三項(xiàng)αgδε(φ)對(duì)應(yīng)加權(quán)面積項(xiàng)，用于加快演化曲線運(yùn)動(dòng)，當(dāng)α為正時(shí)，輪廓向內(nèi)收縮，α為負(fù)時(shí)，輪廓向外擴(kuò)張，從而驅(qū)動(dòng)零水平集曲線向邊界演化。

式（15）中后兩項(xiàng)都涉及控制曲線演化位置的邊緣指示函數(shù)g。由式（11）知，g是嚴(yán)格非負(fù)遞減函數(shù)，它依賴于圖像的灰度，通常選高斯平滑后的梯度信息來(lái)進(jìn)行運(yùn)算，因此g對(duì)噪聲極其敏感。一般情況下，當(dāng)演化曲線位于圖像中邊緣時(shí)，灰度變化較大，梯度較大，g近似為0，曲線停止演化；曲線遠(yuǎn)離邊緣時(shí)，梯度值相對(duì)較小，g不為0，曲線繼續(xù)演化。但是，超聲圖像中噪聲較大，可能導(dǎo)致曲線演化停留在當(dāng)前噪聲點(diǎn)，而不繼續(xù)演化，致使曲線無(wú)法演化到目標(biāo)邊界。而且，若圖像存在弱邊緣，曲線演化到邊界附近，g的值仍充分大，曲線還有較大演化速度，致使曲線越過邊界繼續(xù)演化，造成邊界泄露，定位不準(zhǔn)確。故本研究提出一個(gè)新的邊緣指示函數(shù)。

2.3 改進(jìn)邊緣指示函數(shù)

鑒于甲狀腺超聲圖像的噪聲大且邊界模糊，提出綜合考慮噪聲和邊界強(qiáng)弱的邊緣指示函數(shù)，定義為：

其中，ρ為控制曲線收斂速率的參數(shù)，θ為控制噪聲敏感度的參數(shù)。ρ＞0，θ＞0，且均為常量。經(jīng)實(shí)驗(yàn)分析，分割甲狀腺超聲圖像時(shí)，ρ、θ最優(yōu)取值分別為：ρ取1～100，步長(zhǎng)是1；θ取0.1～2，步長(zhǎng)0.1。實(shí)際應(yīng)用中，ρ和θ相互作用，具體取值應(yīng)根據(jù)各個(gè)圖像特征實(shí)時(shí)調(diào)整。

若圖像目標(biāo)邊界較強(qiáng)，則取較小ρ值；反之，取較大ρ值，以加快gI收斂到0的速率，使其在弱邊界也能快速收斂。

若圖像噪聲較小，可取較大θ值來(lái)加速演化；若圖像噪聲較大，則取較小θ值，使演化曲線可以跳過噪聲點(diǎn)繼續(xù)演化，直到停止在目標(biāo)邊界。

3 實(shí)驗(yàn)結(jié)果及分析

本研究分別對(duì)六幅受噪聲污染程度和邊界清晰度均不同的甲狀腺超聲圖像進(jìn)行分割實(shí)驗(yàn)，以驗(yàn)證算法的有效性。實(shí)驗(yàn)環(huán)境為Intel（R）Core（TM）i7-6700 CPU 3.40 GHz，RAM：8 GB。操作系統(tǒng)：Windows 7，軟件環(huán)境：MATLAB 2014a。實(shí)驗(yàn)圖像均來(lái)源于南京同仁醫(yī)院。

六幅實(shí)驗(yàn)圖均使用三種不同邊緣指示函數(shù)，并對(duì)比迭代次數(shù)和運(yùn)行時(shí)間。除改進(jìn)的gI外，另兩種邊緣指示函數(shù)分別為L(zhǎng)i在DRLSE中使用的和劉哲等[17]改進(jìn)的也依賴于圖像灰度信息，對(duì)噪聲非常敏感。

圖3～圖8是本實(shí)驗(yàn)選用的六幅甲狀腺超聲圖像分割結(jié)果。

圖3 噪聲較大邊界不連續(xù)的甲狀腺超聲圖像的分割Fig 3 Segmentation of thyroid ultrasound images with large noise and discontinuous boundaries

圖4 噪聲大邊界較模糊的甲狀腺超聲圖像的分割Fig 4 Segmentation of thyroid ultrasound images with large noise and blurred boundaries

圖5 噪聲較小邊界較弱的甲狀腺超聲圖像的分割Fig 5 Segmentation of thyroid ultrasound image with small noise and weak boundary

圖6 噪聲較小邊界較模糊的甲狀腺超聲圖像的分割Fig 6 Segmentation of thyroid ultrasound images with small noise and blurred boundaries

圖7 噪聲小邊界較清晰的甲狀腺超聲圖像的分割Fig 7 Segmentation of thyroid ultrasound images with small noise and clear boundaries

圖8 噪聲小邊界清晰的甲狀腺超聲圖像的分割Fig 8 Segmentation of thyroid ultrasound images with small noise and clear boundary

對(duì)比圖3～圖8中的（a）、（b）圖可知，對(duì)（a）圖進(jìn)行雙邊濾波得到（b）圖，此過程有效降低了斑點(diǎn)噪聲，保護(hù)了甲狀腺區(qū)域的邊界信息。

本研究改進(jìn)gI的參數(shù)設(shè)置見表1（分別對(duì)應(yīng)圖3～圖8中（e）圖）。其中，圖3（a）、圖4（a）噪聲污染嚴(yán)重、甲狀腺邊界較強(qiáng)，故實(shí)驗(yàn)采用較小θ值（0.8和0.7）和較小的 ρ值（20和23）。圖5（a）、圖6（a）甲狀腺邊界較弱、噪聲斑點(diǎn)也較弱，故選取了較大ρ值（50和80）和 θ值（1.2和 1.4）。圖 7（a）、圖 8（a）是本實(shí)驗(yàn)中甲狀腺邊界最強(qiáng)的兩幅圖，故選用本研究最小ρ值（15和18）。

表1 改進(jìn)邊緣指示函數(shù)gI參數(shù)設(shè)置Table 1 Parameter setting of improved edge indicator function gI

在保證其他參數(shù)設(shè)置均不變的情況下，實(shí)驗(yàn)比較了使用g、g I、gL的DRLSE模型在分割甲狀腺圖像時(shí)所需迭代次數(shù)和運(yùn)行時(shí)間（對(duì)應(yīng)圖3～圖8中（c）、（d）、（e）），見表2。

表2 分別使用g、gL和gI函數(shù)時(shí)的迭代次數(shù)和運(yùn)行時(shí)間Table 2 The number of iterations and running time when using g、gL and gI functions respectively

分析表2，使用本研究改進(jìn)gI的DRLSE模型分割甲狀腺，明顯減少了曲線演化迭代次數(shù)和運(yùn)行時(shí)間，極大提高了曲線演化速率和分割效率。

本研究分割結(jié)果的準(zhǔn)確性，采用Shattuck等[20]提出的骰子相似系數(shù)（dice similarity coefficient，DSC）來(lái)表示，其表達(dá)式為（SSEG是算法分割結(jié)果；SGT是手動(dòng)分割結(jié)果；N(S)表示分割區(qū)域的面積）PDSC的數(shù)值越接近1，則表示分割精度越高，分割結(jié)果越準(zhǔn)確。本研究中圖3～圖 8（c）、（d）、（e）圖分割精度見表 3。

表3 分別使用g、gL和gI函數(shù)時(shí)的分割精度Table 3 Segmentation accuracy using g，gL and gI functions respectively

由表3知，DRLSE模型分別采用g、gL和gI，均能達(dá)到較高分割精度，表明均能大致分割出超聲圖像中的甲狀腺區(qū)域。

綜上，在保證圖像分割精度的同時(shí)，使用本研究方法分割超聲圖像中的甲狀腺區(qū)域，可以明顯減少迭代次數(shù)和運(yùn)行時(shí)間，有效地提高了曲線演化速率和分割效率。

4 結(jié)論

本研究針對(duì)甲狀腺超聲圖像受噪聲污染嚴(yán)重等特點(diǎn)，首先使用雙邊濾波對(duì)其進(jìn)行預(yù)處理，降低了甲狀腺及其周圍組織區(qū)域的復(fù)雜度。然后對(duì)Li提出的DRLSE模型中的邊緣指示函數(shù)進(jìn)行了改進(jìn)，減少了水平集演化過程中的迭代次數(shù)和運(yùn)行時(shí)間，在保證分割精度的同時(shí)，有效提高了分割效率。

但是，改進(jìn)邊緣指示函數(shù)中參數(shù)ρ和θ的值需要人為設(shè)置，缺乏自適應(yīng)性，一定程度上影響了甲狀腺超聲圖像分割效率。因此，如何更好的、更合理的設(shè)置好參數(shù)，提高分割效率，是需要繼續(xù)研究的課題。