随着医学成像技术及图像处理技术的发展，许多疾病的分析、诊断和治疗都将医学图像作为定性和定量分析组织性能的重要依据，同时如图像的去噪、分割、配准、分析和重建等图像处理方法也成为医学图像处理的主要技术。其中，医学图像的分割是医学图像处理中的经典问题，也是医学图像处理的关键技术。然而医学图像在采集过程中受到光照和射频线圈等的影响，同一组织常具有不同灰度，且伴有强噪声。这些缺陷使得基于灰度均匀假设的常用分割方法效果不理想甚至失败，进而影响医学诊断。因此，灰度不均匀带噪图像的精确分割成为医学图像处理领域的重点和难点。

过去几十年里，基于活动轮廓模型图像分割方法广泛应用于图像分割^[1-8]，并取得较好结果。现有的活动轮廓模型大致分为两类：基于边缘的和基于区域的；基于边缘的活动模型如最初的GAC(geodesic active contours)^[2]，利用图像的梯度引导曲线演化。该类方法没有进行图像灰度均匀性假设，可以分割灰度不均匀图像，但对轮廓初始化和噪声相当敏感，计算开销也过于昂贵，限制了该类方法的使用。基于区域的活动轮廓模型如最经典的C-V^[3]，利用灰度等区域描述符拟合每个区域来引导轮廓曲线的向边界靠近，该类方法在处理弱边界、噪声图像时比基于边缘的方法具有更理想的性能，而且初始化不敏感。但是该方法是在灰度均匀的假设下进行的，因此处理不了灰度不均匀图像。医学图像对比度低，质量差，为实现医学图像的精确分割就要求分割算法要同时具有如下特性。图像灰度的充分描述；具有较强的抗噪性；可克服灰度不均匀因素。

针对灰度不均图像的分割问题，国内外的学者提出了许多基于局部区域信息的活动轮廓模型^[3-11]。2002年Vese和Chan^[3]提出分段平滑的PS(piecewise smooth)模型，具有一定的处理灰度不均匀图像的能力。但PS模型的计算量巨大，限制了它的实际应用。2007年Li等人^[9]提出基于局部二值拟合能量活动轮廓模型，该模型利用一定形式的核函数选定某像素点附近的局部区域，用两个常数分别估计轮廓曲线内外局部区域内的灰度分布，选用均值作为局部区域描述子驱动轮廓曲线向目标边界靠近从而完成分割。因为使用了局部信息，LBF(local binary fitting)能够分割灰度不均匀的图像。最近，一些基于LBF模型的算法相继提出^[12-13]。这类方法在一定程度上克服了灰度不均匀性，但因为仅将灰度均值作为描述子，不足以为精确的分割提供信息。特别是在图像均值相同而纹理不同的情况下，往往导致分割失败。2009年Wang等人^[10]提出了局部高斯分布拟合能量的活动轮廓模型，该模型在贝叶斯框架下使用均值和方差两个描述子，在一定程度上解决了上述问题，但在强噪声下效果并不理想。

针对上述情况，提出一种基于多个描述子局部熵、均值和局部熵的活动轮廓模型，在能量泛函中引入灰度不均匀因子克服灰度不均现象，该算法简记为MDAC。实验结果表明本文方法能够分割具有灰度不均匀现象的低对比度医学图像，完成偏移场的估计，在较强噪声下获得较高的分割精度。

Li等人^[14]根据的灰度不均匀图像的局部灰度聚类特点，引入偏移因子$ b $，提出了局部灰度聚类能量，定义了能量泛函

$ \begin{array}{*{20}{c}} {{E^{{\rm{LIC}}}}\left( {\phi ,b,{z_1},{z_2}} \right) = }\\ {\nu \int_\Omega {\frac{1}{2}{{\left( {\left| {\nabla \phi } \right| - 1} \right)}^2}{\rm{d}}x} + \mu \int_\Omega {{\delta _\varepsilon }\left( \phi \right)\left| {\nabla \phi } \right|{\rm{d}}x} + }\\ {\int {\left[ {\int {K\left( {x - y} \right){{\left( {I\left( y \right) - b\left( x \right){z_1}} \right)}^2}{H_\varepsilon }\left( {\phi \left( y \right)} \right){\rm{d}}y} } \right]{\rm{d}}x} + }\\ {\int {\left[ {\int {K\left( {x - y} \right){{\left( {I\left( y \right) - b\left( x \right){z_2}} \right)}^2}} \times } \right.} }\\ {\left. {\left( {1 - {H_\varepsilon }\left( {\phi \left( y \right)} \right)} \right){\rm{d}}y} \right]{\rm{d}}x} \end{array} $

