CN111462099A

CN111462099A - 一种基于快速积分图监测的图像细胞区域定位方法

Info

Publication number: CN111462099A
Application number: CN202010262698.XA
Authority: CN
Inventors: 何昆仑; 赵晓静; 石金龙; 贾志龙; 贾倩; 张诗慧; 刘春祥; 穆欣; 阎岩
Original assignee: Beijing Boco Inter Telecom Technology Co ltd; Chinese PLA General Hospital
Current assignee: Beijing Boco Inter Telecom Technology Co ltd; Chinese PLA General Hospital
Priority date: 2020-04-05
Filing date: 2020-04-05
Publication date: 2020-07-28
Anticipated expiration: 2040-04-05
Also published as: CN111462099B

Abstract

本发明公开了一种基于快速积分图监测的图像细胞区域定位方法，该方法包括：将图像转换为灰度图像；对灰度图像进行去噪滤波，生成滤波后图像；对滤波后图像进行边缘检测，得到边缘图像A；对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，所述自适应阈值根据各像素块区域内的灰度确定；对二值化图像进行开运算滤波，将非细胞部分滤除，并将滤除后的细胞部分进行闭运算填充，获取到图像C；利用一个遍历框分别在图像A、C中进行遍历，计算框在不同位置的密度值，并根据该密度值确定细胞所在图像区域。通过本发明，能够准确的识别出细胞区域。

Description

一种基于快速积分图监测的图像细胞区域定位方法

技术领域

本发明涉及医学图像处理技术领域，尤其涉及一种细胞区域定位或确定方法及装置。

背景技术

目前在细胞图像的分割方法中主要分为传统方法与机器学习(深度学习)方法，在实际应用中采用深度学习速度快，精度高，但是深度学习的效果取决于样本质量，在某些图片质量不佳的情况下需要结合传统方法，传统方法大致采用滤波、全局分割、特征提取等步骤，但是全局分割对于目标与背景相差不大的情况效果不理想，会将很多背景中的杂块也当成目标分割出，影响分割效果。

发明内容

鉴于上述问题，提出了本发明以便提供一种克服上述问题或者至少部分地解决上述问题的技术方案。因此，本发明的一个方面，提供了一种细胞区域定位或确定方法，其特征在于，该方法包括：

将图像转换为灰度图像；

对灰度图像进行去噪滤波，生成滤波后图像；

对滤波后图像进行边缘检测，得到边缘图像A；

对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，所述自适应阈值根据各像素块区域内的灰度确定；

对二值化图像进行开运算滤波，将非细胞部分滤除，并将滤除后的细胞部分进行闭运算填充，获取到图像C；

利用一个遍历框分别在图像A、C中进行遍历，计算框在不同位置的密度值，并根据该密度值确定细胞所在图像区域。

可选的，对滤波后图像进行边缘检测为采用canny算子进行边缘检测。

可选的，对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，包括：

计算当前像素区块内的灰度均值；

根据各像素块的所述灰度均值设定不同的阈值区间；

根据所述灰度均值、阈值区间确定当前像素块内的二值阈值。

可选的，利用积分图计算框在不同位置的密度值。

可选的，根据该密度值确定细胞所在图像区域，包括：筛选满足预定匹配度的目标框所框的区域，确定细胞所在图像区域。

本发明还提供一种细胞区域确定装置，该装置包括：

图像转换单元，用于将图像转换为灰度图像；

滤波单元，用于对灰度图像进行去噪滤波，生成滤波后图像；

边缘检测单元，用于对滤波后图像进行边缘检测，得到边缘图像A；

二值化单元，用于对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，所述自适应阈值根据各像素块区域内的灰度确定；

开运算滤波单元，用于对二值化图像进行开运算滤波，将非细胞部分滤除，并将滤除后的细胞部分进行闭运算填充，获取到图像C；

遍历单元，用于利用一个遍历框分别在图像A、C中进行遍历，计算框在不同位置的密度值，并根据该密度值确定细胞所在图像区域。

可选的，所述边缘检测单元对滤波后图像采用canny算子进行边缘检测。

可选的，所述二值化单元包括：

灰度计算模块，用于计算当前像素区块内的灰度均值；

阈值区间设定模块，用于根据各像素块的所述灰度均值设定不同的阈值区间；

阈值确定模块，用于根据所述灰度均值、阈值区间确定当前像素块内的二值阈值。

可选的，利用积分图计算框在不同位置的密度值。

本申请实施例中提供的技术方案，至少具有如下技术效果或优点：

本文提出的方法可定位出细胞的部分，在此基础上进行分割，可避免分割出背景杂块，提高分割正确率。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述技术方案和其目的、特征和优点能够更明显易懂，以下特举本发明的具体实施方式。

附图说明

通过阅读下文优选实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：

图1示出了本发明提出的细胞区域确定方法的流程图；

图2示出了积分图像示意图；

图3示出了积分示意图；

具体实施方式

下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。

本发明提供一种细胞区域确定方法，如图1所示，该方法包括：

S1.将图像转换为灰度图像；

S2.对灰度图像进行去噪滤波，生成滤波后图像；

S31.对滤波后图像进行边缘检测，得到边缘图像A；

S32.对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，所述自适应阈值根据各像素块区域内的灰度确定；

S33.对二值化图像进行开运算滤波，将非细胞部分滤除，并将滤除后的细胞部分进行闭运算填充，获取到图像C；

S4.利用一个遍历框分别在图像A、C中进行遍历，计算框在不同位置的密度值，并根据该密度值确定细胞所在图像区域。

本发明先定位出细胞部分，在此基础上进行分割，相对于现有技术中使用的全局分割方式，能够避免分割出背景杂块，极大地提高分割正确率。在本发明中，将去噪滤波后的图像分成了两份，其中一份采用自适应阈值进行局部二值化，得到二值化图像B，自适应阈值求得的方法过程如下：

计算当前像素区块内的灰度均值；

根据各像素块的所述灰度均值设定不同的阈值区间；

先通过求当前像素领域内的灰度均值，根据均值确定阈值系数(根据不同的灰度级设定不同的阈值区间)，因此每一块领域内阈值系数均与当前领域内的灰度有关。此领域与图像块的概念相同，为以当前点为中心的n*n正方形，n为经验值，根据当前图像大小及图像中目标物体所占像素数而定，理论上16*16、32*32，64*64均可，为一经验值(可以设定为当前图像大小的1/8或1/16)，也可设定为多个值集合，进行多次取值。

利用灰度均值与阈值系数的乘积来确定当前的二值阈值，此种方法中，阈值系数跟随当前领域内的灰度均值进行调整，较传统的固定系数的方式可以很好的处理有复杂背景、阴影、曝光等情况的图像。

对上述得到的二值化图像B采用复合形态学开运算滤波，采用开运算将B的非细胞小块进行滤波，然后进行闭运算使细胞部分“充实”得到图像C。开运算是通过先对图像腐蚀再膨胀实现，在图像处理上应用很多。腐蚀可以使目标区域范围“变小”，其实质造成图像的边界收缩，可以用来消除小且无意义的目标物。在腐蚀过程中，用数结构F腐蚀图像E，需要注意的是结构F中需要定义一个原点，而且F的移动的过程与卷积核移动的过程一致，同卷积核与图像有重叠之后再计算一样。当F的原点平移到图像E的像元(x,y)时，如果F在(x,y)处，完全被包含在图像E重叠的区域(也就是F中为1的元素位置上对应的E图像值全部也为1)，则将输出图像对应的像元(x,y)赋值为1，否则赋值为0。

F依顺序在E上移动(和卷积核在图像上移动一样，然后在E的覆盖域上进行形态学运算)，当其覆盖F的区域为[1,1；1,1]或者[1,0；1,1]时，(也就是F中‘1’是覆盖区域的子集)对应输出图像的位置才会为1。

而膨胀会使目标区域范围“变大”，将于目标区域接触的背景点合并到该目标物中，使目标边界向外部扩张。作用就是可以用来填补目标区域中某些空洞以及消除包含在目标区域中的小颗粒噪声。

