CN113269756B - 基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 - Google Patents
基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 Download PDFInfo
- Publication number
- CN113269756B CN113269756B CN202110588676.7A CN202110588676A CN113269756B CN 113269756 B CN113269756 B CN 113269756B CN 202110588676 A CN202110588676 A CN 202110588676A CN 113269756 B CN113269756 B CN 113269756B
- Authority
- CN
- China
- Prior art keywords
- image
- scale
- segmentation
- blood vessel
- particle swarm
- 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.)
- Active
Links
- 230000011218 segmentation Effects 0.000 title claims abstract description 93
- 238000000034 method Methods 0.000 title claims abstract description 78
- 239000002245 particle Substances 0.000 title claims abstract description 76
- 238000001914 filtration Methods 0.000 title claims abstract description 50
- 210000001210 retinal vessel Anatomy 0.000 title claims description 30
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 86
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 32
- 238000005457 optimization Methods 0.000 claims abstract description 29
- 210000001525 retina Anatomy 0.000 claims abstract description 13
- 238000012545 processing Methods 0.000 claims abstract description 9
- 230000006870 function Effects 0.000 claims description 19
- 230000000694 effects Effects 0.000 claims description 18
- 238000000605 extraction Methods 0.000 claims description 14
- 230000002792 vascular Effects 0.000 claims description 14
- 238000003709 image segmentation Methods 0.000 claims description 12
- 238000012805 post-processing Methods 0.000 claims description 10
- 238000007781 pre-processing Methods 0.000 claims description 8
- 230000002207 retinal effect Effects 0.000 claims description 8
- 238000011156 evaluation Methods 0.000 claims description 7
- 238000002474 experimental method Methods 0.000 claims description 7
- 239000013598 vector Substances 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000004044 response Effects 0.000 claims 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 claims description 2
- 241000764238 Isis Species 0.000 claims 1
- 230000035945 sensitivity Effects 0.000 abstract description 10
- 238000005286 illumination Methods 0.000 abstract description 6
- 230000004256 retinal image Effects 0.000 description 7
- 238000012360 testing method Methods 0.000 description 5
- 238000012549 training Methods 0.000 description 5
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 201000010099 disease Diseases 0.000 description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000000877 morphologic effect Effects 0.000 description 2
- 230000001105 regulatory effect Effects 0.000 description 2
- 238000012706 support-vector machine Methods 0.000 description 2
- 208000024172 Cardiovascular disease Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 206010012601 diabetes mellitus Diseases 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- 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
- 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
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
-
- 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/10024—Color 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/20081—Training; Learning
-
- 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/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30041—Eye; Retina; Ophthalmic
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Mathematical Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Computing Systems (AREA)
- Medical Informatics (AREA)
- Quality & Reliability (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Multimedia (AREA)
- Image Processing (AREA)
- Eye Examination Apparatus (AREA)
Abstract
本发明公开了一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置。首先运用MSR算法对图像亮度进行处理，改善眼底图像光照不均匀现象；其次运用多尺度高斯匹配滤波以增强目标血管与背景对比度；最后，运用粒子群优化OTSU多阈值对图像进行分割。在DRIVE数据集进行分割，其准确率、灵敏度、特异性分别为：95.72％、77.98％、97.29％。分割结果证明，该方法准确度较高，可以分割出更多细小血管。
Description
技术领域
本发明涉及图像识别领域，尤其涉及一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法。
背景技术
视网膜血管图像分割是医疗图像研究中的一个重要课题，能够有效地辅助医生进行快速的心血管疾病、糖尿病等疾病的临床诊断和治疗。近年来，许多学者都在研究视网膜血管图像分割方面的研究内容并已经取得一定的成果。但是，由于视网膜图像的复杂性和图像采集过程中噪声和光照因素的影响，进行精确的视网膜血管图像分割仍然是一项具有挑战性的任务[1]。目前，国内外提出的视网膜血管分割方法主要分为基于监督学习的分割方法与无监督的分割方法。
在有监督分割方法方面，Staal等[2]提出基于KNN分类器自动视网膜血管分割方法，而后有学者陆续提出基于改进支持向量机(SVM)的分割方法[3]，基于U-net网络的视网膜血管分割方法[4]，基于KCN网络的视网膜血管分割方法[5]，基于CNN网络的分割方法[6]，基于局部回归的视网膜血管多标签分类方法[7]等。
在无监督传统分割方法方面，可以分为基于滤波器、基于聚类、基于跟踪、基于形态学、基于阈值的视网膜血管分割方法。Chaudhuri等[8]提出一种二维匹配滤波方法对视网膜血管进行分割，Odstrcilik等[9]提出改进匹配滤波视网膜分割方法，Gabor[10]首次提出了Gabor滤波器理论，之后Daugman[11]以该理论为基础设计出二维Gabor滤波器，Remco[12]提出Cake滤波器，文献[13]提出基于方向分数和Frangi滤波器的视网膜血管分割算法；基于模糊C均值聚类算法(FCM)的血管分割方法[14]，基于K-means和FCM的血管分割方法[15]；Tolias等[16]最早提出视网膜血管自动分割方法，文献[17]提出基于血管跟踪的方法提取出视网膜血管部分；Fraz等[18]提出基于形态学顶帽变换的分割方法，文献[19]提出一种基于验证的多阈值探测自适应局部阈值法并将其应用在视网膜血管分割中，汪维华等[20]提出改进的形态学与Otsu相结合的视网膜血管分割方法。
基于有监督的分割方法在训练时需要对应的血管图共同参与训练，通常需要从训练图像中选取较为关键的特征作为向量进行训练[1]，需要大量先验数据进行建模，但目前国内外公开视网膜图像较少，存在局限性。非监督类分割方法整体分割性能良好，但是易受噪声影响，细小血管容易丢失，分割精度不高。综上，本发明提出基于多尺度匹配滤波与粒子群优化的多阈值分割相结合的方法以分割出更多细小血管，提高分割准确率。
发明内容
1、本发明的目的
为解决视网膜血管分割时容易出现细小血管丢失断裂，分割效果不佳，影响疾病的诊断等问题，而提出了一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法及装置。
2、本发明所采用的技术方案
本发明公开了一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，具体为：
图像预处理，通过对原始图像进行多尺度MSR处理，调整图像亮度不均并减少噪声，提取灰度值分布均匀的绿色通道进行后处理；
图像对比度增强，通过大尺度高斯匹配滤波提取血管主要轮廓特征，小尺度提取血管细节特征，中尺度提取包括与背景对比度不强的轮廓特征与细节特征，将三个尺度下匹配滤波结果进行叠加，增强目标血管与背景对比度，获取多尺度下的血管特征；
图像多阈值分割，运用粒子群优化OTSU三阈值对图像预分割，运用OTSU单阈值将图像转化为二值图像；
图像后处理，最后运用图像后处理对图像操作包括断点连接、去噪、去除边缘轮廓。
更进一步，图像预处理，通过多尺度MSR处理，选取绿色通道图像后续处理：
多尺度MSR定义为：
其中Si(x,y)表示原始图像，F表示高斯函数，w表示每个尺度进行高斯卷积的权重，N表示尺度数量，R表示图像在对数域输出。
更进一步，提取血管特征：
采用二维匹配滤波方法对视网膜血管进行提取，二维高斯滤波器内核表达式为：
K(x,y)＝-exp(-X2/2σ2),|y|<2/L (2)
其中，L为假定血管具有固定方向的段的长度，由此设定高斯核的长度；σ为血管的尺度，由此设定高斯核的尺度；由于血管的多个生长方向，旋转高斯核使其与各个方向上血管的进行卷积，高斯核由0至180度每隔15度旋转一次，θ＝0°，15°，…，180°，共12个方向，保留各个像素响应最大值，旋转矩阵为：
更进一步，匹配滤波：
选择多尺度的匹配滤波器对血管图像的特征进行提取，选用大尺度对图像进行滤波时，提取粗血管，即当σ1＝1.9时可以有效提取血管的主轮廓特征；增加一个介于大小尺度之间的中间尺度对图像进行滤波，当σ2＝0.5时可以得到同时包括与背景对比度不强的轮廓特征与细节特征的图像；选用小尺度σ3＝0.13对图像进行滤波时细小血管检测。
更进一步，血管特征信息融合：
将大中小三个尺度进行匹配滤波结果叠加能有效增强目标血管与背景对比度，数学表达式如下：
I＝ω1I1+ω2I2+ω3I3 (4)
w1+w2+w3＝1 (5)
其中I为叠加后图像，I1、I2、I3分别为σ1＝1.9、σ2＝0.5、σ3＝0.13匹配滤波提取血管特征信息图。
更进一步，采用Otsu算法进行图像分割：
设给定图像的像素用L个灰度级表示[1,2,...,L]；第i级的像素数用ni表示，总像素数用N＝n1+n2+...+nL表示；灰度值为i的像素存在概率记为pi：
用t阈值将给定图像分为c0和c1两个区域；c0表示像素级别[1,2...,K]，c1表示级别为[K+1,...,L]的像素；该区域的概率和平均灰度值分别由式(7)和(8)给出，原始图像的总均值水平uT由式(9)给出；
公式(10)表示的以下两个关系进行验证：
w0u0+w1u1＝uT,w0+w1＝1 (10)
Otsu方法的目标函数定义为式(11)：
σ2 B＝w0(u0-uT)2+w1(u1-uT)2 (11)
当满足时σ2 B(t)＝Argmaxσ2 B，t取得最优值；
将OTSU单阈值扩展到多阈值，类间方差最大值为，
满足式(12)时，取得最佳阈值组合(t1，t2，...tm)。
更进一步，粒子群优化算法：
每个粒子相当于解空间中的一个待定解，粒子群的规模设为M；每个粒子i都有两个基本特征:当前的位置Xi＝(xi1,xi2,...xid)和当前速度Vi＝(vi1,vi2,...vid)；实践中评价解的好坏程度由所求优化问题的适应度函数决定；粒子i每次迭代中自身所经历的最佳位置记为Pbesti＝(pi1,pi2,...pid)，粒子群中所经历的最佳位置记为gbest＝(g1,g2,...gd)，粒子通过上述两个最佳位置动态更新自己的速度和位置，迭代公式如下[22]:
Vk id＝wVk-1 id+c1r1(pbestid-xk-1 id)+c2r2 (13)
xk id＝xk-1 id+vk-1 id (14)
其中vk id表示第k次迭代粒子i飞行速度矢量的第d维分量，xk id表示第k次迭代粒子i位置矢量的第d维分量，c1、c2为加速度常数，调节学习的步长，r1、r2两个随机函数，取值范围[0，1]，增加搜索随机性，ω>0，调节对解空间的搜索范围。
更进一步，粒子群优化OTSU图像分割：
由于经匹配滤波后的视网膜血管图像，含有背景，目标血管，噪声，运用一维OTSU单阈值，将其分为背景与目标两类，分割效果不佳；本发明采取多阈值对图像进行预分割，将图像分为多组目标与背景，但这需要在全灰度范围内搜索一个最佳门限组合，耗时较多，实际应用困难，为简化计算，提高运行速度，采用粒子群优化算法来搜索最佳阈值组合，式(12)便是粒子群优化算法的适应度函数；
经多次实验发现，当分割阈值数为3时，能达到较好的预分割效果；粒子群优化OTSU图像分割算法具体步骤如下：
Step1：初始化迭代次数M，种群规模N，维度D；
Step2：初始化粒子的速度和位置，以及个体极值pBesti和全局极值gBest；Step3：将式(12)作为粒子的适应度函数，计算每个粒子的适应度值，以此更新个体极值pBesti和全局极值gBest；
Step4：根据式(13)、(14)来更新粒子的速度和位置；
Step5：判断是否满足迭代停止条件，是则算法结束；否则转到Step3，继续进行迭代循环，最终找到最优解。
本发明公开了一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法的分割装置，包括存储器和处理器，存储器存储有计算机程序，所述处理器执行所述计算机程序时实现所述的方法步骤。
本发明公开了一种计算机可读存储介质，其上存储有计算机程序，所述的计算机程序被处理器执行时实现所述的方法步骤。
3、本发明所采用的有益效果
本发明运用MSR算法对原始视网膜图像进行亮度调节去除光照不均的影响，运用多尺度匹配滤波以增强目标血管与背景对比度，最后经粒子群优化OTSU多阈值对视网膜血管分割。对DRIVE数据集中的20幅测试集图像进行分割及评估，评估所获得的平均准确率为95.72％，平均灵敏度为77.98％，平均特异性为97.58％，保证特异性同时，准确率、灵敏度较高，本发明整体分割性能良好。
附图说明
图1为血管分割流程图；
图2为预处理结果图，(a)原始彩色图像、(b)MSR亮度调节、(c)绿色通道图像；
图3为视网膜血管信息提取图；(a)σ1＝1.9、(b)σ2＝0.5、(c)σ3＝0.13；
图4为多尺度与单尺度血管提取细节对比图；
图5为单阈值与多阈值分割结果对比图；
图6为分割结果与金标准对比图。
具体实施方式
下面结合本发明实例中的附图，对本发明实例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域技术人员在没有做创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
下面将结合附图对本发明实例作进一步地详细描述。
实施例1
基于多尺度匹配滤波与粒子群优化视网膜分割方法主要分为四个阶段：图像预处理、图像对比度增强、图像多阈值分割、图像后处理。首先，通过对原始图像进行多尺度Retinex(MSR)处理，调整图像亮度不均并减少噪声，提取灰度值分布均匀的绿色通道便于后处理；其次，通过大尺度高斯匹配滤波提取血管主要轮廓特征，小尺度提取血管细节特征，中尺度提取包括与背景对比度不强的轮廓特征与细节特征，将三个尺度下匹配滤波结果进行叠加，增强目标血管与背景对比度，获取多尺度下的血管特征；然后运用粒子群优化OTSU三阈值对图像预分割，运用OTSU单阈值将图像转化为二值图像；最后运用图像后处理对图像进行断点连接、去噪、去除边缘轮廓等操作。
本发明实验是在Windows 10系统，2.30GHz处理器，3.8GB运行内存，Matlab2016a编程环境下进行。使用DRIVE[23]数据库，DRIVE数据库是大小为565×585的彩色视网膜图像，该数据集分为两个子集：测试数据库子集和训练数据库子集，各含有20幅图像。
图像预处理
针对视网膜图像采集时受光照影响，图像亮度不均衡，使得目标血管与背景对比度低，对后期分割不利，本发明采用Jobson DJ，Rahman Z，Woodell GA等人[21]提出的多尺度Retinex(MSR)算法对视网膜图像进行处理，去除光照影响，进行亮度调节。
多尺度Retinex(MSR)定义为：
其中Si(x,y)表示原始图像，F表示高斯函数，w表示每个尺度进行高斯卷积的权重，N表示尺度数量，R表示图像在对数域输出。
通过图像红绿蓝三通道提取实验，发现绿色通道图像亮度均衡，目标血管与背景对比度高，灰度分布较均匀，红色通道目标与背景对比度低，蓝色通道噪声较大，故选取绿色通道图像以便后续处理。原始图像预处理结果如图2所示。
血管特征提取
本文采用Chaudhuri等[8]提出的二维匹配滤波方法对视网膜血管进行提取，二维高斯滤波器内核表达式为：
K(x,y)＝-exp(-X2/2σ2),|y|<2/L (2)
其中，L为假定血管具有固定方向的段的长度，由此设定高斯核的长度；σ为血管的尺度，由此设定高斯核的尺度；由于血管的多个生长方向，应旋转高斯核使其与各个方向上血管的进行卷积，高斯核由0至180度每隔15度旋转一次(θ＝0°，15°，…，180°)共12个方向，保留各个像素响应最大值，旋转矩阵为：
匹配滤波
由于视网膜血管生长方向不一，血管粗细长短不一，单一尺度下难以准确提取血管特征信息，为此，本文选择多尺度的匹配滤波器对血管图像的特征进行提取。选用大尺度对图像进行滤波时，主要提取的是粗血管，细小血管与背景融为一体，难以提取，经多次试验得到，当σ1＝1.9时可以有效提取血管的主轮廓特征；为了提取更多细小血管，增加一个介于大小尺度之间的中间尺度对图像进行滤波，当σ2＝0.5时可以得到同时包括与背景对比度不强的轮廓特征与细节特征的图像；选用小尺度对图像进行滤波时细小血管被检测出，当σ3＝0.13时可以有效提取血管的细节特征。尺度分别为1.9、0.5、0.13时的血管信息提取效果如图3所示。
血管特征信息融合
将大中小三个尺度进行匹配滤波结果叠加能有效增强目标血管与背景对比度，数学表达式如下：
I＝ω1I1+ω2I2+ω3I3 (4)
w1+w2+w3＝1 (5)
其中I为叠加后图像，I1、I2、I3分别为σ1＝1.9、σ2＝0.5、σ3＝0.13匹配滤波提取血管特征信息图，多尺度血管提取细节对比图如图4所示。
图4将多尺度与单尺度血管提取的结果进行了对比，其中，第一行分别为单尺度σ1＝1.9，σ3＝0.13和多尺度的提取效果，第二行分别对红色区域的细节信息进行放大显示。由图4可以看出，选择单一尺度对视网膜血管进行特征提取时，往往得不到较好的提取效果。如图4.(a)所示，σ1＝1.9时大量细小血管丢失，部分主血管也丢失，使血管结构不完整；图4.(b)中σ3＝0.13时少量细小血管丢失，存在血管断裂，血管连通性差，噪声被增强。而选择多尺度对图像特性进行提取并融合后，如图4.(c)和图4.(f)所示，在保留血管完整性同时，能有效提取较多细小血管，较小尺度特征提取而言减少了噪声影响，提取效果较好。
图像分割
Otsu算法最早提出于1979年，它通过使被分割的类的类方差最大来选择最优阈值。设给定图像的像素用L个灰度级表示[1,2,...,L]。第i级的像素数用ni表示，总像素数用N＝n1+n2+...+nL表示。灰度值为i的像素存在概率记为pi：
用t阈值将给定图像分为c0和c1两个区域。c0表示像素级别[1,2...,K]，c1表示级别为[K+1,...,L]的像素。该区域的概率和平均灰度值分别由式(7)和(8)给出，原始图像的总均值水平uT由式(9)给出。
公式(10)表示的以下两个关系可以很容易验证：
w0u0+w1u1＝uT,w0+w1＝1 (10)
Otsu方法的目标函数可以定义为式(11)：
σ2 B＝w0(u0-uT)2+w1(u1-uT)2 (11)
当满足时σ2 B(t)＝Argmaxσ2 B，t取得最优值。
将OTSU单阈值扩展到多阈值，类间方差最大值为，
满足式(12)时，取得最佳阈值组合(t1，t2，...tm)。
粒子群优化算法
在粒子群优化算法中，每个粒子相当于解空间中的一个待定解，粒子群的规模设为M。每个粒子i都有两个基本特征:当前的位置Xi＝(xi1,xi2,...xid)和当前速度Vi＝(vi1,vi2,...vid)。实践中评价解的好坏程度由所求优化问题的适应度函数决定。粒子i每次迭代中自身所经历的最佳位置记为Pbesti＝(pi1,pi2,...pid)，粒子群中所经历的最佳位置记为gbest＝(g1,g2,...gd)，粒子通过上述两个最佳位置动态更新自己的速度和位置，迭代公式如下[22]:
Vk id＝wVk-1 id+c1r1(pbestid-xk-1 id)+c2r2 (13)
xk id＝xk-1 id+vk-1 id (14)
其中vk id表示第k次迭代粒子i飞行速度矢量的第d维分量，xk id表示第k次迭代粒子i位置矢量的第d维分量，c1、c2为加速度常数，调节学习的步长，r1、r2两个随机函数，取值范围[0，1]，以增加搜索随机性，ω>0，调节对解空间的搜索范围。
粒子群优化OTSU图像分割
由于经匹配滤波后的视网膜血管图像，含有背景，目标血管，噪声，运用一维OTSU单阈值，将其分为背景与目标两类，分割效果不佳。本发明采取多阈值对图像进行预分割，将图像分为多组目标与背景，但这需要在全灰度范围内搜索一个最佳门限组合，耗时较多，实际应用困难，为简化计算，提高运行速度，采用粒子群优化算法来搜索最佳阈值组合，式(12)便是粒子群优化算法的适应度函数。
经多次实验发现，当分割阈值数为3时，能达到较好的预分割效果。粒子群优化OTSU图像分割算法具体步骤如下：
Step1：初始化迭代次数M，种群规模N，维度D；
Step2：初始化粒子的速度和位置，以及个体极值pBesti和全局极值gBest；
Step3：将式(12)作为粒子的适应度函数，计算每个粒子的适应度值，以此更新个体极值pBesti和全局极值gBest；
Step4：根据式(13)、(14)来更新粒子的速度和位置；
Step5：判断是否满足迭代停止条件，是则算法结束；否则转到Step3，继续进行迭代循环，最终找到最优解。
多尺度匹配滤波融合后的视网膜血管图像含有背景，目标血管，噪声，运用一维OTSU单阈值，将其分为背景与目标两类，分割效果不佳。本文采取多阈值对图像进行预分割，将图像分为多组目标与背景，但这需要在全灰度范围内搜索一个最佳阈值组合，耗时较多。为简化计算，提高运行速度，本文采用粒子群优化算法来搜索最佳阈值组合，式(12)便是粒子群优化算法的适应度函数。实验发现，分割阈值数为3时，能达到较好的预分割效果。多阈值与单阈值分割结果见图5。图5是在Drive数据集的测试集中随机选取的一张图片，在多尺度匹配框架下运用粒子群优化OTSU的单阈值与三阈值分割结果对比。其中图5.(a)是多尺度匹配滤波结果图，图5.(b)与图5.(c)分别为粒子群优化OTSU单阈值与三阈值分割结果图；图5.(d-f)分别为粒子群优化OTSU单阈值与三阈值分割细节放大图。在多尺度匹配滤波框架下，如图5.(e)所示，单阈值分割较多细小血管丢失，主血管存在断裂现象与如图5.(f)所示，相较于单阈值分割，三阈值分割可以分割出更多细小血管，且血管具有较好的连通性。
分割评价指标
视网膜图像血管分割就是将图像中的每个像素分为血管像素和非血管像素。像素点的统计有4种情况:真阳性(TP)表示正确分割为血管的点数；假阳性(FP)表示被错误分割的血管点数；真阴性(TN)表示正确分割为背景的点数；假阴性(FN)表示被错误分割的背景点数.可以通过像素点分类结果计算出准确率(Acc)、灵敏度(Sn)、特异度(Sp)，准确率表示分割正确的像素点数占整体像素点数的比率，灵敏度表示分割正确的血管点数占像素点数的比率，特异度表示分割正确的背景点数占像素点数的比率[12]，衡量指标计算公式如下:
分割结果与分析
为了保证分割结果以及评估指标计算的准确性，需要用一维OTSU单阈值对上节所获取的视网膜血管图像进行再分割，以去除非黑白的灰度值，将图像转换为二值图像。二值化的图像含有大量噪声，部分血管存在断裂现象，有误分割的眼底照相机的视场边缘存在。针对这些问题对图像进行后处理，首先采用中值滤波对图像进行初步去噪，同时可连接断裂血管，其次运用形态学处理使连通域面积小于一定值，此时可以去除大块噪声。最后提取图像掩膜，求取原始多尺度匹配滤波图像与掩膜的差值图像，将差值图像二值化，再对二值图进行膨胀，用分割后的视网膜血管图与膨胀边缘图相减，最终去除轮廓边缘。
本文利用Drive数据集对所提出算法进行测试，并将分割结果与第二专家的金标准图像、Chaudhuri[8]分割图像进行比较，视网膜血管分割效果对比图见图6。图6是在DRIVE测试集中随机选择的三幅视网膜图像，分别应用两种方法之后的分割效果对比。其中，第一列为原始彩色眼底视网膜图像，第二列为第二专家金标准图像，第三列为Chaudhuri分割结果，第四列为本文多尺度框架下分割结果。对比结果证明，文所提出的多尺度多阈值分割算法在保证血管结构完整性同时，可以分割出较多细小血管与第二位专家手工分割结果较相近。
分割性能
本发明分割金标准采用Drive数据集第二专家的手动分割结果，采用3个衡量指标来评价分割算法性能，并与文献[4][6][8][19]性能指标作比较，如表1所示：
表1 Drive数据库不同方法性能指标比较
表1中文献[19]提出的方法特异性最高，略高于本文算法的1.68％，但是本文算法的准确率与特异性分别比文献[19]高了1.90％与2.11％；相较于其他几种方法，本文所提出的方法获得的灵敏度(Sn)最高，较其它方法中最高的灵敏度也高出5.48％；相较于其他几种方法，本文提出的方法算法所得准确率(Acc)也最高，较其它方法中最高的准确率也高出0.52％。性能结果表明，本文分割方法在保证特异性同时，准确率、灵敏度均优于其它方法，所以发明整体分割性能较好。
为解决视网膜血管分割时容易出现细小血管丢失断裂，分割效果不佳，影响疾病的诊断等问题，本发明运用MSR算法对原始视网膜图像进行亮度调节去除光照不均的影响，运用多尺度匹配滤波以增强目标血管与背景对比度，最后经粒子群优化OTSU多阈值对视网膜血管分割。对DRIVE数据集中的20幅测试集图像进行分割及评估，评估所获得的平均准确率为95.72％，平均灵敏度为77.98％，平均特异性为97.58％，保证特异性同时，准确率、灵敏度较高，本方法整体分割性能良好。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。
Claims (7)
1.一种基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于：
图像预处理，通过对原始图像进行多尺度MSR处理，调整图像亮度不均并减少噪声，提取灰度值分布均匀的绿色通道进行后处理；
图像对比度增强，通过大尺度高斯匹配滤波提取血管主要轮廓特征，小尺度提取血管细节特征，中尺度提取包括与背景对比度不强的轮廓特征与细节特征，将三个尺度下匹配滤波结果进行叠加，增强目标血管与背景对比度，获取多尺度下的血管特征；
图像多阈值分割，运用粒子群优化OTSU三阈值对图像预分割，运用OTSU单阈值将图像转化为二值图像；
图像后处理，最后运用图像后处理对图像操作包括断点连接、去噪、去除边缘轮廓；
提取血管特征：
采用二维匹配滤波方法对视网膜血管进行提取，二维高斯滤波器内核表达式为：
（2）
其中，为假定血管具有固定方向的段的长度，由此设定高斯核的长度；/>为血管的尺度，由此设定高斯核的尺度；由于血管的多个生长方向，旋转高斯核使其与各个方向上血管的进行卷积，高斯核由0至180度每隔15度旋转一次，θ = 0°，15° ，…，180°，共12个方向，保留各个像素响应最大值，旋转矩阵为：
（3）；
匹配滤波：
选择多尺度的匹配滤波器对血管图像的特征进行提取，选用大尺度对图像进行滤波时，提取粗血管，即当时可以有效提取血管的主轮廓特征；增加一个介于大小尺度之间的中间尺度对图像进行滤波，当/>时可以得到同时包括与背景对比度不强的轮廓特征与细节特征的图像；选用小尺度/>对图像进行滤波时细小血管检测。
2.根据权利要求1所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于，图像预处理，通过多尺度MSR处理，选取绿色通道图像后续处理：
多尺度MSR定义为：
（1）
其中表示原始图像，F表示高斯函数，w表示每个尺度进行高斯卷积的权重，N表示尺度数量，R表示图像在对数域输出。
3.根据权利要求1所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于，血管特征信息融合：
将大中小三个尺度进行匹配滤波结果叠加能有效增强目标血管与背景对比度，数学表达式如下：
(4)
(5)
其中为叠加后图像，/>、/>、/>分别为/>、/>、/>匹配滤波提取血管特征信息图。
4.根据权利要求1所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于，采用Otsu算法进行图像分割：
设给定图像的像素用个灰度级表示/>；第/>级的像素数用/>表示，总像素数用表示；灰度值为/>的像素存在概率记为/>：
（6）
用t阈值将给定图像分为和/>两个区域；/>表示像素级别/>，/>表示级别为的像素；该区域的概率和平均灰度值分别由式（7）和（8）给出，原始图像的总均值水平/>由式（9）给出；
（7）
（8）
（9）
公式（10）表示的以下两个关系进行验证：
（10）
Otsu方法的目标函数定义为式(11)：
（11）
当满足时，t取得最优值；
Maximize （12）
满足式（12）时，取得最佳阈值组合（t1，t2，...tm）。
5.根据权利要求1所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于，粒子群优化算法：
每个粒子相当于解空间中的一个待定解，粒子群的规模设为M；每个粒子都有两个基本特征:当前的位置/>和当前速度 />；实践中评价解的好坏程度由所求优化问题的适应度函数决定；粒子i每次迭代中自身所经历的最佳位置记为，粒子群中所经历的最佳位置记为 />，粒子通过上述两个最佳位置动态更新自己的速度和位置，迭代公式如下:
（13）
（14）
其中表示第/>次迭代粒子/>飞行速度矢量的第/>维分量，/>表示第k次迭代粒子/>位置矢量的第/>维分量，/>为加速度常数，调节学习的步长，/>两个随机函数，取值范围[0，1]，增加搜索随机性，/>>0，调节对解空间的搜索范围。
6.根据权利要求5所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法，其特征在于，粒子群优化OTSU图像分割：
由于经匹配滤波后的视网膜血管图像，含有背景，目标血管，噪声，运用一维OTSU单阈值，将其分为背景与目标两类，分割效果不佳；本发明采取多阈值对图像进行预分割，将图像分为多组目标与背景，但这需要在全灰度范围内搜索一个最佳门限组合，耗时较多，实际应用困难，为简化计算，提高运行速度，采用粒子群优化算法来搜索最佳阈值组合，式（12）便是粒子群优化算法的适应度函数 ；
经多次实验发现，当分割阈值数为3时，能达到较好的预分割效果；粒子群优化OTSU图像分割算法具体步骤如下：
Step1：初始化迭代次数M，种群规模N，维度D；
Step2：初始化粒子的速度和位置，以及个体极值pBesti 和全局极值 gBest；Step3：将式（12）作为粒子的适应度函数，计算每个粒子的适应度值，以此更新个体极值pBesti 和全局极值 gBest；
Step4：根据式（13）、（14）来更新粒子的速度和位置；
Step5：判断是否满足迭代停止条件，是则算法结束；否则转到 Step3，继续进行迭代循环，最终找到最优解。
7.一种如权利要求1-6任一所述的基于多尺度匹配滤波与粒子群优化视网膜血管分割方法的分割装置，包括存储器和处理器，存储器存储有计算机程序，其特征在于；所述处理器执行所述计算机程序时实现如权利要求1-4任一所述的方法步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110588676.7A CN113269756B (zh) | 2021-05-28 | 2021-05-28 | 基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110588676.7A CN113269756B (zh) | 2021-05-28 | 2021-05-28 | 基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113269756A CN113269756A (zh) | 2021-08-17 |
| CN113269756B true CN113269756B (zh) | 2024-02-27 |
Family
ID=77233487
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110588676.7A Active CN113269756B (zh) | 2021-05-28 | 2021-05-28 | 基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113269756B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114677497A (zh) * | 2022-05-18 | 2022-06-28 | 浙江大华技术股份有限公司 | 一种图像处理方法及装置 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008091568A2 (en) * | 2007-01-23 | 2008-07-31 | Dtherapeutics, Llc | Applications of scaling laws of tree structures |
| CN102496156A (zh) * | 2011-11-17 | 2012-06-13 | 西安电子科技大学 | 基于协同量子粒子群算法的医学图像分割方法 |
| CN104899862A (zh) * | 2015-04-01 | 2015-09-09 | 武汉工程大学 | 基于全局与局部阈值的视网膜血管分割算法 |
| CN104899876A (zh) * | 2015-05-18 | 2015-09-09 | 天津工业大学 | 一种基于自适应高斯差分的眼底图像血管分割方法 |
| WO2016177337A1 (en) * | 2015-05-05 | 2016-11-10 | Shanghai United Imaging Healthcare Co., Ltd. | System and method for image segmentation |
| CN106407917A (zh) * | 2016-09-05 | 2017-02-15 | 山东大学 | 基于动态尺度分配的视网膜血管提取方法及系统 |
| CN106651846A (zh) * | 2016-12-20 | 2017-05-10 | 中南大学湘雅医院 | 视网膜血管图像的分割方法 |
| CN106683061A (zh) * | 2017-01-05 | 2017-05-17 | 南京觅踪电子科技有限公司 | 一种基于修正的多尺度retinex算法对医疗图像进行增强的方法 |
| CN111351801A (zh) * | 2020-03-11 | 2020-06-30 | 春光线缆有限公司 | 一种电线电缆缺陷检测系统 |
| CN112053363A (zh) * | 2020-08-19 | 2020-12-08 | 苏州超云生命智能产业研究院有限公司 | 视网膜血管分割方法、装置及模型构建方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7769513B2 (en) * | 2002-09-03 | 2010-08-03 | Automotive Technologies International, Inc. | Image processing for vehicular applications applying edge detection technique |
| US20060251325A1 (en) * | 2004-11-08 | 2006-11-09 | Charles Florin | Particle filter based vessel segmentation |
-
2021
- 2021-05-28 CN CN202110588676.7A patent/CN113269756B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008091568A2 (en) * | 2007-01-23 | 2008-07-31 | Dtherapeutics, Llc | Applications of scaling laws of tree structures |
| CN102496156A (zh) * | 2011-11-17 | 2012-06-13 | 西安电子科技大学 | 基于协同量子粒子群算法的医学图像分割方法 |
| CN104899862A (zh) * | 2015-04-01 | 2015-09-09 | 武汉工程大学 | 基于全局与局部阈值的视网膜血管分割算法 |
| WO2016177337A1 (en) * | 2015-05-05 | 2016-11-10 | Shanghai United Imaging Healthcare Co., Ltd. | System and method for image segmentation |
| CN104899876A (zh) * | 2015-05-18 | 2015-09-09 | 天津工业大学 | 一种基于自适应高斯差分的眼底图像血管分割方法 |
| CN106407917A (zh) * | 2016-09-05 | 2017-02-15 | 山东大学 | 基于动态尺度分配的视网膜血管提取方法及系统 |
| CN106651846A (zh) * | 2016-12-20 | 2017-05-10 | 中南大学湘雅医院 | 视网膜血管图像的分割方法 |
| CN106683061A (zh) * | 2017-01-05 | 2017-05-17 | 南京觅踪电子科技有限公司 | 一种基于修正的多尺度retinex算法对医疗图像进行增强的方法 |
| CN111351801A (zh) * | 2020-03-11 | 2020-06-30 | 春光线缆有限公司 | 一种电线电缆缺陷检测系统 |
| CN112053363A (zh) * | 2020-08-19 | 2020-12-08 | 苏州超云生命智能产业研究院有限公司 | 视网膜血管分割方法、装置及模型构建方法 |
Non-Patent Citations (2)
| Title |
|---|
| 基于粒子群优化算法的最大类间方差多阈值图像分割;周晓伟;葛永慧;;测绘科学(第02期);90-91+124 * |
| 粒子群优化的多阈值图像自分割算法;马培培;胡敏;;微计算机信息(第29期);正文1-3小节 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113269756A (zh) | 2021-08-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sheng et al. | Retinal vessel segmentation using minimum spanning superpixel tree detector | |
| CN110276356B (zh) | 基于r-cnn的眼底图像微动脉瘤识别方法 | |
| Singh et al. | Image processing based automatic diagnosis of glaucoma using wavelet features of segmented optic disc from fundus image | |
| Rehman et al. | Multi-parametric optic disc segmentation using superpixel based feature classification | |
| Memari et al. | Supervised retinal vessel segmentation from color fundus images based on matched filtering and AdaBoost classifier | |
| Kovács et al. | A self-calibrating approach for the segmentation of retinal vessels by template matching and contour reconstruction | |
| Wisaeng et al. | Exudates detection using morphology mean shift algorithm in retinal images | |
| Saha et al. | Spatial shape constrained fuzzy c-means (FCM) clustering for nucleus segmentation in pap smear images | |
| Hashemzadeh et al. | Retinal blood vessel extraction employing effective image features and combination of supervised and unsupervised machine learning methods | |
| CN108986106A (zh) | 面向青光眼临床诊断的视网膜血管自动分割方法 | |
| Khowaja et al. | A framework for retinal vessel segmentation from fundus images using hybrid feature set and hierarchical classification | |
| Harangi et al. | Automatic exudate detection using active contour model and regionwise classification | |
| Mahapatra et al. | A field of experts model for optic cup and disc segmentation from retinal fundus images | |
| Rekhi et al. | Automated classification of exudates from digital fundus images | |
| Giachetti et al. | Multiresolution localization and segmentation of the optical disc in fundus images using inpainted background and vessel information | |
| Wang et al. | Segmenting retinal vessels with revised top-bottom-hat transformation and flattening of minimum circumscribed ellipse | |
| Elbalaoui et al. | Automatic detection of blood vessel in retinal images | |
| Dharmawan et al. | A new optic disc segmentation method using a modified Dolph-Chebyshev matched filter | |
| Karkuzhali et al. | Robust intensity variation and inverse surface adaptive thresholding techniques for detection of optic disc and exudates in retinal fundus images | |
| CN113269756B (zh) | 基于多尺度匹配滤波与粒子群优化视网膜血管分割方法、装置 | |
| Mamilla et al. | Extraction of microaneurysms and hemorrhages from digital retinal images | |
| Fang et al. | Automatic segmentation of hard exudates in fundus images based on boosted soft segmentation | |
| Lyu et al. | HRED-net: high-resolution encoder-decoder network for fine-grained image segmentation | |
| Jana et al. | A semi-supervised approach for automatic detection and segmentation of optic disc from retinal fundus image | |
| Karmawat et al. | Glaucoma detection using fuzzy C-means optic cup segmentation and feature classification |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |