基于混合轮廓模型的细胞核仁和细胞膜的分割方法

专利名称：基于混合轮廓模型的细胞核仁和细胞膜的分割方法
技术领域：
本发明涉及图像处理、生物医学、计算机视觉、计算方法、尤其 是生物细胞显微图像的分割方法。
背景技术：
这些年来，生物学发明和其转换为临床应用治疗迅速发展。病理 学家使用病人的组织切片病理图像，并在显微镜下检测。当检测这类 图像时，病理学家常常用组织图像检测出的细胞结构变化和细胞核仁 在细胞质中的比例变化来评估病变程度。自动检测出细胞的轮廓和其 核仁的方法显得尤为重要。
在细胞图像处理领域，已经有了大量的研究成果。但同时也有些 挑战存在于这些研究中。噪声污染是出现频率最高的一类问题。这类 问题产生于组织细胞的染色过程，因此在处理染色颜料时经常会出现 染色分布不均的问题。除了噪声，阴影效果经常出现在图像中使得细 胞图像表现的不平滑。
在细胞分析前，比如细胞形态和细胞行为，需要准确的细胞分割。 主动轮廓模型成为成功的模型之一。基于主动轮廓模型的技术具有更 好地估计细胞形态的潜力。现存的主动轮廓模型可以划分为两类基 于边界的轮廓模型和基于区域的轮廓模型。 一方面，基于边界的模型 直接使用灰度梯度信息驱使轮廓线朝着目标边界运动。因此这类模型 在处理弱边界物体对象时表现的不好。这是因为细胞图像由于在细胞 膜附近的低对比度从而显现出边界模糊性。另一方面，基于区域的模型通过某类区域描述子，旨在区别出不同灰度的区域。它驱使着轮廓 的运动，从某种程度上说，并且对初始轮廓线位置也不敏感。这种模 型相对于前一种模型更适合细胞分割。
T. Chan和L. Vese在2001年出版的IEEE图像处理学报第10巻 266-277页中发表的文章Active contours without edges (无需边界的主 动轮廓)提出了一种流行的基于区域的主动轮廓模型。这种模型已经 成功用于图像分割。但是图像中感兴趣区域在统计上经常表现的不平 滑，因此这种模型并不能直接适用于细胞图像。Li等人在2007年IEEE 计算机视觉和模式识别会议上1-7页中Implicit active contours driven by local binary fitting energy(局部二值拟合能量驱动的隐式主动轮廓模 型)提到的灰度非齐次性，这是在低信噪比图像中最大的问题，经常 出现在组织切片图和医学图像中。