用结构元素F膨胀图像E，将结构元素F的原点平移到图像像元(x,y)位置。如果E在图像像元(x,y)处与F的交集不为空(也就是F中为1的元素位置上对应E的图像值至少有一个为1)，则输出图像对应的像元(x,y)赋值为1，否则赋值为0。

无论腐蚀还是膨胀，都是把结构元素像卷积操作那样，在图像上平移，结构元素中的原点就相当于卷积核的核中心，结果也是存储在核中心对应位置的元素上。只不过腐蚀是结构元素被完全包含在其所覆盖的区域，而膨胀时结构元素与其所覆盖的区域有交集即可。

通过开运算能够将二值化图像B的非细胞小块进行滤波，然后进行闭运算能够使细胞部分“充实”得到图像C。闭运算也是由腐蚀和膨胀结合形成，闭运算就是先膨胀再腐蚀，闭运算在图像处理上已经有所应用，在此不再进行赘述。通过闭运算能够对通过开运算得到的细胞部分进行充实，得到图像C。

作为一种优选实施方式，对滤波后图像进行边缘检测为采用canny算子进行边缘检测。本来，Canny边缘检测是从不同视觉对象中提取有用的结构信息并大大减少要处理的数据量的一种技术，目前已广泛应用于各种计算机视觉系统。但在其他图像处理领域应用还很少。在本发明中，需要获取细胞精确的边缘，边缘检测要达到如下技术效果：以低的错误率检测边缘，也即意味着需要尽可能准确的捕获图像中尽可能多的边缘；检测到的边缘应精确定位在真实边缘的中心；图像中给定的边缘应只被标记一次，并且在可能的情况下，图像的噪声不应产生假的边缘。Canny边缘检测算法可以分为以下5个步骤：

1)使用高斯滤波器，以平滑图像，滤除噪声。

2)计算图像中每个像素点的梯度强度和方向。

3)应用非极大值(Non-Maximum Suppression)抑制，以消除边缘检测带来的杂散响应。

4)应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。

5)通过抑制孤立的弱边缘最终完成边缘检测。

在1)中，为了尽可能减少噪声对边缘检测结果的影响，所以必须滤除噪声以防止由噪声引起的错误检测。为了平滑图像，使用高斯滤波器与图像进行卷积，该步骤将平滑图像，以减少边缘检测器上明显的噪声影响。大小为(2k+1)x(2k+1)的高斯滤波器核的生成方程式由下式给出：

下面是一个sigma＝1.4，尺寸为3x3的高斯卷积核的例子(需要注意归一化)：

若图像中一个3x3的窗口为A，要滤波的像素点为e，则经过高斯滤波之后，像素点e的亮度值为：

其中*为卷积符号，sum表示矩阵中所有元素相加求和。

高斯卷积核大小的选择将影响Canny检测器的性能。尺寸越大，检测器对噪声的敏感度越低，但是边缘检测的定位误差也将略有增加。一般5x5是一个比较不错的trade off。

在2)计算梯度强度和方向中，图像中的边缘可以指向各个方向，因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。边缘检测的算子(如Roberts，Prewitt，Sobel等)返回水平Gx和垂直Gy方向的一阶导数值，由此便可以确定像素点的梯度G和方向theta。

其中G为梯度强度，theta表示梯度方向，arctan为反正切函数。下面以Sobel算子为例讲述如何计算梯度强度和方向。

x和y方向的Sobel算子分别为：

其中Sx表示x方向的Sobel算子，用于检测y方向的边缘；Sy表示y方向的Sobel算子，用于检测x方向的边缘(边缘方向和梯度方向垂直)。

若图像中一个3x3的窗口为A，要计算梯度的像素点为e，则和Sobel算子进行卷积之后，像素点e在x和y方向的梯度值分别为：

其中*为卷积符号，sum表示矩阵中所有元素相加求和。根据公式(3-2)便可以计算出像素点e的梯度和方向。

在3)中非极大值抑制

非极大值抑制是一种边缘稀疏技术，非极大值抑制的作用在于“瘦”边。对图像进行梯度计算后，仅仅基于梯度值提取的边缘仍然很模糊。为了图像中给定的边缘应只被标记一次，对边缘有且应当只有一个准确的响应。而非极大值抑制则可以帮助将局部最大值之外的所有梯度值抑制为0，对梯度图像中每个像素进行非极大值抑制的算法是：

1]将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。

2]如果当前像素的梯度强度与另外两个像素相比最大，则该像素点保留为边缘点，否则该像素点将被抑制。

通常为了更加精确的计算，在跨越梯度方向的两个相邻像素之间使用线性插值来得到要比较的像素梯度。像素点P的梯度方向为theta，则像素点P1和P2的梯度线性插值为：

tan(θ)＝G_y/G_x

G_p1＝(1-tan(θ))×E+tan(θ)×NE

G_p2＝(1-tan(θ))×W+tan(θ)×SW

因此非极大值抑制的伪代码描写如下：

ifG_p≥G_p1andG_p≥G_p2

G_pmay be an edge

else

G_p should be sup pressed

需要注意的是，如何标志方向并不重要，重要的是梯度方向的计算要和梯度算子的选取保持一致。

关于4)双阈值检测，在施加非极大值抑制之后，剩余的像素可以更准确地表示图像中的实际边缘。然而，仍然存在由于噪声和颜色变化引起的一些边缘像素。为了解决这些杂散响应，必须用弱梯度值过滤边缘像素，并保留具有高梯度值的边缘像素，可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值，则将其标记为强边缘像素；如果边缘像素的梯度值小于高阈值并且大于低阈值，则将其标记为弱边缘像素；如果边缘像素的梯度值小于低阈值，则会被抑制。阈值的选择取决于给定输入图像的内容。

双阈值检测的伪代码描写如下：

if G_p≥HighThreshold

G_p is an strong edge

else if G_p≥LowThreshold

G_p is an weak edge

else

G_p should be sup pressed

关于5)抑制孤立低阈值点，到目前为止，被划分为强边缘的像素点已经被确定为边缘，因为它们是从图像中的真实边缘中提取出来的。然而，对于弱边缘像素，将会有一些争论，因为这些像素可以从真实边缘提取也可以是因噪声或颜色变化引起的。为了获得准确的结果，应该抑制由后者引起的弱边缘。通常，由真实边缘引起的弱边缘像素将连接到强边缘像素，而噪声响应未连接。为了跟踪边缘连接，通过查看弱边缘像素及其8个邻域像素，只要其中一个为强边缘像素，则该弱边缘点就可以保留为真实的边缘。

抑制孤立边缘点的伪代码描述如下：

if G_p＝＝LowThreshold and G_p connected to a strong edge pixel

G_p is an strong edge

else

G_p should be sup pressed

通过以上5个步骤即可完成基于Canny算法的边缘提取。

本发明发现，在细胞区域的边缘图像中边缘密度与二值化图像中的二值密度在一定的范围内，采用一个遍历框分别在A与C中遍历图像，采用快速积分图计算框在不同位置的密度值，求得A与C中各自的“目标框”。利用积分图计算框在不同位置的密度值。

图像是由一系列的离散像素点组成，因此图像的积分其实就是求和.图像积分图中每个点的值是原图像中该点左上角的所有像素值之和.

首先建立一个数组A作为积分图像，其宽高与原图像相等.然后对这个数组赋值，每个点存储的是该点与图像原点所构成的矩形中所有像素的和：

SAT(x，y)＝∑I(xi，yi) (1)

其中I(x，y)表示图像(x，y)位置的像素值。

积分图像可以采用增量的方式计算：

SAT(x,y)＝SAT(x,y-1)+SAT(x-1,y)-SAT(x-1,y-1)+I(x,y) (2)

初始边界：SAT(-1,y)＝SAT(x,-1)＝SAT(-1,-1)＝0

SAT(x,y-1)+SAT(x,y-1)后，有一部分重合的区域，即SAT(x-1,y-1)，所以需减掉，最后还需要将当前坐标(x,y)的像素值I(x,y)I(x,y)包含进来。

定义了积分图的概念，就可以很方便的计算任意区域内的像素和，如图2所示。

积分图数组初始化之后,我们就得到了一张积分图，如图3所示:

点1的积分SAT1＝Sum(Ra)，

点2的积分SAT2＝Sum(Ra)+Sum(Rb)，

点3的积分SAT3＝Sum(Ra)+Sum(Rc)，

点4的积分SAT4＝Sum(Ra)+Sum(Rb)+Sum(Rc)+Sum(Rd)

那么为了计算某个矩形像素和，比如区域Rd内所有点的像素值之和(积分)可以表示为：

Sum(Rd)＝SAT1+SAT4-SAT2-SAT3 (3)

所以无论矩形的尺寸大小，只需查找积分图像4次就可以快速计算任意矩形内像素值的和,即算法复杂度为O(4)。在本发明中在获取到细胞区域后，对每个细胞区域计算其区域内像素和(积分)，再构建区域密度值，因此对一个灰度图而言，事先将其积分图构建好，当需要计算灰度图某个区域内所有像素点的像素值之和的时候，利用积分图，通过查表运算，可以迅速得到图像A，C中的各自的“目标框”。

根据该密度值确定细胞所在图像区域，包括：筛选满足预定匹配度的目标框所框的区域，确定细胞所在图像区域。分别计算A中目标框与C中目标框的IOU,＝(A交C)/(A并C)),设置阈值(经验值例如0.7)，大于0.7的两个框位置进行合并,合并的原则如下

设A的左上角坐标x1A,y1A,右下角坐标x2A,y2A

C的左上角坐标x1C,y1C,右下角坐标x2C,y2C

合并框左上角坐标(x1,y1),右下角坐标(x2,y2)

(x1,y1)＝(min(x1A,x1C),min(y1A,y1C))

(x2,y2)＝(max(x2A,x2C),max(y2A,y2C))

本发明还提供一种细胞区域确定装置，该装置包括：

图像转换单元，用于将图像转换为灰度图像；

可选的，所述二值化单元包括：

灰度计算模块，用于计算当前像素区块内的灰度均值；

可选的，利用积分图计算框在不同位置的密度值。

在此处所提供的说明书中，说明了大量具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

类似地，应当理解，为了精简本发明并帮助理解各个发明方面中的一个或多个，在上面对本发明的示例性实施例的描述中，本发明的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本发明要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，发明方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本发明的单独实施例。

应该注意的是上述实施例对本发明进行说明而不是对本发明进行限制，并且本领域技术人员在不脱离所附权利要求的范围的情况下可设计出替换实施例。在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。

Claims

1.一种细胞区域定位方法，其特征在于，该方法包括：

将图像转换为灰度图像；

对灰度图像进行去噪滤波，生成滤波后图像；

对滤波后图像进行边缘检测，得到边缘图像A；

2.根据权利要求1所述的方法，其特征还在于，对滤波后图像进行边缘检测为采用canny算子进行边缘检测。

3.根据权利要求1所述的方法，其特征还在于，对滤波后图像利用自适应阈值进行二值化，得到二值化图像B，包括：

计算当前像素区块内的灰度均值；

根据各像素块的所述灰度均值设定不同的阈值区间；

4.根据权利要求1所述的方法，其特征还在于，利用积分图计算框在不同位置的密度值。

5.根据权利要求1所述的方法，其特征还在于，根据该密度值确定细胞所在图像区域，包括：筛选满足预定匹配度的目标框所框的区域，确定细胞所在图像区域。

6.一种细胞区域定位装置，其特征在于，该装置包括：

图像转换单元，用于将图像转换为灰度图像；

滤波单元，用于对灰度图像进行去噪滤波，生成滤波后图像；边缘检测单元，用于对滤波后图像进行边缘检测，得到边缘图像A；

7.根据权利要求6所述的装置，其特征还在于，所述边缘检测单元对滤波后图像采用canny算子进行边缘检测。

8.根据权利要求6所述的装置，其特征还在于，所述二值化单元包括：灰度计算模块，用于计算当前像素区块内的灰度均值；

9.根据权利要求6所述的装置，其特征还在于，利用积分图计算框在不同位置的密度值。

10.根据权利要求6所述的装置，其特征还在于，根据该密度值确定细胞所在图像区域，包括：筛选满足预定匹配度的目标框所框的区域，确定细胞所在图像区域。