朱家明，李祥健，徐婷宜

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

0 引言

随着科技的不断进步，医学影像成为医生后期诊断的重要依据。常用于医学领域的影像学成像技术有：计算机断层扫描成像技术(CT)、核磁共振成像(MRI)、单光子辐射断层成像(SPECT)、正电子放射断层成像技术(PET)、光学断层纤维成像、核磁共振显微成像以及正电子断层显微成像等。其中MRI图像由于对含水分的软组织比较敏感，可以用于对脑部等组织进行解剖结构分析检测的场合[1]。然而MRI图像的成像质量也存在不足：受到设备、环境以及射频信号不均匀的影响，图像经常存在噪声和灰度不均的现象[2]。

医学图像分割是医学图像处理的一个关键技术[3]，是服务于临床医生的计算机辅助诊断的重要组成部分。其任务是从医学图像中提取感兴趣的目标，即解剖组织。国内外学者提出了很多图像分割方法，如阈值法、区域生长法；结合特定理论工具的方法，如模式识别法、人工神经网络法等；基于模糊分割的方法，如小波变换法；基于遗传算法的方法以及水平集方法等[4]。其中水平集方法是解决曲线演化问题的一种经典方法[5-7]。水平集方法最初由Osher和Sethian[8]提出，用于解决遵循热力学方程的火苗外形变化过程。Mumford和Shah[9]于1989年提出了M-S模型。在此基础上，Chan和Vese[10-11]分别于2001年和2002年提出了单水平集CV模型和多水平集CV模型。Li和Xu[12]于2010年提出了距离正则化水平集模型。

1 结合小波变换与中值滤波去噪

1.1 小波变换去噪

有效抑制噪声是进行图像分割的前提条件。小波变换[13]以其较好的时频特性在图像去噪领域受到了广泛的关注。小波变换方法有2个关键步骤，分别是阈值函数的选取和阈值的选取。

阈值函数直接影响到信号去噪的效果。常用的阈值函数包括：硬阈值法、软阈值法及模平方法等[14]。

设ω为小波系数，wλ为加入阈值后小波系数的大小，λ为阈值。

硬阈值函数如下所示：

(1)

软阈值函数如下所示：

(2)

模平方函数如下所示：

(3)

硬阈值法去除高斯噪声效果很好，但是图像灰度变化较大的邻域出现了振荡，即振铃现象。软阈值法去除高斯噪声效果较好，但是图像边缘出现了模糊。模平方法效果约等于硬阈值法与软阈值法的折中。

针对上述阈值函数的不足，本文提出S型分段阈值函数，其特点是保持函数值的连续与光滑。在保持较好的去噪性能前提下，既避免了硬阈值函数的振铃现象，又避免了软阈值函数的边缘模糊现象。

S型分段阈值函数如下：

(4)

函数图像如图1所示。

图1 S型阈值函数图像Fig.1 S-type threshold function

阈值的选取对于去噪效果同样很关键。本文采用Donoho和Johnstorn[15]提出的统一阈值：

(5)

式中，δ为噪声标准方差，N为信号的尺寸或长度。

1.2 中值滤波去噪

中值滤波去噪法是一种对图像局部进行平均平滑的技术[16]，属于一种非线性滤波。中值滤波的基本原理是将邻域中的所有像素值按灰度级进行排列，取队列的中间值为输出像素值。它处理图像的效果主要取决于2个要素：邻域的空间范围和模板所覆盖的像素数。时间序列x(t)的中值滤波可以表示如下：

y(t)=MED{x(t-k),…,x(t),…,x(t+k)}，

式中，中值滤波器的长度为2k+1。改变滤波器的长度或者模板可改变滤波器的性能。

2 改进的双水平集分割方法

2.1 单水平集模型描述

单水平集是依据水平集函数φ的取值将图像划分为目标和背景两相。其适用于分割单个目标和灰度均匀目标的任务。图2为单水平集模型原理图。可见水平集函数φ将图像分割成2个互不重合的部分。

图2 单水平集区域划分图Fig.2 Single level set area division diagram

曲线拟合能量泛函定义为：

E(c1,c2,c)=μ·Length(c)+λ1∬Ω1(u(x,y)-c1)2dxdy+

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

(6)

式中，μ,λ1,λ2为正常数，表示各能量项的权重系数，采用水平集φ后能量泛函为：

E(c1,c2,φ)=μ·Length(φ)+λ1∬Ω1H(φ)(u(x,y)-c1)2dxdy+λ2∬Ω2(1-H(φ))(u(x,y)-c2)2dxdy，

