2020年，全球男性前列腺癌新发病例数为141万，位居男性癌症发病率第4；中国男性前列腺癌新发病例数为12万，位居中国男性癌症发病率第6^[1].前列腺癌对男性健康构成严重威胁.临床上常使用穿刺和指诊被来检查前列腺^[2]，但这些方法会给患者带来一定的创伤和不适.随着计算机辅助诊断技术的迅速发展，多序列磁共振成像（magnetic resonance imaging，MRI）在前列腺疾病的诊断中发挥着越来越重要的作用，而感兴趣区域的分割^[3]是计算机辅助诊断的重要组成部分^{[4, 5]}.

目前已有的有效的分割模型包括传统的机器学习模型^[6-9]、图割模型^{[10, 11]}和活动轮廓模型^[12-17].深度学习也被应用于前列腺区域分割领域^{[18, 19]}.Vafaie等^[7]利用马尔可夫随机场构建概率评估模型，对前列腺磁共振图像中的像素进行分类，并估计像素所属的区域，但该模型计算量大、分割效果不稳定.目前，在图像分割中被广泛应用的活动轮廓模型可分为两类：基于边缘的模型^[12-14]和基于区域的模型^[15-17].Li等^[12]提出了一种基于边缘的水平集分割模型，该模型加入了距离正则化，解决了水平集演化过程中的周期初始化问题，加快了演化速度，然而当边界较弱时，效果并不理想.Li等^[15]利用图像的局部尺度信息构造能量函数，但该模型在分割灰度不均匀的图像时，效果不理想.前列腺的中央腺靠近膀胱处存在一定的弱边界，图像分割时往往会忽略前列腺中央腺的顶部.朱泽华等^[14]在基于改进的距离正则化水平集（distance regularization level set evolution，DRLSE）分割模型中，提出了两个指示函数函数，分别分割前列腺强边界和弱边界.Milletari等^[19]提出了一种基于体素的全卷积三维分割V-net神经网络模型，虽然在分割结果上取得了一些进展，但在边界较弱的情况下，分割效果较差.

在前文所列举的这些活动轮廓模型中，文献[12]是完全基于图像梯度信息进行分割，文献[15]是基于图像的灰度信息进行分割.但是前列腺磁共振图像比较复杂，单纯以一种信息去构建能量函数，难以精确分割前列腺区域.因此，本文提出了一种将图像的灰度信息和梯度信息进行融合的前列腺区域分割模型，其原理是在DRLSE的基础上，进行原始边缘指示函数的优化，并进行局部灰度信息的计算，通过提取前列腺的粗轮廓作为初始进化轮廓，并融合边缘信息和局部灰度信息构造新的能量函数来获得前列腺的精细轮廓.

本文所提的前列腺区域分割模型包括4个步骤：1）读取T₂加权磁共振图像；2）采用半径为5个像素的8邻域模板，计算局部二值模式（local binary pattern，LBP）特征图^[20]；3）利用改进的DRLSE模型对LBP特征图进行粗分割，得到前列腺粗轮廓；4）以粗轮廓为初始状态，采用新的水平集能量函数进行迭代，得到前列腺的精细轮廓.

LBP是由Ojala等^[20]首次提出的一种计算局部图像模式的算子，用来描述图像局部纹理特征.基本LBP算法利用邻域模板，以中心像素的灰度值为阈值，将邻域中其他像素的灰度值与中心点进行比较：如果邻域中像素的灰度值大于中心点，则该像素的位置标记为“1”；否则，标记为“0”.在运算结束时，得到一个二进制序列，然后以十进制数为中心点.

但每个相邻点与正方形邻域中心点之间的距离不一致，无法很好地表达中心点与相邻点之间的灰度关系.为了克服这个缺点，本文使用了一个圆形邻域LBP：

(1)${\rm{LBP}}_{P, R}^{} = \sum\limits_{n = 0}^{P - 1} {s({I_n} - {I_c})} {2^P}$

