2. 湖南科技大学 计算机科学与工程学院, 湖南 湘潭 411100;
3. 中南大学 生物医学工程研究所, 湖南 长沙 410083;
4. 胜利石油管理局滨南医院, 山东 滨州 256600
, ZHAO Yu-qian1,3 
, ZENG Ye-zhan3 , HUANG Zhong-chao3 , ZHANG Bing-kui3,4 , ZOU Bei-ji1
2. School of Computer Science and Engineering, Hunan University of Science and Technology, Xiangtan 411100, China;
3. Institute of Biomedical Engineering, Central South University, Changsha 410083, China;
4. Binnan Hospital of Shengli Petroleum Administration Bureau, Binzhou 256600, China
细胞计数是各种病理检查的常规手段, 传统方法通过显微镜观察来人工检验细胞数目.长时间对着显微镜, 人眼容易疲劳, 错检率高, 且费时费力.研究基于机器视觉的显微细胞图像自动分割与识别方法具有重要意义.显微图像中细胞疏密程度、排列状态不一以及细胞间的接触和重叠等因素, 都给细胞图像的自动分割带来了很大困难.
现有的细胞图像自动分割方法主要包括形态学方法、分水岭算法、活动轮廓模型[1]、椭圆拟合[2]和凹点检测法[3]等.其中, 形态学方法通常通过选用特定的结构元素对图像进行腐蚀操作来分离细胞, 由于不同成像系统得到的细胞往往大小不一, 很难选择一个通用的结构元素来进行处理, 结构元素过小则不容易将重叠细胞分开, 结构元素过大又容易丢失小细胞[4].分水岭算法是一种基于拓扑理论的数学形态学的分割方法, 优点是速度快, 能够检测重叠细胞, 但它存在一个明显的缺点, 即对图像灰度变化极敏感, 当细胞区域的灰度分布不均时, 该方法极易造成过分割[5-6].活动轮廓模型的基本思想是首先定义轮廓曲线的能量函数和初始位置, 然后使轮廓曲线向能量减小的方向演化, 直至到达目标的边界位置[7-8].椭圆拟合和凹点检测法都不依赖于显微图像的灰度信息, 仅通过分析细胞轮廓的形状和凹凸性来达到分离重叠细胞的目的[9-10].例如, Bai等[11]提出基于凹点检测和椭圆拟合的重叠细胞分割方法, 主要包括轮廓预处理和椭圆拟合两个部分.在轮廓预处理中, 首先运用多边形近似提取能够有效地描述细胞轮廓的特征点, 然后通过计算切线夹角进行凹点检测, 识别出重叠细胞并初步定位重叠细胞的重叠位置.在椭圆拟合阶段, 首先根据检测到的凹点, 将重叠细胞轮廓分成多个边缘片段, 然后对所有片段分别进行椭圆拟合, 并根据预先设定的规则对拟合结果进行判断, 分离其中的重叠细胞.该方法中, 凹点检测的可行性建立在假设图像中细胞的大小和形状都基本一致的基础之上, 因此鲁棒性较差.此外, Wang等[12]提出运用形状分类来分割重叠细胞的方法, 首先运用训练好的形状分类器来判断当前区域是否属于重叠细胞区域.如果是, 则根据区域轮廓特征提取分离点对, 并利用原始图像的梯度信息寻找最优分割路径, 对当前重叠细胞区域进行分割.该方法对于过度重叠的细胞易产生过分割, 并且在寻找分离点对和最优分割路径的过程中, 算法复杂度高, 耗时长.
本文针对显微图像中细胞的形状特征, 提出基于SVM和椭圆拟合的重叠细胞自动分割方法.通过提取细胞连通区域的多个形状特征, 对该区域进行判断.对于被判定为重叠细胞的区域应用瓶颈检测寻找分离点对, 运用一种改进的椭圆拟合算法对分离边缘进行修正, 形成新的封闭区域, 最后对新生成的封闭区域进行循环检测, 直至所有重叠细胞分割完毕.
1 算法描述如图 1所示为本文算法的流程图.流程主要包括特征提取、SVM分类、分离点对检测和边缘修正等步骤.首先, 对输入图像进行二值化, 获取细胞区域；然后提取每一个细胞连通域的多个形状特征, 输入已训练的SVM分类器进行判断.若该连通区域被判定为重叠细胞, 则对该区域轮廓进行瓶颈检测获取分离点对, 并应用椭圆拟合对分离的细胞边缘进行修正, 得到新的封闭区域.对新生成的封闭区域重复上述特征提取与判断过程, 直到所有的连通区域均不存在重叠细胞, 算法结束.
|
图 1 细胞分割算法流程 Fig. 1 Flowchart of cell segmentation |
重叠细胞与单个细胞通常在形状上具有显著差异.本文通过提取细胞连通区域的形状特征对该区域进行判定.在进行特征提取之前, 首先对输入图像进行二值化获取细胞区域, 具体过程如图 2所示.首先将显微细胞图像由RGB色彩空间转化为灰度图像, 如图 2(a)所示；然后应用Otsu算法对灰度图像进行二值化, 提取细胞前景区域, 并对二值图像进行形态学孔洞填充, 得到图 2(b)所示的完整细胞区域.
|
图 2 细胞区域提取 Fig. 2 Cell region extraction |
对于每个细胞连通区域, 提取该区域特征并输入已训练好的SVM进行分类.SVM是一种有监督的机器学习方法, 通过寻求结构化风险最小来提高学习机泛化能力, 实现经验风险和置信范围的最小化, 从而达到在统计样本量较小的情况下, 获得良好统计规律的目的.它可以有效地解决模式分类问题和非线性回归问题, 本文使用的是SVM的分类功能.分类结果的好坏很大程度上取决于特征提取的优劣, 提取多个能够有效区分重叠与单个细胞的形状特征, 包括固性、凸性、离心率、半径的方差以及区域的面积, 且分别定义如下.
1) 固性(solidity).
| $ {\rm{solidity}} = \frac{{{S_{\rm{a}}}}}{{{S_{\rm{b}}}}}. $ | (1) |
式中：Sa为当前区域的面积, Sb为包含该区域的最小凸多边形的面积, 如图 3所示.Solidity取值为(0, 1], 由于单个细胞一般为凸区域, 计算得到的Solidity非常接近于1, 而重叠细胞一般为非凸区域, 它对应的Solidity明显小于1.
|
图 3 固性计算示意图 Fig. 3 Illustration of computation of solidity value |
2) 凸性(convexity).
| $ {\rm{convexity}} = \frac{{{c_{\rm{a}}}}}{{{c_{\rm{b}}}}}. $ | (2) |
式中：ca为当前区域轮廓的周长；cb为包含该区域的最小凸多边形的周长, 即图 3中白色区域和虚曲线所示区域的周长.convexity为大于等于1的正数, 由于重叠细胞的凹度明显大于单个细胞, 由重叠细胞计算得到的convexity明显大于单个细胞.
3) 离心率(eccentricity).
| $ {\rm{eccentricity}} = \frac{L}{W}. $ | (3) |
eccentricity表示与当前区域具有相同标准二阶中心矩的椭圆的离心率, 其中L和W分别为该椭圆的半焦距和长半轴的长度.单个细胞通常接近于圆形结构, 而重叠细胞形状不规则、多细胞聚合往往呈长条形, 因此由重叠细胞计算得到的离心率明显大于单个细胞.
4) 半径的方差(=var_radius).
| $ {\rm{var}}\_{\rm{radius}} = \frac{1}{n}\sum\limits_{i = 1}^n {{{({r_i}-\bar r)}^2}} . $ | (4) |
式中：ri为当前区域边缘像素点到该区域质心的距离, 即半径；r为半径的均值；n为边缘像素点个数.与重叠细胞相比, 单个细胞的形状规则、边缘像素点到质心的距离基本一致, 因此由单个细胞计算得到的var_radius明显小于重叠细胞.
此外, 当前区域中的像素总个数, 即面积, 被作为特征之一用于区分重叠与单个细胞.通常, 对于同一幅图像中的同种细胞, 重叠细胞的面积明显大于单个细胞的面积.
1.2 分离点对检测用于分离重叠细胞的分离点对往往出现在细胞轮廓的狭窄部位[12], 即类似瓶颈的位置, 如图 4中的点pa和pb所示.若当前区域被判定为重叠细胞区域, 则对该区域进行瓶颈检测, 可以获得分离点对.对于细胞轮廓上任意的两个像素点pa和pb, 它们之间的瓶颈率e(pa, pb)定义为
|
图 4 瓶颈检测示意图 Fig. 4 Illustration of bottleneck detection |
| $ e\left( {{p_{\rm{a}}}, {p_{\rm{b}}}} \right) = \frac{{{\rm{dist}}\left( {{p_{\rm{a}}}, {p_{\rm{b}}}} \right)}}{{{\rm{min}}\left\{ {{\rm{length}}\left( {{p_{\rm{a}}}, {p_{\rm{b}}}} \right), {\rm{length}}\left( {{p_{\rm{b}}}, {p_{\rm{a}}}} \right)} \right\}}} $ | (5) |
式中：dist(pa, pb)为点pa与pb之间的欧氏距离, length(pa, pb)为点pa沿着细胞边缘顺时针行进至点pb的路径长度, min{·}表示取两条路径中长度较小的值.如图 4所示, length(pa, pb)为两点间上方边缘的长度, length(pb, pa)为两点间下方边缘的长度.计算细胞轮廓上所有像素点两两之间的瓶颈率, 其中瓶颈率最小的边缘点对(pa*, pb*)为分离点对：
| $ \left( {p_{\rm{a}}^*, p_{\rm{b}}^*} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{{p_{\rm{a}}}, {p_{\rm{b}}}} e({p_{\rm{a}}}, {p_{\rm{b}}}). $ | (6) |
若细胞轮廓上的像素点个数为N, 则瓶颈率的计算次数为CN2.为了减少计算次数、提高算法效率, 首先对当前细胞轮廓进行多边形近似[11], 提取轮廓特征点, 再通过计算特征点两两之间的瓶颈率达到检测分离点对的目的.如图 5所示为多边形近似结果示例, 其中, 图 5(a)为重叠细胞的原始轮廓, 图 5(b)为对原始轮廓进行多边形近似得到的特征点.可以看出, 多边形近似能够有效地反映细胞区域的大致轮廓, 且保留轮廓中的一些拐点和凹点.瓶颈检测通常是定位细胞轮廓的凹陷部位, 因此运用多边形近似得到的特征点进行分离点对检测将大大减少计算次数、提高运算效率, 且不影响分割精度.
|
图 5 多边形近似 Fig. 5 Application of polygonal approximation |
为了完整有效地分割某一连通区域的所有重叠细胞, 在检测到分离点对后, 须对分离边缘进行修正.边缘修正的目的是寻找适当的曲线连接分离点对, 将重叠细胞区域分离成两个新的封闭区域, 再对新生成区域进行循环检测, 直至完整分离出其中所有的重叠细胞.大多数现有的细胞分割算法都直接运用直线或其他最短路径准则连接分离点对[9, 12].对于重叠细胞的迭代分割, 这些方法可能导致某些分离的单个细胞被误判为重叠细胞, 从而产生过分割.基于显微图像中细胞的圆形或椭圆形结构, 提出新的基于椭圆拟合的分离边缘自动修正方法.如图 6所示为边缘修正示意图, 图 6(a)为某一重叠细胞的原始轮廓, 图 6(b)、(c)为分离细胞的边缘修正结果, 空心点之间的椭圆曲线为修补的边缘.边缘修正的具体步骤如下.
|
图 6 边缘修正示意图 Fig. 6 Schematic of edge modification |
1) 在分离点处的边缘上各选取一段局部边缘.如图 6(b)、(c)所示, 空心点(为了显示得更清晰, 图中进行了放大)为分离点对, 星形点为局部边缘像素.
2) 对选取的局部边缘应用最小二乘算法[13]进行椭圆拟合.为了确保该椭圆经过分离点对(x1, y1)和(x2, y2), 应用式(7) 对拟合椭圆进行约束.
| $ \left. \begin{array}{l} a{x^2}_1 + b{x_1}{y_1} + c{y^2}_1 + d{x_1} + e{y_1} + f = 0, \\ a{x^2}_2 + b{x_2}{y_2} + c{y^2}_2 + d{x_2} + e{y_2} + f = 0. \end{array} \right\} $ | (7) |
3) 将拟合的椭圆在分离点对处断开, 分成两段椭圆弧线, 并分别与分离边缘形成一个封闭区域；运用式(1) 分别计算两个封闭区域的solidity, 取solidity较大的区域对应的椭圆弧线作为修补的边缘, 如图 6(b)、(c)的实椭圆曲线所示.
2 实验结果与分析为了验证该算法的有效性, 从美国社会血液学数据库[14]中挑选了40幅血液细胞图像作为实验数据集, 其中20幅为训练集, 20幅为测试集.如图 2(a)所示为血液细胞图像样本.可以看出, 细胞密度较高, 存在多细胞粘连和重叠的现象.在实验过程中, 首先从训练图像中分别选取500个单个细胞区域和500个重叠细胞区域作为正、负样本集, 并分别提取每个区域的5个形状特征, 即固性、凸性、离心率、面积以及半径的方差, 构建一个1 000×5维的特征向量, 输入SVM进行训练.对于测试图像, 提取每个细胞连通域的上述5个形状特征, 输入训练好的SVM进行判断.SVM分类器的核函数为线性核, 且误差惩罚参数C设置为1.
如图 7所示为运用不同算法对图 2(a)进行分割的结果, 图 7(a)~(f)分别为文献[4]、[5]、[9]、[12]、[15]和本文方法的分割结果.可以看出, 图 7(a)~(e)对于重叠程度较严重的细胞都存在一定程度的欠分割, 对于从重叠细胞簇中分离得到的形状不规则的单个细胞容易出现过分割, 本文方法能够较准确地分离重叠细胞, 且可以有效地避免单个细胞的过分割.
|
图 7 图 2(a)所示的血液细胞图像分割结果比较 Fig. 7 Segmentation result comparison of blood cell image shown in Fig. 2 (a) |
为了更详细阐述该方法的优势, 图 8展示了几组具有代表性、比较难分割的重叠细胞的分割结果比较.图中, 第一行为5幅原始细胞图像, 这些图像中包括双细胞和多细胞重叠的情况；第(b)~(f)行分别为文献[4]、[5]、[9]、[12]和[15]的分割结果；最后一行为本文方法的分割结果.由于文献[4]是应用形态学腐蚀和重构进行分割, 当细胞重叠程度较大时, 会出现欠分割现象(第(b)行的第1、第3~5列).文献[5]是一种基于分水岭的分割方法, 由于在分割前加入了梯度检测, 有一定的过分割抑制作用, 但当重叠细胞间无明显边界时, 采用该方法不能有效地分割(第(c)行第3列).Farhan等[9]首先在边缘上寻找凹点, 然后进行配对, 配对过程可能会出现误配(第(d)行的第1、2列), 且对于凹陷不明显的细胞不能有效地分割(第(d)行的第3、5列).与本文类似, Wang等[12]通过运用形状分类器对重叠区域进行判断, 不同的是, 在确定分离点对后, 该方法根据原始图像的梯度信息, 运用Dijkstra算法寻找最优分割路径分离重叠细胞.对于细胞重叠程度较大的区域, 在成功分离其中的某一单个细胞之后, 会在单个细胞上形成不规则边缘, 从而导致单个细胞在迭代分割过程中被误判为重叠细胞, 造成过分割(第(e)行的第3、4列).Liao等[15]利用细胞轮廓的瓶颈率和椭圆拟合误差对重叠细胞进行识别和分离, 当重叠细胞形状类似椭圆结构时, 该方法失效(第(f)行第5列).本文方法通过提取连通区域的多个形状特征输入SVM进行判断, 区分单个细胞与重叠细胞, 即使重叠细胞的形状类似椭圆结构, 也可以通过连通域的离心率、面积以及半径的方差等特征对该区域进行识别(第(g)行第5列).在确定分离点对后, 针对细胞的形状特征, 提出运用椭圆拟合对分离细胞的边缘进行修正.修正后的细胞区域反映了重叠在一起的细胞的真实形状, 可以在循环检测阶段有效地抑制欠分割和过分割(第(g)行第1~5列).
|
图 8 重叠细胞分割结果比较 Fig. 8 Segmentation result comparison of overlapping cells |
为了进一步证实该方法的有效性, 应用分割准确性AC对分割结果进行评价：
| $ {A_{\rm{C}}} = \frac{{{N_{{\rm{seg}}}}}}{{{N_{{\rm{seg}}}} + {N_{{\rm{spl}}}} + {N_{{\rm{mer}}}} + {N_{{\rm{add}}}} + {N_{{\rm{mis}}}}}}. $ | (8) |
式中：Nspl和Nmer分别为过分割和欠分割的细胞数目, Nadd为将背景误分割为细胞的数目, Nmis为将细胞被误分割为背景的数目, Nseg为正确分割的细胞数目.
如表 1所示为不同算法分割结果的定量统计分析.可以看出, 本文算法相较其他算法具有较高的分割准确率.文献[4]的分割结果中, Nmer和Nmis较高, 说明采用该方法不能有效地分开重叠细胞, 且由于形态学腐蚀作用, 使得较小细胞容易丢失.文献[5]的分割准确性较高, 但存在一定的过分割和欠分割现象.文献[9]方法的Nmer较小, Nspl较大, 说明采用该方法能够有效地抑制欠分割, 但过分割现象相对严重.这是因为二值化后的细胞边缘不平整, 存在大量局部凹点, 而该方法不能进行有效地凹点匹配与判断.Wang等[12]在确定分离点对后, 利用细胞间的梯度信息分离重叠细胞, 因此往往容易在分离的单个细胞上形成不规则边缘, 导致在循环检测中被误判为重叠细胞, 从而使得分割结果出现较严重的过分割现象, 产生较高的Nspl.Liao等[15]首先根据细胞轮廓的瓶颈率寻找候选分离点对, 然后通过细胞轮廓各部分的拟合误差进一步对分离点对进行判断.若细胞重叠程度较严重或重叠细胞类似椭圆结构, 则该方法将无法检测, 从而产生较高的Nmer.本文方法通过SVM分类器判断重叠细胞, 运用椭圆拟合进行边缘修正, 可以有效地避免欠分割与过分割, 具有较高的分割准确率.
|
表 1 不同算法的分割准确性比较 Table 1 Segmentation accuracy comparison of differentmethods |
表 2给出不同算法对测试图像中细胞粘连与否判断的敏感度Sen和特异度Spe, 即分别将重叠细胞和单个细胞正确识别并分割的百分率.可以看出, 本文方法对于显微图像中细胞粘连与否的判断敏感度和特异度分别达到了93.9%和91.3%, 均高于文献[4, 5, 9]和[12], 表明本文方法与这些方法相比, 在识别重叠细胞以及抑制单个细胞过分割方面具有明显的优势.本文方法依据细胞连通域的形状特征对重叠和单个细胞进行判断, 因此, 当单个细胞形状不规则时, 本文方法易对该细胞产生过分割.Liao等[15]在重叠与单个细胞识别过程中, 设计了一系列的既定规则, 可以有效避免单个细胞的过分割, 因此由文献[15]方法计算得到的特异度稍高于本文方法.此外, 文献[15]方法对于细胞粘连与否判断的敏感度低于本文方法, 表明采用本文方法能够更有效地识别重叠细胞.
|
|
表 2 不同算法对于细胞粘连与否判断的敏感度和特异度 Table 2 Sensitivity and specificity of different methods foroverlapping cell recognition |
本文算法中时间复杂度最高的部分为瓶颈检测, 时间复杂度为O(n2), 当轮廓像素点数增加时, 计算次数快速增长.为了减少计算次数, 提高算法效率, 在进行瓶颈检测前, 运用多边形近似提取轮廓特征点.图 8所示的5幅图像在采用多边形近似(M+APP)和不采用多边形近似(M-APP)两种情况下算法运行时间的平均值和标准差(Intel双核i3-3.29 GHz CPU, 4.0 Gb RAM的PC机, Windows XP环境下Matlab 7.11编程), 分别为(5.97±2.75) s和(3.89±1.62) s.可以看出, 采用多边形近似可以有效地提高算法效率, 减少运算时间.
3 结语本文提出新的基于SVM和椭圆拟合的细胞图像自动分割方法.根据细胞连通区域的多个形状特征, 利用SVM识别单个与重叠细胞；对于重叠细胞, 运用瓶颈检测获取分离点对；将多边形近似应用于瓶颈检测, 减少计算时间, 提高算法效率；对于分割后的边缘, 提出基于椭圆拟合的边缘自动修正方法.结果表明, 本文算法在细胞凹点不明显的情况下能够有效地分开重叠细胞, 且能够有效地抑制过分割现象, 分割准确率达到91.6%, 明显高于其他分割方法.
| [1] | 张明慧, 卢振泰, 张娟, 等. 基于多图谱活动轮廓模型的脑部图像分割[J]. 计算机学报, 2016, 39(7): 1490–1500. DOI:10.11897/SP.J.1016.2016.01490 |
| [2] | CHO M. Performance comparison of two ellipse fitting-based cell separation algorithms[J]. Journal of the Royal Society of Medicine, 2015, 13(3): 215–219. |
| [3] | WEI L. Concave regularisation-based non-rigid feature point matching algorithm[J]. Electronics Letters, 2015, 51(8): 621–623. DOI:10.1049/el.2014.3807 |
| [4] | HUANG J D. An improved algorithm of overlapping cell division[C]//2010 International Conference on Intelligent Computing and Integrated Systems. Guilin: IEEE, 2010: 687-691. |
| [5] | LIM H N, MASHOR M Y, HASSAN R. White blood cell segmentation for acute leukemia bone marrow images[C]// 2012 International Conference on Biomedical Engineering. Penang, Malaysia: IEEE, 2012: 357-361. |
| [6] |
王鑫, 胡洋洋, 杨惠中. 基于迭代腐蚀的粘连细胞图像分割研究[J].
南京理工大学学报:自然科学版, 2016, 40(3): 285–289.
WANG Xin, HU Yang-yang, YANG Hui-zhong. Segmentation of adherent cell image based on iterative erosion[J]. Journal of Nanjing University of Science and Technology, 2016, 40(3): 285–289. |
| [7] | WU P, YI J, ZHAO G, et al. Active contour-based cell segmentation during freezing and its application in cryopreservation[J]. IEEE Transactions on Biomedical Engineering, 2015, 62(1): 284–295. DOI:10.1109/TBME.2014.2350011 |
| [8] | KAUR S, SAHAMBI J S. Curvelet initialized level set cell segmentation for touching cells in low contrast images[J]. Computerized Medical Imaging and Graphics, 2016, 49: 46–57. DOI:10.1016/j.compmedimag.2016.01.002 |
| [9] | FARHAN M, YLIHARJA O, NIEMISTO A. An improved clump splitting method for convex objects[C]//Proceedings of the 7th International Workshop on Computational Systems Biology. Luxembourg: TICSP, 2010, 35-38. |
| [10] | LATORRE A, ALONSONANCLARES L, MUELAS S, et al. Segmentation of neuronal nuclei based on clump splitting and a two-step binarization of images[J]. Expert Systems with Applications, 2013, 40(16): 6521–6530. DOI:10.1016/j.eswa.2013.06.010 |
| [11] | BAI X, SUN C, ZHOU F. Splitting touching cells based on concave points and ellipse fitting[J]. Pattern Recognition, 2009, 42(11): 2434–2446. DOI:10.1016/j.patcog.2009.04.003 |
| [12] | WANG H, ZHANG H, RAY N. Clump splitting via bottleneck detection and shape classification[J]. Pattern Recognition, 2012, 45(7): 2780–2787. DOI:10.1016/j.patcog.2011.12.020 |
| [13] | FITZGIBBON A W, PILU M, FISHER R B. Direct least squares fitting of ellipses[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1999, 21(5): 476–480. DOI:10.1109/34.765658 |
| [14] | The American society hematology[EB/OL]. 2010-06-01. http://imagebank.hematology.org. |
| [15] | LIAO M, ZHAO Y Q, LI X H, et al. Automatic segmentation for cell images based on bottleneck detection and ellipse fitting[J]. Neurocomputing, 2016, 173(P3): 615–622. |