结合新型边界扩展函数的水平集医学图像分割方法

朱家明，许 诺，张倚天，包 清，杭 俊

(扬州大学 信息工程学院，江苏 扬州 225127)

0 引言

随着科技的发展，很多医生把医学影像作为疾病诊断的重要依据之一。由于核磁共振成像(MRI)图像具有对含水分的软组织较为敏感的特性，经常用于对脑部组织进行解剖结构的分析与检测[1]。但是MRI图像的成像也有缺点,例如由于受到设备、环境和不均匀的射频信号影响，图像大多存在噪声和灰度不均的问题[2]。

在医学图像处理领域，图像分割是一个重要的技术[3]，其任务是从图像中提取有诊断价值的目标，是用于医生临床诊断的计算机辅助系统的重要组成部分。对于图像分割，各国学者提出了多种方法，例如：阈值法、模式识别法、区域生长法、神经网络法、小波变换法、模糊分割法、水平集方法[4]以及遗传算法等。其中水平集方法是曲线演化问题的一种有效的求解方法[5-7]，其最初由Osher和Sethian[8]提出，用于分析火苗的外形变化过程。在1989年，Mumford和Shah[9]提出了M-S模型。Chan[10]于2001年提出了单水平集CV模型，Vese[11]于2002年提出了多水平集CV模型，Li和Xu[12]于2010年提出了距离正则化水平集模型。为了解决MRI图像的分割问题，本文提出了一种结合阈值法预分割的水平集分割方法。考虑到MRI图像具有的高斯噪声，先对其作低频滤波处理。理论分析和实验结果都验证了该方法的有效性。

1 频域滤波去噪

1.1 高斯噪声与频域滤波

作为图像分割的前提条件，要求对噪声先作有效抑制。高斯噪声是医学图像中广泛存在的一种噪声，通常是由于电子元件的敏感性以及各种外界的电磁干扰造成的。消除高斯噪声常用的方法有空间滤波与频域滤波。因为噪声信号一般都是高频信号，当采用低通滤波时，选择合适的截止频率就能够滤除大部分噪声，同时保留图像的边缘和细节。而采用低通滤波的前提就是先对图像进行傅里叶变换，使之转换成频域模型；然后针对频域模型设计低通滤波器进行滤波；最后对滤波后的频域模型进行傅里叶反变换，得到去噪的图像。

计算机处理的都是离散化的系统，一幅灰度图像可以表示为二维的离散函数，离散函数的二维傅里叶变换与傅里叶反变换公式分别为：

(1)

(2)

1.2 低通滤波

低通滤波方法有巴特沃斯滤波、切比雪夫滤波以及理想低通滤波等。其中理想低通滤波器属于假想的低通滤波器，对于高频(即高于截止频率)信号完全截止，而对于低频信号完全无失真传输。由于理想低通滤波去噪的效果要优于其他低通滤波，因此本文采用理想低通滤波去除噪声，如图1所示。

(a) FFT频谱

(b) FFT频谱

(c) ButterworthLPF频谱

(d) 理想LPF频谱

(e) 原始图像

(f) 噪声图像

(g) ButterworthLPF图像

(h) 理想LPF图像

2 改进的水平集分割方法

2.1 阈值法的预分割

用阈值法先预分割，即设定水平集函数初值。水平集分割法步骤：先定一个初值(初始轮廓)，再按照演化方程进行演化,即不断更新水平集函数取值。相应改变了零水平集的轮廓曲线，完成了图像分割。当达到设定标准，就停止演化并输出结果。因为水平集的初值会对分割结果有一定的影响，所以选取合适的初值是有必要的。本文采用阈值法预分割(定水平集初值)，用整幅图像灰度区间的黄金分割数作为阈值，即取权值α=0.382。阈值的计算公式如下：

Ith=Ilow+0.382(Ihigh-Ilow)，

(3)

式中，Ith为阈值，Ilow为整幅图像的灰度值的下限，Ihigh为灰度值的上限。其中默认目标灰度值比背景灰度值小，反之如果目标灰度值比背景灰度值大，则取α=0.618，计算公式如下：

Ith=Ilow+0.618(Ihigh-Ilow)。

(4)

阈值法确定的初始轮廓图如图2所示。

(a) α=0.382

(b) α=0.618

(c) α=0.85

2.2 水平集模型描述

水平集方法是依据水平集函数φ的取值将图像划分为目标和背景两相。水平集函数定义为符号距离函数。即图上一点(x,y)与目标轮廓曲线C的距离。表达式如下：

(5)

由式(5)可见，水平集函数取值为零的点的集合(零水平集)即为目标的边缘轮廓。水平集函数原理如图3所示，从图中可见水平集函数φ将整个图像分割成互不重合的两个部分。

