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

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

1 引 言

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

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

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

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

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

2 本研究算法

2.1 BF模型

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

實際應(yīng)用中，在圖像非邊緣區(qū)域，像素差值較小，σr變大，此時，空間域方差σ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）對應(yīng)函數(shù)）。其中，D函數(shù)表示區(qū)域的位置關(guān)系，R函數(shù)體現(xiàn)像素的灰度關(guān)系。當(dāng)輸入有噪聲和邊緣的圖像，通過核函數(shù)W 運算后，就可以得到降低噪聲并增強(qiáng)邊緣信息的圖像。

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

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

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

2.2 DRLSE模型

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

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

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

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

其中，λ＞0，α∈R。Lg(φ)是以g為權(quán)重的加權(quán)長度項，Ag(φ ) 是以g為權(quán)重的加權(quán)面積項。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。

能量泛函的最小化可由以下梯度流來求解：

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

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

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

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

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

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

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

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

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

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

圖3～圖8是本實驗選用的六幅甲狀腺超聲圖像分割結(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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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