(1)

式中，${z_1} $，$ {z_2} $为均值，$ ϕ $为水平集函数，K(·)为核函数。

LIC模型能够在完成图像分割同时，实现图像的偏移场校正。但因仅使用灰度均值拟合图像，在强噪声下不能满足医学图像的精确分割。

Wang等人^[10]根据最大后验概率准则，将灰度概率密度函数引入到能量泛函中，并使用了高斯核函数，LGDF能量泛函定义为

$ \begin{array}{*{20}{c}} {{E^{{\rm{LGDF}}}}\left( {\phi ,{u_1},{u_2},{\sigma _1},{\sigma _2}} \right) = }\\ {\nu \int_\Omega {\frac{1}{2}{{\left( {\left| {\nabla \phi } \right| - 1} \right)}^2}{\rm{d}}x} + \mu \int_\Omega {{\delta _\varepsilon }\left( \phi \right)\left| {\nabla \phi } \right|{\rm{d}}x} - }\\ {\int {\left[ {\int {K\left( {x - y} \right){p_{1,x}}\left( {I\left( y \right)} \right){H_\varepsilon }\left( {\phi \left( y \right)} \right){\rm{d}}y} } \right]{\rm{d}}x} - }\\ {\int {\left. {\left[ {\int {K\left( {x - y} \right){p_{2,x}}\left( {I\left( y \right)} \right)\left( {1 - {H_\varepsilon }\left( {\phi \left( y \right)} \right)} \right){\rm{d}}y} } \right]} \right]{\rm{d}}x} } \end{array} $

(2)

式中，$ {p_{i, x}}(I(y))$为子邻域$ \left\{ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i} \cap {\mathit{\boldsymbol{O}}_x}} \right\}_{_{i = 1}}^{^N} $中像素的灰度概率密度函数。

该模型使用均值和方差两个区域描述子拟合图像灰度分布，在分割精度及应用范围上有所改善，但在分割精度要求高的医学图像时效果仍不理想，且没有实现偏移场的估计。

MDAC模型利用局部熵完善图像灰度信息，提高精度，改善抗噪性能。引入LIC的灰度不均匀图像模型，在完成分割的同时，实现偏移场的估计。

根据Shannon信息论的定义^[15]，概率系统由$ n $个事件组成，第$ i $个事件出现的概率为$ p_i $，则该事件具备的信息量为

$ I = - {\log _2}\left( {{p_i}} \right)\;\;\;\;\left( {0 \le {p_i} \le 1} \right) $

(3)

系统全部事件含信息量的期望定义为熵

$ H = - \sum\limits_{i = 1}^n {\left( {{p_i}{{\log }_2}\left( {{p_i}} \right)} \right)} $

(4)

Shiozaki提出用SHANNON熵的定义来表示图像的熵^[16]。一幅大小为M×N图像的熵定义为

$ H = - \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {{p_{ij}}{{\log }_2}\left( {{p_{ij}}} \right)/{{\log }_2}\left( {MN} \right)} } $

(5)

式中，$ {p_{ij}} $是图像的灰度分布。目前图像熵没有统一的定义，关于灰度分布的选取通常有两种形式。一种是图像灰度直方图统计，表现了像素灰度在灰度级上的分布情况，适用于图像的阈值分割方法。另一种是像素灰度的归一化，即单个像素灰度值对整幅图像灰度和的比值，适用于水平集分割方法。故选择后者，定义式为

$ {p_{ij}} = f\left( {i,j} \right)/\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {f\left( {i,j} \right)} } $

(6)

式中，$ f(i, j) $是像素的灰度。

图像的不同区域具有不同的信息量，为了研究图像的不同部分，需计算相应的熵值。定义一个空间上连续的区域$ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x} \subset \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $，则该区域的局部熵为

