一种基于多相互斥水平集的黏连细胞图像分割方法

摘要

本发明公开了一种基于多相互斥水平集的黏连细胞图像分割方法。通过对显微细胞图中的类圆细胞进行Hough变换圆检测，获取所有细胞的圆心，以每个细胞的圆心自动设置一个小圆作为初始曲线，通过图像梯度信息的引导，采用多水平集函数同时演化的方式，在互斥能量的作用下对黏连细胞区域进行独立演化，且能保证不同封闭曲线间的互斥，从而达到分割黏连细胞图像的效果，解决传统分割方法难以分割黏连拥挤的细胞群区域的问题。本发明具有处理步骤少，所分割的细胞轮廓线自然，可用于类圆细胞图的分割和计数，测量精度高，适应能力强的优点。

说明书

技术领域

本发明属于医学图像处理技术领域，更具体地说，涉及一种基于多相互斥水平集的黏连 细胞图像分割方法。

背景技术

在计算机辅助分析乳腺组织切片的显微图像方法中，精确分割是重要的环节，它涉及细 胞计数和细胞形态分析的准确性。调查发现借助于自动检测设备进行乳腺组织切片显微图像 的病理分析，如果细胞之间拥挤或黏连，一般的图像分割技术常常会将两个黏连的细胞分割 到一个细胞区域中，从而导致检测区域中细胞的计数产生误差，因此，针对黏连细胞群，彼 此独立分割呈现，成为细胞显微图像分析的技术难点。

目前，常见的显微细胞图像分割技术主要是利用细胞图像特有的边缘、区域、颜色等特 征信息对细胞图像进行分割。基于边缘的分割方法，如Laplacian、Canny算子等，但是其分 割结果易受到噪声和纹理的影响，且分割的细胞边界常存在不连续性。基于区域信息的分水 岭算法是细胞分割邻域运用最为普遍的方法之一，学者们通过在分割目标的内部设置标记点， 而后基于这些标记点进行分水岭分割，从而能较好的抑制噪声及复杂背景的干扰，且能获得 封闭的细胞边界，然而此类方法对黏连细胞的分割效果仍不理想，且细胞轮廓外形不自然。

当前基于水平集理论的图像分割算法往往能获得较为自然的轮廓，其将可变形曲线作为 零水平集嵌入到高一维度的水平集函数中，通过图像统计信息的引导，使得零水平集不断向 物体轮廓演化，以达到分割物体的目的。其中经典的Chan-Vese(CV)模型是一种基于区域 信息的模型，其能够自动进行拓扑结构的变化，分割出图像的内部轮廓，但是其无法处理强 度分布不均匀的图像，且容易受到细胞复杂背景环境的干扰。李春明等提出了基于距离正则 化的水平集演化(DistanceRegularizedLevelSetEvolution，DRLSE)模型，其利用改进的惩 罚项和图像的梯度信息引导曲线演化，并且利用从分割目标的内部向外演化的方式，极大的 减少了背景的干扰，取得了较好的分割效果，但该模型依旧存在将两个同质黏连的细胞区域 演化成同一区域的现象。

经检索，关于黏连细胞分割方法的技术方案已有公开，但多存在处理步骤多，方法繁琐 的不足，且不能有效地解决细胞间拥挤和黏连的问题。如中国专利申请号201510253573.X， 申请日为2015年5月18日，发明创造名称为：基于多尺度生长与双策略去粘连模型的乳腺 细胞分割方法；该申请案先将乳腺组织图像转换为灰度图像；然后增强对比度；接着利用小 波分解进行细胞定位；然后进行多尺度区域生长；接着通过投票选择实现细胞区域的初步分 割；然后判断分割出的区域是否存在细胞粘连；如果不存在细胞粘连，则为单个细胞区域， 输出分割结果；如果存在细胞粘连，则为粘连细胞区域，进行粘连细胞分割；然后利用形态 学腐蚀-膨胀运算和角点检测分割算法构建的双策略去粘连模型进行粘连细胞分割，直至所 有细胞分割结束。该申请案能够抑制乳腺组织切片图像的复杂背景对乳腺细胞分割的影响， 在一定程度上提高对粘连细胞分割线的识别精度，达到提高粘连细胞分割精度的目的。但该 申请案中如果结构元半径过小，则不能将粘连区域分割开；反之，则会导致区域失真。此外， 该申请案提到的角点检测存在易受噪声干扰导致黏连区域判断失败等问题，且其处理步骤甚 多。所以，仍需一种适应性强、处理步骤简便的黏连细胞分割方法。

