李建華 王藝衡* 楊來(lái)俠 李素麗 張聚良

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

0 引 言

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

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

1 基本算法簡(jiǎn)介

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

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

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

(1)

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

(2)

div是計(jì)算散度，dp的定義如下：

(3)

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

(4)

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

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

(5)

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

(6)

1.2 Canny算子簡(jiǎn)介

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

2 算法設(shè)計(jì)

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

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

本文算法的處理過(guò)程如下：

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

(7)

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

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

(8)

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

Step2gG(x,y)是進(jìn)行高斯平滑后的圖像，對(duì)其求一階導(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)

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

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)

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

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

(12)

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

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

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

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

3 實(shí) 驗(yàn)

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

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

3.1 分割精度對(duì)比實(shí)驗(yàn)

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

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

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

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

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

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

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

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

3.2 算法的抗噪聲性能對(duì)比實(shí)驗(yàn)

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

觀察圖5(a)可見(jiàn)右下角外輪廓出現(xiàn)了欠分割的情況。與圖3(a)相比較，添加噪聲后的圖像右下角欠分割現(xiàn)象更加明顯。從圖5(b)可以看出本文算法良好地完成了植入物外輪廓的分割任務(wù)，與圖3(b)并沒(méi)有太大差異。圖6(a)的左邊部分出現(xiàn)了比圖4(a)更嚴(yán)重的邊緣泄漏。而圖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。

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

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

4 結(jié) 語(yǔ)

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

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

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