(7)

式中，H(*)为Hesviside函数。

曲线演化的目标是使得上述能量泛函取得极小值。根据多元函数求极值的理论，求泛函对c1的一阶偏导数，并令它为0，可得c1，同理可得c2。公式如下：

(8)

根据求函数极值的梯度下降法，令水平集函数φ的导数等于能量泛函的负梯度得到水平集函数的演化方程如下：

(9)

2.2 双水平集模型描述

为了将水平集方法推广到分割多目标与灰度不均目标的任务，Chan和Vese提出了多相CV模型，即可通过n条水平集函数将图像分割成2n个互不重合的部分。

图3为四相CV模型原理图。水平集函数φ1和φ2将图像分割成Ω11，Ω12，Ω21，Ω22四个区域，它们之间互不重合。

图3 四相CV区域划分图Fig.3 Division diagram of four phase CV region

定义能量函数如下：

(10)

式中，c11，c12，c21，c22分别对应Ω11，Ω12，Ω21，Ω22各区域的灰度均值，它们的计算方法如下：

(11)

考虑到演化过程中出现的重新初始化[15]的问题，在原来的能量函数中增加惩罚项p(φ)，表达式如下：

(12)

式中，ν为能量惩罚项的系数。

新的能量函数更新如下：

Ep=E(c,φ)+p(φ)。

(13)

针对传统的水平集函数演化方程在演化后期容易出现演化速度放缓停滞等问题，提出了改进的演化方程。增加一组加速因子K1,K2作为乘法因子，自适应地调整加快演化速度。

(14)

图4 加速因子Ki的函数图像Fig.4 Function image of acceleration factor Ki

其中，

(15)

根据Euler-Lagrange方程和梯度下降法，最终的水平集函数φ1,φ2的演化方程为：

(16)

(17)

2.3 算法流程

① 首先加载原始图像；

② 其次用小波变换去除高斯噪声；

③ 然后用中值滤波去除椒盐噪声；

④ 最后用多水平集快速分割算法对图像进行分割；

⑤ 输出分割效果图。

3 实验结果与分析

本文实验环境：软件 Windows10, MATLAB 7.0;硬件Huawei MateBook X Pro, CPU CORE i7-8550U 1.80 GHz，RAM 8.0 G。

为了比较中值滤波和本文算法对高斯噪声图像的去噪效果，在脑部MRI图像中加入高斯噪声。分别用中值滤波和本文算法进行去噪，由图5(b)和图5(c)可以看出，中值滤波只能去除部分高斯噪声，而本文算法能够完全去除噪声。

图5 中值滤波与本文算法去噪效果对比图Fig. 5 Comparison of denoising effect between median filter and the algorithm in this paper

为了比较小波变换和本文算法对椒盐噪声的去噪效果，分别用两种算法对含有椒盐噪声的图6(a)进行去噪，通过对比分割效果图可以看出本文算法可以完全去除图像中的椒盐噪声。需要说明：如果椒盐噪声密度很高,那么一次中值滤波可能无法消除全部噪声。此时应该进行两次或多次中值滤波。

图6 小波变换与本文算法去噪效果对比图Fig.6 Comparison of denoising effect between wavelet transform and the algorithm in this paper

图7 传统水平集分割效果图Fig.7 Traditional level set segmentation renderings

图8 本文模型分割效果图Fig.8 Model segmentation renderings in this paper

通过对比图7和图8可知，对于传统的双水平集算法，当迭代到30次时，分割的效果只显示出部分轮廓；当迭代到100次时，分割的效果基本理想，显示出比较准确的轮廓。而对于快速水平集算法，当迭代到第30次时，分割的效果已经达到理想效果。继续迭代到100次时，分割的结果也没有明显变化。实验结果验证了快速算法的有效性。

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

(18)

JS值越大表明表示分割效果越好，表1为本文算法与传统水平集算法的JS指标，可看出本文算法的分割效果更好。

表1 2种算法的JS指标
Tab.1 JS index of two algorithms

图像组织分割算法传统水平集本文算法脑部MRI图像1白质0.780.86灰质0.730.85

4 结束语

本文结合小波变换与中值滤波2种去噪算法的特点，对MRI图像进行去噪预处理，兼顾了较好的去噪效果和保留图像细节的目标。在传统水平集方法中增加了自适应的加速因子以加快图像分割速度，并且通过在能量函数中增加惩罚项解决曲线在演化过程中重新初始化的问题。实验表明，本文的方法可以快速有效地分割含有高噪声、灰度不均的MRI图像。