图3 水平集区域划分图Fig.3 Level set area division diagram

水平集模型主要包括两大类模型：一类是基于边缘信息的模型，如Snake模型；另一类是基于区域信息的模型，如M-S模型、C-V模型和变分水平集等。本文考虑结合边缘信息和区域信息的模型。

2.3 边界扩展函数

本文设计一种新型的边界扩展函数，并将它结合到传统的能量函数中，用于解决目标边界划分细节不准确的问题。基本设想如下：在目标轮廓线附近，当目标外部灰度接近目标内部的灰度时，就使得目标轮廓线扩张，把外面这部分包含到目标内部来。该边界扩展函数设计如下：

(1-H(φ))dxdy，

(6)

式中，u(x,y)表示目标外部某一点的灰度，cδ表示该点附近的目标局部的平均灰度。

。

(7)

水平集函数φ的演化方程为：

(8)

2.4 能量函数与演化方程

主动轮廓模型的能量泛函定义为：

E(C)=∮Ωg(|∇I(C(s))|)ds。

(9)

对式(9)引入Hesviside函数，结合格林公式，可改写成：

Es(φ)=∬Ωg(x,y)‖∇H(φ)‖dxdy，

(10)

式中，φ表示水平集函数，g(x,y)为反映边缘特征的函数，表达式如下：

p≥1。

(11)

曲线演化的目标是使能量泛函取得极小值。根据E-L方程和梯度下降法，解得演化方程为：

(12)

C-V模型的能量泛函定义为：

E(c)=μ·length(C)+

λ1∬Ω1(u(x,y)-c1)2dxdy+

λ2∬Ω2(u(x,y)-c2)2dxdy，

(13)

式中，μ、λ1、λ2为正的常数，u(x,y)为图像灰度函数。引入Hesviside函数后能量泛函为：

Ecv(φ)=μ∬Ω‖∇H(φ)‖dxdy+

λ1∬ΩH(φ)(u(x,y)-c1)2dxdy+

λ2∬Ω(1-H(φ))(u(x,y)-c2)2dxdy，

(14)

对应的演化方程为：

(15)

将泛函(14)对c1求偏导，令它等于0可求得c1。同理可求得c2,公式如下：

(16)

(17)

为了避免在演化过程中出现重新初始化的问题[15]，在能量泛函(14)中增加惩罚项Ep(φ):

(18)

式中，ν为正的常数。综上所述，能量函数定义为：

E(φ)=Ecv(φ)+Ep(φ)+Es(φ)+Eex(φ)=

μ∬Ω‖∇H(φ)‖dxdy+

λ1∬ΩH(φ)(u(x,y)-c1)2dxdy+

λ2∬Ω(1-H(φ))(u(x,y)-c2)2dxdy+

α∬Ωg(x,y)‖∇H(φ)‖dxdy+

。(19)

对应的演化方程为：

λ1(u-c1)2+λ2(u-c2)2+

(20)

式中，Δ表示拉普拉斯算子。

2.5 算法流程

算法流程为：① 加载原始图像；② 用理想低频滤波消除噪声；③ 用阈值法进行预分割；④ 用水平集方法进行图像分割；⑤ 输出分割结果图。

3 实验结果与分析

实验环境：硬件Huawei MateBook X Pro,RAM 8.0 G，CPU CORE i7-8550U 1.80 GHz; 软件 Windows10中文版,MATLAB 7.0。选取一幅脑部MRI图像进行去噪预处理，迭代次数为200次。参数设置为：α=3,μ=5,v=1.5,λ1=λ1=10-6,采样时间间隔Δt=0.02。分别采用C-V算法和本文算法来分割去噪后的MRI图像，分割结果如图4所示。由图可见，本文的算法分割的效果优于C-V算法。

(a) 原始图像

(b) C-V算法

(c) LSE算法

采用Jaccard Similarity(JS)指标来对分割结果进行评价[13]：

(21)

指标值越大代表分割效果越好，表1为C-V算法与本文算法的JS指标，数据显示本文算法的分割效果优于C-V算法。

表1 两种算法的JS指标

4 结束语

本文提出一种结合新型边界扩展函数的水平集图像分割方法。首先用理想低频滤波对图像进行去噪处理，同时兼顾去噪效果和保留细节两个目标，然后用阈值法确定目标的初始轮廓。提出新型边界扩展函数，并整合到能量函数中。在其中增加惩罚项，避免了在演化过程中的重新初始化。实验结果表明，该方法能够准确地分割含有高斯噪声的灰度不均的MRI图像。