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

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

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

0 引 言

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

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

1 基本算法簡介

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

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

(5)

δε是Heaviside函數，其定義如下：

(6)

1.2 Canny算子簡介

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

2 算法設計

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

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

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

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

(7)

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

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

(8)

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

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

(9)

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

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

(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算子得到的邊緣代入水平集的邊緣指示函數中。

邊緣指示函數g解析式為：

(11)

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

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

(12)

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

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

Step6將水平集函數輸出，得到最終的分割結果。

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

3 實 驗

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

表1 實驗環境與圖像相關參數

3.1 分割精度對比實驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結 語

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

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

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