$ \begin{array}{*{20}{c}} {H\left( {x,{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} \right) = - \frac{1}{{{{\log }_2}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} \right|}} \times }\\ {\int_{{\Omega _x}} {p\left( {y,{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} \right)\left( {p\left( {y,{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} \right) - 1} \right){\rm{d}}y} } \end{array} $

(7)

式中，$ p(y, {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}) $是灰度分布函数。

图像的局部熵反映图像信息丰富程度的特性，适合于表征图像分割区域的同质性，结合其他描述子对图像灰度分布进行联合拟合，可以提高图像描述的精准度。局部熵的取值是由窗口内多个像素点共同决定的，单个像素的灰度值对熵的影响不大，因而对单个噪声点也不敏感，局部熵具有一定的滤波作用。另外，灰度概率分布的归一化处理也具有平滑噪声作用。医学图像通常具有较强的噪声，利用局部熵作为描述子可以有效的改善分割算法的抗噪性能。此外，局部熵对光照稳定，具有亮度不变性，引入局部熵使得分割模型对偏移场具有良好的鲁棒性。

MDAC模型采用图像均值、方差和局部熵作为描述子，分割同时完成偏移场的估计，在贝叶斯框架下定义模型的数据拟合项。一点的数据能量为

$ {E_x} = {H_x} \cdot \sum\limits_{i = 1}^N {{\lambda _i}\int_{{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i} \cap {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} { - {{\log }_2}\left[ {{p_{i,x}}\left( {I\left( y \right)} \right)} \right]{\rm{d}}y} } $

(8)

式中，为$ {\lambda _i} $为调节系数，$ {H_x} $为$ x $点邻域的局部熵，$ {p_{i, x}}(I(y)) $为子邻域$ \left\{ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i} \cap {\mathit{\boldsymbol{O}}_x}} \right\}_{_{i = 1}}^{^N} $中像素的灰度概率密度函数。这里认为概率密度函数服从高斯分布，均值为$z_i $，方差为$ {\sigma _i}(x) $，即

$ \begin{array}{*{20}{c}} {{p_{i,x}}\left( {I\left( y \right)} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} {\sigma _i}\left( x \right)}} \times }\\ {\exp \left( { - \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2{\sigma _i}{{\left( x \right)}^2}}}} \right)} \end{array} $

(9)

一点的数据能量整理为

$ {E_x} = {H_x} \cdot \sum\limits_{i = 1}^N {{\lambda _i}\int_{{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i} \cap {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x}} {\left[ \begin{array}{l} {\log _2}{\sigma _i}\left( x \right) + \\ \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2\sigma _i^2\left( x \right)}} \end{array} \right]{\rm{d}}y} } $

(10)

为方便数值计算，引入如下窗口函数$ k(u) $代替式(10)中$ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_x} $限定的像素$ x $的邻域，即

$ k\left( u \right) = \left\{ \begin{array}{l} a\;\;\;\left| u \right| \le \rho \\ 0\;\;\;\left| u \right| > \rho \end{array} \right. $

(11)

式中，$ a $为正常数保证$ \smallint k\left( u \right) = 1 $。一点数据能量为

$ \begin{array}{*{20}{c}} {{E_x} = {H_x} \cdot }\\ {\sum\limits_{i = 1}^N {{\lambda _i}\int_{{\Omega _i}} {k\left( {x - y} \right) \cdot \left[ \begin{array}{l} {\log _2}{\sigma _i}\left( x \right) + \\ \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2\sigma _i^2\left( x \right)}} \end{array} \right]{\rm{d}}y} } } \end{array} $

(12)

最小化局部能量$ E_x $可以得到最优的局部分割，为了分割整幅图像，对$ E_x $在整个图像域内积分来定义模型的数据拟合项

$ \begin{array}{*{20}{c}} {{E_d} = \int_\mathit{\Omega } {{H_x}} \cdot }\\ {\sum\limits_{i = 1}^N {{\lambda _i}\int_{{\Omega _i}} {k\left( {x - y} \right) \cdot \left[ \begin{array}{l} {\log _2}{\sigma _i}\left( x \right) + \\ \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2\sigma _i^2\left( x \right)}} \end{array} \right]{\rm{d}}y{\rm{d}}x} } } \end{array} $

(13)

式中，利用局部熵$ H_x $量化像素邻域内的灰度变化，并作为能量项的权重。当$ x $是邻域内概率密度很小的噪点时，$ H_x $的值较小，对整个能量项影响很小；反之，则加快轮廓曲线向目标边界的演化速度。

在水平集方法中，进行两相医学图像的分割时，利用一个水平集函数$ ϕ $将图像分成前景和背景两个区域，故将式(13)中的积分上限N取2。本法对应水平集为

