引言

近年来，随着人们生活水平的提高，前列腺癌的发病率和死亡率越来越高^[1]，是位居欧美国家男性发病率第一位、死亡率第二位的恶性肿瘤.2018年，中国国家癌症中心^[2]发布的我国肿瘤登记地区表明我国前列腺发病率也有明显的增长趋势，上海、北京、广州发达城市的发病率已经达到了0.01%.临床经验表明，前列腺癌如果能尽早发现并及时治疗，会大大提高患者的生存率.医学影像在前列腺癌诊断和治疗中发挥着重要作用.其中磁共振图像对组织分辨率高，具有多方向、多参数成像方式，不仅能显示人体解剖结构及其病理变化的信息，而且还能提供有关器官的功能信息.而图像分割^[3]是基于磁共振图像的前列腺癌早期诊断和治疗的基础^{[4, 5]}.

几何活动轮廓模型是图像分割领域的重要研究方向之一^[6-13].2010年，Langerak等^[14]利用标记融合方法进行了前列腺磁共振图像分割研究.该方法首先将目标图像与训练图像库中的多幅图像进行配准，然后将配准后的训练图像进行手工分割，配准的质量直接影响分割结果.2010年，Gao等^[15]利用形状配准方法进行了前列腺磁共振图像分割研究.该方法基于面配准，而点相关信息的缺乏导致基于面配准算法必须建立在迭代搜索的基础上，占用机时长.2010年，Martin等^[16]利用概率统计方法进行了前列腺磁共振图像分割研究.该方法在多个统计形状模型中选择最优模型，最优模型的选取直接影响图像的分割结果.2013年，Khalvati等^[17]对双向配准方法进行了研究.该方法需要选取特征点，但特征点的选取较为困难，且直接影响图像分割结果.综上所述，由于前列腺磁共振图像的周围组织结构复杂，图像比较模糊而且噪声比较大，其分割一直都是一个挑战.

本文基于改进边缘的距离正则化水平集演化（Distance Regularized Level Set Evolution，DRLSE）^[18]模型，对前列腺磁共振图像分阶段分割法进行了研究.相比其他水平集模型（如CV^[19]、LBF^[20]等），该方法提高了分割精度，而且操作简单，可以直接对现有图像进行处理.本文还并采用Dice相似性系数（Dice similarity coefficient，DSC）和Hausdorff距离对其分割效果进行评估.

1 改进的DRLSE模型

本文提出的改进DRLSE模型是在传统DRLSE方法^[21]的基础上使用两个指示函数，另外添加局部灰度信息的方法.

1.1 传统的DRLSE方法

初始化函数-取点连线法在前列腺磁共振图像内部任意取N个点，使其首尾连接成为一个封闭的区域，设为${S_N}$.该方法不需要设置坐标范围，可直接在图像上选点来确定初始轮廓.理论上，初始化轮廓越接近期望，形状分割精度越高，初始化水平集函数如(1)式所示：