发明内容
为了克服已有的细胞图像分割方法的精度低、分割效果差的不足， 本发明提供一种精度高、分割效果良好的基于混合轮廓模型的细胞核 仁和细胞膜的分割方法。
本发明解决其技术问题所采用的技术方案是
一种基于混合轮廓模型的细胞核仁和细胞膜的分割方法，所述分 割方法包括以下步骤
1)、将待分割的细胞图像依照混合轮廓模型建立能量函数，形式如
￡￡W,Cl,c2,y;,/2) = ￡MW,Cl,C2,/;,/2)+JpW)+JtW) (11);
其中，i^^,q,c2,/;,/2)表示基于区域模型和局部二值拟合模型的
能量函数，其计算公式为
￡MW,Cl,c2,y;,/2) = (i-Ag)￡MFW,y;,/2)<formula>formula see original document page 6</formula>
其中4是策略权重参数， 五，w,y;,A)如公式(6)中定义，五g((^,c2)如 公式(5)中定义
其中H是Heaviside函数，刑为符号距离函数，p;代表了 和^ 的比值，&代表了^2和^的比值；
该策略权重参数4(x,力定义如下 '
(l + alV/(x,力l) 其中IV/(x,力l表示图像/中的灰度梯度信
(8);
a是一正常数；
抑制水平集函数^偏离标准距离符号函数来规则化水平集函数， 该项定义如下，
(9);<formula>formula see original document page 6</formula>
(10);
尸(^= J5(iw"')卜1)血 长度约束项用来平滑零水平集轮廓， 湖=
拟合函数A和6如下所示，
<formula>formula see original document page 6</formula>
常数C;和。；
2)、采用最速下降法最小化公式(11)，如下
<formula>formula see original document page 6</formula>
(14)
其中^， ^， v，和^是常数，外)是Dirac函数；q和s如下:《(x卜J"人'。(j' —x)j/(x)-- 乂(.,Oi:力'， / = 1， 2.
通过离散化偏微分方程(14)得到迭代公式，水平集函数-在—o时算出 细胞的细胞核仁和细胞膜的边界。
本发明的技术构思为依据区域相互间差别最大化原则，这种模 型能够成功地探测局部模糊物体边界和分割出主要不同区域。此模型 定义的能量函数包括了两项主要项局部拟合项和全局拟合项。局部 拟合项产生强大的力吸引边界并使轮廓停在目标边界上，尽管物体边 界不十分清晰甚至有些模糊。全局拟合项在区域差别最大化原则下确 保了曲线能够提取出图像中的主要不同部分。除此之外，利用图像灰 度梯度信息引入了一种策略权重参数，它使得上述两项拟合项作用在 一起，从而构成了混合轮廓模型。这种参数模型引入了一类由局部和 全局信息作为驱动的混合力。
本发明的有益效果主要表现在精度高、分割效果良好。


图1是一种分割后的结果示意图。
图2是另一种分割后的结果示意图。
图3是将图像中嗜中性粒细胞核区域中面积小于该阈值的部分填
补掉的结果示意图。
具体实施例方式
下面结合附图对本发明作进一步描述。
一种基于混合轮廓模型的细胞核仁和细胞膜的分割方法，所述分 割方法包括以下步骤
1)、将待分割的细胞图像依照混合轮廓模型建立能量函数，形式如
￡，,Cl,c2,./;,/2) = ￡MW，Cl,C2，/;,/2)+/^)+z#) (11); 其中，r'W,c,，c"/,/。表示基于区域模型和局部二值拟合模型的<formula>formula see original document page 8</formula>
(7);
其中^是策略权重参数，￡，w,y;,/2)如公式(6)中定义 公式(5)中定义
其中是Heaviside函数，0w为符号距离函数，A代表了 A;和爿 的比值，&代表了 4和J的比值； ' 该策略权重参数/lg(x,;;)定义如下
<formula>formula see original document page 8</formula> (8);
其中 IW(x，力l表不图f象/中白勺灰度梯度<5， a 是一正常数；
抑制水平集函数^偏离标准距离符号函数来规则化水平集函数， 该项定义如下，
<formula>formula see original document page 8</formula> (9);
长度约束项用来平滑零水平集轮廓， ，=f )*
拟合函数和^如下所示，
脸
(12)
<formula>formula see original document page 8</formula>
常数禾口 。；
2)、采用最速下降法最小化公式(11)，如下
<formula>formula see original document page 8</formula>
(14)其中^，^， v，和〃是常数，《O是DiraC函数；A和e,如下
sW= JX^-x)i/(x)-y;(力i2如 i = 1， 2.
通过离散化偏微分方程(14)得到迭代公式，水平集函数^在—o时
算出细胞的细胞核仁和细胞膜的边界。
Chan和Vese提到的基于区域的主动轮廓模型，通过最小化简化 的Mumford-Shah泛函，能将一种图像分割为两类集合。这个模型基
本思想如下。假设图像区域Qc:W,并且/: Q^9l为给定的图像。
Mumford和Shah将图像分割问题考虑为寻找一个最佳的边界轮廓C， 将图像区域分割为近似于分段灰度常量区域w,和M。 。 C则表示边 界。于是在Chan-Vese模型中全局数据拟合项定义如下
E^o^cg- l(/—c)2血办十(/—C2)2血办 (1)
其中q和5分别代表在轮廓线C外和轮廓线C内的区域，q和c2是 轮廓线外和轮廓线内两个拟合常量。
此模型认为同在一个区域内的像素具有最大的相似性，从而弥补 了边界检测子的缺点。当轮廓线准确的捕获物体边界时，两项拟合项 使得能量函数拟合值最小。在每个分割区域内，聚合的像素均值分别 等于。和q 。因此带参数c;和q的拟合项是在内部区域齐次性准 则下，驱动轮廓曲线运动。
在图像分割中由于区域差别性是指导原则，因此不同区域间的差
别性应被看作是模型的驱动，如下所示
这类基于区域的主动轮廓模型能量函数是根据区域间差别最大化而定
义的。最小化在公式(2)中定义的能量函数等同于最大化区域间差别。 公式(2)定义了全局指导项。本文作者修改该指导项如下，￡ = ,-c2)2 、" 中^代表了整个图像区域。轮廓线C的曲线演化能量函数定义如下
￡G(Cl,CXCl-C2)(^ + ^%v (4) 其中^ = |^办和^ = |^办项分别定义了轮廓线外部和内部的区域面
积，iV代表轮廓线C的外法线方向。当轮廓线分割面积过小时，为了 让上述能量项(4)有个平稳的值，而不是表现出震荡发散的值，也就是.
说避免(4)式中分母项趋于零，将五g项修改如下
五cW,c', c2) = J(c! - c2)(/(x) - qp2 - c2a)//W(x))血 (5 ).
其中是Heaviside函数，^)为符号距离函数，/ 7代表了 jc/和J 的比值，A代表了^2和j的比值。
虽然在公式(5)中定义的能量函数f考虑区域间的差别性作为曲 线演化的规则，但是该区域模型没有办法处理在成像过程中产生的灰 度非齐次性和阴影效果。因为基于区域的模型将全局差异作为驱动力， 而忽视了局部细节信息，因此该模型需要一个局部拟合能量项来改进 模型在上述情况中的表现。
提出了局部二值拟合模型。其中被称为核函数的这个重要参数被 引入用来规定局部区域的范围。另外，力和/2是两个随空间变化的 拟合函数项，用来逼近局部灰度。在LBF模型中，这个局部数据拟合 项定义如下
4 f[-力i /(力-乂o) i2 + ^ j"[ &—力i /(力-/2(x) |2 (i - (6)
其中H是Heaviside函数，^是高斯核函数，当W增加时^w减小 并趋向于零。函数/(^H十算了点x附近的拟合程度，点x可以看作是局部区域的中心点。该LBF模型较Chan-Vese模型最大的优势就是它 比Chan-Vese模型更准确地捕捉到目标局部细节。 本文提出的模型引入全局和局部灰度信息项，充分利用了公式(5)中提 到的基于区域模型和局部二值拟合模型的优点。整个模型能量函数定 义如下，
<formula>formula see original document page 11</formula>
其中^是策略权重参数， &SFW，,,/2)如公式(6)中定义， E。W,CpC2)如
公式(5)中定义。
该策略权重参数足定义如下
<formula>formula see original document page 11</formula>
其中，IV/(x,力l表示图像/中的灰度梯度信，a是一疋常数。
这个策略权重参数实现了混合模型中的权衡因子思想。在混合能 量约束项的作用下，曲线的驱动力表现的具有策略性。当图像局部区 域趋于平坦或远离目标边界，IW(x,力l的值相对较小，因此公式(8)的右 边部分的值接近于一。混合模型被引导成最大化区域间差别模型。当 点(3C，^落边界所在的灰度非齐次区域内时，IW(x,州的值将会很大，因 此4(x,力的值趋近于零。混合模型具有提取局部边界的潜在能力。它结 合了局部拟合模型和全局区域差别模型的优点，并利用灰度梯度信息 引出其策略参数。
最小化能量函数(7)，意味着解偏微分方程。本文采用梯度最速下 降法，给出该偏微分方程的数值解。
另夕卜，为了得到更精确的曲线演化结果，本文还引入了文献[6]中
提到的一项拟合项，通过抑制水平集函数w扁离标准距离符号函数来规则化水平集函数。该项定义如下，
，=]"会(1，)卜1)2血 (9)
除此之外，本文引入在Chan-Vese模型中使用的长度约束项用来平滑 零水平集轮廓
iW)=J|w/ (10)
至此，本文介绍的整个模型定义如下
e￡W,Ci,C2，/;,/2) = ￡mW,ci,c2,/;,/2)+7^)+jlw) (11) 拟合函数力和/2如下所示，
C)/) (12)
你》
/2W =
(13)
h(州
常数q和q的定义和文献[13]中相同。
为了最小化公式(ll)中定义的关于-的能量函数，本文给出了采用
最速下降法的解决方案，如下
+(i - w⑨[,v, - V2]+w(,(2)
I，
.+Mvwv(!》i (14)
I，
其中^ 4, v,和^是常数，外)是Dimc函数.q和《如下
e,W=JX(y—x)|/(x)—y:O0l2^,/=1,2. (15) 通过离散化偏微分方程(14)得到迭代公式。水平集函数-在^o时给出 了目标的边界。更进一步，在表达式a^,力中，作者给出^ = ^2=1,
"60, 等于三。高斯核函数&中的a等于二。这些参数根据不同的
图像设置为不同的常数。
本实施例提出的模型可以检测嗜中性粒细胞的弱边界，同时准确地将核仁从细胞质中区别出来。全局驱动力通过利用最大化全局区域 差别确保了模型能够获取目标的边界。
采用淋巴肿瘤细胞作为实验图片，分割现象可以通过后续的形态 学处理消除。首先，将^0区域设置为零，并将^0区域设置为一，而 —O区域则是物体的边界。如此一来，可以得到如图l和图2所示的 二值图像。
其次，对于图1给出后续操作。设置某一阈值，将图像中嗜中性 粒细胞核区域中面积小于该阈值的部分填补掉，结果如图3所示。
权利要求
1、一种基于混合轮廓模型的细胞核仁和细胞膜的分割方法，所述分割方法包括以下步骤1)、将待分割的细胞图像依照混合轮廓模型建立能量函数，形式如EE(φ，c1，c2，f1，f2)＝EM(φ，c1，c2，f1，f2)+P(φ)+L(φ)(11)；其中，EM(φ，c1，c2，f1，f2)表示基于区域模型和局部二值拟合模型的能量函数，其计算公式为EM(φ，c1，c2，f1，f2)＝(1-λg)ELBF(φ，f1，f2)+λg(x，y)EG(φ，c1，c2)(7)；其中λg是策略权重参数，ELBF(φ，f1，f2)如公式(6)中定义，EG(φ，c1，c2)如公式(5)中定义EG(φ，c1，c2)＝∫(c1-c2)(I(x)-c1p2-c2p1)H(φ(x))dx(5)；其中H是Heaviside函数，φ(·)为符号距离函数，p1代表了Ac1和A的比值，p2代表了Ac2和A的比值；该策略权重参数λg(x，y)定义如下<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msub> <mi>&lambda;</mi> <mi>g</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><mo>(</mo><mn>1</mn><mo>+</mo><mi>&alpha;</mi><mo>|</mo><mo>&dtri;</mo><mi>I</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><mo>|</mo><mo>)</mo> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math></maths>其中 id="icf0002" file="A2009101015590002C2.tif" wi="13" he="3" top= "175" left = "35" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>表示图像I中的灰度梯度信，α是一正常数；抑制水平集函数φ偏离标准距离符号函数来规则化水平集函数，该项定义如下，<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><mi>P</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mo>=</mo><mo>&Integral;</mo><mfrac> <mn>1</mn> <mn>2</mn></mfrac><msup> <mrow><mo>(</mo><mo>|</mo><mo>&dtri;</mo><mi>&phi;</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>|</mo><mo>-</mo><mn>1</mn><mo>)</mo> </mrow> <mn>2</mn></msup><mi>dx</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math></maths>长度约束项用来平滑零水平集轮廓，<maths id="math0003" num="0003" ><math><![CDATA[ <mrow><mi>L</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mo>=</mo><mo>&Integral;</mo><mo>|</mo><mo>&dtri;</mo><mi>H</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mrow><mo>(</mo><mi>x</mi><mo>)</mo> </mrow> <mo>)</mo></mrow><mo>|</mo><mi>dx</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math></maths>拟合函数f1和f2如下所示，<maths id="math0004" num="0004" ><math><![CDATA[ <mrow><msub> <mi>f</mi> <mn>1</mn></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mrow><msub> <mi>K</mi> <mi>&sigma;</mi></msub><mo>*</mo><mrow> <mo>(</mo> <mi>H</mi> <mrow><mo>(</mo><mi>&phi;</mi><mo>)</mo> </mrow> <mi>I</mi> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>K</mi> <mi>&sigma;</mi></msub><mo>*</mo><mi>H</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow> </mrow></mfrac><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo></mrow> </mrow>]]></math></maths><maths id="math0005" num="0005" ><math><![CDATA[ <mrow><msub> <mi>f</mi> <mn>2</mn></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mrow><msub> <mi>K</mi> <mi>&sigma;</mi></msub><mo>*</mo><mo>[</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mrow><mo>(</mo><mi>&phi;</mi><mo>)</mo> </mrow> <mo>)</mo></mrow><mi>I</mi><mo>]</mo> </mrow> <mrow><msub> <mi>K</mi> <mi>&sigma;</mi></msub><mo>*</mo><mo>[</mo><mn>1</mn><mo>-</mo><mi>H</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mo>]</mo> </mrow></mfrac><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>常数c1和c2；2)、采用最速下降法最小化公式(11)，如下<maths id="math0006" num="0006" ><math><![CDATA[ <mrow><mfrac> <mrow><mo>&PartialD;</mo><mi>&phi;</mi> </mrow> <mrow><mo>&PartialD;</mo><mi>t</mi> </mrow></mfrac><mo>=</mo><msub> <mi>&lambda;</mi> <mi>g</mi></msub><mi>&delta;</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mrow> <mo>(</mo> <msub><mi>c</mi><mn>1</mn> </msub> <mo>-</mo> <msub><mi>c</mi><mn>2</mn> </msub> <mo>)</mo></mrow><mrow> <mo>(</mo> <mi>I</mi> <mo>-</mo> <msub><mi>c</mi><mn>1</mn> </msub> <msub><mi>p</mi><mn>2</mn> </msub> <mo>-</mo> <msub><mi>c</mi><mn>2</mn> </msub> <msub><mi>p</mi><mn>1</mn> </msub> <mo>)</mo></mrow> </mrow>]]></math></maths><maths id="math0007" num="0007" ><math><![CDATA[ <mrow><mo>+</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub><mi>&lambda;</mi><mi>g</mi> </msub> <mo>)</mo></mrow><mi>&delta;</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mo>[</mo><msub> <mi>&lambda;</mi> <mn>1</mn></msub><msub> <mi>e</mi> <mn>1</mn></msub><mo>-</mo><msub> <mi>&lambda;</mi> <mn>2</mn></msub><msub> <mi>e</mi> <mn>2</mn></msub><mo>]</mo><mo>+</mo><mi>v&delta;</mi><mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo></mrow><mi>div</mi><mrow> <mo>(</mo> <mfrac><mrow> <mo>&dtri;</mo> <mi>&phi;</mi></mrow><mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo></mrow> </mfrac> <mo>)</mo></mrow> </mrow>]]></math></maths><maths id="math0008" num="0008" ><math><![CDATA[ <mrow><mo>+</mo><mi>&mu;</mi><mrow> <mo>(</mo> <msup><mo>&dtri;</mo><mn>2</mn> </msup> <mi>&phi;</mi> <mo>-</mo> <mi>div</mi> <mrow><mo>(</mo><mfrac> <mrow><mo>&dtri;</mo><mi>&phi;</mi> </mrow> <mrow><mo>|</mo><mo>&dtri;</mo><mi>&phi;</mi><mo>|</mo> </mrow></mfrac><mo>)</mo> </mrow> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中λ1，λ2，v，和μ是常数，δ(·)是Dirac函数；e1和e2如下ei(x)＝∫Kσ(y-x)|I(x)-fi(y)|2dy，i＝1，2.(15)通过离散化偏微分方程(14)得到迭代公式，水平集函数φ在φ＝0时算出细胞的细胞核仁和细胞膜的边界。
全文摘要
一种基于混合轮廓模型的细胞核仁和细胞膜的分割方法，所述分割方法包括以下步骤1)将待分割的细胞图像依照混合轮廓模型建立能量函数，形式如E^E(φ，c₁，c₂，f₁，f₂)＝E^M(φ，c₁，c₂，f₁，f₂)+P(φ)+L(φ)；2)构造并使用最速下降法最小化公式，并通过离散化偏微分方程的迭代方法，利用水平集函数φ在φ＝0时算出细胞的细胞核仁和细胞膜的边界。本发明提供一种精度高、分割效果良好的基于混合轮廓模型的细胞核仁和细胞膜的分割方法。
文档编号G06T5/00GK101661614SQ200910101559
公开日2010年3月3日 申请日期2009年8月10日 优先权日2009年8月10日
发明者姚春燕, 秋 管, 赵明珠, 敏 陈, 陈胜勇 申请人:浙江工业大学