眼底图像血管增强与分割--（1）匹配滤波算法原理及实现

       视网膜血管分割是眼科计算机辅助诊断和大规模疾病筛查系统的基础，当眼器官发生视觉疾病的时候，网膜血管的直径、颜色和弯曲程度等就会出现异常，科医生可以据此作出诊断。医疗上对血管图像进行分割常用的方法有：基于血管跟踪的方法、基于匹配滤波的方法、基于形态学处理的方法、基于形变模型的方法和基于机器学习的方法等 。本篇将介绍基于匹配滤波算法。

匹配滤波器(matched filter)算法

原理部分：

      匹配滤波算法最早见于论文【1】《Detection of  Blood Vessels in Retinal Images Using Two-Dimensional Matched Filters

》中，后来又有许多针对该算法的不足进行改进的算法出现。

      匹配滤波器是指经过滤波后，使得滤波器输出端的信号瞬时功率与噪声平均功率的比值（即是信噪比SNR）最大，当有用信号与噪声同时进入滤波器时，有用信号在某一瞬间出现尖峰值，而噪声信号受到抑制。

      论文【1】对眼底视网膜上的血管特点进行分析，基于这些特点设计匹配滤波器，主要参考特点：血管的宽度在较小的范围内变动，血管壁的两条线是近似平行的，血管有方向，血管横切线上的灰度曲线如下，近似满足高斯分布：



        要理解论文中的公式推导，则需要对匹配滤波有更深层次的理解和掌握。首先，匹配滤波器也是一种在某一准则条件下的最优设计，其对应的最优的准则是输出信噪比（SNR）最大，而且前提条件是噪声信号为白噪声。论文得到的匹配滤波器的表达式为：H(f)=S*(f)，即匹配滤波器的频率响应是输入信号频率响应的共轭。从物理直观上理解：

（1）从幅频特性来看，匹配滤波器和输入信号的幅频特性一致。在信号越强的频率点，滤波器的放大倍数也越大；在信号越弱的频率点，滤波器的放大倍数越小。这就是信号处理中的“马太效应”，匹配滤波器是让有用信号尽可能通过。白噪声的功率谱在各个频率点都一样，这种情况下，有用信号尽可能通过，噪声则相反。

（2）从相频特性上看，匹配滤波器的相频特性和输入信号则相反。通过匹配滤波器后，信号的相位为0，正好能实现信号时域上的相干叠加。而噪声的相位是随机的，只能实现非相干叠加。这样在时域上保证了输出信噪比的最大。

       在信号的幅频特性与相频特性中，幅频特性表征了频率特性，频特性表征了时间特性。匹配滤波器无论是从时域还是从频域，都充分保证了信号尽可能大的通过，噪声尽可能小的通过。

       回到论文【1】中的描述，将血管想象成一小段一小段的平行区域的组合，设定长度为L，宽度为3σ，然后用一个高斯曲线来模拟血管横切面上的灰度曲线， 匹配滤波器如下：



原文中讲到：“When the concept of matched filter is extended to two- dimensional images, it must be appreciated that a vessel may be oriented at any angle θ (0 <=θ<=pi ). The matched filter s ( t ) will have its peak response only when
it is aligned at an angle θ+- pi / 2 . ” （不是很理解为什么，有知道的请告知？）

由于眼底中血管是有方向的，并且血管的方向是随机伸展的，而滤波器在垂直于血管的方向上出现峰值。所以设计滤波器从0度到180度每15度旋转一次，得到12个方向上的滤波器，然后分别进行卷积，选取最大响应的一个作为最终的响应输出。

旋转矩阵为：


，

算法实现中用到的核函数及掩码设置参见原文，怕翻的不精准（捂脸~）


，

算法实现：

参考1中的代码试过效果不是很好，可能是参数需要调试，这里附录下：

%% matched filter2，测试程序
clc,close all,clear all;
img=imread('../image/1.jpg');
if length(size(img))==3
img=rgb2gray(img);
end
[g,bg]=matchedFilter2(img);
function [g,bg]=matchedFilter2(f)

f=double(f);
subplot(1,2,1);
imshow(uint8(f));
title('origin image');
% mean filter
f=medfilt2(f,[5,5]);
% f=medfilt2(f,[21 1]);
% f=medfilt2(f,[1,7]);
% 参数
os=12; % 角度的个数
sigma=2;
tim=3;
L=9;
t=180; % 全局阈值,需要多次尝试

thetas=0:(os-1);
thetas=thetas.*(180/os);
N1=-tim*sigma:tim*sigma;
N1=-exp(-(N1.^2)/(2*sigma*sigma));
N=repmat(N1,[2*floor(L/2)+1,1]);
r2=floor(L/2);
c2=floor(tim*sigma);
[m,n]=size(f);
RNs=cell(1,os); % rotated kernals
MFRs=cell(1,os); % filtered images
g1=f;

% matched filter
for i=1:os
theta=thetas(i);
RN=imrotate(N,theta);
%去掉多余的0行和零列
RN=RN(:,any(RN));
RN=RN(any(RN'),:);
meanN=mean2(RN);
RN=RN-meanN;
RNs{1,i}=RN;
MFRs{1,i}=imfilter(f,RN,'conv','symmetric');
end

% get the max response
g=MFRs{1,1};
for j=2:os
g=max(g,MFRs{1,j});
end
bg=g<t;

subplot(1,2,2);
imshow(bg);
title('filterd image');
end
上述代码中的三个参数sigma,L,T的选择需要针对不同的应用去尝试，如果是视网膜血管图像，参数可以设置为2,9,3。全局阈值t的选择，可以多试几次找最好的。测试结果如下：



我重点测试了参考3的算法实现，测试代码如下，整个工程文件见后面链接：

%% 测试函数，匹配滤波算法血管分割
clc,close all,clear all;
sigma=2;
yLength=10;
direction_number=12;
img=imread('../image/1.bmp');
if length(size(img))==3
img=rgb2gray(img);
end
MF = MatchFilter(img, sigma, yLength,direction_number);
mask=[0 0 0 0 0;
0 1 1 1 0;
0 1 1 1 0;
0 1 1 1 0;
0 0 0 0 0;];
MF(mask==0) = 0;
MF = normalize(double(MF));
% Adding to features
features = MF;
subplot(1,2,1);
imshow(img);title('origin image');
subplot(1,2,2);
imshow(features);title('filted image');

效果图如下：





全部工程文件下载链接

算法优缺点：

滤波之后图像中的血管结构得到增强，但是该方法可能会丢失血管分叉点和细小的血管；

参考：
http://www.cnblogs.com/naniJser/archive/2012/12/09/2810551.html http://blog.csdn.net/chiooo/article/details/43564069
http://blog.csdn.net/u013288466/article/details/66473284