$ \begin{array}{*{20}{c}} {{E_d}\left( {\phi ,b,Z,{\rm{\Theta }}} \right) = }\\ {\int_\Omega {{H_x} \cdot \left( {\sum\limits_{i = 1}^2 {{\lambda _i}\int_{{\Omega _i}} {k\left( {x - y} \right) \cdot } } } \right.} }\\ {\left. {\left[ \begin{array}{l} {\log _2}{\sigma _i}\left( x \right) + \\ \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2\sigma _i^2\left( x \right)}} \end{array} \right]{M_i}{\rm{d}}y} \right){\rm{d}}x} \end{array} $

(14)

式中，$ {M_1}\left( \phi \right) = H(\phi ) $，$ {M_2}\left( \phi \right) = 1-H(\phi ) $。

在实际应用中，Heaviside函数采用光滑版本

$ {H_\varepsilon }\left( z \right) = \frac{1}{2}\left( {1 + \frac{2}{{\rm{ \mathsf{ π} }}}\arctan \left( {\frac{z}{\varepsilon }} \right)} \right) $

(15)

其导数为

$ {\delta _\varepsilon }\left( z \right) = \frac{1}{{\rm{ \mathsf{ π} }}} \cdot \frac{\varepsilon }{{{\varepsilon ^2} + {z^2}}} $

(16)

MDAC模型的完整水平集公式定义为

$ \begin{array}{*{20}{c}} {E\left( {\phi ,b,Z,{\rm{\Theta }}} \right) = {E_d}\left( {\phi ,b,Z,{\rm{\Theta }}} \right) + }\\ {\mu R\left( \phi \right) + \upsilon L\left( \phi \right)} \end{array} $

(17)

式中，第2项为水平集正则项，在演化过程中保证水平集的稳定性，即保持为符号距离函数

$ R\left( \phi \right) = \int_\Omega {\frac{1}{2}{{\left( {\left| {\nabla \phi } \right| - 1} \right)}^2}{\rm{d}}x} $

(18)

第3项为水平集长度项，保证轮廓的光滑

$ L\left( \phi \right) = \int_\Omega {{\delta _\varepsilon }\left( \phi \right)\left| {\nabla \phi } \right|{\rm{d}}x} $

(19)

使用标准梯度下降流方法实现能量的最小化。能量泛函中有多个变量，在求解一个变量时分别固定其他变量，得到求解公式为

$ {z_i} = \frac{{\int {\left( {k\frac{b}{{\sigma _i^2}}} \right)I\left( y \right){M_i}{\rm{d}}y} }}{{\int {\left( {k\frac{b}{{\sigma _i^2}}} \right){M_i}{\rm{d}}y} }};i = 1,2 $

(20)

式中，${\sigma _i} $为轮廓线内外方差，$ M_i $为成员函数，$ k $为核函数。

