基于Canny算子和距離正則化水平集的乳腺植入物圖像分割算法

李建華 王藝衡* 楊來俠 李素麗 張聚良

1(西安科技大學 陜西 西安 710054) 2(第四軍醫(yī)大學西京醫(yī)院 陜西 西安 710032)

0 引 言

乳房一期重建手術(shù)患者在植入假體后，醫(yī)生通常會對患者的恢復(fù)情況進行跟蹤，以評估患者的恢復(fù)情況以及獲取科研數(shù)據(jù)。通常是對患者的植入?yún)^(qū)域利用CT或者核磁共振圖像進行三維重建。而圖像分割是三維重建的基礎(chǔ)，因此，圖像分割的精確度會直接影響三維建模的精度[1]。由于傳統(tǒng)圖像分割算法對噪聲敏感，因此在分割時會出現(xiàn)很多虛假的邊緣。無法得到滿意的分割效果。國內(nèi)外學者為了改善傳統(tǒng)圖像分割算法的缺陷，進行了很多的研究。Li等[2]將動態(tài)輪廓模型結(jié)合Snake模型在HSV空間里使分割曲線收斂到了相應(yīng)區(qū)域，但由于只采用了圖像的梯度信息來作為Snake模型的外部能量項，因此這種方法很難收斂到最佳的邊緣。陸東瑩等[3]利用了Level Set算法在演化時能量平均變化的特性，提出了在低對比度環(huán)境下用于醫(yī)學分割的Level Set方法，使得算法的魯棒性和自動性得到增強。Li等[4]在Level Set算法中引入了距離正則化項，消除了由于保留水平集方程符號距離特性而引起的不良邊緣效應(yīng)。使得水平集函數(shù)不需要在演化的過程中反復(fù)地重新初始化水平集函數(shù)，大大提高了Level Set算法的實用性，但在實際使用中依然存在算法無法收斂到最佳邊緣的缺點。

針對Canny算子邊緣易出現(xiàn)斷裂以及距離正則化水平集算法(DRLSE)不易收斂到最佳邊緣等問題，本文提出基于Canny算子和DRLSE算法的乳腺植入物圖像分割算法。算法結(jié)合了Canny算法定位邊界精確和DRLSE空間連續(xù)演化的思想，可形成連續(xù)且準確的體內(nèi)植入物分割邊界。

1 基本算法簡介

1.1 距離正則化水平集的基本原理

DRLSE是在基于外部能量的水平集方程里面添加了距離正則化項實現(xiàn)對水平集函數(shù)的形狀進行控制[5-6]。

E(Φ)=μRp(Φ)+Eext(Φ)

(1)

式中：Rp(Φ)是正則化項，用來控制水平集方程的符號距離屬性[7]；μ>0是正則化項的參數(shù)；Eext(Φ)是外部能量項。計算正則化項首先需要求Rp(Φ)的Gateaux導(dǎo)數(shù)：

(2)

div是計算散度，dp的定義如下：

(3)

式中：p是一個潛在方程(或者也稱為能量密度方程)。

(4)

為了滿足水平集函數(shù)的條件，dp需要在s=1和s=0處有最小點。因為p的作用，當dp(|▽Φ|)為正的時候，能量方程是前向擴散，|▽Φ|增加；當dp(|▽Φ|)為負的時候，能量方程是后向擴散，|▽Φ|減少。這種擴散稱為前進和后退(FAB)擴散[7]。正是由于這種前后向擴散方式，使得|▽Φ|能夠得到適當?shù)恼{(diào)整，讓它逼近p(s)兩個最小值中的一個，使水平集函數(shù)可以保持期望的形狀，消除了不良邊緣的影響[8]。

將距離正則化項代入到水平集方程中，于是可以得到距離正則化水平集能量方程的最小化梯度求解方法：

(5)

δε是Heaviside函數(shù)，其定義如下：

(6)

1.2 Canny算子簡介

Canny邊緣檢測算法是由Canny在1986年提出的一種多級邊緣檢測算法[11]，具有比較好的信噪比和較高的邊緣檢測的精度。Canny在提出Canny算法時，還提出了邊緣檢測的三項準則：(1) 邊緣檢測時的錯誤率要盡可能的低：邊緣檢測算法要求能夠精確地、盡可能多地找出圖像中的邊緣，同時盡可能地減少誤檢和漏檢；(2) 最佳的定位能力：檢測出的邊緣點應(yīng)該精確地定位在邊緣的中心；(3) 足夠低的響應(yīng)：圖像的邊緣區(qū)域不可以被多次標記，同時不應(yīng)該由于噪音而產(chǎn)生虛假的邊緣[12]。Canny算法首次采用數(shù)學的方法表示了上述判據(jù)，同時利用最優(yōu)化數(shù)值方法獲得了對于特定邊緣最優(yōu)的檢測模板。

2 算法設(shè)計

本文提出的基于Canny算子和DRLSE算法的乳腺植入物圖像分割算法的具體流程如圖1所示。

Canny邊緣檢測算法針對的是圖像的一維邊緣，其用于檢測階躍邊緣的最佳模板形狀類似于高斯函數(shù)的一階微分。由于二維的高斯函數(shù)具有圓對稱和可分解等性質(zhì)，使得高斯函數(shù)在圖像的任意方向的方向?qū)?shù)和卷積可以很輕松地計算出來。于是將高斯函數(shù)的一階微分用作次最優(yōu)的檢查算子來對圖像進行邊緣檢測更符合實際需求[13]。

本文算法的處理過程如下：

Step1對乳腺植入物核磁共振圖像進行高斯濾波，高斯濾波函數(shù)在二維情況下的表達式如下：

(7)

設(shè)原圖像為f(x,y)，將它做高斯平滑后的圖像是：

gG(x,y)=f(x,y)*G(x,y,σ)

(8)

本文采用代碼“cv2.GaussianBlur(Img,ksize,Sigma)”對圖像進行高斯平滑。

Step2gG(x,y)是進行高斯平滑后的圖像，對其求一階導(dǎo)數(shù)，然后表示為梯度向量的形式：

(9)

式中：Gx(x,y,σ)和Gy(x,y,σ)是式(6)在x方向和y方向的一階偏導(dǎo)數(shù)。

本文采用代碼“Iy,Ix=numpy.gradient(Img)”獲得圖像在X和Y方向的一階導(dǎo)數(shù)。

(10)

本部分對求解圖像梯度大小與方向角的核心代碼如下：

M=numpy.sqrt(numpy.square(Ix)+numpy.square(Iy))

#求圖像梯度大小

theta[i,j]=math.atan(dx[i,j]/(dy[i,j]+0.000000001))

#求圖像方向角

Step4將使用Canny算子得到的邊緣代入水平集的邊緣指示函數(shù)中。

邊緣指示函數(shù)g解析式為：

(11)

式中：▽是梯度算子；*是卷積運算；Gσ是標準方差為σ的二維高斯濾波器，用于去除圖像中的噪聲。I(x,y)是原圖的灰度圖。

將gxy(x,y)引入邊緣指示函數(shù)。

(12)

式中：gxy(x,y)是圖像邊緣梯度圖，gxy(x,y)在圖像的邊緣區(qū)灰度會顯著地大于平滑區(qū)的灰度。當邊緣指示函數(shù)的坐標在gxy(x,y)的邊緣時，|Gσ*gxy(x,y)|的梯度會進一步增加，g(x,y)的值快速趨近于0；當邊緣指示函數(shù)的坐標在gxy(x,y)的平滑區(qū)域時，|Gσ*gxy(x,y)|的梯度會進一步減小，g(x,y)的值快速趨近于1。通過判斷g(x,y)的值是與0接近還是與1接近便可以控制零水平集函數(shù)逐漸貼近圖像的輪廓[14]。

Step5將邊緣指示函數(shù)g(x,y)代入式(5)中，待用戶選擇好初始零水平集、擴散方向和迭代次數(shù)后，程序按照迭代次數(shù)循環(huán)更新水平集函數(shù)。

Step6將水平集函數(shù)輸出，得到最終的分割結(jié)果。