$ {\phi }_{0}(x, y)\left\{ {\begin{array}{*{20}{l}}-{c}_{0}, \ \ \ \ \ \ {\phi }_{0}(x, y)\in {S}_{N}\\ {c}_{0}, \ \ \ \ \ \ \ \ \ \ 其它\end{array}} \right.$

(1)

为了解决水平集函数在演化过程中需要不断初始化为符号距离函数的问题，Li提出了DRLSE方法^[21]，在传统的活动轮廓模型中的能量泛函上再添加一项内部能量泛函，该方法可以来纠正水平集函数在演化过程中与距离函数的偏差，从而实现水平集函数在演化过程中无需周期性地重新初始化水平集函数的要求，其能量泛函如(2)式所示：

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

(2)

其中${R_p}(\phi )$是水平集的距离调整项，${\varepsilon _{ext}}(\phi )$是轮廓的驱动能量项.驱动水平集函数曲线向前列腺的边界方向运动.驱动能量项有很多形式，但根据前列腺图像的特点，其应满足条件：当水平集曲线达到前列腺边界时，能量达到最小值.

内部能量泛函${R_p}(\phi )$表示为(3)式：

${R_p}(\phi ) = \mu \int_\mathit{\Omega} {p(|\nabla \phi |){\rm{d}}x} $

(3)

(3) 式为函数ϕ的内部能量泛函（Ω是图像的区域），用于消除水平集函数和符号距离函数之间的偏差，使其始终保持符号距离函数的属性；$\mu $为大于0的常数，是权重参数；∇为拉普拉斯算子；p为能量密度函数，$p:[0, \infty ) \to R$.为了避免边界效应，能量密度函数如(4)式所示：

$p(s)\left\{ {\begin{array}{*{20}{l}} \frac{1}{{{{(2\pi )}^2}}}[1 - \cos (2\pi s)], \ \ \ s \leqslant 1 \\ \frac{{{{(s - 1)}^2}}}{2}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ s > 1 \\ \end{array}} \right.$

(4)

${\varepsilon _{ext}}(\phi )$为外部能量函数，仅仅依赖于图像特征，用于驱使水平集函数向期望位置演化，并且达到期望位置时能量最小，其能量方程如(5)式所示：

${\varepsilon _{ext}}(\phi ) = \lambda {L_g}(\phi ) + \alpha {A_g}(\phi ) = \lambda \int_\mathit{\Omega} {g{\delta _\varepsilon }(\phi )(|\nabla \phi |){\rm{d}}x + \alpha } \int_\mathit{\Omega} {g{H_\varepsilon }( - \phi ){\rm{d}}x} $

(5)

其中ϕ为水平集函数.${L_g}(\phi )$和${A_g}(\phi )$分别表示水平集的长度和水平集所围成的面积.λ、α均为权重系数.g为边界指示函数.${\delta _\varepsilon }(\phi )$和${H_\varepsilon }( - \phi )$分别为一维的Dirac函数和Heaviside函数.Heaviside函数如(6)式所示：

$ {H_\varepsilon } = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{2}\left[ {1 + \frac{x}{\varepsilon } + \frac{1}{{\pi}}\cos \left( {\frac{{{\pi}x}}{\varepsilon }} \right)} \right], }&{\left| x \right| \leqslant \varepsilon } \\ {0, }&{x > \varepsilon } \\ {0, }&{x < - \varepsilon } \end{array}} \right. $

(6)

$\varepsilon$表示函数宽度，Dirac函数为${H_\varepsilon }( - \phi )$的导数，如(7)式所示：

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

(7)

1.2 改进的DRLSE方法 1.2.1 改进边界指示函数

g为边界指示函数.本文首先观察了在传统DRLSE下使用两个边界函数进行前列腺图像分割的情况.前一阶段使用${g_1}(I)$，将曲线演化到前列腺靠近尿道侧及左右两侧的边界，后一阶段使用${g_2}(I)$，将曲线演化到前列腺靠近膀胱一侧的边界.由于前列腺靠近尿道以及左右两侧的边界比较清晰、靠近膀胱一侧比较模糊，因此本研究在不同阶段采用不同的指示函数，以达到最优的边界分割结果.${g_1}(I)$如(8)式所示：

$ {g_1}(x, y) = \frac{1}{{1 + |\nabla {G_{{\sigma _1}}}I(x, y){|^2}}} $

(8)

(8)式中$\nabla $表示空间梯度算子；$I(x, y)$为图像的灰度函数，${G_\sigma }$表示标准方差为σ的二维高斯函数滤波器，作用是消除噪声和平滑图像.当边缘指示函数到期望边界区域时，图像梯度的值$|\nabla {G_{{\sigma _1}}}I(x, y)|$趋大，而$g(x, y)$较小，且趋近0，迫使整个外部能量函数趋于0，使外部能量函数不再演化；相反，远离图像中的期望边缘时，图像的梯度值$|\nabla {G_{{\sigma _1}}}I(x, y)|$相对较小，因而$g(x, y)$的值趋近1，曲线还会不断进行演化，由此可以根据函数$g(x, y)$的值是趋近0还是趋近1来判断对应像素$I(x, y)$是否位于边界.

${g_2}(I)$表示如(9)式所示：

${g_2}(x, y) = \rho {{\rm{e}}^{ - \frac{{|\nabla {G_{{\sigma _2}}}I(x, y){|^2}}}{{2{\sigma _2}^2}}}}$

(9)

ρ表示权重参数.

对于磁共振图像，特别是对前列腺磁共振图像来说，由于前列腺靠近膀胱一侧的边界比较模糊，对比度比较低，传统DRLSE模型中的停止函数无法准确定位弱边界^[22]，容易导致过分割.改进的径向基高斯核函数为边缘停止函数^[23]，该函数满足单调性，自变量趋向于∞时，函数值趋近0；自变量趋向于0时，函数值趋近于ρ，使曲线在前列腺弱边界处停止演化.通过最小化泛函能量函数$\varepsilon (\phi )$，将水平集函数[(2)式]转化为偏微分方程，如(10)式所示：