$ \sigma _i^2\left( x \right) = \frac{{\int {k\left( {x - y} \right){{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}{M_i}{\rm{d}}y} }}{{\int {k\left( {x - y} \right){M_i}{\rm{d}}y} }} $

(21)

偏移场$ b(x) $为

$ b\left( x \right) = \frac{{\sum\limits_{i = 1}^2 {\int_\Omega {k\left( {x - y} \right)\frac{{{z_i}}}{{\sigma _i^2}}I\left( y \right){M_i}{\rm{d}}y} } }}{{\sum\limits_{i = 1}^2 {\int_\Omega {k\left( {x - y} \right)\frac{{{z_i}}}{{\sigma _i^2}}{M_i}{\rm{d}}y} } }} $

(22)

固定$ b $, $Z $, Θ，得到水平集演化方程

$ \begin{array}{*{20}{c}} {\frac{{\partial \phi }}{{\partial t}} = \mu \left( {\Delta \phi - {\rm{div}}\left( {\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right)} \right) + }\\ {v{\delta _\varepsilon }\left( \phi \right){\rm{div}}\left( {\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right) - {\delta _\varepsilon }\left( \phi \right)\left( {{\lambda _1}{e_1} - {\lambda _2}{e_2}} \right)} \end{array} $

(23)

式中，$ e_i $为

$ {e_i} = \int {{H_x} \cdot k\left( {x - y} \right)\left[ \begin{array}{l} {\log _2}{\sigma _i}\left( x \right) + \\ \frac{{{{\left( {I\left( y \right) - b\left( x \right){z_i}} \right)}^2}}}{{2{\sigma _i}\left( x \right)}} \end{array} \right]{\rm{d}}y} $

(24)

选取了合成图像和医学图像对MDAC的有效性进行测试。针对分割精度及抗噪性能，与经典的LIC和LGDF算法进行了对比实验。实验在PC机上MATLAB2013环境下进行，配置为Inter CPU 3.20，4 GB内存，windows7操作系统。实验中MDAC方法的实验参数选取为$ μ=1.0 $，$ ρ=21 $，$ r=4.5$，$υ=0.05×255×255 $，$ Δt=0.1 $，$ {\lambda _1} $和$ {\lambda _2} $依据图像和分割效果在值1附近做微调。

为了验证算法能够很好地分割带噪、灰度不均匀图像，且可以同时完成灰度纠偏，设计了以下实验来验证算法的有效性。

选取带有灰度不均的人工合成图像和医学心脏图像，在加噪的情况下利用本算法对其进行分割，实验结果如图 1、图 2所示。其中，人工合成图像$ {\lambda _1} = {\lambda _2} = 1 $，医学心脏图像中${\lambda _1} = 1 $，$ {\lambda _2} = 1.5$。

很显然，在明显噪声的情况下，本文算法可以成功地分割图像，同时完成灰度纠偏。与图 1(a)相比，纠偏图像图 1(d)中的每一部分的灰度更加均匀。

图 2(a)所示为两幅原始图像的直方图，图 2(b)是两幅纠偏图像的直方图。由纠偏图像直方图可见，图像中的两相内容很好地体现为界限清晰的两个峰。而原始图像的直方图中由于灰度不均匀的存在，多出了一个峰，且各峰之间分界不如纠偏后图像的分解效果好。

为了测试算法的分割精度和抗噪性能，利用本文算法和经典的LIC和LGDF算法分割相同的医学图像，对实验结果进行比较。本实验选取的是血管和脑部医学影像，实验中，血管影像加入了均值为零方差为0.001的高斯白噪声，${\lambda _1} = {\lambda _2} = 1 $，脑MR影像加入了均值为零方差为0.005的高斯白噪声，${\lambda _1} = {\lambda _2} = 1 $。实验结果如图 3所示。从图 3(b)可以看出，由于噪声的干扰LIC方法有很多错误的划分。图 3(c)中LGDF的分割结果要强于LIC，但仍然存在着很多错误。图 3(d)为MDAC方法的分割结果，可见MDAC方法具有最好的分割结果，可以正确分割出图像前景和背景。

为了定量评估算法的性能，使用相似系数DSC作为分割精度的指标。分割结果越接近真实结果，DSC的值越大，各方法分割精度如表 1所示。

由表 1可见，LIC算法的精度最低，LGDF精度略高于LIC，本文提出的MDAC精度明显高于前两者。

本文算法利用熵、均值和方差作为灰度描述子，并在模型中引入偏移场因子。与LGDF相比，既增加了描述子熵，用来提高灰度描述精准度，又补充了偏移场的全局信息，因此，性能优于LGDF。本文算法与LIC相比，增加了方差和熵两个灰度描述子，大幅提高了局部灰度的描述精准度，这对具有灰度不均带噪图像的分割是至关重要的，并且，熵不仅充实了灰度分布的描述，还具有较好的抗噪性能，因此，本文MDAC算法的分割精度比LIC的分割精度提高了30%多。

值得一提的是，在通常情况下LIC的性能是优于LGDF的，但针对本文所分割的低对比度医学图像，LGDF因引入方差作为描述子，导致分割精度高于LIC。

本文提出了一种针对医学图像分割的多描述子活动轮廓模型MDAC，解决了医学图像带有灰度不均和受噪声干扰的问题。图像的局部熵可以反映图像信息丰富程度，适合于表征图像分割区域的同质，可以提高灰度不均的医学图像的描述精准度。MDAC算法在均值和方差的基础上引入局部熵，提高了图像灰度分布描述的准确性，从而提高了分割精度。局部熵的取值由窗口内多个像素点共同决定的，单个像素的灰度值对熵的影响不大，因而对单个噪声点也不敏感，局部熵具有一定的滤波作用。局部熵的引入使得MDAC算法具有较强的抗噪性能。实验结果表明，算法可以同时实现图像的分割和图像纠偏，而且算法同经典的LIC和LGDF相比具有更高的分割精度和更强的抗噪性能。未来的工作将针对脑部医学影像的复杂性和组织结构的特点展开多相分割的研究，同时提高算法精度和抗噪性。