CN114648511A - 一种大肠杆菌轮廓的精确提取及识别方法 - Google Patents
一种大肠杆菌轮廓的精确提取及识别方法 Download PDFInfo
- Publication number
- CN114648511A CN114648511A CN202210313888.9A CN202210313888A CN114648511A CN 114648511 A CN114648511 A CN 114648511A CN 202210313888 A CN202210313888 A CN 202210313888A CN 114648511 A CN114648511 A CN 114648511A
- Authority
- CN
- China
- Prior art keywords
- image
- escherichia coli
- point
- contour
- pixel
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 241000588724 Escherichia coli Species 0.000 title claims abstract description 184
- 238000000034 method Methods 0.000 title claims abstract description 79
- 238000000605 extraction Methods 0.000 title claims description 14
- 230000008569 process Effects 0.000 claims abstract description 45
- 238000012545 processing Methods 0.000 claims abstract description 27
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 13
- 230000000877 morphologic effect Effects 0.000 claims description 64
- 230000011218 segmentation Effects 0.000 claims description 59
- 239000013598 vector Substances 0.000 claims description 30
- 238000004364 calculation method Methods 0.000 claims description 18
- 238000001914 filtration Methods 0.000 claims description 18
- 230000005764 inhibitory process Effects 0.000 claims description 18
- 230000008859 change Effects 0.000 claims description 17
- 238000003705 background correction Methods 0.000 claims description 16
- 238000004458 analytical method Methods 0.000 claims description 15
- 230000002146 bilateral effect Effects 0.000 claims description 13
- 238000011049 filling Methods 0.000 claims description 13
- 230000009466 transformation Effects 0.000 claims description 11
- 238000001514 detection method Methods 0.000 claims description 10
- 239000012535 impurity Substances 0.000 claims description 9
- 239000003242 anti bacterial agent Substances 0.000 claims description 7
- 229940088710 antibiotic agent Drugs 0.000 claims description 7
- 230000003115 biocidal effect Effects 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 5
- 239000000853 adhesive Substances 0.000 claims description 4
- 230000001070 adhesive effect Effects 0.000 claims description 4
- 229940079593 drug Drugs 0.000 claims description 4
- 239000003814 drug Substances 0.000 claims description 4
- 230000002401 inhibitory effect Effects 0.000 claims description 4
- 238000004891 communication Methods 0.000 claims description 3
- 230000009977 dual effect Effects 0.000 claims description 3
- 238000003708 edge detection Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000005530 etching Methods 0.000 claims description 2
- 238000013461 design Methods 0.000 abstract description 8
- 238000002474 experimental method Methods 0.000 abstract description 8
- 238000010191 image analysis Methods 0.000 abstract description 7
- 238000005516 engineering process Methods 0.000 abstract description 5
- 230000002829 reductive effect Effects 0.000 abstract description 5
- 238000012896 Statistical algorithm Methods 0.000 abstract description 2
- 230000009471 action Effects 0.000 description 5
- 244000005700 microbiome Species 0.000 description 4
- 230000036961 partial effect Effects 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 239000011248 coating agent Substances 0.000 description 3
- 238000000576 coating method Methods 0.000 description 3
- 239000003086 colorant Substances 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000000813 microbial effect Effects 0.000 description 3
- 238000002360 preparation method Methods 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 230000001580 bacterial effect Effects 0.000 description 2
- 230000001413 cellular effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 206010059866 Drug resistance Diseases 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013527 convolutional neural network Methods 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 230000002906 microbiologic effect Effects 0.000 description 1
- 238000000399 optical microscopy Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/045—Combinations of networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4007—Scaling of whole images or parts thereof, e.g. expanding or contracting based on interpolation, e.g. bilinear interpolation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/12—Edge-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10056—Microscopic image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
- G06T2207/20028—Bilateral filtering
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Software Systems (AREA)
- Evolutionary Computation (AREA)
- Computational Linguistics (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Geometry (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种大肠杆菌轮廓的精确提取及识别方法；本发明依靠计算机视觉的图像处理技术，结合大肠杆菌的几何特性，针对性地设计了一套图像处理算法，深度结合了实验所得出的数据，提高不同场景下系统的个性化和适配性。并利用统计学的算法，从培养皿复杂的大肠杆菌形态图像中将大肠杆菌形态的轮廓特征提取出来，从而方便了使用者统计信息，在大量的图像分析中减少了使用者的操作负担，加快实验的进度。此外，本发明设计一种轮廓亚像素重建算法，在不改变硬件的前提下，提高图像的分辨率，使系统的轮廓变得更加细致，从而能够精确地提取轮廓。同时本发明采用自适应双阈值二值化来处理二值化过程中的手动调参，使系统更加智能和高效。
Description
技术领域
本发明属于生物医学图像处理领域，主要涉及一种大肠杆菌轮廓的精确提取及识别方法。
背景技术
细胞形态和亚细胞空间结构对微生物细胞的功能有着重大影响，因此通过对微生物形态的研究具有非常重要的理论意义和实际应用价值。且长期以来，光学显微镜一直是表征和理解微生物世界的重要技术，广泛用于测量活细胞和群落内部的空间结构，这通常会导致大量的或难以手动分析的图像。在过去，微生物的标本检查等工作是由人工完成的，存在检测效率低、精度有限等一系列的问题，当需要分析数百，数千甚至更多的图像时，手动分析变得根本不可能。而利用计算机对细胞图像进行自动判读可以有效地避免这些问题，因而也成为国内外学者研究的一个热点。随着基于计算机图像分析的快速发展，使得传统图像处理技术和卷积神经网络都能够量化单个细胞和群落的大量属性。由于神经网络需要大量的样本数据和训练时间且泛化能力弱。所以目前大多数应用于微生物学的图像分析工具都依赖于传统的图像分析技术，传统的图像处理技术构成了目前大多数微生物学图像分析工具的基础。然而在传统微生物图像处理领域中存在微生物形态轮廓提取不够精确，轮廓特征较少导致识别的准确率低的问题。形态轮廓和特征提取对实验分析至关重要，但也是最困难、最具有挑战性的步骤，这是一个非常重要的环节，处理结果的好坏直接关系到后续识别的效果。由于细胞图像的多样性和复杂性，目前还没有一种完全通用的处理方法可以实现对所有细胞图像的分析。通常要与该领域知识充分结合适配，才能更为有效地解决问题。
在生物医学领域的研究实验中，为了保证结果的准确性，通常需要大量的重复实验，而在这个过程中就会产生大量的细胞形态图像，手动分析必然会带来低效和误判。所以具有一套该专业领域的图像处理算法对细胞形态图像分析至关重要。
鉴于此，将图像处理技术应用到微生物领域，针对性的设计了一套图像处理算法，解决了形态的完整性、形态分辨率的精度、充分多的形态特征。三者构成了大肠杆菌精确轮廓提取及识别算法的充要条件，实现了大肠杆菌精确轮廓提取及识别，实现了实验中能够精确快速的对抗生素作用下大肠杆菌的抑制率的分析，从而有利于提高科研工作者的效率和准确性。
发明内容
本发明的目的在于提供一种大肠杆菌轮廓的精确提取及识别方法。以解决现有技术中黏连分割效果不佳、细菌形态轮廓分辨率精度不够、分类识别判断条件不充分、提取和识别准确度不够，以及因图像处理技术的专业性而导致在分析实验数据时不能充分结合适配的技术问题。
该大肠杆菌轮廓的精确提取及识别方法，具体步骤如下：
步骤1、采集被识别的大肠杆菌形态图像；
步骤2、对步骤1所得的大肠杆菌形态图像进行预处理，得到预处理图像。
步骤3、对步骤2所得的预处理图像进行自适应双阈值二值化，得到二值图像。
步骤4、对步骤3所得的二值图像进行区域填充，得到填充图像。
步骤5、对步骤4的填充图像进行两次形态学开运算，得到菌落开运算图像。
步骤6、对步骤5所得的开运算图像进行黏连分割，得到黏连分割图像。
步骤7、对步骤6所得的黏连分割图像进行轮廓绘制并进行计数，得到轮廓图像。
步骤8、将步骤7所得的轮廓图像进行亚像素重建，得到轮廓亚像素重建图像。
步骤9、对步骤8所得的轮廓亚像素重建图像进行轮廓分割，得到轮廓分割图像。
步骤10、提取大肠杆菌的轮廓特征。
步骤11、根据步骤10得到的大肠杆菌的轮廓特征，进行大肠杆菌形态分析。
作为优选，被测的大肠杆菌培养皿中投加有抗生素；每隔预设时长重复执行一次步骤1至11，得到在不同时间段下，大肠杆菌受到抗生素的抑制效果。
作为优选，步骤1的具体过程为：
1-1.向大肠杆菌培养皿中添加抗生素。
1-2.用透射电子显微镜对培养皿中的大肠杆菌进行形态图像采集，得到三通道的大肠杆菌形态图像。
1-3.将步骤1-2所得的大肠杆菌形态图像传输到计算机上。
作为优选，步骤2中所述预处理的具体过程为：
2-1.对步骤1所得的大肠杆菌形态图像使用分量权重法进行灰度化，其中三通道形态图像得到单通道形态灰度图像。
2-2.对步骤2-1所得的形态灰度图像进行背景校正，得到形态背景校正图像，背景校正的具体过程为：
2-2-1.选取3*3的卷积核模板，该卷积核以中心位置作为坐标原点。
2-2-2.使用步骤2-2-1所得的3*3卷积核对步骤2-1所得的大肠杆菌形态灰度图像做卷积，得到轻模糊处理图像o(x,y)。
2-2-3.选取25*25的卷积核模板，对步骤2-1所得的大肠杆菌形态灰度图像进行高斯滤波，得到强模糊处理图像b(x,y)。并求取强模糊处理图像b(x,y)的平均灰度M。
2-2-4.对步骤2-1所得的大肠杆菌灰度图像进行背景校正，得到形态背景校正图像f(x,y)如下：
f(x,y)＝o(x,y)-b(x,y)+M
2-3.对步骤2-2所得的形态背景校正图像进行双边滤波，得到大肠杆菌形态双边滤波灰度图像。双边滤波方法具体为将空间域核与颜色域核叠加得到双核。用双核对被识别图像进行滤波。
2-4.对步骤2-3所得的大肠杆菌形态双边滤波灰度图像进行直方图均衡化增强，得到大肠杆菌形态灰度增强图像。
2-5.对步骤2-4所得的大肠杆菌形态灰度增强图像进行形态学梯度变换。得到大肠杆菌形态学梯度变换图像g如下：
式(2)中，f表示大肠杆菌形态灰度增强图像，B表示3*3的交叉形内核，
表示膨胀操作，
表示腐蚀操作。大肠杆菌形态学梯度变换图像g即为步骤2得到的预处理图像。
作为优选，步骤3中所述自适应双阈值二值化的具体过程为：
3-1.对步骤2所得的预处理图像进行灰度梯度计算，梯度计算包括计算梯度方向与幅值，结合Sobel算子，计算出四个方向上的梯度变化强度，再通过线性插值计算两个对角线方向上的梯度变化强度。
3-2.将当前像素梯度变化强度与沿正负、对角线梯度方向上的相邻像素的梯度强度进行比较，若一个像素点为比较中的最大值，则保留该像素为边缘点，像素值设为255。若不是最大，则对其进行抑制，不将其作为边缘点。像素值设为0。
3-3.通过Otsu算法，将梯度图像非真实边缘部分划分为伪边缘部分和非边缘部分。将迭代分割运算中取得最大类间方差时的分割阈值作为最优分割阈值，抑制孤立的弱边缘，完成边缘检测。类间方差的表达式如下：
其中，w0为目标点数占总图像比例，u0为其平均灰度；w1为背景点数占图像比例，u1为其平均灰度；u为图像总平均灰度；g′为类间方差。
3-4.取类间方差值取最大时对应的阈值作为最优分割阈值，并将其定义为低阈值；将低阈值的2.5倍定义为高阈值。梯度强度低于低阈值的像素点被抑制，不作为边缘点，像素值设为0；高于高阈值的像素点被定义为强边缘，保留为边缘点，像素值设为255；处于高低阈值之间的定义为弱边缘；若弱边缘的邻域像素中有至少一个为强边缘像素点，则将该弱边缘保留为强边缘，像素值设为255；否则，该弱边缘不作为边缘点，像素值设为0。
作为优选，步骤4中进行区域填充的具体过程为：
4-1.对步骤3所得的二值图像中的连通域之外的像素值全部赋值为0，再将所得图像进行按位取反处理。
4-2.对步骤4-1所得的图像与步骤3所得的二值图像进行按位取或处理，得到填充图像。
作为优选，步骤6中所述黏连分割的具体过程为：
6-1.对黏连分割图像进行八连通链码标记。链码标记具体过程为：
式(4)中，j为链码序号；k为链码步进长度；ζj,k为链码第j个点的链码值；(xj,yj)、(xj-k,yj-k)分别表示链码第j、j-k个点坐标。
6-1-1.以起始点开始，并记录为起始点坐标，以8连通的方式从编码为1的方向开始搜索下一个边界像素点；之后，记录方向编码。
6-1-2.从找到的像素开始重复步骤6-1-1，得出一组链码。若该组链码闭合，则在编码完成后返回到起始点；否则，省略起始点坐标。
6-1-3.对步骤6-1-2所得的链码进行归一化，得归一化链码。将链码看作n位自然数，将链码按一个方向循环，使其构成的n位自然数最小，形成起点唯一的链码。
6-2.确定连通域边界上的凹点。
6-2-1.对步骤6-1所得的图像边界链码，再根据向量叉乘运算公式进行叉乘运算，所得向量叉乘运算结果存放在集合Se中。定义一个阈值Th。若向量叉乘运算结果大于阈值Th，则该两个向量对应的像素点满足凹点特性条件。
向量叉乘运算公式如下：
s＝a×b＝|a||b|sin(θ)
其中，s为向量叉乘运算结果；a表示以链码的起始坐标为起点，以步长k为距离所做的向量；b表示以向量a的终点作为起点，以步长k为距离所做的向量。θ为向量a到向量b的夹角。
6-2-2.对步骤6-2-1所得的集合Se中满足凹点特性条件的像素点进行标记，得到多个凹点待定集合；每个凹点待定集合均包括一系列位于凹点周围的点；再将每个凹点待定集合中绝对值最大的那个点作为该集合的凹点，得到图像中所有的凹点。
6-3.对步骤6-2所得的连通域边界上的凹点进行匹配分割，凹点匹配分割的具体过程为：
6-3-1.将凹点分类为内凹点和外凹点，内凹点表示沿着图像的内边缘依次寻找得到的凹点；外凹点表示沿着图像的外边缘依次寻找得到的凹点。
6-3-2.若连通域边界上的外凹点个数小于1，则不进行分割。
6-3-3.若连通域边界上的外凹点个数为2，则通过将该两个外凹点直线相连的方式进行分割。
6-3-4.若连通域边界上的外凹点个数大于2。则将各外凹点两两之间配对；若两个外凹点满足欧式距离最小且外凹点的分割线在同方向上，则该两个外凹点直接配对。将各组配对的外凹点直线相连，完成分割。
作为优选，步骤8中进行轮廓亚像素重建的具体步骤为：
8-1.对步骤7所得的轮廓图像进行灰度化。
8-2.灰度化后的轮廓图像进行梯度检测，梯度检测包含幅值和方向，设置8个3*3方向模板，用8个方向模板在图像上进行卷积运算，得到8个梯度图像。模板的方向就是梯度方向和X轴的夹角。
梯度检测使用的梯度图像变换公式如下：
R(i,j)＝∣f(i+1,j+1)+2f(i,j+1)+f(i-1,j+1)-f(i-1,j-1)-2f(i,j-1)-f(i+1,j-1)∣+∣f(i-1,j-1)+2f(i-1,j)+f(i-1,j+1)-f(i+1,j+1)-2f(i+1,j)-f(i+1,j-1)∣
其中，R(i,j)为坐标(i,j)上的梯度幅值。f(i-1,j-1)、f(i,j-1)、f(i+1,j-1)、f(i-1,j)、f(i+1,j)、f(i-1,j+1)、f(i,j+1)、f(i+1,j+1)分别为坐标(i,j)左上、上、右上、左、右、左下、下、右下八个方向上的灰度值。
8-3.通过梯度方向上的变化规律，在3*3的邻域中根据梯度的方向比较中心点与该方向相邻点的幅值，幅值大的为边缘点，进行非极大值抑制。得到只有一个像素宽度的单像素边缘。
8-4.对步骤8-3所得的单像素边缘进行多项式插值，得到亚像素边缘。多项式插值的公式如下：
其中，(Xc，Yc)表示亚像素坐标；(m，n)表示像素级坐标；W为相邻像素点到边缘点的距离；θ为梯度方向与X轴正向夹角；R0为边缘点(m,n)的灰度梯度的模，R-1、R1分别是在梯度方向上与边缘点相邻的两像素点的梯度幅值。
8-5.对步骤8-4得到的边缘亚像素坐标进行存储，创建两个通道的图像，第一个通道表示x，第二个通道表示y，x通道存储像素值，y通道存储的是亚像素坐标。
8-6.对步骤8-5所得的各个亚像素轮廓分别进行绘制连接。
作为优选，步骤9中进行轮廓分割的具体流程为：
9-1.对步骤8所得的轮廓亚像素重建图像进行链码标记。
9-2.对步骤9-1中的链码进行一阶差分链码计算以及一阶差分链码的归一化。
9-3.依据拐点在切线处不连续，曲率连续的特点，对步骤8所得的轮廓亚像素重建图像进行拐点的提取。
9-4.对步骤8所得的轮廓图像进行切点的提取。
9-5.根据步骤9-3、步骤9-4将拐点和切点全部提取出来，将拐点之间用线段连接，切点之间用最小二乘法进行椭圆拟合，并统计线段和椭圆弧的数量。得到轮廓分割为直线段、椭圆弧段的图像。
9-6.对步骤9-5所得的直线段和椭圆弧段进行绘制。
步骤10中，提取轮廓特征的具体过程为：
10-1.取连通域个数作为大肠杆菌数量。
10-2.取单个连通域像素总数，作为对应的单个大肠杆菌面积。
10-3.取单个连通域的边界像素总数，作为对应的单个大肠杆菌周长。
10-4.取单个连通域的长轴短轴比，作为对应的单个大肠杆菌长短比。
10-5.取单个连通域的区域面积与π乘中心点到所有轮廓像素距离的最大值平方之比，作为对应的单个大肠杆菌圆形度。
10-6.统计单个大肠杆菌的轮廓直线和椭圆弧数量。
步骤11中进行大肠杆菌形态分析的具体流程为：
11-1.若大肠杆菌的面积大于100或周长大于40，则判断大肠杆菌的长短比是否大于3且小于10。
11-2.若否，则判定为杂质，若是，则判断大肠杆菌的轮廓中直线数量是否大于1且椭圆数量是否大于2。
11-3.若是，将其判定为被抑制的大肠杆菌；若否，将其判定为已耐药的大肠杆菌。
11-4.若大肠杆菌的面积小于100且周长小于40，则判断大肠杆菌的圆形度是否大于0.5。
11-5.若是，则判定为被抑制的大肠杆菌；若否，则判定为杂质。
11-6.取被抑制的大肠杆菌与大肠杆菌的总数量的比值，作为大肠杆菌的抑制率。
作为优选，步骤9-2的具体流程为：
9-2-1.计算一阶差分链码δj,k如下：
其中，j为链码序号；k为链码步进长度；δj,k为链码第j个点的一阶差分值；(xj-k,yj-k)、(xj,yj)、(xj+k,yj+k)分别表示链码第j-k、j、j+k个点坐标。
9-2-2.对步骤9-2-1所得的一阶差分链码进行归一化处理，把一阶差分链码看作由各方向数所构成的自然数，再按一个方向循环，使其构成自然数最小。
步骤9-3的具体流程为：
9-3-1.计算差分链码偏离平均值的绝对值
用该绝对值
建立高斯分布，并选取阈值。将偏离均值最大的一部分点作为分段点。绝对值
的表达式如下：
其中，j为链码序号；k为链码步进长度；δj,k为链码第j个点的一阶差分值；μδ表示链码均值。
9-3-2.设定阈值T；将大于阈值T且为邻域最大值的分段点作为拐点。阈值T的表达式如下：
T＝z1-2ασ
其中，z为步骤9-3-1中建立的二维高斯分布的横坐标轴，1-2α为步骤9-3-1中建立的二维高斯分布的置信区间；z1-2α为步骤9-3-1中建立的二维高斯分布的置信区间横坐标值；σ为步骤9-3-1中建立的二维高斯分布的标准方差。
步骤9-4的具体流程为：
9-4-1.对步骤8所得的轮廓图像进行多边形相似逼近，轮廓变为多折线段，得到平滑分裂点。计算轮廓和近似函数之间的最大误差E∞(u1,ui+1)。若达到分段数量或者分裂点最大误差低于设定阈值时停止分裂。最大误差E∞(u1,ui+1)的表达式如下：
其中，ei,j是区间i的第j个点与其近似函数之间的欧几里德距离；(u1,ui+1)是区间i的两个端点。
9-4-2.对步骤9-4-1所得到的分裂点进行分段近似椭圆拟合，选取第一个分裂点为起点，根据相邻轮廓段的三个分裂点不可能处在同一直线上这一前提下，将起点与相邻的两个分裂点作为集合进行近似椭圆拟合，并计算该近似椭圆与该分裂点集合对应轮廓段之间的最大误差。
9-4-3.若最大误差小于该分裂点集合对应轮廓段与多边形相似逼近线段的最大距离，则将该分裂点集合对应轮廓段就标记为合并处理对象，保存标记该分裂点集合，将下一个分裂点作为起点。执行步骤9-4-2的操作。否则，将下一个分裂点加入该分裂点集合中，将分裂点集合进行近似椭圆拟合，计算该近似椭圆与该分裂点集合对应轮廓段之间的最大误差。
9-4-4.重复执行步骤9-4-3，经过多次迭代直到没有合并处理的对象为止。
9-4-5.将标记的每个分裂点集合的起点和终点视为切点。
本发明具有的有益效果是：
1、本发明依靠计算机视觉的图像处理技术，结合大肠杆菌的几何特性(由直线段和圆弧段组成)，针对性地设计了一套图像处理算法，深度结合了实验所得出的数据，提高不同场景下系统的个性化和适配性。并利用统计学的算法，从培养皿复杂的大肠杆菌形态图像中将大肠杆菌形态的轮廓特征提取出来，从而方便了使用者统计信息，在大量的图像分析中减少了使用者的操作负担，加快实验的进度。
2、本发明设计了一种对大肠杆菌形态轮廓做分割的方法，将其轮廓分割为线段、椭圆弧两种形状，作为轮廓的一种特征，增加了大肠杆菌的形态轮廓特征，提高了识别率。结合大肠杆菌抑制前和抑制后的几何特性，即在抑制前为椭圆形，抑制后呈现类似圆形态，通过判断大肠杆菌所包含的椭圆弧和线段的数量以及面积、周长、圆形度、长短比等特征，可以准确度的对被抑制的大肠杆菌、耐药性的大肠杆菌、杂质做区分判断。通常一些不规则的杂质很容易满足形态面积、周长、长短比特征，此时增加的轮廓分割为直线段和椭圆弧这一特征之后，可以有效地避免此类误判。将其应用到大肠杆菌抑制前、抑制后、杂质识别中。增加其判断条件，减少系统误判，提高了识别的准确度，进而提高了抗生素作用下大肠杆菌抑制率的准确度。
3、本发明设计黏连分割、轮廓亚像素、轮廓分割的图像算法从形态的完整性、分辨率、几何特性三个角度去设计优化轮廓提取的精度和系统识别的准度。对于图像的完整性保证了大肠杆菌轮廓的几何特征，为后续识别做铺垫；图像的分辨率保证了大肠杆菌轮廓提取的精确度。提高了轮廓几何特性准度；大肠杆菌的几何特性增加了大肠杆菌的轮廓特征，保证在识别判断中减少系统的误判。三者结合大大地增加了系统的精度和准度。
4、本发明设计针对黏连大肠杆菌的分割算法，通过精确分割来保证大肠杆菌的完整性，进而提高形态轮廓的特征。为了更精确、更快速地分割粘连区域，该算法采用链码标记和夹角运算以及凹点配对分割，实现了分割的准确性；利用多线程的形式，实现并发，相对于循环执行在效率方面得到了改善，提升了运行速率，实现了快速性的要求。同时降低了硬件资源要求。
5、本发明设计一种轮廓亚像素重建算法，在不改变硬件的前提下，提高图像的分辨率，使系统的轮廓变得更加细致，从而能够精确地提取轮廓。同时本发明采用自适应双阈值二值化来处理二值化过程中的手动调参，使系统更加智能和高效。
附图说明
图1为本发明的实施步骤框图；
图2为本发明的大肠杆菌形态分析流程图；
图3为本发明所得的大肠杆菌形态图像的预处理图像；
图4为本发明所得的大肠杆菌形态图像的二值图像；
图5为本发明所得的大肠杆菌形态图像的区域填充图像；
图6为本发明所得的大肠杆菌形态图像的开运算图像；
图7为本发明所得的大肠杆菌形态图像的黏连分割图像；
图8为本发明所得的大肠杆菌形态图像的像素级轮廓图像；
图9为本发明所得的的细菌形态图像的亚像素级轮廓图像；
图10为本发明所得的大肠杆菌形态图像的亚像素级轮廓局部放大对比图；
图11为本发明所得的大肠杆菌形态图像的亚像素级轮廓基元分段图像；
图12为本发明所得的大肠杆菌形态图像的亚像素级轮廓基元分段局部放大对比图。
具体实施方式
以下结合附图对本发明作进一步说明。
如图1所示，一种大肠杆菌轮廓的精确提取及识别方法，用于提取图像中大肠杆菌的轮廓和数量；依据所得的大肠杆菌形态轮廓的变化能够反映出其在抗生素作用下的抑制效果，来进行抗生素作用下大肠杆菌抑制率的分析实验。
如图1所示，该提取及识别方法的具体步骤如下：
步骤1、通过透射电子显微镜获取抗生素作用下的大肠杆菌形态图像。
1-4.向大肠杆菌培养皿中添加低浓度(1至15mg·kg-1)的抗生素，等待大约30分钟。
1-5.用透射电子显微镜对培养皿中的大肠杆菌进行形态图像采集，得到三通道的大肠杆菌形态图像。
1-6.将步骤1-2所得的大肠杆菌形态图像传输到计算机上。
步骤2、对步骤1所得的大肠杆菌形态图像进行预处理，得到预处理的图像，如图3所示。预处理的具体过程为：
2-1.根据下面公式(1)，对步骤1所得的大肠杆菌形态图像进行灰度化，其中三通道形态图像得到单通道形态灰度图像。
Gray(i,j)＝R(i,j)*0.299+G(i,j)*0.587+B(i,j)*0.114 式(1)
式(1)中，Gray(i,j)为灰度化后的灰度图像中第i行和第j列的像素点的灰度值，R(i,j)，G(i,j)，B(i,j)分别为原彩色图像的中第i行和第j列的像素点的红色，绿色，蓝色的分量值，0.299，0.587，0.114分别为各像素点三分量的权值。
2-2.对步骤2-1所得的形态灰度图像进行背景校正，得到形态背景校正图像，背景校正的具体过程为：
2-2-1.选取3*3的卷积核模板，该卷积核以中心位置作为坐标原点，将卷积核各点坐标位置代入下面二维高斯函数公式(2)中，得到3*3卷积核模板系数。
式(2)，中σ是标准方差，决定了高斯滤波器的宽度。(x,y)表示坐标，G(x,y)表示坐标(x,y)处的灰度值，e表示自然对数函数的底数。
2-2-2.使用步骤2-2-1所得的3*3卷积核模板系数对步骤2-1所得的大肠杆菌形态灰度图像做卷积(模板滑过的地方的像素值和模板中对应的值相乘，最后全部加总)，得到高斯滤波图像，作为轻模糊处理图像。得到的图像灰度用o(x,y)表示。
2-2-3.选取25*25的卷积核模板，根据步骤2-2-1、步骤2-2-2做高斯滤波，得到高斯滤波图像，作为强模糊处理图像。得到的图像灰度用b(x,y)表示，再求其平均灰度，并用M表示。
2-2-4.根据下面公式(3)，对步骤2-1所得的大肠杆菌灰度图像进行背景校正，得到大肠杆菌背景校正的灰度图像。
f(x,y)＝o(x,y)-b(x,y)+M 式(3)
2-3.对步骤2-2所得的大肠杆菌背景校正的灰度图像进行双边滤波，突现边缘细节。得到大肠杆菌形态双边滤波灰度图像。
2-3-1.非线性双边滤波结合图像的邻近度和像素值相似度，滤除噪声的同时可以保留原图的边缘信息。双边滤波方法将空间域核与颜色域核叠加得到双核。用双核对被识别图像进行滤波。
2-4.对步骤2-3所得的大肠杆菌形态双边滤波灰度图像进行直方图均衡化增强，通过改变图像的直方图来改变图像中各像素的灰度，得到大肠杆菌形态灰度增强图像。
2-5.对步骤2-4所得的大肠杆菌形态灰度增强图像和下面灰度级形态学梯度公式(4)，进行形态学梯度变换。得到大肠杆菌形态学梯度变换图像。
式(4)中，g表示梯度图像，f表示灰度图像，B表示3*3的交叉形内核，
表示膨胀，
表示腐蚀。
步骤3、对步骤2所得的预处理图像进行自适应双阈值二值化，得到二值图像，如图4所示，自适应双阈值二值化的具体过程为：
3-1.对步骤2所得的预处理图像进行灰度梯度计算，梯度计算包括计算梯度方向与幅值，结合Sobel算子，计算出四个方向上的梯度变化强度，再通过线性插值计算两个对角线方向上的梯度变化强度。
3-2.将当前像素梯度变化强度与沿正负、对角线梯度方向上的相邻像素的梯度强度进行比较，若为最大(即为极值)，则保留该像素为边缘点，像素值设为255。若不是最大，则对其进行抑制，不将其作为边缘点。像素值设为0。
3-3.通过Otsu法，将梯度图像非真实边缘部分，分为伪边缘部分和非边缘部分。当迭代分割运算取得最大类间方差时的分割阈值即为最优分割阈值。抑制孤立的弱边缘完成边缘检测。根据下面公式(5)求最大类间方差：
式(5)中，w0为目标点数占总图像比例，u0为其平均灰度；w1为背景点数占图像比例，u1为其平均灰度；u为图像总平均灰度；g′为类间方差。
3-4.类间方差值取最大时对应的阈值为最优分割阈值，将其定义为低阈值，低阈值的2.5倍定义为高阈值。梯度强度低于低阈值的像素点被抑制，不作为边缘点，像素值设为0；高于高阈值的像素点被定义为强边缘，保留为边缘点，像素值设为255；处于高低阈值之间的定义为弱边缘，只要其邻域像素中有一个为强边缘像素点，则该弱边缘可以保留为强边缘，即真实边缘，像素值设为255，否则不作为边缘点，像素值设为0。
步骤4、对步骤3所得的二值图像进行区域填充，得到填充图像，如图5所示。区域填充的具体过程为：
4-1.对步骤3所得的二值图像中的连通域之外的像素值全部赋值为0，再将像素值按位取反。
4-2.对步骤4-1所得的图像和步骤3所得的二值图像按位或，得到填充图像。
步骤5、对步骤4的填充图像进行两次形态学开运算，得到菌落开运算图像，如图6所示。
5-1.对步骤4的填充图像先进行两次形态学腐蚀，再进行两次形态学膨胀。
步骤6、对步骤5所得的开运算图像进行黏连分割，得到黏连分割图像，如图7所示。黏连分割的具体过程为：
6-1.根据下面公式(6)，进行八连通链码标记。链码标记具体过程为：
式(6)中，j为链码序号；k为链码步进长度；ζj,k为链码第j个点的链码值；(xj,yj)、(xj-k,yj-k)分别表示链码第j、j-k个点坐标。
6-1-1.以起始点开始，并记录为起始点坐标，以8连通的方式从编码为1的方向(45度)开始搜索下一个边界像素点，找到之后记录方向编码(0到7)。
6-1-2.从找到的像素开始重复步骤6-1-1，得出一组链码。如果是闭合的，编码完成后返回到起始点，否则可以省略起始点坐标。
6-1-3.对步骤6-1-2所得的链码进行归一化，得归一化链码。将链码看作n位自然数，将链码按一个方向循环，使其构成的n位自然数最小，形成起点唯一的链码。
6-2.根据下面向量叉乘运算公式(7)，确定连通域边界上的凹点。
a×b＝|a||b|sin(θ) 式(7)
式(7)中，a表示以链码的起始坐标开始，以距离步长k所做的向量，其步长k视具体情况而定，在本算法中选取k＝1，b表示以向量a的终点作为起点，以同步长k所做的向量。其中θ为向量a到向量b的夹角，叉乘方向满足右手螺旋准则。
6-2-1.对步骤6-1所得的图像边界链码，再根据公式(7)进行叉乘运算，得到结果存放在集合Se中。定义一个阈值Th。若两个向量叉乘运算结果大于阈值Th，则这两个向量位于凹点周围，反之则位于非凹点周围，将此定义为凹点特性条件。
6-2-2.对步骤6-2-1所得的集合Se中满足凹点特性条件的结果进行标记，得到一系列位于凹点周围的点的多个集合，再将每个集合中绝对值最大的那个点作为该集合的凹点。对其它点集合进行同样的处理，得到图像中所有的凹点。
6-3.对步骤6-2所得的连通域边界上的凹点进行匹配分割，凹点匹配分割的具体过程为：
6-3-1.将凹点分类为内凹点和外凹点，内凹点就是沿着图像的内边缘依次寻找得到的凹点，外凹点就是沿着图像的外边缘依次寻找得到的凹点。
6-3-2.若连通域边界上的外凹点个数小于1，不进行分割。
6-3-3.若连通域边界上的外凹点个数为2，两个外凹点用四个像素宽的直线相连分割。
6-3-4.若连通域边界上的外凹点个数大于2。外凹点与外凹点配对，若满足欧式距离最小且外凹点的分割线在同方向上，直接配对；若不满足，则进行下一个外凹点配对。
6-3-5.用四个像素宽的直线连接所有配对上的凹点实现分割。
步骤7、对步骤6所得的黏连分割图像进行轮廓绘制并进行计数，且使用随机颜色进行显示，得到轮廓图像。如图8所示。
步骤8、将步骤7所得的轮廓图像进行亚像素重建，得到轮廓亚像素重建图像。如图9所示，局部放大如图10所示。轮廓亚像素重建的具体步骤为：
8-1.对步骤7所得的轮廓图像进行灰度化，过程根据步骤2中的步骤2-1流程一样。
8-2.根据下面梯度图像变换公式(8)，进行梯度检测，梯度检测包含幅值和方向，设置8个3*3方向模板，用8个方向模板在检测图像上进行卷积运算(模板滑过的地方的像素值和模板中对应的值相乘，最后全部加总),得到8个梯度图像。模板的方向就是和X轴的夹角。
式(8)中，R(i,j)为坐标(i,j)上的梯度幅值，其代表梯度图像。f(i-1,j-1)、f(i,j-1)、f(i+1,j-1)、f(i-1,j)、f(i+1,j)、f(i-1,j+1)、f(i,j+1)、f(i+1,j+1)分别为坐标(i,j)左上、上、右上、左、右、左下、下、右下八个方向上的灰度值。
8-3.定义梯度幅值最大的方向为梯度方向，由于8个梯度图像中都是幅值比较大的梯度方向，边缘可能有多个像素表示，而梯度幅值图中最亮的点才是边缘，即梯度幅值最大值处的点就是边缘的整数坐标。通过梯度方向上的变化规律，在3*3的邻域中根据梯度的方向比较中心点与该方向相邻点的幅值，幅值大的为边缘点，进行非极大值抑制。得到只有一个像素宽度的单像素边缘。
8-4.根据下面公式(9)，对步骤8-3所得的单像素边缘进行多项式插值，得到亚像素边缘。
式(9)中，(Xc，Yc)表示亚像素坐标；(m，n)表示像素级坐标；W为相邻像素点到边缘点的距离(其中梯度方向为0°,90°,180°,270°时等于1,其他角度为
)；θ为梯度方向与X轴正向夹角；R0为边缘点(m,n)的灰度梯度的模，R-1、R1分别是在梯度方向上与边缘点相邻的两像素点的梯度幅值。
8-5.对步骤8-4得到的边缘亚像素坐标进行存储，创建两个通道的图像，第一个通道表示x,第二个通道表示y,x通道存储像素值，y通道存储的是亚像素坐标。该方式保持了和图像相关联的状态。
8-6.对步骤8-5所得的亚像素轮廓使用随机颜色进行绘制连接。
步骤9、对步骤8所得的轮廓亚像素重建图像进行轮廓分割，得到轮廓分割图像。如图11所示、局部放大如图12所示。轮廓分割的具体流程为：
9-1.对步骤8所得的轮廓亚像素重建图像进行链码标记，过程根据步骤6中的步骤6-1流程一样。
9-2.对步骤9-1中的链码进行一阶差分链码计算以及一阶差分链码的归一化，具体流程为：
9-2-1.根据下面公式(10)计算一阶差分链码，其由相邻两个元素方向变化(逆时针方向)得到。即前一个数字变化到后一个数字需要经过的步数。
式(10)中，j为链码序号；k为链码步进长度；δj,k为链码第j个点的一阶差分值；(xj-k,yj-k)、(xj,yj)、(xj+k,yj+k)分别表示链码第j-k、j、j+k个点坐标。
9-2-2.对步骤9-2-1所得的一阶差分链码进行归一化处理，把一阶差分链码看作由各方向数所构成的自然数，再按一个方向循环，使其构成自然数最小。
9-3.依据拐点在切线处不连续，曲率连续的特点，对步骤8所得的轮廓图像进行拐点的提取。拐点提取的具体流程为：
9-3-1.根据步骤9-2-1中的公式(10)，得一阶差分链码，也表示轮廓切线的变化率。根据下面公式(11)，为差分链码减去链码均值的绝对值，表示偏离平均值的绝对值。用绝对值建立高斯分布，在某个置信范围内选取阈值。将偏离均值最大的一部分点作为分段点。
式(11)中，j为链码序号；k为链码步进长度；
为差分链码减去链码均值的绝对值，代表偏离平均值的绝对值；δj,k为链码第j个点的一阶差分值；μδ表示链码均值。
9-3-2.根据下面公式(12)，得阈值。若符合大于阈值且为邻域最大值的分段点视为拐点。
T＝z1-2ασ 式(12)
式(12)中，T为所求阈值；z为步骤9-3-1中建立的二维高斯分布的横坐标轴，1-2α为步骤9-3-1中建立的二维高斯分布的置信区间；z1-2α为步骤9-3-1中建立的二维高斯分布的置信区间横坐标值；σ为步骤9-3-1中建立的二维高斯分布的标准方差。
9-4.依据切点在切线处连续，曲率不连续的特点。对步骤8所得的轮廓图像进行切点的提取。切点提取的具体流程为：
9-4-1.对步骤8所得的轮廓图像进行多边形相似逼近，轮廓变为多折线段，得到平滑分裂点。根据下面误差检测分裂点的标准公式(13)，得轮廓和近似函数之间的最大误差。若达到分段数量或者分裂点最大误差低于设定阈值时停止分裂。
式(13)中，E∞(u1,ui+1)为逼近误差，其中ei,j是区间i的第j个点与其近似函数之间的欧几里德距离，并且(u1,ui+1)是区间的两个端点。
9-4-2.对步骤9-4-1所得到的分裂点进行分段近似椭圆拟合，选取第一个分裂点为起点，根据相邻轮廓段的三个分裂点不可能处在同一直线上这一前提下，将起点与相邻的两个分裂点作为集合进行近似椭圆拟合，根据公式(13)计算该近似椭圆与该分裂点集合对应轮廓段之间的最大偏差值。
9-4-3.若最大偏差值比该分裂点集合对应轮廓段与多边形相似逼近线段的最大距离要小，将该分裂点集合对应轮廓段就标记为合并处理对象，保存标记该分裂点集合，将下一个分裂点作为起点。执行步骤9-4-2的操作。
9-4-4.否则将下一个分裂点加入该分裂点集合中，将分裂点集合进行近似椭圆拟合，根据公式(13)计算该近似椭圆与该分裂点集合对应轮廓段之间的最大偏差值。
9-4-5.重复执行步骤9-4-3、步骤9-4-4，经过多次迭代直到没有合并处理的对象为止。
9-4-6.将步骤9-4-5所标记的每个分裂点集合的起点和终点视为切点。
9-5.根据步骤9-3、步骤9-4将拐点和切点全部提取出来，将拐点之间用线段连接，切点之间用最小二乘法进行椭圆拟合，并统计线段和椭圆弧的数量。得到轮廓分割为直线段、椭圆弧段的图像。
9-6.对步骤9-5所得的直线段和椭圆弧段使用随机颜色进行绘制。
步骤10、提取轮廓特征，包括数量、面积、周长、长短比、圆形度、轮廓线段和椭圆弧的统计。提取轮廓特征的具体过程为：
10-1.大肠杆菌数量：连通域个数。
10-2.单个大肠杆菌面积：连通域像素总数。
10-3.单个大肠杆菌周长：边界像素总数。
10-4.单个大肠杆菌长短比：长轴短轴比。
10-5.单个大肠杆菌圆形度：区域面积与π乘中心点到所有轮廓像素距离的最大值平方之比。
10-6.单个大肠杆菌轮廓直线和椭圆弧的统计：根据步骤9中的步骤9-5，得单个大肠杆菌轮廓的线段和椭圆弧的数量。
步骤11、对步骤10所得到的大肠杆菌的轮廓特征进行大肠杆菌形态分析，包括抗生素对大肠杆菌的抑制率的分析，如图2所示。大肠杆菌形态分析的具体流程为：
11-1.若大肠杆菌的面积大于100或周长大于40，则判断大肠杆菌的长短比是否大于3且小于10。
11-2.若否，则判定为杂质，若是，则判断大肠杆菌的轮廓中直线数量是否大于1且椭圆数量是否大于2。
11-3.若是，将其判定为被抑制的大肠杆菌；若否，将其判定为已耐药的大肠杆菌。
11-4.若大肠杆菌的面积小于100且周长小于40，则判断大肠杆菌的圆形度是否大于0.5。
11-5.若是，则判定为被抑制的大肠杆菌；若否，则判定为杂质。
11-6.循环遍历整个轮廓，执行步骤11-1到步骤11-5之间的流程。将被抑制的大肠杆菌比上大肠杆菌的总数量定义为大肠杆菌的抑制率，得大肠杆菌的抑制率。
11-7.每隔30分钟重复执行一次步骤1至11，从而得到在不同时间段下抗生素对大肠杆菌的抑制效果。进而可以推出大肠杆菌的耐药性演化。且根据每个大肠杆菌在不同时间段的信息，可以计算出大肠杆菌的生长速度、状况，统计出整体大肠杆菌的生长特征。
Claims (10)
1.一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：包括以下步骤：
步骤1、采集被识别的大肠杆菌形态图像；
步骤2、对步骤1所得的大肠杆菌形态图像进行预处理，得到预处理图像；
步骤3、对步骤2所得的预处理图像进行自适应双阈值二值化，得到二值图像；
步骤4、对步骤3所得的二值图像进行区域填充，得到填充图像；
步骤5、对步骤4的填充图像进行两次形态学开运算，得到菌落开运算图像；
步骤6、对步骤5所得的开运算图像进行黏连分割，得到黏连分割图像；
步骤7、对步骤6所得的黏连分割图像进行轮廓绘制并进行计数，得到轮廓图像；
步骤8、将步骤7所得的轮廓图像进行亚像素重建，得到轮廓亚像素重建图像；
步骤9、对步骤8所得的轮廓亚像素重建图像进行轮廓分割，得到轮廓分割图像；
步骤10、提取大肠杆菌的轮廓特征；
根据步骤10得到的大肠杆菌的轮廓特征，进行大肠杆菌形态分析。
2.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：被测的大肠杆菌培养皿中投加有抗生素；每隔预设时长重复执行一次步骤1至11，得到在不同时间段下，大肠杆菌受到抗生素的抑制效果。
3.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤2中所述预处理的具体过程为：
2-1.对步骤1所得的大肠杆菌形态图像使用分量权重法进行灰度化，其中三通道形态图像得到单通道形态灰度图像；
2-2.对步骤2-1所得的形态灰度图像进行背景校正，得到形态背景校正图像，背景校正的具体过程为：
2-2-1.选取3*3的卷积核模板，该卷积核以中心位置作为坐标原点；
2-2-2.使用步骤2-2-1所得的3*3卷积核对步骤2-1所得的大肠杆菌形态灰度图像做卷积，得到轻模糊处理图像o(x,y)；
2-2-3.选取25*25的卷积核模板，对步骤2-1所得的大肠杆菌形态灰度图像进行高斯滤波，得到强模糊处理图像b(x,y)；并求取强模糊处理图像b(x,y)的平均灰度M；
2-2-4.对步骤2-1所得的大肠杆菌灰度图像进行背景校正，得到形态背景校正图像f(x,y)如下：
f(x,y)＝o(x,y)-b(x,y)+M
2-3.对步骤2-2所得的形态背景校正图像进行双边滤波，得到大肠杆菌形态双边滤波灰度图像；双边滤波方法具体为将空间域核与颜色域核叠加得到双核；用双核对被识别图像进行滤波；
2-4.对步骤2-3所得的大肠杆菌形态双边滤波灰度图像进行直方图均衡化增强，得到大肠杆菌形态灰度增强图像；
2-5.对步骤2-4所得的大肠杆菌形态灰度增强图像进行形态学梯度变换；得到大肠杆菌形态学梯度变换图像g如下：
式(2)中，f表示大肠杆菌形态灰度增强图像，B表示3*3的交叉形内核，
表示膨胀操作，
表示腐蚀操作；大肠杆菌形态学梯度变换图像g即为步骤2得到的预处理图像。
4.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤3中所述自适应双阈值二值化的具体过程为：
3-1.对步骤2所得的预处理图像进行灰度梯度计算，梯度计算包括计算梯度方向与幅值，结合Sobel算子，计算出四个方向上的梯度变化强度，再通过线性插值计算两个对角线方向上的梯度变化强度；
3-2.将当前像素梯度变化强度与沿正负、对角线梯度方向上的相邻像素的梯度强度进行比较，若一个像素点为比较中的最大值，则保留该像素为边缘点，像素值设为255；若不是最大，则对其进行抑制，不将其作为边缘点；像素值设为0；
3-3.通过Otsu算法，将梯度图像非真实边缘部分划分为伪边缘部分和非边缘部分；将迭代分割运算中取得最大类间方差时的分割阈值作为最优分割阈值，抑制孤立的弱边缘，完成边缘检测；类间方差的表达式如下：
其中，w0为目标点数占总图像比例，u0为其平均灰度；w1为背景点数占图像比例，u1为其平均灰度；u为图像总平均灰度；g′为类间方差；
取类间方差值取最大时对应的阈值作为最优分割阈值，并将其定义为低阈值；将低阈值的2.5倍定义为高阈值；梯度强度低于低阈值的像素点被抑制，不作为边缘点，像素值设为0；高于高阈值的像素点被定义为强边缘，保留为边缘点，像素值设为255；处于高低阈值之间的定义为弱边缘；若弱边缘的邻域像素中有至少一个为强边缘像素点，则将该弱边缘保留为强边缘，像素值设为255；否则，该弱边缘不作为边缘点，像素值设为0。
5.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤4中进行区域填充的具体过程为：
4-1.对步骤3所得的二值图像中的连通域之外的像素值全部赋值为0，再将所得图像进行按位取反处理；
对步骤4-1所得的图像与步骤3所得的二值图像进行按位取或处理，得到填充图像。
6.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤6中所述黏连分割的具体过程为：
6-1.对黏连分割图像进行八连通链码标记；链码标记具体过程为：
式(4)中，j为链码序号；k为链码步进长度；ζj,k为链码第j个点的链码值；(xj,yj)、(xj-k,yj-k)分别表示链码第j、j-k个点坐标；
6-1-1.以起始点开始，并记录为起始点坐标，以8连通的方式从编码为1的方向开始搜索下一个边界像素点；之后，记录方向编码；
6-1-2.从找到的像素开始重复步骤6-1-1，得出一组链码；若该组链码闭合，则在编码完成后返回到起始点；否则，省略起始点坐标；
6-1-3.对步骤6-1-2所得的链码进行归一化，得归一化链码；将链码看作n位自然数，将链码按一个方向循环，使其构成的n位自然数最小，形成起点唯一的链码；
6-2.确定连通域边界上的凹点；
6-2-1.对步骤6-1所得的图像边界链码，再根据向量叉乘运算公式进行叉乘运算，所得向量叉乘运算结果存放在集合Se中；定义一个阈值Th；若向量叉乘运算结果大于阈值Th，则该两个向量对应的像素点满足凹点特性条件；
向量叉乘运算公式如下：
s＝a×b＝|a||b|sin(θ)
其中，s为向量叉乘运算结果；a表示以链码的起始坐标为起点，以步长k为距离所做的向量；b表示以向量a的终点作为起点，以步长k为距离所做的向量；θ为向量a到向量b的夹角；
6-2-2.对步骤6-2-1所得的集合Se中满足凹点特性条件的像素点进行标记，得到多个凹点待定集合；每个凹点待定集合均包括一系列位于凹点周围的点；再将每个凹点待定集合中绝对值最大的那个点作为该集合的凹点，得到图像中所有的凹点；
6-3.对步骤6-2所得的连通域边界上的凹点进行匹配分割，凹点匹配分割的具体过程为：
6-3-1.将凹点分类为内凹点和外凹点，内凹点表示沿着图像的内边缘依次寻找得到的凹点；外凹点表示沿着图像的外边缘依次寻找得到的凹点；
6-3-2.若连通域边界上的外凹点个数小于1，则不进行分割；
6-3-3.若连通域边界上的外凹点个数为2，则通过将该两个外凹点直线相连的方式进行分割；
若连通域边界上的外凹点个数大于2；则将各外凹点两两之间配对；若两个外凹点满足欧式距离最小且外凹点的分割线在同方向上，则该两个外凹点直接配对；将各组配对的外凹点直线相连，完成分割。
7.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤8中进行轮廓亚像素重建的具体步骤为：
8-1.对步骤7所得的轮廓图像进行灰度化；
8-2.灰度化后的轮廓图像进行梯度检测，梯度检测包含幅值和方向，设置8个3*3方向模板，用8个方向模板在图像上进行卷积运算，得到8个梯度图像；模板的方向就是梯度方向和X轴的夹角；
梯度检测使用的梯度图像变换公式如下：
R(i,j)＝∣f(i+1,j+1)+2f(i,j+1)+f(i-1,j+1)-f(i-1,j-1)-2f(i,j-1)-f(i+1,j-1)∣+∣f(i-1,j-1)+2f(i-1,j)+f(i-1,j+1)-f(i+1,j+1)-2f(i+1,j)-f(i+1,j-1)∣
其中，R(i,j)为坐标(i,j)上的梯度幅值；f(i-1,j-1)、f(i,j-1)、f(i+1,j-1)、f(i-1,j)、f(i+1,j)、f(i-1,j+1)、f(i,j+1)、f(i+1,j+1)分别为坐标(i,j)左上、上、右上、左、右、左下、下、右下八个方向上的灰度值；
8-3.通过梯度方向上的变化规律，在3*3的邻域中根据梯度的方向比较中心点与该方向相邻点的幅值，幅值大的为边缘点，进行非极大值抑制；得到只有一个像素宽度的单像素边缘；
8-4.对步骤8-3所得的单像素边缘进行多项式插值，得到亚像素边缘；多项式插值的公式如下：
其中，(Xc，Yc)表示亚像素坐标；(m，n)表示像素级坐标；W为相邻像素点到边缘点的距离；θ为梯度方向与X轴正向夹角；R0为边缘点(m,n)的灰度梯度的模，R-1、R1分别是在梯度方向上与边缘点相邻的两像素点的梯度幅值；
8-5.对步骤8-4得到的边缘亚像素坐标进行存储，创建两个通道的图像，第一个通道表示x，第二个通道表示y，x通道存储像素值，y通道存储的是亚像素坐标；
对步骤8-5所得的各个亚像素轮廓分别进行绘制连接。
8.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤9中进行轮廓分割的具体流程为：
9-1.对步骤8所得的轮廓亚像素重建图像进行链码标记；
9-2.对步骤9-1中的链码进行一阶差分链码计算以及一阶差分链码的归一化；
9-3.依据拐点在切线处不连续，曲率连续的特点，对步骤8所得的轮廓亚像素重建图像进行拐点的提取；
9-4.对步骤8所得的轮廓图像进行切点的提取；
9-5.根据步骤9-3、步骤9-4将拐点和切点全部提取出来，将拐点之间用线段连接，切点之间用最小二乘法进行椭圆拟合，并统计线段和椭圆弧的数量；得到轮廓分割为直线段、椭圆弧段的图像；
9-6.对步骤9-5所得的直线段和椭圆弧段进行绘制；
步骤10中，提取轮廓特征的具体过程为：
10-1.取连通域个数作为大肠杆菌数量；
10-2.取单个连通域像素总数，作为对应的单个大肠杆菌面积；
10-3.取单个连通域的边界像素总数，作为对应的单个大肠杆菌周长；
10-4.取单个连通域的长轴短轴比，作为对应的单个大肠杆菌长短比；
10-5.取单个连通域的区域面积与π乘中心点到所有轮廓像素距离的最大值平方之比，作为对应的单个大肠杆菌圆形度；
10-6.统计单个大肠杆菌的轮廓直线和椭圆弧数量；
步骤11中进行大肠杆菌形态分析的具体流程为：
11-1.若大肠杆菌的面积大于100或周长大于40，则判断大肠杆菌的长短比是否大于3且小于10；
11-2.若否，则判定为杂质，若是，则判断大肠杆菌的轮廓中直线数量是否大于1且椭圆数量是否大于2；
11-3.若是，将其判定为被抑制的大肠杆菌；若否，将其判定为已耐药的大肠杆菌；
11-4.若大肠杆菌的面积小于100且周长小于40，则判断大肠杆菌的圆形度是否大于0.5；
11-5.若是，则判定为被抑制的大肠杆菌；若否，则判定为杂质；
取被抑制的大肠杆菌与大肠杆菌的总数量的比值，作为大肠杆菌的抑制率。
9.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤9-2的具体流程为：
9-2-1.计算一阶差分链码δj,k如下：
其中，j为链码序号；k为链码步进长度；δj,k为链码第j个点的一阶差分值；(xj-k,yj-k)、(xj,yj)、(xj+k,yj+k)分别表示链码第j-k、j、j+k个点坐标；
对步骤9-2-1所得的一阶差分链码进行归一化处理，把一阶差分链码看作由各方向数所构成的自然数，再按一个方向循环，使其构成自然数最小。
10.根据权利要求1所述的一种大肠杆菌轮廓的精确提取及识别方法，其特征在于：步骤9-3的具体流程为：
9-3-1.计算差分链码偏离平均值的绝对值
用该绝对值
建立高斯分布，并选取阈值；将偏离均值最大的一部分点作为分段点；绝对值
的表达式如下：
其中，j为链码序号；k为链码步进长度；δj,k为链码第j个点的一阶差分值；μδ表示链码均值；
9-3-2.设定阈值T；将大于阈值T且为邻域最大值的分段点作为拐点；阈值T的表达式如下：
T＝z1-2ασ
其中，z为步骤9-3-1中建立的二维高斯分布的横坐标轴，1-2α为步骤9-3-1中建立的二维高斯分布的置信区间；z1-2α为步骤9-3-1中建立的二维高斯分布的置信区间横坐标值；σ为步骤9-3-1中建立的二维高斯分布的标准方差；
步骤9-4的具体流程为：
9-4-1.对步骤8所得的轮廓图像进行多边形相似逼近，轮廓变为多折线段，得到平滑分裂点；计算轮廓和近似函数之间的最大误差E∞(u1,ui+1)；若达到分段数量或者分裂点最大误差低于设定阈值时停止分裂；最大误差E∞(u1,ui+1)的表达式如下：
其中，ei,j是区间i的第j个点与其近似函数之间的欧几里德距离；(u1,ui+1)是区间i的两个端点；
9-4-2.对步骤9-4-1所得到的分裂点进行分段近似椭圆拟合，选取第一个分裂点为起点，根据相邻轮廓段的三个分裂点不可能处在同一直线上这一前提下，将起点与相邻的两个分裂点作为集合进行近似椭圆拟合，并计算该近似椭圆与该分裂点集合对应轮廓段之间的最大误差；
9-4-3.若最大误差小于该分裂点集合对应轮廓段与多边形相似逼近线段的最大距离，则将该分裂点集合对应轮廓段就标记为合并处理对象，保存标记该分裂点集合，将下一个分裂点作为起点；执行步骤9-4-2的操作；否则，将下一个分裂点加入该分裂点集合中，将分裂点集合进行近似椭圆拟合，计算该近似椭圆与该分裂点集合对应轮廓段之间的最大误差；
9-4-4.重复执行步骤9-4-3，经过多次迭代直到没有合并处理的对象为止；
将标记的每个分裂点集合的起点和终点视为切点。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210313888.9A CN114648511A (zh) | 2022-03-28 | 2022-03-28 | 一种大肠杆菌轮廓的精确提取及识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210313888.9A CN114648511A (zh) | 2022-03-28 | 2022-03-28 | 一种大肠杆菌轮廓的精确提取及识别方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114648511A true CN114648511A (zh) | 2022-06-21 |
Family
ID=81995539
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210313888.9A Pending CN114648511A (zh) | 2022-03-28 | 2022-03-28 | 一种大肠杆菌轮廓的精确提取及识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114648511A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115082441A (zh) * | 2022-07-22 | 2022-09-20 | 山东微山湖酒业有限公司 | 一种基于计算机视觉的酿酒蒸馏工艺中甑料平铺方法 |
| CN116030460A (zh) * | 2023-02-22 | 2023-04-28 | 青岛海关技术中心 | 一种基于图像特征的细菌显微图像识别方法 |
| CN116704209A (zh) * | 2023-08-08 | 2023-09-05 | 山东顺发重工有限公司 | 一种法兰轮廓快速提取方法及系统 |
| CN117495851A (zh) * | 2023-12-29 | 2024-02-02 | 陕西中医药大学 | 一种基于图像轮廓处理的水环境微生物检测方法 |
-
2022
- 2022-03-28 CN CN202210313888.9A patent/CN114648511A/zh active Pending
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115082441A (zh) * | 2022-07-22 | 2022-09-20 | 山东微山湖酒业有限公司 | 一种基于计算机视觉的酿酒蒸馏工艺中甑料平铺方法 |
| CN115082441B (zh) * | 2022-07-22 | 2022-11-11 | 山东微山湖酒业有限公司 | 一种基于计算机视觉的酿酒蒸馏工艺中甑料平铺方法 |
| CN116030460A (zh) * | 2023-02-22 | 2023-04-28 | 青岛海关技术中心 | 一种基于图像特征的细菌显微图像识别方法 |
| CN116704209A (zh) * | 2023-08-08 | 2023-09-05 | 山东顺发重工有限公司 | 一种法兰轮廓快速提取方法及系统 |
| CN116704209B (zh) * | 2023-08-08 | 2023-10-17 | 山东顺发重工有限公司 | 一种法兰轮廓快速提取方法及系统 |
| CN117495851A (zh) * | 2023-12-29 | 2024-02-02 | 陕西中医药大学 | 一种基于图像轮廓处理的水环境微生物检测方法 |
| CN117495851B (zh) * | 2023-12-29 | 2024-04-05 | 陕西中医药大学 | 一种基于图像轮廓处理的水环境微生物检测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114648511A (zh) | 一种大肠杆菌轮廓的精确提取及识别方法 | |
| CN105069389B (zh) | 二维码分块解码方法及系统 | |
| CN114723701A (zh) | 基于计算机视觉的齿轮缺陷检测方法和系统 | |
| CN111340824B (zh) | 一种基于数据挖掘的图像特征分割方法 | |
| CN111178173A (zh) | 一种目标菌落生长特征识别方法 | |
| CN108615239B (zh) | 基于阈值技术和灰度投影的舌图像分割方法 | |
| CN106503704B (zh) | 一种自然场景中圆形交通标志定位方法 | |
| CN102346850A (zh) | 复杂金属背景下DataMatrix条码区域定位方法 | |
| CN108376403B (zh) | 基于霍夫圆变换的网格菌落图像分割方法 | |
| CN114219805B (zh) | 一种玻璃缺陷智能检测方法 | |
| CN104268602A (zh) | 一种基于二进制特征匹配的遮挡工件识别方法及装置 | |
| CN109035274A (zh) | 基于背景估计与u型卷积神经网络的文档图像二值化方法 | |
| CN109903282B (zh) | 一种细胞计数方法、系统、装置和存储介质 | |
| CN112464942B (zh) | 基于计算机视觉的重叠烟叶智能化分级方法 | |
| CN109389033B (zh) | 一种新型的瞳孔快速定位方法 | |
| CN106096491B (zh) | 一种眼底彩色照相图像中微动脉瘤自动化识别方法 | |
| CN115330758B (zh) | 基于去噪处理的焊接质量检测方法 | |
| CN110472479A (zh) | 一种基于surf特征点提取和局部lbp编码的指静脉识别方法 | |
| CN115457004A (zh) | 基于计算机视觉的锌膏的智能检测方法 | |
| CN113902713B (zh) | 一种基于图像分割算法的多菌落体计数方法 | |
| CN109359604B (zh) | 面向巡检机器人的阴影干扰下仪表识别方法 | |
| CN109872326B (zh) | 基于深度强化网络跳跃连接的轮廓检测方法 | |
| CN116934763B (zh) | 一种基于视觉特征的医用橡胶胶塞缺陷检测方法 | |
| CN112215865A (zh) | 一种荧光显微图像下的微液滴自动检测方法 | |
| CN109949334B (zh) | 基于深度强化网络残差连接的轮廓检测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |