基于优化水平集的细胞图像阈值分割算法

张瑞华

(江汉大学 物理与信息工程学院，武汉 430000)

阈值法[1-5]由于其计算量小、简单有效而成为图像分割中使用较广泛的一种方法.何国生等[1]通过分析细胞灰度图像的分布特性，应用新的自适应鉴别特征方法来分割干细胞，该方法对椭圆类含噪细胞图像的效果较好，但对其他形状细胞的分割结果不太理想.夏欣等[6]利用最大类间方差法进行迭代处理，自适应地搜寻最佳分割阈值，该算法对带噪、对比度低的粘连细胞图像的分割效果一般.杨辉华等[7]提出一种基于水平集和凹点区域检测的细胞分割算法，与传统水平集相比，该方法在准确度上有所提高，但在实时性上有待加强.

传统水平集方法为C-V模型[6]和Shi快速水平集法[7].前者由于不采用目标边界的梯度算子，而是基于区域灰度信息，因而对弱噪声图像的分割效果较好，但由于计算量大，在实时性上难以保证[8]；后者为采用双链表方式的窄带模型，该算法实现简单，但速度函数中的阈值无法自适应调整[9].本文将Bayes分类器和KNN分类器引入到水平集外部速度函数的设计中，两种分类器轮流作用，无需设定阈值便能产生水平集驱动力.同时，算法将Shi模型的双链表和C-V模型的全局分割相结合，以加快曲线演化.

本文方法分为2个步骤，首先通过优化水平集算法对细胞图像进行边缘检测，然后将目标与背景的类内平均距离引入到OTSU阈值法的阈值选择函数中，对细胞边缘图像进行分割.试验结果表明，本文算法相较与传统水平集和阈值法，能更准确的实现对带噪、对比度低的复杂粘连细胞图像的分割.

1 优化水平集算法

1.1 算法原理

在Shi快速水平集算法中背景Ωb和目标Ω0是基于边缘C的不断演化而分割的，其本质是曲线C对图像像素点x的二分类，即若x∈Lout/Lin，判断x∈Ω0/Ωb.其中Lin、Lout分别为Shi算法的内部链表和外部链表.基于此，将Bayes分类器和KNN分类器引入到水平集外部速度函数的设计中.其中驱动力h(x)为：

(1)

式中，I(x),ε,cb,c0分别是像素x灰度值，某正数，背景灰度均值和目标灰度均值.p(Ω0)和p(I(x)|Ω0)为目标区域的先验概率和条件概率密度；p(Ωb)和p(I(x)|Ωb)为背景区域的先验概率和条件概率密度.

(2)

(3)

式中：Num[φ(x)<0]为目标类像素点总数，Num[φ(x)>0]为背景类像素点总数.但当l0(x)=lb(x)，即|I(x)-c0|=|I(x)-cb|时分类失败.因此，将Bayes分类器g(x)引入到h(x)的设计中.

如上所述，q0(x)=qb(x)时Bayes分类器失效，l0(x)=lb(x)时KNN分类器失效，而本文将Bayes分类器和KNN分类器同时引入到水平集外部速度函数的设计中，两种分类器轮流作用，当一种分类器失效时，另一种分类器将产生水平集驱动力以促进曲线继续演化.只有当∀x∈Lout(或∀x∈Lin)且q0(x)=qb(x),l0(x)=lb(x)成立时，演化终止.

为了简化计算，这里以Ω0和Ωb内像素点的个数来近似等效先验概率p(Ω0)和p(Ωb),即p(Ω0)=∑(φ(x)<0)/Num(x∈Ω)，p(Ωb)=∑(φ(x)>0)/Num(x∈Ω)；用Ω0和Ωb内像素点个数的统计直方图来近似等效类条件概率密度p(I(x)|Ω0)和p(I(x)|Ωb),即p(I(x)|Ω0)=Num(I(x))/∑φ(x)<0,p(I(x)|Ωb)=Num(I(x))/∑φ(x)>0.同时采用加权均值滤波器G(x)以提高h(x)的抗噪性能.

