CN107146211A - 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 - Google Patents
基于线扩散函数和双边滤波的视网膜血管图像降噪方法 Download PDFInfo
- Publication number
- CN107146211A CN107146211A CN201710428689.1A CN201710428689A CN107146211A CN 107146211 A CN107146211 A CN 107146211A CN 201710428689 A CN201710428689 A CN 201710428689A CN 107146211 A CN107146211 A CN 107146211A
- Authority
- CN
- China
- Prior art keywords
- mrow
- convolution kernel
- pixel
- msub
- filtering
- 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
- 238000001914 filtration Methods 0.000 title claims abstract description 83
- 238000000034 method Methods 0.000 title claims abstract description 57
- 230000002146 bilateral effect Effects 0.000 title claims abstract description 41
- 230000002792 vascular Effects 0.000 title claims abstract description 26
- 230000002207 retinal effect Effects 0.000 title claims abstract description 25
- 230000004256 retinal image Effects 0.000 claims abstract description 22
- 230000006870 function Effects 0.000 claims description 48
- 210000004204 blood vessel Anatomy 0.000 claims description 47
- 230000008569 process Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000010276 construction Methods 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 5
- 238000003860 storage Methods 0.000 claims description 5
- 238000010606 normalization Methods 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 abstract description 21
- 238000009826 distribution Methods 0.000 abstract description 8
- 230000000694 effects Effects 0.000 abstract description 7
- 230000009441 vascular protection Effects 0.000 abstract description 3
- 238000001514 detection method Methods 0.000 description 8
- 238000012545 processing Methods 0.000 description 6
- 239000008280 blood Substances 0.000 description 5
- 210000004369 blood Anatomy 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 210000001210 retinal vessel Anatomy 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 230000011218 segmentation Effects 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000004224 protection Effects 0.000 description 2
- 238000011946 reduction process Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 208000024172 Cardiovascular disease Diseases 0.000 description 1
- 239000006002 Pepper Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 208000026106 cerebrovascular disease Diseases 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 230000006378 damage Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- 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
-
- 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/20172—Image enhancement details
- G06T2207/20192—Edge enhancement; Edge preservation
-
- 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
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于线扩散函数和双边滤波的视网膜血管图像降噪方法，包括：利用线扩散函数为视网膜图像构建具有不同方向和尺度的类高斯卷积核，求以该像素为中心的方形邻域与m*n个卷积核的内积，获取与最大内积对应的卷积核的方向和尺度，利用获取到的卷积核的方向、尺度以及线扩散函数，为该像素计算最优空间卷积核，最后基于双边滤波的框架，将每个像素的最优空间卷积核与其灰度卷积核相结合，生成新的滤波卷积核，利用新的卷积核对视网膜图像进行滤波。本发明的技术方案保证了滤波的权重分配的准确性，提高了血管保护的效果，并且能够极大提高图像降噪算法的运行速度。
Description
技术领域
本发明涉及医学图像处理领域，尤其涉及一种基于线扩散函数和双边滤波的视网膜血管图像降噪的方法。
背景技术
视网膜图像的血管分析是诊断人体各种眼科疾病和心脑血管疾病的重要依据。近年来，随着社会医疗水平的提高，视网膜血管图像数据不断增长，为了弥补医生人工观察和经验诊断的效率低、主观性强等问题，借助计算机来自动分析视网膜血管图像逐渐成为了临床诊断的必要手段之一。然而，视网膜图像中的噪声给视网膜血管的自动分析带来了困难，具体表现为破坏了血管的形态结构，降低了血管探测和分割等自动分析算法的精度。所以图像降噪成为视网膜图像自动分析和诊断过程中的第一个步骤，也是最基础和关键的一步。
目前常用的图像降噪方法可以分为基于高斯滤波的方法和基于边缘保留的方法。基于高斯滤波的方法利用高斯函数对图像进行平滑处理，具有计算速度快、易于实现等优点，但是，由于高斯滤波无法保留图像的边缘结构，容易造成图像特征信息严重丢失。基于边缘保留的方法能够在降噪的过程中根据图像本身的灰度或梯度信息来保护图像的边缘特征，是目前运用最为广泛的降噪方法，然而，该方法往往只能用于处理一些简单的边缘结构，对于视网膜图像中对比度较低、局部形状不固定的毛细血管结构很难适用。
线扩散函数模型起源于1969年Rossmann等人对光学成像理论的研究，最初用来模拟一条线形的光源以高斯横截面的形式向两侧扩散。1989年Chaudhuri等人发现图像中血管的横截面结构可以被模拟成具有高斯横截面的线扩散函数，并将线扩散函数模型成功应用于探测视网膜图像中复杂的血管结构，提出了著名的匹配滤波方法用于探测视网膜图像中复杂的血管结构。
双边滤波是由Tomasi于1998年提出的一个新的具有保留边缘效果的图像滤波方法，由空间卷积核和灰度幅度卷积核两部分构成，其重要思想是利用像素之间的空间距离和灰度值的差异来分配平滑滤波的权重，通过参数调节到一个最佳权衡状态，在对图像进行降噪的同时极大程度上保留了图像的边缘信息。作为目前最流行的图像保边降噪方法，双边滤波已经在多个领域得到了广泛的应用。然而双边滤波缺少相应的函数来表示复杂的血管结构特征，因此在降噪后，视网膜图像中的较细的对比度较低的血管会丢失。
因此，目前需要本领域技术人员迫切解决的一个技术问题是：如何在降噪的同时保留血管结构，提高视网膜血管分析的准确度。
发明内容
为了解决上述问题，本发明结合了线扩散函数和双边滤波的优点，建立一种新的图像滤波模型，弥补了基于传统的空间卷积核的双边滤波方法在保留对比度较低的毛细血管结构方面的不足，并采用该模型实现视网膜血管图像分割之前的预处理。利用最优线扩散函数来探测血管结构，能够将图像滤波的权重在空间上自动根据血管的局部方向和尺度进行分配；基于双边滤波的框架能够将像素的空间关系与灰度关系相结合，进一步保证了滤波的权重分配的准确性，提高了血管保护的效果。
为了实现上述目的，本发明采用如下技术方案：
一种基于线扩散函数和双边滤波的视网膜血管图像降噪方法，包括以下步骤：
步骤1：选取一张视网膜图像，为图像中的血管定义m个尺度和n个方向，利用线扩散函数为该图像构建m*n个具有不同方向和尺度的类高斯卷积核；
步骤2：对于视网膜图像中的每个像素，分别求以该像素为中心的方形邻域与m*n个卷积核的内积，获取与最大内积对应的卷积核的方向和尺度；
步骤3：对于视网膜图像中的每个像素，利用获取到的卷积核的方向、尺度以及线扩散函数，为该像素计算最优空间卷积核；
步骤4：基于双边滤波的框架，将每个像素的最优空间卷积核与其灰度卷积核相结合，生成新的滤波卷积核，利用新的卷积核对视网膜图像进行滤波。
所述步骤1中用于构建m*n个类高斯卷积核的线扩散函数为：
且i∈(1,2,3...n)；j∈(1,2,3...m)；
其中K′i,j(x,y)表示在第i个方向和第j个尺度下的卷积核中，坐标为(x,y)的像素的灰度值，通过调节高斯函数的方差σ可以将每个方向的卷积核变换出n个尺度，最终得到m*n个类高斯卷积核；公式中的u计算公式为：
且
其中，k表示局部方形卷积核的长度，θi表示在第i个方向上的血管与y轴的夹角，[x,y]表示当血管方向与y轴夹角为0°时，类高斯卷积核中任意一点p的坐标，[u,v]表示点p经过卷积核旋转θi角后的坐标。
所述步骤2中获取方向和尺度的公式为：
且i,t∈(1,2,3...n)；j,s∈(1,2,3...m)
其中，p表示图像中的某个像素，F(p)是以p为中心的方形邻域矩阵，K′t,s是从m*n个卷积核中选出的与F(p)中的血管方向和尺度最匹配的卷积核。
所述步骤3中最优空间卷积核的计算公式为：
W′s(dvq)＝|K′t,s(x，y)|
其中，W′s代表最优空间卷积核，q代表以点p为中心的邻域中的任意一点；dvq定义像素的空间距离。
所述步骤4中生成新的滤波卷积核的公式为：
且
其中，I(q)和D(p)分别表示像素p的输入和输出值，W(·)是基本的高斯函数，dvq和fpq定义像素的空间距离和欧式距离，kp是归一化项；W′s、Wr分别是最优空间卷积核和灰度卷积核，W’sWr表示二者进行卷积运算。
进一步地，对双边滤波框架的优化内容包括以下至少一项：
(1)将步骤1获取到的大小和方向分别用1，2，3…n这样的数字来代替，并用矩阵来存储，将问题转化为计算机能够识别的语言；
(2)步骤4中最优空间卷积核与灰度卷积核相结合生成新的滤波核对图像滤波，该计算过程利用引导图像滤波的方法对其进行优化；
(3)利用频域转换的方法将卷积运算转化为频域里的乘积运算。
本发明的有益效果：
1、在处理效果上，本发明将基于线扩散函数的血管探测与基于双边滤波的图像滤波相结合，在对含有血管结构的图像进行降噪的同时，极大程度的保留了图像中的血管结构，包括对比度较低的细血管结构。
2、在运算速度上，首先，本方法中的线扩散函数空间卷积核是线性的，且能够独立于图像灰度进行计算，因此计算速度较快；其次，在滤波的过程中可以采用“引导图像”、频域转换等优化方法对其进行优化，进一步提高了运算速度；
3、在适用性和扩展性上，由于双边滤波框架能够适用于大部分情景，本方法又能够突出图像中的血管结构，因此适用于大部分含有血管结构的医学图像；此外，双边滤波的灰度卷积核可以进一步扩展成梯度、纹理等特征，因此本方法具有很好的扩展性。
附图说明
图1是本发明视网膜血管图像降噪方法的流程图；
图2是利用本发明的方法探测局部血管方向和尺度的一个实例；
图3是传统的基于高斯函数的空间卷积核权重分布示意图；
图4是本发明中基于线扩散函数的最优空间卷积核权重分布示意图；
图5是本发明中的新的卷积核生成过程示意图；
图6展示了本法明在手工制作的图像上的处理的一个实例；
图7展示了本发明在真实的视网膜血管图像上的处理的一个实例；
图8展示了本发明对某个经典血管分割算法的预处理的一个实例。
具体实施方式
术语解释：
(1)线扩散函数：评价成像系统成像质量的一个重要参数,线扩散函数的模拟分析与验证对成像仪的研制至关重要。
(2)双边滤波:一种可以对图像进行保边降噪的滤波器。
(3)图像降噪:英文名称是Image Denoising.是图像处理中的专业术语。现实中的数字图像在数字化和传输过程中常受到成像设备与外部环境噪声干扰等影响，称为含噪图像或噪声图像。减少数字图像中噪声的过程称为图像降噪,有时候又称为图像去噪。
(4)平滑滤波：低频增强的空间域滤波技术。它的目的有两类：一类是模糊；另一类是消除噪音。空间域的平滑滤波一般采用简单平均法进行，就是求邻近像元点的平均亮度值。邻域的大小与平滑的效果直接相关，邻域越大平滑的效果越好，但邻域过大，平滑会使边缘信息损失的越大，从而使输出的图像变得模糊，因此需合理选择邻域的大小。
下面结合附图与实施例对本发明作进一步说明。
实施例1：
图1是本发明视网膜血管图像降噪方法的流程图。
在整个方法实施之前需要输入和定义的关参数为：待降噪的视网膜血管图像I，卷积核的长k，以及空间卷积核和灰度差异卷积核的方差σs、σr，以及尺度和方向的数量m,n。经大量实验测试证明，对于大部分视网膜血管图像的最佳参数为：k＝9；σs＝2；σr＝0.04；m＝4；n＝12。
该方法可以总结为血管探测和图像滤波两部分：首先利用线扩散函数对视网膜图像中的血管结构进行探测，选出与图像局部血管结构相一致的最佳卷积核，再将该卷积核与灰度幅度卷积核相结合，构成新的图像滤波方法对图像进行降噪。
具体包括以下几个步骤：
步骤1：选取一张视网膜图像，为图像中的血管定义m个尺度和n个方向，利用线扩散函数公式为该图像构建m*n个具有不同方向和尺度的类高斯卷积核；
步骤1.1：设置图像中尺度和方向的总数m和n，输入视网膜血管图像I；
步骤2.2：利用公式(1)和(2)计算m个不同方向的类高斯卷积核，通过调节高斯函数的方差σ可以将每个方向的卷积核变换出n个尺度，最终得到m*n个类高斯卷积核。
其中K′i,j(x,y)表示在第i个方向和第j个尺度下的卷积核中，坐标为(x,y)的像素的灰度值。若令k表示局部方形卷积核的长度，θi表示在第i个方向上的血管与y轴的夹角，则公式中的u可以由坐标(x,y)利用下公式计算得出：
其中，[x,y]表示当血管方向与y轴夹角为0°(即血管的方向与y轴平行)时类高斯卷积核中任意一点p的坐标，则[u,v]表示点p经过卷积核旋转θi角后的坐标。
步骤2：对于视网膜图像中的每个像素，分别求以该像素为中心的方形邻域与m*n个卷积核的内积，比较大小选出最大内积，获取与最大内积对应的卷积核的方向和尺度，作为该像素的局部血管结构特征；
步骤2.1：对于图像中的每一个像素p和以p为中心的方形邻域F(p)
do
{
1.求出F(p)与m*n个类高斯卷积核的内积，即Doti,j＝K′i,j·F(p)；
2.比较求出的m*n个内积的大小，从中选出最大的内积Dot,并记录下该内积对应的滤波核大小和方向t、尺度s；
3.记录并更新像素的坐标。
}While(待处理的像素坐标值超出了图像的最大索引)；
其中，令t，s分别表示获取到的某个像素局部血管的方向和尺度，则获取t，s的公式为：
其中p表示图像中的某个像素，F(p)是以p为中心的方形邻域矩阵，K′t,s是从m*n个卷积核中选出的与F(p)中的血管方向和尺度最匹配的滤波核；
步骤2.2:输出并保存图像中每一个像素邻域对应的t和s。
步骤(1)和(2)所述的血管结构探测的过程如图2所示。
步骤3：对于视网膜图像中的每个像素，利用获取到的血管方向、尺度以及线扩散函数，为该像素计算一个最优空间卷积核；
所述步骤(3)中的最优空间核是本发明的核心。令Ws'代表最优空间卷积核，q代表以点p为中心的邻域中的任意一点，则q点的值W′s(dvq)可以由公式(3)中的K′t,s得出：
W′s(dvq)＝|K′t,s(x,y)|公式(4)
图3展示了传统空间卷积核权重分布的三维立体图，x、y轴确定的平面代表权重分布的位置，z轴表示相应位置上的权重值，权重值Ws由二维高斯函数公式计算得出：
其中dpq表示在二维平面中，像素点q到点p的距离，由该公式可以看出，传统的空间卷积核中每一点的权重值与该点到中心点的距离呈反比。
图4给出了一个本发明的步骤3中计算得到的最优线扩散函数空间卷积核的权重分布的示意图(3维立体图)，图中的卷积核利用公式(4)计算得出，其中与y轴的角度被设置为0，用来处理与y轴平行方向的血管。注意该空间核与图2中的传统的高斯函数的卷积核的区别，其权重值的大小与像素点到中心线的距离呈反比。
步骤4：基于双边滤波的框架，将每个像素的最优空间卷积核与其灰度卷积核相结合，生成新的滤波卷积核，利用新的卷积核中的权重来计算该像素的加权平均值。
最优空间卷积核与灰度卷积核相结合，生成新的滤波核的过程中的权重分布变化的二维和三维示意图。利用结合生成的卷积核对图像中进行平滑滤波处理，实质是利用局部血管的空间结构和灰度关系对像素进行加权平均。图像滤波降噪的具体过程如下：
步骤4.1：设置图像中尺度和方向的总数m和n，初始化各变量：输入视网膜血管图像I，高斯函数方差σs、σr，卷积核的长k。
步骤4.2：对于图像中的每一个像素p及以p为中心的方形邻域F(p)
do
{
1.利用步骤(2)中血管探测的结果t和s以及公式(1)、(2)、(3)和(4)，
计算该邻域中基于线扩散函数的最优空间卷积核，即：
2.利用传统的双边滤波的方法计算灰度卷积核，灰度卷积核公式为：
且fpq＝(I(p)-I(q))2
其中I(p)和I(q)分别是p、q两点的灰度值。
3.利用公式(4)将W′s和Wr结合起来，与F(p)求点积，并除以归一化项。
其中，I(q)和D(p)分别表示像素p的输入和输出值，W(·)是基本的高斯函数，dvq和fpq定义像素的空间距离和欧式距离，kp是归一化项；W′s、Wr分别是所述步骤(4)中的最优空间卷积核和灰度卷积核。
4.记录并更新像素的坐标。
}While(处理的像素坐标值超出了图像的最大索引)；
图像滤波过程中的最优空间核W′s和Wr相结合的方式是矩阵的卷积运算，可以根据卷积运算定理：
I(x,y)*w(x,y)<＝>F(u,v)H(u,v)公式(7)
其中I(x,y)、w(x,y)分别是图像像素的值和卷积核中的权重值，在频域中分别对应F(u,v)和H(u,v)。将卷积放在频域中处理能够使运算速度得到显著提高。
由公式(5)可以看出，本发明的核心----最优空间卷积核Ws'值的实质上是由滤波核中某一像素点到血管的中心直线的距离决定，而传统的高斯空间卷积核中的值则是由像素点到点的距离决定。由于血管的局部具有像素个数少、管状、曲率较小、方向不定等复杂的特征，且其局部几何形状可以近似为一条直线，因此相比传统的双边滤波，本发明能够在降噪的同时保留血管结构，提高视网膜血管分析的准确度。
在实施过程中为了提高算法的运行速度，本发明对传统的双边滤波框架基做了进一步优化，优化的内容主要包括：
(1)将步骤1获取到的大小和方向分别用1，2，3…n这样的数字来代替，并用矩阵来存储，将问题转化为计算机能够识别的语言。
(2)步骤4中最优空间卷积核与灰度卷积核相结合生成新的滤波核对图像滤波，其实质是求加权平均，该计算过程可以利用著名的“引导图像滤波”的方法对其进行优化。
(3)由于图像滤波采用的基本方式是卷积运算，可以利用频域转换的方法将卷积运算转化为频域里的乘积运算，进一步提高程序的运行速度。
本实施例提出了一种基于最优线扩散函数双边滤波的视网膜血管图像降噪方法。该方法可以总结为血管探测和图像滤波两部分：首先利用线扩散函数对视网膜图像中的血管结构进行探测，选出与图像局部血管结构相一致的卷积核，再将该卷积核与灰度幅度卷积核相结合，构成新的图像滤波方法对图像进行平滑滤波。其中，利用最优线扩散函数来探测血管结构，能够将图像滤波的权重在空间上自动根据血管的局部方向和尺度进行分配；基于双边滤波的框架能够将像素的空间关系与灰度关系相结合，进一步保证了滤波的权重分配的准确性，提高了血管保护的效果。此外，基于线扩散函数的空间核在计算上不需要遍历图像提取像素值，计算量大大减少，能够极大提高图像降噪算法的运行速度。
实施例2：
本实施例提供了一个本发明在手工制作的血管图像数据上的实例，如图6所示：
图6(a)是利用在Windows7操作系统下，利用“画图”工具绘制出的手工图像数据，在图像中用颜色暗于背景的直线来模拟视网膜图像的血管结构，血管和背景的像素值分别为125和180，并用Matlab图像处理软件在图像中添加了少量的高斯或椒盐噪声。由图7(b)和(c)可以看出，传统的双边滤波对图像噪声虽然又抑制作用，然而当血管与背景的对比度较低并且较细时，双边滤波极易破坏血管结构。采用本发明的最优线扩散函数双边滤波能够在很好的降噪的同时，保留并加深血管与背景的对比度，对血管结构起到较好的增强作用。
实施例3：
本实施例提供了一个本发明在真实的视网膜血管图像数据上的实例，如图7所示：
图7(a)是一张真实的视网膜血管图像，该图像是从著名的视网膜血管(数据库官方网站：http://www.isi.uu.nl/Research/Databases/DRIVE/)图像数据库DRIVE里面随机挑选的一张图像并取该图像的绿色通道，(b)、(c)、(d)分别是利用高斯滤波、双边滤波、基于最优线扩散函数的双边滤波的降噪处理结果。从图上可以看出，高斯滤波丢失大量的血管结构，相比之下，双边滤波能够保留血管结构，但图中较细的对比度较低的血管结构也遭到了破坏；本发明中的最优线扩散函数双边滤波的方法能够极大程度的保留不同尺度的血管，同时去除图像中的不均匀的噪声灰度。
实施例4：
本实施例提供了一个本发明在真实的视网膜血管图像数据上的实例，如图8所示：
图8(a)是图7(a)的彩色图像；图8(b)是该图像的绿色通道；图8(c)是该视网膜图像对应的手工分割结果，用于作为图像分割算法评估的标准；图8(d)是未经过降噪处理直接利用著名的血管探测算法Frangi对血管进行探测的结果；图8(e)是利用传统的双边滤波对其进行降噪然后探测的结果；图8(f)是利用本发明的方法对图像进行降噪，再用Frangi方法进行探测的结果。通过对比可以发现：在血管探测之前对图像进行降噪，有助于抑制图像噪声，然而即使是目前比较流行的双边滤波方法，在降噪后也破坏一些细的血管。而本方法相比另外两种，能够极大程度的保留血管结构，同时降低图像的噪声，提高后续血管探测的准确度。
本领域技术人员应该明白，上述的本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。
Claims (6)
1.一种基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，包括以下步骤：
步骤1：选取一张视网膜图像，为图像中的血管定义m个尺度和n个方向，利用线扩散函数为该图像构建m*n个具有不同方向和尺度的类高斯卷积核；
步骤2：对于视网膜图像中的每个像素，分别求以该像素为中心的方形邻域与m*n个卷积核的内积，获取与最大内积对应的卷积核的方向和尺度；
步骤3：对于视网膜图像中的每个像素，利用获取到的卷积核的方向、尺度以及线扩散函数，为该像素计算最优空间卷积核；
步骤4：基于双边滤波的框架，将每个像素的最优空间卷积核与其灰度卷积核相结合，生成新的滤波卷积核，利用新的卷积核对视网膜图像进行滤波。
2.如权利要求1所述的基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，所述步骤1中用于构建m*n个类高斯卷积核的线扩散函数为：
<mrow>
<msubsup>
<mi>K</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
</mrow>
<mo>,</mo>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<msup>
<mi>l</mi>
<mrow>
<mo>(</mo>
<mo>-</mo>
<mfrac>
<msup>
<mi>u</mi>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>j</mi>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</msup>
</mrow>
且i∈(1,2,3...n)；j∈(1,2,3...m)；
其中K′i,j(x,y)表示在第i个方向和第j个尺度下的卷积核中，坐标为(x,y)的像素的灰度值，通过调节高斯函数的方差σ可以将每个方向的卷积核变换出n个尺度，最终得到m*n个类高斯卷积核；公式中的u计算公式为：
<mrow>
<mo>&lsqb;</mo>
<mi>u</mi>
<mo>,</mo>
<mi>v</mi>
<mo>&rsqb;</mo>
<mo>=</mo>
<mo>&lsqb;</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>&rsqb;</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>cos&theta;</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>sin&theta;</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>sin&theta;</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>cos&theta;</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
且
其中，k表示局部方形卷积核的长度，θi表示在第i个方向上的血管与y轴的夹角，[x,y]表示当血管方向与y轴夹角为0°时类高斯卷积核中任意一点p的坐标，则[u,v]表示点p经过卷积核旋转θi角后的坐标。
3.如权利要求1所述的基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，所述步骤2中获取方向和尺度的公式为：
<mrow>
<msubsup>
<mi>K</mi>
<mrow>
<mi>t</mi>
<mo>,</mo>
<mi>s</mi>
</mrow>
<mo>,</mo>
</msubsup>
<mo>=</mo>
<munder>
<mrow>
<mi>arg</mi>
<mi>max</mi>
</mrow>
<msubsup>
<mi>K</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
</mrow>
<mo>,</mo>
</msubsup>
</munder>
<mrow>
<mo>(</mo>
<msubsup>
<mi>K</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
</mrow>
<mo>,</mo>
</msubsup>
<mo>&CenterDot;</mo>
<mi>F</mi>
<mo>(</mo>
<mi>p</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
且i,t∈(1,2,3...n)；j,s∈(1,2,3...m)
其中，p表示图像中的某个像素，F(p)是以p为中心的方形邻域矩阵，K′t,s是从m*n个卷积核中选出的与F(p)中的血管方向和尺度最匹配的卷积核。
4.如权利要求3所述的基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，所述步骤3中最优空间卷积核的计算公式为：
W′s(dvq)＝|K′t,s(x,y)|
其中，W′s代表最优空间卷积核，q代表以点p为中心的邻域中的任意一点；dvq定义像素的空间距离。
5.如权利要求4所述的基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，所述步骤4中生成新的滤波卷积核的公式为：
<mrow>
<mi>D</mi>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msubsup>
<mi>k</mi>
<mi>p</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>q</mi>
<mo>&Element;</mo>
<msub>
<mi>R</mi>
<mi>p</mi>
</msub>
</mrow>
</munder>
<msubsup>
<mi>W</mi>
<mi>s</mi>
<mo>,</mo>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mrow>
<mi>v</mi>
<mi>q</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<msub>
<mi>W</mi>
<mi>r</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>f</mi>
<mrow>
<mi>p</mi>
<mi>q</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
</mrow>
且
其中，I(q)和D(p)分别表示像素p的输入和输出值，W(·)是基本的高斯函数，dvq和fpq定义像素的空间距离和欧式距离，kp是归一化项；W′s、Wr分别是最优空间卷积核和灰度卷积核，W’sWr表示二者进行卷积运算。
6.如权利要求1所述的基于线扩散函数和双边滤波的视网膜血管图像降噪方法，其特征是，其中，对双边滤波框架的优化内容包括以下至少一项：
(1)将步骤1获取到的大小和方向分别用1，2，3…n这样的数字来代替，并用矩阵来存储，将问题转化为计算机能够识别的语言；
(2)步骤4中最优空间卷积核与灰度卷积核相结合生成新的滤波核对图像滤波，该计算过程利用引导图像滤波的方法对其进行优化；
(3)利用频域转换的方法将卷积运算转化为频域里的乘积运算。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710428689.1A CN107146211A (zh) | 2017-06-08 | 2017-06-08 | 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710428689.1A CN107146211A (zh) | 2017-06-08 | 2017-06-08 | 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107146211A true CN107146211A (zh) | 2017-09-08 |
Family
ID=59781107
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710428689.1A Pending CN107146211A (zh) | 2017-06-08 | 2017-06-08 | 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107146211A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112819733A (zh) * | 2021-01-29 | 2021-05-18 | 成都国科微电子有限公司 | 一种定向双边图像滤波方法及装置 |
| CN113189634A (zh) * | 2021-03-02 | 2021-07-30 | 四川新先达测控技术有限公司 | 一种类高斯成形方法 |
| CN114882596A (zh) * | 2022-07-08 | 2022-08-09 | 深圳市信润富联数字科技有限公司 | 行为预警方法、装置、电子设备及存储介质 |
| CN117041531A (zh) * | 2023-09-04 | 2023-11-10 | 无锡维凯科技有限公司 | 一种基于图像质量评估的手机摄像头聚焦检测方法和系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101763627A (zh) * | 2008-12-12 | 2010-06-30 | 新奥特（北京）视频技术有限公司 | 一种高斯模糊的实现方法和装置 |
| CN102226738A (zh) * | 2011-03-25 | 2011-10-26 | 宁波大学 | 一种红外玻璃非均匀性检测装置及检测方法 |
| CN104574374A (zh) * | 2014-12-23 | 2015-04-29 | 苏州大学 | 视网膜浆液性色素上皮层脱离的自动分割方法 |
| CN106651753A (zh) * | 2016-09-28 | 2017-05-10 | 沈阳东软医疗系统有限公司 | 提高ct图像显示效果的方法及装置 |
-
2017
- 2017-06-08 CN CN201710428689.1A patent/CN107146211A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101763627A (zh) * | 2008-12-12 | 2010-06-30 | 新奥特（北京）视频技术有限公司 | 一种高斯模糊的实现方法和装置 |
| CN102226738A (zh) * | 2011-03-25 | 2011-10-26 | 宁波大学 | 一种红外玻璃非均匀性检测装置及检测方法 |
| CN104574374A (zh) * | 2014-12-23 | 2015-04-29 | 苏州大学 | 视网膜浆液性色素上皮层脱离的自动分割方法 |
| CN106651753A (zh) * | 2016-09-28 | 2017-05-10 | 沈阳东软医疗系统有限公司 | 提高ct图像显示效果的方法及装置 |
Non-Patent Citations (1)
| Title |
|---|
| YUNLONG HE等: "Retinal Image Denoising via Bilateral Filter with a Spatial Kernel of Optimally Oriented Line Spread Function", 《COMPUTATIONAL AND MATHEMATICAL METHODS IN MEDICINE》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112819733A (zh) * | 2021-01-29 | 2021-05-18 | 成都国科微电子有限公司 | 一种定向双边图像滤波方法及装置 |
| CN112819733B (zh) * | 2021-01-29 | 2024-04-16 | 成都国科微电子有限公司 | 一种定向双边图像滤波方法及装置 |
| CN113189634A (zh) * | 2021-03-02 | 2021-07-30 | 四川新先达测控技术有限公司 | 一种类高斯成形方法 |
| CN114882596A (zh) * | 2022-07-08 | 2022-08-09 | 深圳市信润富联数字科技有限公司 | 行为预警方法、装置、电子设备及存储介质 |
| CN114882596B (zh) * | 2022-07-08 | 2022-11-15 | 深圳市信润富联数字科技有限公司 | 行为预警方法、装置、电子设备及存储介质 |
| CN117041531A (zh) * | 2023-09-04 | 2023-11-10 | 无锡维凯科技有限公司 | 一种基于图像质量评估的手机摄像头聚焦检测方法和系统 |
| CN117041531B (zh) * | 2023-09-04 | 2024-03-15 | 无锡维凯科技有限公司 | 一种基于图像质量评估的手机摄像头聚焦检测方法和系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107154023B (zh) | 基于生成对抗网络和亚像素卷积的人脸超分辨率重建方法 | |
| CN111524106B (zh) | 颅骨骨折检测和模型训练方法、装置、设备和存储介质 | |
| CN107492071A (zh) | 医学图像处理方法及设备 | |
| Tian et al. | Multi-path convolutional neural network in fundus segmentation of blood vessels | |
| Zhao et al. | Saliency driven vasculature segmentation with infinite perimeter active contour model | |
| CN109886179A (zh) | 基于Mask-RCNN的子宫颈细胞涂片的图像分割方法和系统 | |
| CN106339998A (zh) | 基于对比度金字塔变换的多聚焦图像融合方法 | |
| Li et al. | Robust retinal image enhancement via dual-tree complex wavelet transform and morphology-based method | |
| CN108764342B (zh) | 一种对于眼底图中视盘和视杯的语义分割方法 | |
| CN107146211A (zh) | 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 | |
| CN111275686B (zh) | 用于人工神经网络训练的医学图像数据的生成方法及装置 | |
| Wang et al. | Multimodal medical image fusion based on Gabor representation combination of multi-CNN and fuzzy neural network | |
| CN106504221B (zh) | 基于四元数小波变换上下文结构的医学图像融合方法 | |
| CN106709431A (zh) | 虹膜识别方法及装置 | |
| CN115375711A (zh) | 基于多尺度融合的全局上下文关注网络的图像分割方法 | |
| CN112150564A (zh) | 基于深度卷积神经网络医学图像融合算法 | |
| CN111724345A (zh) | 可自适应调节感受野大小的肺炎图片检定装置与方法 | |
| CN108305268A (zh) | 一种图像分割方法及装置 | |
| CN113313728B (zh) | 一种颅内动脉分割方法及系统 | |
| Zhang et al. | Medical image fusion based on quasi-cross bilateral filtering | |
| CN113362360B (zh) | 基于流体速度场的超声颈动脉斑块分割方法 | |
| Wang et al. | Automatic measurement of fetal head circumference using a novel GCN-assisted deep convolutional network | |
| Tuyet et al. | A Deep Bottleneck U-Net Combined with Saliency Map for Classifying Diabetic Retinopathy in Fundus Images. | |
| Zhao et al. | A multi-module medical image fusion method based on non-subsampled shear wave transformation and convolutional neural network | |
| Yue et al. | Generative adversarial network combined with SE-ResNet and dilated inception block for segmenting retinal vessels |
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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20170908 |
|
| RJ01 | Rejection of invention patent application after publication |