发明内容

1.发明要解决的技术问题

本发明的目的在于解决现有乳腺切片图像黏连区域难以分割的问题，为提高细胞图自动 分析结果的可靠性，提供了一种基于多相互斥水平集的黏连细胞图像分割方法；本发明主要 针对细胞具有类圆特性的乳腺切片细胞图，通过Hough圆检测确定细胞图中所有细胞的圆心， 并通过选定最小细胞内的初始圆，设置所有细胞的初始曲线，通过水平集能量泛函的演化， 使得轮廓向黏连细胞边界运动，通过设定互斥能量项确保相邻的轮廓不过界。

2.技术方案

为达到上述目的，本发明提供的技术方案为：

本发明的一种基于多相互斥水平集的黏连细胞图像分割方法，其步骤为：

步骤一、读取细胞图像，并对其进行二值化及阈值处理，获得代表细胞区域的近似圆形 斑块图；

步骤二、利用Hough圆检测方法对斑块图进行检测，获取图像中全部近似圆形斑块的类 圆圆心；

步骤三、以小于最小细胞的CHT检测圆半径作为所有细胞区域的初始圆半径，依据步骤 二所得每个细胞类圆圆心坐标设置多个水平集的圆形初始曲线；

步骤四、对多相互斥水平集能量泛函采用梯度下降法进行求解得到对应的偏微分方程；

步骤五、将步骤三所得初始曲线作为φ⁰代入到偏微分方程进行求解，从而引导图像中多 个初始曲线同时向目标边界演化；

步骤六、重复步骤五，继而得到细胞图像的最终边界轮廓。

更进一步地，步骤二中为滤掉半径过大或过小带来的干扰，必须通过设置合适的检测圆 的半径范围，获取所有类圆细胞斑块圆的圆心。

更进一步地，步骤三中为了避免最小细胞的CHT检测圆作为水平集初始轮廓在实际细胞 区域之外，必须通过以最小细胞的CHT检测圆半径的0.5～0.9倍作为所有细胞区域的初始圆 半径。

更进一步地，步骤四中多相互斥水平集模型将互斥能量项ε^R(φ_i)与DRLSE模型以多相水 平集的方式进行结合，对于水平集函数φ_i，所述多相互斥水平集能量泛函为：

$\left(\begin{array}{c}{E}_{MRLSE}={ϵ}_{ϵ}^{DRLSE}\left({\phi }_i\right)+\frac{\omega }{2}{ϵ}_{ϵ}^R\left({\phi }_i\right)\\ ={\mathrm{\mu R}}_p\left({\phi }_i\right)+{\mathrm{\lambda L}}_g\left({\phi }_i\right)+{\mathrm{\alpha A}}_g\left({\phi }_i\right)+\frac{\omega }{2}{ϵ}_{ϵ}^R\left({\phi }_i\right)\end{array}\right)$

式中，R_p(φ)为距离规则能量项，L_g(φ)为曲线长度项，A_g(φ)为面积项，ε^R(φ_i)为互斥能量 项，μ、λ、α、ω均为加权系数；

其中：

$R_p\left(\phi \right)=\underset{\Omega }{\int }p\left(\right|▿\phi \left|\right)dxdy$

$L_g\left(\phi \right)=\underset{\Omega }{\int }{\mathrm{g\delta }}_{ϵ}\left(\phi \right)|▿\phi |dxdy$

$A_g\left(\phi \right)=\underset{\Omega }{\int }{\mathrm{gH}}_{ϵ}(-\phi )dxdy$

${ϵ}_{ϵ}^R\left({\phi }_i\right)=\underset{\Omega }{\int }\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{(1-H_{ϵ}\left({\phi }_j\right))}^{2}dxdy$

上式中，Ω为图像全域，g为边界检测函数，H_ε为近似光滑函数，δ_ε为Hε的导函数； 多相互斥水平集能量泛函可进一步表示为：

$E_{MRLSE}=\mu \underset{\Omega }{\int }p\left(\left|▿{\phi }_i\right|\right)dxdy+\lambda \underset{\Omega }{\int }{\mathrm{g\delta }}_{ϵ}\left({\phi }_i\right)dxdy+\alpha \underset{\Omega }{\int }{\mathrm{gH}}_{ϵ}\left(-{\phi }_i\right)dxdy+\frac{\omega }{2}\underset{\Omega }{\int }\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{\left(1-H_{ϵ}\left({\phi }_j\right)\right)}^{2}dxdy.$

更进一步地，步骤四利用梯度下降法对步骤四所述的多相互斥水平集能量泛函进行求解， 得到如下偏微分方程：

$\frac{\partial {\phi }_i}{\partial t}=\mu div\left(p^{\prime }\left(\left|▿{\phi }_i\right|\right)\frac{▿{\phi }_i}{\left|▿{\phi }_i\right|}\right)+{\mathrm{\lambda \delta }}_{ϵ}\left({\phi }_i\right)div\left(g\frac{▿{\phi }_i}{▿{\phi }_i}\right)+{\mathrm{\alpha g\delta }}_{ϵ}\left({\phi }_i\right)+{\mathrm{\omega \delta }}_{ϵ}\left({\phi }_j\right)\underset{j=1,j\ne i}{\overset{N}{\Sigma }}\left(1-H_{ϵ}\left({\phi }_j\right)\right)$

其中div为散度算子，变量j为演化优先级参数。

3.有益效果

采用本发明提供的技术方案，与已有的公知技术相比，具有如下显著效果：

(1)本发明的一种基于多相互斥水平集的黏连细胞图像分割方法，以图像二值化常规处 理方法，得到代表细胞区域类圆斑块，通过Hough变换圆检测，获取类圆细胞斑块的圆半径 和圆心位置，能够统计出圆心的个数，解决细胞图中细胞精确计数的关键问题；

(2)本发明的一种基于多相互斥水平集的黏连细胞图像分割方法，以所有圆心为坐标设 置半径统一的圆形初始曲线，进一步提高了多相水平集初始曲线的设置效率，有利于后继水 平集曲线演化；

(3)本发明的一种基于多相互斥水平集的黏连细胞图像分割方法，在DRLSE模型基础 上增加互斥能量项有助于水平集曲线演化到黏连细胞区域边界处不过界，围绕各自细胞的圆 心自动完成彼此分离的轮廓曲线演化，解决了长期以来黏连区域细胞难以分割的问题；

(4)本发明的一种基于多相互斥水平集的黏连细胞图像分割方法，处理步骤少，所分割 的细胞轮廓线自然，可用于类圆细胞图的分割和计数，测量精度高，适应能力强。

附图说明

图1为本发明中多水平集函数重叠的示意图；

图2中的(a)～(f)为本发明中Hough变换圆心检测过程的示例图；其中，图(a)为 原始细胞图像；图(b)为二值化图像；图(c)为面积阈值法获得的图像；图(d)为Hough 变换法获得的圆心检测结果图；图(e)为多处黏连的圆心结果图；图(f)为黏连严重处的 圆心结果图；

图3中的(a)为利用CV模型对细胞图像的分割效果图；图3中的(b)为利用基于标 记的分水岭分割方法对细胞图像的分割效果图；图3中的(c)为本发明中利用基于多相互斥 水平集模型对细胞图像的分割效果图。

具体实施方式

为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。

实施例1

结合附图，本实施例的一种基于多相互斥水平集的黏连细胞图像分割方法，其步骤为：

步骤一、读取细胞图像，将细胞图像二值化，而后根据细胞的面积特征，利用阈值法， 反复实验选取效果较好的经验值对图像进行处理，将面积较小的噪声区域过滤掉，使得细胞 图中所有的细胞个体近似圆形斑块，获得代表细胞区域的近似圆形斑块图。

步骤二、针对获得的近似圆形斑块图，利用基于Hough变换的圆检测方法对图像进行处 理，求解细胞图像中全部近似圆形斑块的类圆圆心，同时通过调节模型中的参数去设置合适 的检测圆的半径范围，过滤掉半径过大圆和半径过小圆的干扰，以提高Hough变换圆检测方 法的准确率，进而获取每个类圆斑块的圆心，统计所得圆心个数，即为细胞图中细胞个数。

步骤三、为了避免最小细胞的CHT检测圆作为水平集初始轮廓在实际细胞区域之外，以 最小细胞的CHT检测圆半径的0.7倍作为所有细胞区域的初始圆半径，依据每个细胞类圆圆 心坐标设置多个水平集的圆形初始曲线。

步骤四、对多相互斥水平集能量泛函进行求解得到对应的偏微分方程，多相互斥水平集 模型将互斥能量项ε^R(φ_i)与DRLSE模型以多相水平集的方式进行结合，其中DRLSE能量项 用于引导轮廓曲线向边界演化，而互斥能量项ε^R(φ_i)用于保证多水平集之间的相互独立，对 于水平集函数φ_i，该模型能量泛函定义如下：

$\begin{array}{c}{E}_{MRLSE}={ϵ}_{ϵ}^{DRLSE}\left({\phi }_i\right)+\frac{\omega }{2}{ϵ}_{ϵ}^R\left({\phi }_i\right)\\ ={\mathrm{\mu R}}_p\left({\phi }_i\right)+{\mathrm{\lambda L}}_g\left({\phi }_i\right)+{\mathrm{\alpha A}}_g\left({\phi }_i\right)+\frac{\omega }{2}{ϵ}_{ϵ}^R\left({\phi }_i\right)\end{array}---\left(1\right)$

其中L_g(φ)为曲线长度项，λ为其加权系数(λ＞0)，可确保水平集函数演化的稳定性，其 求解公式如下：

$L_g\left(\phi \right)=\underset{\Omega }{\int }{\mathrm{g\delta }}_{ϵ}\left(\phi \right)|▿\phi |dxdy---\left(2\right)$

(1)式中α为面积项A_g(φ)的加权系数，通过设置系数的正负控制曲线演化的方向： α＞0，演化曲线由外向内收缩，α＜0，曲线向外扩张，Ω为图像全域。A_g(φ)定义如下；

$A_g\left(\phi \right)=\underset{\Omega }{\int }{\mathrm{gH}}_{ϵ}(-\phi )dxdy---\left(3\right)$

(2)和(3)式中g为边界检测函数，其定义如下:

$g=\frac{1}{1+|▿G_{\sigma }*I{|}^2}---\left(4\right)$

其中G_σ为高斯核函数，▽为梯度算子，I是原始图像。此函数利用图像的梯度信息引导 初始轮廓向目标边界演化。

(3)式中H_ε为Heaviside公式的近似光滑函数，其定义为：

$H_{ϵ}\left(x\right)=\left(\begin{array}{cc}\frac{1}{2}(1+\frac{x}{ϵ}+\frac{1}{\pi }sin(\frac{x}{ϵ})),& \left|x\right|\le ϵ\\ 1,& x>ϵ\\ 0,& x<-ϵ\end{array}\right)---\left(5\right)$

式中ε为常数，H_ε的导函数δ_ε(x)表示如下：

${\delta }_{ϵ}\left(x\right)=H_{ϵ}^{\prime }\left(x\right)=\left(\begin{array}{cc}\frac{1}{2ϵ}\left(1+\cos \left(\frac{\pi x}{ϵ}\right)\right),& \left|x\right|\le ϵ\\ 0,& \left|x\right|>ϵ\end{array}\right)---\left(6\right)$

(1)式中R_p(φ)为距离规则能量项，μ为其加权系数,R_p(φ)被定义如下：

$R_p\left(\phi \right)=\underset{\Omega }{\int }p\left(\right|▿\phi \left|\right)dxdy---\left(7\right)$

(1)式中ε^R(φ_i)为水平集函数φ_i的互斥能量项，ω为其系数，本文将其设置为是为 了求解偏微分方程后的系数一致性。ε^R(φ_i)表达式如下：

${ϵ}_{ϵ}^R\left({\phi }_i\right)=\underset{\Omega }{\int }\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{(1-H_{ϵ}\left({\phi }_j\right))}^{2}dxdy---\left(9\right)$

因此能量泛函(1)可有如下表示：

$E_{MRLSE}=\mu \underset{\Omega }{\int }p\left(\left|▿{\phi }_i\right|\right)dxdy+\lambda \underset{\Omega }{\int }{\mathrm{g\delta }}_{ϵ}\left({\phi }_i\right)dxdy+\alpha \underset{\Omega }{\int }{\mathrm{gH}}_{ϵ}\left(-{\phi }_i\right)dxdy+\frac{\omega }{2}\underset{\Omega }{\int }\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{\left(1-H_{ϵ}\left({\phi }_j\right)\right)}^{2}dxdy---\left(10\right)$

利用梯度下降法对公式(10)的能量泛函进行求解，得到如下偏微分方程：

$\frac{\partial {\phi }_i}{\partial t}=\mu div\left(p^{\prime }\left(\left|▿{\phi }_i\right|\right)\frac{▿{\phi }_i}{\left|▿{\phi }_i\right|}\right)+{\mathrm{\lambda \delta }}_{ϵ}\left({\phi }_i\right)div\left(g\frac{▿{\phi }_i}{▿{\phi }_i}\right)+{\mathrm{\alpha g\delta }}_{ϵ}\left({\phi }_i\right)+{\mathrm{\omega \delta }}_{ϵ}\left({\phi }_j\right)\underset{j=1,j\ne i}{\overset{N}{\Sigma }}\left(1-H_{ϵ}\left({\phi }_j\right)\right)---\left(11\right)$

其中div为散度算子，而变量j∈{1,2,3…i-1,i+1…N}为演化优先级参数，其取值区间决 定了水平集函数之间的演化优先程度。对于偏微分方程(11)，增加条件j≠k，则表示水平集 函数φ_i的演化优先级高于水平集函数φ_k，进而可调整各细胞之间的优先级关系。

为了较好的说明互斥能量项在水平集演化过程中的作用，设置了三个相互重叠的水平集 函数，如图1所示。其中粗黑实线、细黑实线和虚线三条可变形曲线分别代表了水平集函数 φ₁、φ₂、φ₃的演化曲线，且曲线内侧点的水平集函数值皆小于零，外侧皆大于零，曲线上点的 值可视为0，即零水平集函数。

对于水平集函数φ₁，其互斥能量项可表示为：

${ϵ}_{ϵ}^R\left({\phi }_1\right(x,y,t))=\underset{\Omega }{\int }\underset{j=2,3}{\Sigma }{(1-H_{ϵ}({\phi }_j(x,y,t)\left)\right)}^{2}dxdy---\left(12\right)$

根据图1中的区域划分，对水平集函数φ₁中每一点的互斥能量进行求解可得：

${ϵ}_{ϵ}^R\left({\phi }_1\right(x,y,t))\approx \left(\begin{array}{c}0,(x,y)\in {\Omega }_1\\ 1,(x,y)\in {\Omega }_{12}\cup {\Omega }_{13}\cup {\Omega }_2\cup {\Omega }_3\\ 2,(x,y)\in {\Omega }_{123}\cup {\Omega }_{23}\end{array}\right)---\left(13\right)$

式(13)表明，互斥能量项将给予水平集函数φ₁的重叠区域所有点一个正值能量，并 且此能量会根据同区域重叠数量的增加而变大，而此正值能量与曲线内侧点的水平集函数值 小于零的特性互相“排斥”，因此零水平函数的位置会受到该正值能量的影响，从而发生变 化，这样重叠区域内的点将会被演化到曲线外侧，而且其他水平集曲线内部区域的互斥能量 也为正值，这样能够较好的避免水平集函数在之后的演化过程中再次出现重叠的情况。因此 在多相互斥水平集模型中，水平集函数φ_i将在DRLSE能量项与互斥能量项的共同作用下进 行演化，并对重叠区域和其他水平集曲线内部区域进行排斥，从而达到分割黏连细胞的目的。

