基于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算法边缘泄漏的情况，拥有更高的分割精度和抗噪声性能。提升了通过曲线演化方法对乳腺植入物进行分割的实用性。

但是，算法的迭代次数和初始零水平集的选取需要人工参与。一定程度上降低了算法的效率，在面对大量核磁共振切片时短板尤其明显。因此，如何让算法自动选择初始水平集函数和迭代次数，是接下来需要研究的课题。