$ \frac{{\partial \phi }}{{\partial t}} = \mu div({d_p}|\nabla \phi |\nabla \phi ) + \lambda {\delta _\varepsilon }(\phi )div(g\frac{{|\nabla \phi |}}{{\nabla \phi }}) + \alpha g{\delta _\varepsilon }(\phi ) $

(10)

其中∇为拉普拉斯算子；div为散度.

图 1为传统DRLSE下使用不同指示函数的演化结果.由图像分割结果可见，图 1(a)前列腺靠近膀胱一侧梯度不明显，仅使用${g_1}(I)$函数演化出现了过分割现象.图 1(b)仅使用${g_2}(I)$函数演化，曲线停止在弱边界处并且不会导致过分割，但达不到期望的边界处，从而导致欠分割.图 1(c)为结合${g_1}(I)$和${g_2}(I)$分阶段作用演化的结果，曲线既能在强边界停止又能在弱边界停止.因为前列腺图像内部中央腺和外周带存在虚假边界，使用该算法结果未能越过虚假边界，可以证明边界函数对图像分割结果具有重要的影响.尽管水平集方法高效简单，但由于前列腺图像部位特殊、噪声较大、存在虚假边界及靠近膀胱侧模糊等问题，基于传统DRLSE方法的前列腺分割效果很不理想.

1.2.2 增加能量牵制项

鉴于传统DRLSE方法在前列腺轮廓提取方面的局限，本文从局部灰度信息出发对其进行改进.即在传统DRLSE的外部能量泛函中添加新的能量项，本文称之为牵制能量项${C_g}$.它能反映图像演化曲线上各点的位置相对于真实边界的临近程度，可有效驱动零水平集曲线向前列腺真实边界方向运动，避免轮廓演化停止在背景区域的局部梯度极值处.此时，水平集函数ϕ的外部能量函数泛函为(11)式所示：

${E_{ext}}(\phi ) = \lambda {L_g}(\phi ) + \alpha {A_g}(\phi ) + k{C_g}(\phi )$

(11)

其中，k为牵制能量项的权重系数.牵制能量项为(12)式所示：

${C_g} = \int_l {XH} ( - \phi ){\rm{d}}x$

(12)

$X = 1/\sum\nolimits_1^n | {f_j} - {f_i}|$

(13)

其中，n为水平集曲线上选取的点的个数，以选取的点M作曲线的法线，在M点的法线上选取位于曲线两侧的点i和j（图 2），${f_j}$为曲线弧外侧点的像素，${f_i}$为曲线弧内侧点的像素.l为图 2曲线的一周，H为单位阶跃函数.此时改进的能量泛函如(14)式所示：

${E_{ext}}(\phi ) = k{C_g}(\phi )$

(14)

外部能量仅使用牵制能量项泛函进行演化的结果如图 3所示.结果显示，使用能量牵制项可以使曲线向边界靠近，特别是向梯度变化比较大的地方演化，而且可加速曲线的演化.

本文改进的DRLSE方法保留了传统DRLSE的长度项和面积项，为了得到极小值，将(2)式转化为偏微分方程，如(15)式所示：

$\frac{{\partial \phi }}{{\partial t}} = \mu div({d_p}|\nabla \phi |\nabla \phi ) + \lambda {\delta _\varepsilon }(\phi )div(g\frac{{|\nabla \phi |}}{{\nabla \phi }}) + \alpha g{\delta _\varepsilon }(\phi ) + k{C_g}{\delta _\varepsilon }(\phi )$

(15)

由上式可得，改进DRLSE方法不仅利用了图像边界的梯度信息，而且还利用了图像的局部灰度信息，能够拖动零水平集向边界处运动.当水平集的能量泛函最小时，牵制能量项达到最小值，曲线恰好达到边界值，此时，零水平集即为前列腺的真实轮廓.

用改进DRLSE方法对前列腺磁共振图像进行前列腺轮廓的提取，演化结果如图 4所示.

实验结果表明，改进DRLSE方法对前列腺轮廓的提取具有较好的效果，有效地解决了传统的DRLSE方法对前列腺轮廓提取所存在的问题，添加了牵制能量项来驱使曲线越过虚假边界继续演化，同时根据外部能量泛函中的长度和面积最小原则，使零水平集接近梯度极值，最终到达目标区域梯度极值处.

2 材料与方法