从数值求解的角度可将偏微分方程(11)表示为如下迭代形式：

$\frac{{\phi }_i^{n+1}-{\phi }_i^n}{\Delta t}=\mu div\left(p^{\prime }\right(|▿{\phi }_i^n|)\frac{▿{\phi }_i^n}{|▿{\phi }_i^n|})+{\mathrm{\lambda \delta }}_{ϵ}\left({\phi }_i^n\right)div\left(g\frac{▿{\phi }_i^n}{|▿{\phi }_i^n|}\right)+{\mathrm{\alpha g\delta }}_{ϵ}\left({\phi }_i^n\right)+{\mathrm{\omega \delta }}_{ϵ}\left({\phi }_i^n\right)\underset{j=1,j\ne i}{\overset{N}{\Sigma }}(1-H_{ϵ}({\phi }_j))---\left(14\right)$

将多个初始曲线作为零水平集函数带入上式中，进行水平集迭代演化。每次迭代结果 都会得到其对应的演化曲线，从而驱动曲线向目标边界靠近。重复上述步骤当时， 说明演化完全，即曲线到达目标边界位置。

本实施例在DRLSE模型基础上增加互斥能量项有助于水平集曲线演化到黏连细胞区域 边界处不过界，围绕各自细胞的圆心自动完成彼此分离的轮廓曲线演化，解决了长期以来黏 连区域细胞难以分割的问题。

采用含有不同黏连程度的乳腺细胞图像进行实验。首先将细胞图像二值化，而后根据细 胞的面积特征，利用阈值法过滤掉噪声区域，然后利用Hough变换检测方法求取细胞类圆斑 块圆心，限制圆半径范围为5～20像素点。图2中的(d)、(e)、(f)为检测结果，图中“+”号， 为检测出类圆细胞的圆心，实验结果表明，使用本实施例的检测方案具有类圆特性的细胞都 能被检测到圆心，即使在拥挤黏连严重的细胞群区域都具有良好的检测结果，可为水平集演 化模型提供较为可靠的初始曲线定位，同时较为准确地完成了细胞计数工作。

依据圆心坐标，以小于最小细胞区域的圆半径0.7倍作为统一半径参数，设置所有细胞 内的水平集初始轮廓，设置多相互斥水平集模型的参数：ω＝3，μ＝0.01，λ＝5，α＝-3。 而后对黏连细胞图像进行分割，并与CV模型[ChanTF,VeseLA.Activecontourswithout edges[J].IEEETransactionsonImageProcessing,2001,10(2):266-277]，以及基于标记的分水岭 方法[SunY,HeGJ.Segmentationofhigh-resolutionremotesensingimagebasedonmarker-based watershedalgorithm[C]//ProceedingsofInternationalConferenceonFuzzySystemsand KnowledgeDiscovery.WashingtonDC:IEEEComputerSocietyPress,2008,271-276]所得的分 割结果进行对比。其中CV模型的分割结果如图3中的(a)所示，该模型对于背景纹理较为 复杂的图像其分割结果便不够理想，且不具备对黏连细胞进行分割的能力。依据Hough圆检 测方法对细胞前景进行标记，利用基于标记的分水岭方法对细胞图像进行分割结果如图3中 的(b)所示，该方法能消除背景纹理的影响，有效分割出细胞前景图像，然而其并不具备处 理黏连细胞图像的能力。而多相互斥水平集模型分割结果如图3中的(c)所示，不受背景纹 理的影响，具有较好地获取细胞图像的轮廓，且所得轮廓较为圆滑，即使黏连严重的拥挤细 胞图像仍具有良好的处理能力。

表1各方法分割精度对比

由于Hough变换圆检测方法的细胞计数平均精度为93％，为多相互斥水平集模型提供了 可靠的初始曲线，使得本发明对黏连细胞的分割具有较高的准确性。由表1可知，相较于其 他方法，本发明方法分割精度更高，且其分割精度波动较低，能较好地满足细胞检验工作中 对细胞计数和形态分析的要求。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也 只是本发明的实施方式之一，实际的结构并不局限于此。所以，如果本领域的普通技术人员 受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结 构方式及实施例，均应属于本发明的保护范围。