(2)$s(x) = \left\{ {\begin{array}{*{20}{c}} {1, {\rm{ }}x \geqslant 0} \\ {0, {\rm{ }}x < 0} \end{array}} \right.$

其中s为符号函数，R表示圆形邻域的半径，P表示邻域像素的数目，${I_n}$表示中心像素灰度值，${I_c}$表示邻域中其他点的灰度值.本文使用半径为5个像素的圆形邻域模板来构建LBP特征图.

DRLSE是由Li等^[12]在2010年提出的改进的活动轮廓模型，是一种基于边缘信息的水平集分割方法.它通过引入一个新的距离正则化项，在传统活动轮廓模型的能量函数中加入一个内能函数来修正水平集与符号距离函数之间的偏差，从而避免水平集的重新初始化.模型的能量函数如(3)式所示：

(3)$\varepsilon (\phi ) = \mu P(\phi ) + {\varepsilon _{ext}}(\phi )$

其中，$\varepsilon (\phi )$为能量函数；$\mu > 0$为常数；$P(\phi )$为水平集算法中的距离正则化项，为了修正水平集与符号距离函数之间的偏差，避免演化过程中曲线的重新初始化；${\varepsilon _{ext}}(\phi )$是一个外部能量函数，它驱动演化曲线向目标方向移动，用如下公式表示：

(4)${\varepsilon _{ext}}{\rm{(}}\phi {\rm{) = }}\lambda {L_g} + \alpha {A_g}$

(5)${L_g}(\phi ) \triangleq \int_\Omega {g\delta (\phi )} {\rm{ }}\left| {\nabla \phi } \right|{\rm{d}}x$

(6)${A_g}(\phi ) \triangleq \int_\Omega {gH( - \phi )} {\rm{ d}}x$

(7)$g \triangleq \frac{1}{{1 + |\nabla {G_\sigma }^*I{|^{\rm{2}}}}}$

其中，$\phi $为能量泛函，用于消除水平集函数和符号距离函数之间的偏差，使其始终保持符号距离函数的属性；$\nabla $为拉普拉斯算子；Ω表示图像域；${L_g}{\rm{(}}\phi {\rm{)}}$表示长度惩罚项；g为边缘指示函数，它可以防止曲线演化超出目标区域，使目标边界处的能量最小化；${A_g}(\phi )$是区域驱动项，可以用来加速演化；λ、α为常数项；$\nabla {G_\sigma }$表示高斯核算子；$I$为定义域Ω上的一幅图像；x表示图像像素点；δ和H分别为Dirac^[21]和Heaviside^[22]函数，如(8)式和(9)式所示：

(8)${\delta _\varepsilon }{\rm{(}}x{\rm{) = }}\left\{ {\begin{array}{*{20}{c}} {\frac{1}{{2\varepsilon }}\left[ {{\rm{1 - }}\sin (\frac{{{\rm{ \mathsf{ π} }}x}}{\varepsilon })} \right], }&{|x| \leqslant \varepsilon } \\ {0, }&{|x| > \varepsilon } \end{array}} \right.$

(9)${H_\varepsilon }{\rm{(}}x{\rm{) = }}\left\{ {\begin{array}{*{20}{c}} {\frac{1}{2}\left[ {{\rm{1 + }}\frac{x}{\varepsilon }{\rm{ + }}\frac{1}{{\rm{ \mathsf{ π} }}}{\rm{cos(}}\frac{{{\rm{ \mathsf{ π} }}x}}{\varepsilon }{\rm{)}}} \right], }&{\left| x \right| \leqslant \varepsilon } \\ {0, }&{\left| x \right| > \varepsilon } \end{array}} \right.$

其中，ε表示函数宽度，常取1.5.

在原始的DRLSE模型中，g函数在进化过程中用于外部约束.对于灰度分布不均匀的图像，边缘指示函数收敛速度慢，在弱边界容易产生误差.因此，本文提出了一种改进的基于图像阈值的边缘指示函数，并讨论了其减小误差的能力.

新的边缘指示函数g_new利用了Otsu^[23]提出的自适应阈值原理计算阈值：记录T作为目标区域和背景区域的分割阈值；将目标区域中的像素点与图像的比例设置为${w_0}$，将平均灰度设置为${u_0}$，将背景像素数的比例设置为${w_1}$，将平均灰度设置为${u_1}$，则图像的综合平均灰度为$u = {w_0} \times {u_0} + {w_1} \times {u_1}$；当T最大化类间方差值${v^2} = {w_0} \times {({u_0} - u)^2} + {w_1} \times {({u_1} - u)^2}$时，T为最优阈值，简化的类间方差可以用${v^2} = {w_0} \times {w_1} \times {({u_0} - {u_1})^2}$来描述.

(10)${g_{new}} \triangleq \frac{1}{{1 + K^*|\nabla {G_\sigma }^*I{|^{\rm{2}}}}}$

(11)$K = {v^2}$

图 1表示图像梯度分布情况，白色部分数值为“1”，黑色部分数值为“0”.当能量函数在演化过程中碰到黑色时，说明已经到达图像的边缘，即停止演化；碰到白色时，说明还没有到达图像的边缘，需要继续演化.由图 1可以看出，原g函数较易陷入局部极小值，导致曲线的演化速度和进程受到限制，从而使曲线不能收敛到期望的目标边界.利用改进后g_new函数可以保持前列腺区域的边缘，平滑图像的灰度信息，有助于加速曲线的演化.

采用基于g_new函数的改进的DRLSE模型对圆形邻域LBP特征图进行分割，提取了前列腺粗轮廓，如图 2所示，然而，提取的前列腺粗轮廓与真实前列腺轮廓仍有一定差距，有待进一步完善，需要进行细分割.

将粗分割的结果作为精分割的初始水平集轮廓，并利用所提出的能量函数对前列腺区域进行分割.在Li等^[24]的LBF模型中，引入高斯计算器，利用局部灰度信息构造活动轮廓模型，对灰度分布不均匀的前列腺磁共振图像具有一定的应用价值.其能量函数可用(12)式表示：

(12)$F(\phi , {f_1}, {f_2}) = {E_\varepsilon }(\phi , {f_1}, {f_2}) + \mu {R_P}(\phi )$

其中，${R_P}(\phi ) = \frac{{\rm{1}}}{{\rm{2}}}{\int_\Omega {(|\nabla \phi | - 1)} ^2}{\rm{d}}x$, ${R_P}$为距离惩罚项，纠正水平集函数与符号距离函数之间的偏差；μ为正常数；${f_1}(x)$、${f_2}(x)$表示邻域的平均灰度，如(13)式和(14)式所示，

(13)${f_1}(x){\rm{ = }}\frac{{{K_\sigma }^*[{H_\varepsilon }(\phi (x))I(x)]}}{{{K_\sigma }^*{H_\varepsilon }(\phi (x))}}$

(14)${f_{\rm{2}}}(x){\rm{ = }}\frac{{{K_\sigma }^*[({\rm{1 - }}{H_\varepsilon }(\phi (x)))I(x)]}}{{{K_\sigma }^*(1 - {H_\varepsilon }(\phi (x)))}}$

其中${K_\sigma }$是高斯核函数，$H(\phi )$定义与(9)式相同.

DRLSE模型要求图像具有很强的梯度变化.对于边界较弱的灰度图像，分割不够理想^[14].LBF模型将图像的局部灰度信息作为曲线的演化条件，但不计算梯度信息.如果不能在目标边界及时停止演化，往往会导致过度分割^[24].本研究提出一种新的活动轮廓模型来探索解决上述问题的可能性.

2000年，Xu等^[25]提出了一种活动轮廓算法即符号压力（signed pressure force，SPF）的概念，通过统计区域信息改变函数的符号来驱动曲线演化.传统的SPF函数对对比度明显的图像效果较好，但对边界较弱或模糊的图像效果不理想.Zhang等^[26]提出了选择性二值滤波和高斯滤波正则化水平集（selective binary and Gaussian filtering regularized level set，SBGFRLS）模型，利用Chan-Vese（CV）^[27]模型提出的能量函数对SPF函数进行了改进，但是，该模型计算的是图像的全局灰度，由于缺少惩罚约束，对T₂磁共振图像中前列腺区域的分割效果并不理想.

本文利用LBF提出的局部灰度信息，构造了一个新的SPF函数来解决上述问题，其取值范围为$[ - 1, 1]$，具体表示为：

(15)${\rm{Local\_}}spf(I(x)) = \frac{{I(x) - \frac{{{f_1} + {f_2}}}{2}}}{{\max (|I(x) - \frac{{{f_1} + {f_2}}}{2}|)}}, x \in \Omega $

(15) 式中的${f_i}$函数与(13)式和(14)式相同，$I(x)$表示图像中像素点的灰度值.${\rm{Local\_}}spf(I(x))$利用局部区域的灰度信息来驱动曲线的演化，将其应用于灰度不均匀的图像可以得到更好的结果.

前列腺磁共振图像很复杂.由于器官和组织之间的灰度差很小，梯度不明显，而且进化过程中水平集可能超出目标区域的边界.针对上述问题，本研究将原始水平集的能量函数进行优化，构建了一个新的能量函数来约束曲线的演化，可以写成：

(16)$E(\phi ) = {E_e}(\phi ) + \mu {E_{DR}}(\phi ) + \lambda {E_{LR}}(\phi )$

在(16)式中，${E_e}(\phi )$表示外部能量项，${E_{DR}}(\phi )$表示水平集的距离正则化项，${E_{LR}}(\phi )$表示基于SPF的局部区域灰度信息.新的能量函数可以减少重新初始化的次数，不仅可以加快曲线向目标边界演化的速度，而且可以修正水平集能量函数与有符号距离函数之间的偏差.新的能量函数可表达为：

(17)$\begin{gathered} E(\phi ) = v\int_\Omega {{g_{new}}\delta (\phi )|\nabla \phi |{\rm{d}}x} + \alpha \int_\Omega {{g_{new}}{H_\varepsilon }( - \phi )} \, {\rm{d}}x \\ {\rm{ }} + \mu \int_\Omega {\frac{{\rm{1}}}{{\rm{2}}}} {(|\nabla \phi | - 1)^2}{\rm{d}}x + \lambda \int_\Omega {Local\_spf} (I){H_\varepsilon }( - \phi ){\rm{d}}x \\ \end{gathered} $

其中，ν、α、μ、λ为常数项，均为函数的权重系数。将新的边缘指示函数[(10)式]和新的SPF[(15)式]代入(17)式，div表示散度算子，用梯度下降流最小化的公式可以表示为：

(18)$\begin{gathered} \frac{{\partial \phi }}{{\partial t}} = \nu {\delta _\varepsilon }(\phi )div\left( {{g_{new}}\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right) + \alpha {g_{new}}\delta (\phi ) \\ {\rm{ }} + \mu \left[ {\Delta \phi - div\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right] + \lambda Local\_spf(I){\delta _\varepsilon }(\phi ) \\ \end{gathered} $

这项研究使用了203组来自SPIE-AAPM-NCI（国际光学与光子学学会-美国医学物理学家协会-国家癌症研究所）前列腺MR分类挑战赛数据库的T2W序列图像集，每组图像分别来自一位不同的患者.

本实验是在WIN10x64位操作系统，Matlab R2018a下进行的。经过优化，使用的参数为：$\Delta t = 1$（表示以像素为单位进行计算，一次计算一个像素），$v = 5$，$\alpha = - 3$，$\mu = 0.2$，$\lambda = - 0.5$，迭代次数为350次.该模型对参数v和μ不敏感；在203例图像分割实验中，λ数值取-0.5时可以得到较好的分割效果.

本文将医生手工分割结果作为金标准，分别采用95% Hausdorff距离（95% HD）、Dice系数（Dice similarity coefficient）、相对体积差（relative volume difference，RVD）作为指标，对所提模型分割精度进行评价.评价指标均采用文献[28]所提出的公式.

Dice系数用(19)式表示：

(19)$Dice = \frac{{2|X \cap Y|}}{{|X| + |Y|}}$

在(19)式中，$|X|$和$|Y|$分别代表由分割模型和医生手动分割得到的前列腺区域的面积，由像素数表示.其分布范围为$[0, 1]$：0表示两个数据集之间没有交集，说明与医生手动分割结果相差大；1表示两个数据集相同，说明与医生手动分割结果更一致.

RVD通过以下公式定义：

(20)$RVD = (\frac{{|X|}}{{|Y|}} - 1)^*100\% $

其中$|X|$、$|Y|$与(19)式相同.RVD用于评估过分割或欠分割：如果RVD为正，则表示欠分割；否则表示过分割.

HD用于评估两条曲线之间的误差.将医生手工分割及模型分割的前列腺轮廓分别定义为像素集A和像素集B，HD可表示为：

(21)$HD(A, B) = \max (h(A, B), h(B, A))$

(22)$h(A, B) = \mathop {\max }\limits_{x \in A} [\mathop {\min }\limits_{y \in B} d(a, b)]$

其中$d(a, b)$是两条曲线上点之间的欧氏距离.由于正态Hausdorff距离对异常值非常敏感，因此参考文献[28~30]，我们使用非对称Hausdorff距离的第95个百分位，即95% HD，而不是最大值.

图 3显示了两个患者的前列腺轮廓，分别使用本文所提模型和手动分割获得.

本文所提出模型的Dice系数平均值±标准差为0.94±0.01，表明该模型具有较高的精度.RVD平均值±标准差为-1.21%±2.44%，说明该模型有过分割的倾向，但是由于绝对值接近于0，说明总体分割误差不大，且鲁棒性较好.95% HD数值平均值±标准差为6.15±0.66 mm，该参数表示的是金标准轮廓与模型分割轮廓的误差.

文献[7]提出了一种基于马尔可夫随机场的分割模型，文献[12]提出了一种基于DRLSE模型的分割模型.另外，文献[15]中提出了基于LBF模型的分割模型，文献[19]中提出了基于V-net神经网络模型的分割模型.本研究将所提模型与四种现有模型的前列腺分割性能进行了比较.实验结果表明，本文提出的分割模型具有较强的鲁棒性，随机选择一位患者，不同模型的分割结果如图 4所示.

基于文献[19]模型的分割结果如图 4(c)所示，表明前列腺周围区存在过分割现象，这可能是由于V-net在训练中使用的损失函数对分割没有较好的限制，导致过分割的产生.在文献[15]的基础上提出的分割模型，在图像灰度分布均匀的情况下，可以得到较好的分割效果.但对于结构复杂、灰度分布不均匀的前列腺图像，分割效果并不理想，如图 4(d)所示，当目标边界较弱，梯度变化不明显时，分割结果中往往存在过分割现象.前列腺图像的梯度信息是用文献[12]中的分割模型计算得到，但是前列腺图像的灰度分布是不均匀的，总是存在弱边界.如图 4(e)所示，模型的边缘指示函数存在许多局部极小值，限制了曲线的演化，导致无法达到目标的真实边界.而本文所提模型通过计算图像的局部信息，较好地分割出前列腺，使曲线更接近真实前列腺边界.文献[7]分割模型，利用概率模型对图像中的区域像素进行建模，并利用概率模型对像素进行分类和评价，判断它们是否属于同一区域.由于前列腺图像存在较大的灰度和结构差异，使用该模型很难获得较好的分割结果，如图 4(f)所示.不同模型的分割精度如图 5所示.

图 5表明，本文所提模型具有较好的前列腺分割性能.与现有的其他模型相比，该模型与金标准之间的距离差更小、相似度更高，用于前列腺区域分割时，具有较高的精度和鲁棒性.

目前，水平集算法是一种常用的医学图像处理方法.参考文献[11]中的分割模型考虑了先验形状信息.在基于形状信息的分割模型中，形状参数需要通过训练数据获得.然而，在真实的前列腺图像中，图像的形状、大小和纹理往往存在较大个体差异，这给应用带来了一定的限制。

本文提出了一种新的前列腺分割模型，可用于分割不同大小和形状的前列腺.在粗分割阶段，采用一种新的边缘指示函数，以图像的阈值为调整因子，改变图像的梯度分布，减小最小梯度，同时保留前列腺的边界，使曲线快速向目标边界演化.本文提出的新分割模型用于粗分割结果的细化.通过计算局部灰度信息，加入了外部能量驱动项和距离正则化项，克服了水平集周期初始化的缺点.利用高斯核函数对图像进行平滑处理，以减少噪声的影响.基于改进了的边缘指示函数，图像梯度分布均匀，较好地保留了前列腺的边缘.实验结果表明，当前列腺与背景区域的对比度较小时，该模型具有较好的性能.使用我们的新模型均可以获得更好的分割结果，是因为该模型结合了前列腺的梯度信息和局部灰度信息进行曲线演化，并将粗分割的结果作为精分割的初始轮廓，避免了零水平集位置对水平集分割的影响.利用该模型得到的分割结果更接近前列腺真实轮廓.

将本文所提模型的分割结果与金标准进行了比较表明，所提模型可以对不同形状的前列腺进行分割，更重要的是，其分割效果明显优于其他分割模型.因此，该模型具有很强的鲁棒性.

本文提出了一种基于局部灰度信息和梯度信息的主动轮廓模型来更准确地分割前列腺，对LBF模型中基于局部灰度信息的活动轮廓进行了改进，并提出了一种新的基于局部灰度信息的活动轮廓，引入SPF，并将外部能量项和距离正则化项融合应用于前列腺T₂加权磁共振图像序列的分割.首先读取T₂加权磁共振图像，然后在构造LBP特征图的同时，利用改进的DRLSE进行粗分割并快速提取粗轮廓.提取的粗分割轮廓作为细分割的初始轮廓，在一定程度上降低了训练消耗的代价.采用基于局部区域和梯度信息的新模型进行精细分割和调整迭代次数，得到了较好的前列腺分割效果.实验结果也表明该模型具有较高的性能.

由于患者的不同，前列腺的外观存在很大的差异，未来的工作将集中在将该模型采用更大的数据集进行测试.