h′(x)=G(x)⊗h(x).

(4)

综上所述，Fext(x)设计如下：

(5)

其中h′(x)=G(x)⊗h(x)，

(6)

综上所述，本文将Bayes分类器和KNN分类器引入到水平集外部速度函数的设计中，两种分类器轮流作用，无需设定阈值便能产生水平集驱动力.算法将Shi模型的双链表和C-V模型的全局分割相结合，以加快曲线演化.

1.2 算法步骤

本文提出的的优化水平集算法可按以下步骤进行.

1)初始化.设置Shi快速水平集中双链表Lin、Lout和水平集函数φ(x)，根据φ(x)初始化p(Ω0),p(Ωb),p(I(x)|Ω0),p(I(x)|Ωb)，设置曲线演化的最大迭代次数Nmax.

2)Fori=1:Nmaxdo

第1步 (第一个循环，在外部速度函数Fext(x)的作用下曲线的演化过程)

Fori=1:Noutdo

向外演化：根据式(5)计算链表Lout中每点的驱动力h′(x)，进而得到Fext(x)的符号，若Fext(x)>0，设置φ(x)=-1，如果∀y∈N4(x)满足φ(y)=3，则令φ(y)=1.

清除Lin冗余点：检查链表Lin上的每个点x，若满足∀y∈N4(x)，且φ(y)<0，则将点x作为冗余点删除，同时设置φ(x)=-3.

更新驱动力h′(x)：根据φ(x)符号，更新p(Ω0),p(Ωb),p(I(x)|Ω0)，p(I(x)|Ωb)和平均灰度值c0,cb.

向内演化：根据式(5)计算链表Lin中每点的驱动力h′(x)，进而得到Fext(x)的符号，若Fext(x)<0，设置φ(x)=1，如果∀y∈N4(x)满足φ(y)=-3，则令φ(y)=-1.

清除Lout冗余点：检查链表Lout上的每个点x，若满足∀y∈N4(x)，且φ(y)>0，则将点x作为冗余点删除，同时设置φ(x)=3.

更新驱动力h′(x)：更新p(Ω0),p(Ωb),p(I(x)|Ω0)，p(I(x)|Ωb)和平均灰度值c0,cb.

终止条件判断：∀x∈Lout或∀x∈Lin若满足q0(x)=qb(x)和l0(x)=lb(x)，则跳转至第2步；否则，返回第1步.

第2步 (第二个循环，在内部速度函数Fin(x)的作用下曲线的演化过程)

Fori=1:Nindo

循环过程与第1步的基本一致，只需将其中的Fext(x)用Fin(x)代替，并去掉其中的终止条件判断步骤，Fin(x)中的G(x)采用加权均值滤波器.

第3步 迭代停止条件.

∀x∈Lout或∀x∈Lin满足q0(x)=qb(x)且l0(x)=lb(x)或者达到最大迭代次数Nmax，则曲线演化停止，退出主循环，否则返回第1步.

神经元干细胞(neural stem cells，NSC)图像具有边缘模糊、对比度低的特点，对其分别采用本文算法及三种经典水平集算法:Shi算法、C-V算法、Li算法进行边缘检测，结果如图1所示.本文算法的初始轮廓为位于图像中心的30×42像素矩形曲线，参数设置如下：λ0=λb=1,Nmax=80,Nout=20,Nin=3；C-V算法的参数设置为u=0.01×2552,λ0=λb=1；v=0,ε=1,Δt=0.1，Li算法的参数设置为：u=0.04,v=1.5，λ=5，Δt=0.1.

(a) 迭代次数为0 (b) 迭代次数为30 (c) 迭代次数为50 (d) 迭代次数为76

(e) Shi算法 (f) 本文算法 (g) C-V算法 (h) Li算法图1 NSC图像的边缘检测结果Fig.1 Edge detection results of NSC cell image

本文算法的检测过程如图1(a)～(d)，可见当迭代次数为76时已收敛到边缘了.Shi算法中阈值I1、I2设置为[175,252]，由图可知，本文算法和Shi算法结果几乎一致，都准确地检测到边缘.C-V算法的检测结果如图1(g)，虽收敛到边缘但伴随着一些冗余轮廓.Li算法的检测结果如图1(h)，未能收敛到边缘便已停止演化了.

对白细胞(white blood cells,WBC)图像加入(0.05,0.2)高斯噪声，如图2所示.可见，Li算法对噪声较敏感未能收敛到边缘便已停止演化，这是由于该算法是以图像边缘梯度信息作为驱动力.当图像的噪声干扰较大时，基于区域灰度信息的C-V算法会将很多噪声点误判为边缘点，如图2(b)所示.Shi算法实现简单，但速度函数中的阈值无法随噪声自动调整，因而将部分噪声误判为边缘，如图2(c)所示.本文算法在强噪声干扰下仍能准确演化到边缘，如图2(d)所示.

(a) Li算法 (b) C-V算法 (c) Shi算法 (d) 本文算法图2 添加(0.05,0.2)高斯噪声的WBC图像边缘检测结果Fig.2 Edge detection results of the WBC image with (0.05,0.2) Gaussian noise

2 改进的OTSU阈值法

经典OTSU 算法[10]和局部递归OTSU 算法[11-12]都是基于灰度统计的分割算法，它们将背景与目标的类间方差最大值作为阈值选择的依据，而忽略了类内平均距离.基于此，本文在阈值函数的设计中引入目标与背景的类内平均距离：

(7)

其中：w0,w1分别为目标点和背景点占整幅图片的像素点比例，u0、σ0、u1、σ1分别为目标点和背景点的平均灰度和类间方差.

采用优化水平集获得细胞轮廓,如图3(b)所示，用改进OTSU法分割,结果如图3(c)所示，最后通过开、闭[13]获得细胞分割结果，如图3(d)所示.

(a)原始图像阈值法分割结果图3 优化水平集的OTSU 阈值分割过程Fig.3 The OTSU threshold segmentation process based on optimized level set

3 试验结果分析

图4(a)按从左到右的顺序为随机抽取某个NSC细胞图像序列中第2、25、55、85、110、130、140帧图像，图4(b)～(d)给出了三种算法的分割结果.文献[1]对OTSU阈值法进行了改进，引入预测因子k和松弛因子s，参数k,s初始值设置为：k=120,s=10.文献[6]的参数设置为：u=0.03,v=1.4,λ=4.5，Δt=0.1.本文算法的参数设置为：λ0=λb=1,Nmax=80,Nout=20,Nin=3.

图4(b)显示了文献[1]对序列的分割结果，出现了多个粘连处未分割且有丢失细胞的情况.这是由于该算法通过分析细胞灰度图像的分布特性，应用新的自适应鉴别特征方法来分割细胞，该方法对低噪细胞图像的分割效果较好，但对细胞边缘模糊、对比度低的粘连细胞图像的分割结果不太理想.文献[6]的分割结果如图4(c)所示，该算法对圆形细胞的分割较好，但对不规则形状粘连细胞的分割结果差强人意，这是由于该算法是基于圆型先验模型来检测细胞边缘.本文算法的分割结果见图4(d)，除丢失了极个别细胞，所有粘连处均得到正确分割.

(a) 原始图像

(b) 文献[1]阈值法分割结果

(c) 文献[6]水平集算法分割结果

(d) 本文算法分割结果图4 3种算法对NSC序列部分图像的分割对比Fig.4 The results of segmentation and comparison of three algorithms for NSC sequence images

4 结语

本文提出了一种基于优化水平集的细胞图像阈值分割法，可在跟踪细胞边缘的基础上实现对粘连处的分割.多个细胞图像序列的分割实验表明：本文算法相较与传统水平集和阈值法，能更准确的实现对带噪、对比度低的粘连细胞图像的分割.