本文选用国际光学和光子学会、美国物理学家协会和国家癌症研究所联合发起的医疗影像识别竞赛PROSTATEx数据库中190幅前列腺磁共振图像进行实验.

3 实验结果与分析 3.1 传统DRLSE方法和改进DRLSE方法进行前列腺磁共振图像分割的对比

本文选取六个点进行初始轮廓设定，均为多边形（图 5）.

前列腺磁共振图像分割实验在Matlab2016b下运行，其中主要参数的取值：$\mathit{\Delta} t = 1$、$\mu = 0.2$、$\lambda = 5$、$\alpha = - 3$、${\sigma _1} = 3$、${\sigma _2} = 1.5$、$k = 2$.图 6为基于传统的DRLSE方法对前列腺采用两个指示函数分阶段演化的结果，可见演化容易停止在前列腺的虚假边界上.图 7为采用改进的DRLSE方法对前列腺边界进行演化的结果，在不同阶段采用不同的指示函数演化150次，本文先使用${g_1}(I)$演化100次到达边界清晰处，再通过${g_2}(I)$演化50次达到弱边界处.其中白色线为专家手动分割结果，红色线为采用算法分割得到的结果.结果表明使用本文提出的改进DRLSE算法得到的结果和专家手动分割结果更加接近，证明了改进方法的有效性.

3.2 Hausdorff距离和DSC分析

文献^[24]提出将Hausdorff距离用于评估两条曲线之间的误差.本文将使用算法分割区域定义为像素集A，将专家手动分割的标准区域定义为像素集B.Hausdorff距离定义如(16)式和(17)式所示：

$H(A, B) = \max (h(A, B), h(B, A))$

(16)

$h(A, B) = {\max _{a \in A}}{\min _{b \in B}}||a - b||$

(17)

$||a - b||$表示专家手动分割曲线和使用算法自动分割曲线上相对应点之间的距离.通过(17)式求范数，通过(16)式得到最终的距离.使用本文提出的改进DRLSE方法对190幅前列腺磁共振图像进行分割的平均Hausdorff距离为(2.58±1.22) mm，远远小于使用传统DRLSE方法得到的平均Hausdorff距离[(9.39±5.03) mm].

文献^[25]提出使用DSC对使用的分割方法的精确性进行定量评估.DSC定义如(18)式所示：

$DSC({A_s}, {A_m}) = 2\frac{{|{A_s} \cap {A_m}|}}{{|{A_s}| + |{A_m}|}} \times 100\% $

(18)

其中，${A_s}$和${A_m}$分别是使用算法自动分割和专家手动分割得到的区域像素.DSC越接近100%，表明结果越好.针对图 6和图 7显示的前列腺病例的第9~13层扫描层，利用传统的和改进的DRLSE方法对前列腺磁共振图像不同扫描层轮廓提取的DSC和Hausdorff距离如表 1所示，使用改进DRLSE方法的DSC均值达到96%，远高于传统的DRLSE方法.

图 8为使用改进的DRLSE方法获得的190幅前列腺磁共振图像的DSC分布直方图，结果表明DSC分布集中在95%~97%区间，DSC均值为96%，可靠性较高.

表 2列举了文献^[14-17]中使用其他方法进行前列腺磁共振图像分割获得的DSC，除双向匹配方法外，DSC值均在85%左右，使用双向匹配方法获得的DSC为92%，均小于本文提出的方法.这表明本文提出的改进DRLSE方法优于其它方法，对前列腺磁共振图像的分割效果较为理想、精度较高.

4 总结

本文提出了一种基于改进DRLSE的前列腺磁共振图像分割方法.该方法除保留了传统DRLSE的梯度信息外，还添加了牵制能量项，并且使用了两种指示函数分阶段进行演化.由于磁共振图像比较模糊且噪声比较大，本方法在轮廓提取的过程中引入了高斯核函数，通过卷积平滑前列腺图像减少噪声，使得边界约束函数更好地促使分割算法收敛于边界，而且能越过虚假边界最终到达到真正的边界.实验结果表明，使用改进DRLSE方法对前列腺磁共振图像的分割精度明显高于传统DRLSE和文献中使用的其它方法，DSC均值可达到96%，从而证明了改进DRLSE方法对于前列腺磁共振图像分割的有效性和精确性，对前列腺疾病的临床诊断和治疗都具有非常重要的意义.

本文方法仅基于二维磁共振图像研究前列腺的分割，由于三维图像存在计算复杂等问题，对三维图像还在探索中，未来将在三维图像上展开研究.