由于Canny算子的結(jié)果比原圖擁有更明顯、更精確的邊緣信息[15]，并且在水平集函數(shù)接近圖像邊緣的時候，|Gσ*gxy(x,y)|的梯度遠高于|Gσ*I(x,y)|，改善了DRLSE算法出現(xiàn)邊緣泄漏的情況。因此該算法擁有更高的分割精度以及抗干擾性能。

3 實 驗

本文實驗使用的圖像是由第四軍醫(yī)大學西京醫(yī)院甲乳外科提供。使用核磁共振系統(tǒng)增強掃描多名2016年3月至2019年9月間乳腺切除一期重建患者得到的乳腺植入物圖片序列。采用Python平臺開發(fā)，在聯(lián)想Intel(R)Core(TM)i5- 6200U CPU，內(nèi)存8 GB的PC機上運行本文算法，實驗環(huán)境與圖像相關(guān)參數(shù)如表1所示。

表1 實驗環(huán)境與圖像相關(guān)參數(shù)

3.1 分割精度對比實驗

實驗中的參數(shù)選擇如下：初始水平集的大小會影響該算法迭代的次數(shù)，但并不影響其最終分割的效果。因此根據(jù)與植入物橫截面積的相對大小在植入物內(nèi)部任意區(qū)域采用手工勾畫的方式選擇了一個10×10像素的矩形作為初始零水平集。參考DRLSE算法在醫(yī)學圖像分割中的應(yīng)用[5]，當a設(shè)置過大時，邊緣泄漏比較嚴重，過小時，容易使演化提前停止。a值為正，水平集函數(shù)向內(nèi)擴散，為負則向外擴散。將外部能量項中a設(shè)置為-9；分別將迭代次數(shù)設(shè)置為50、100、150對算法進行預(yù)實驗，當?shù)螖?shù)為50時，水平集函數(shù)未能達到植入物邊緣算法便已停止，100次迭代與150次迭代算法都收斂到相同的邊緣，于是將水平集算法的迭代次數(shù)設(shè)置為100；Canny算子的二維高斯函數(shù)方差根據(jù)OpenCV中Canny算法的源代碼選擇為0.6；DRLSE算法中邊緣指示函數(shù)的高斯平滑窗口過大會使Canny算子處理后的梯度圖像變得平滑而導(dǎo)致邊緣泄露的情況增加，而過小則會使水平集算法的抗噪聲性能下降，因此將高斯濾波核設(shè)置為3×3，方差設(shè)置為1.4，然后分別使用DRLSE算法與本文算法對圖像序列進行分割，對比兩個算法的分割結(jié)果。

圖1顯示的是基于Canny算子的DRLSE乳腺植入物分割處理結(jié)果。圖2(a)是原始圖像，圖2(b)白色方框是人工勾畫的初始零水平集，圖2(c)的白色線是分割的結(jié)果。

選擇患者王某進行乳房重建后的圖像集里第13幅圖像，分別采用DRLSE算法和基于Canny算子的DRLSE算法進行分割的對比結(jié)果見圖3(a)和圖3(b)。

對患者盛某進行乳房重建后的圖像集里第9幅圖像，分別采用DRLSE算法和基于Canny算子的DRLSE算法進行分割的對比結(jié)果見圖4(a)和圖4(b)。

從視覺上可以看出在圖3(a)中采用DRLSE算法進行分割的植入物在左下角出現(xiàn)了邊緣泄漏的情況，并且在右下角出現(xiàn)了欠分割的情況。而圖3(b)則良好地完成了分割任務(wù)。圖4(a)植入物的左邊出現(xiàn)了明顯的邊緣泄漏，而采用本文算法得到的圖4(b)精確地將植入物的邊緣分割了出來。

表2 王某采用DRLSE與本文算法的分割精度

表3 盛某采用DRLSE與本文算法的分割精度

上述實驗表明本文算法在對乳腺植入物進行分割時比直接使用DRLSE算法更為精確，水平集函數(shù)更不容易出現(xiàn)邊界泄漏的情況。通過將患者王某與盛某的分割結(jié)果進行對比，顯示本文算法在對不同患者和不同的植入物形狀進行分割時，均比直接使用DRLSE效果更優(yōu)。

3.2 算法的抗噪聲性能對比實驗

圖5(a)-圖5(b)與圖6(a)-圖6(b)對比了DRLSE算法與本文算法的抗噪聲性能。分別使用圖3與圖4的原圖像加入均值為0、方差為0.01的高斯噪聲[17]。水平集算法迭代次數(shù)通過預(yù)實驗確定為200次；由于高斯噪聲的加入，為降低噪聲影響，將Canny算子的方差增大為1；其他參數(shù)保持不變。分別使用DRLSE算法與本文算法對加噪圖片進行分割，實驗結(jié)果如圖5、圖6所示。

觀察圖5(a)可見右下角外輪廓出現(xiàn)了欠分割的情況。與圖3(a)相比較，添加噪聲后的圖像右下角欠分割現(xiàn)象更加明顯。從圖5(b)可以看出本文算法良好地完成了植入物外輪廓的分割任務(wù)，與圖3(b)并沒有太大差異。圖6(a)的左邊部分出現(xiàn)了比圖4(a)更嚴重的邊緣泄漏。而圖6(b)與圖4(b)相比較則在分割區(qū)域右部出現(xiàn)了少量的邊緣泄漏的情況，并且在分割區(qū)域內(nèi)部出現(xiàn)了少量的孔洞。

骰子相似系數(shù)顯示出了本文算法在噪聲條件下的分割性能：王某的圖片在DRLSE算法下PDSC=0.912 1，本文算法下，PDSC=0.942 4。盛某的圖片在DRLSE算法下PDSC=0.763 1，本文算法下PDSC=0.810 7。

通過視覺觀察兩種算法在噪聲下的分割情況和骰子相似系數(shù)所顯示出來的分割精度數(shù)據(jù)，在對圖像增加噪聲后，算法的分割精度有所下降。但本文算法的骰子相似系數(shù)下降幅度并不大，并且依然可以獲得一個令人滿意的分割結(jié)果。可知本算法可以在一定程度上降低噪聲對分割產(chǎn)生的影響。

綜上，使用本文算法對乳腺植入物進行分割可以明顯提高分割的精度，減少邊緣泄漏的發(fā)生并提高算法的抗噪聲性能。

4 結(jié) 語

目前，國內(nèi)外醫(yī)院在對乳房重建手術(shù)患者的恢復(fù)情況進行評估以及獲得相應(yīng)科研數(shù)據(jù)時，多采用的是定性的觀察。如果需要定量分析，便要依靠經(jīng)驗豐富的醫(yī)生利用Mimics軟件對患者的核磁共振圖片進行手動閾值分割。對于切片數(shù)量較多的核磁共振圖片集工作量大且效率低下。基于此，本文通過將Canny算子與DRLSE算法結(jié)合，實現(xiàn)了半自動化的植入物分割，降低了醫(yī)生的工作量、提高了分割效率并且一定程度上消除了由于醫(yī)生經(jīng)驗豐富程度不同帶來的分割差異，在醫(yī)學圖像分割領(lǐng)域有一定的應(yīng)用前景。基于該算法的軟件現(xiàn)已經(jīng)在空軍軍醫(yī)大學西京醫(yī)院開始使用。

本文提出一種基于Canny算子和距離正則化水平集的乳腺植入物圖像分割算法。并且在算法實現(xiàn)過程中推導(dǎo)了基于Canny算子的DRLSE方程的數(shù)學解析式。通過對乳腺切除一期重建患者核磁共振圖像的PCL植入物分割對比實驗顯示，本文算法可以改善DRLSE算法邊緣泄漏的情況，擁有更高的分割精度和抗噪聲性能。提升了通過曲線演化方法對乳腺植入物進行分割的實用性。

但是，算法的迭代次數(shù)和初始零水平集的選取需要人工參與。一定程度上降低了算法的效率，在面對大量核磁共振切片時短板尤其明顯。因此，如何讓算法自動選擇初始水平集函數(shù)和迭代次數(shù)，是接下來需要研究的課題。