二维匹配滤波器

Question

我想根据Chaudhuri等人的论文“使用二维匹配滤波器检测视网膜图像中的血管”实现用于血管提取的二维匹配滤波器，IEEE Trans。医学影像，1989 (there's a PDF on the author's web site)。

简单地说，血管的横截面具有高斯分布，因此我想使用高斯匹配滤波器来提高信噪比。这样的内核可以在数学上表示为：

K(x,y) = -exp(-x^2/2*sigma^2)  for |x|<3*sigma,  |y|<L/2

L这里是固定方向的容器的长度。实验性的sigma=1.5和L = 7。

这部分的MATLAB代码是：

s = 1.5; %sigma
t = -3*s:3*s;
theta=0:15:165; %different rotations
%one dimensional kernel
x = 1/sqrt(6*s)*exp(-t.^2/(2*s.^2));

L=7;
%two dimensional gaussian kernel
x2 = repmat(x,L,1);

考虑此滤波器对属于背景视网膜的像素的响应。假设背景具有恒定的强度和零均值加性高斯白噪声，则滤波器输出的期望值理想情况下应为零。因此，通过从函数本身减去s(t)的平均值来修改卷积核。内核的平均值被确定为：m = Sum(K(x,y))/(number of points)。

因此，此算法中使用的卷积掩码由：K(x, y) = K(x,y) - m给出。

我的MATLAB代码：

m = sum(x2(:))/(size(x2,1)*size(x2,2));
x2 = x2-m; 

血管可以以任何角度0<theta<180定位，当它以theta+- 90对齐时，匹配的滤波器响应最大(截面分布是高斯的，而不是血管本身)。

因此，我们需要以15度的增量旋转匹配滤波器12次。

我的MATLAB代码附在这里，但我没有得到理想的结果。任何帮助都是非常感谢的。

%apply rotated matched filter on image
r = {};
for k = 1:12
    x3=imrotate(x2,theta(k),'crop');%figure;imagesc(x3);colormap gray;   
    r{k}=conv2(img,x3);
end
w=[];h = zeros(584,565);
for i = 1:565
    for j = 1:584
        for k = 1:32
            w= [w ,r{k}(j,i)];
        end
        h(j,i)=max(abs(w));
        w = [];
    end
end
%show result
figure('Name','after matched filter');imagesc(h);colormap gray

对于旋转，我使用了imrotate，这对我来说似乎更合理，但在本文中它是不同的:假设p=[x,y]是内核中的一个离散点。为了计算旋转核中的系数，我们使用了[u,v] = p*Rotation_Matrix。

Rotation_Matrix=[cos(theta),sin(theta);-sin(theta),cos(theta)]

内核是：

K(x,y) = -exp(-u^2/2*s^2)

但是新的内核不再具有高斯形状。使用imrotate可以保留高斯形状。那么使用旋转矩阵的好处是什么呢？

输入图像为：

​

​

输出：

​

​

匹配滤波有助于提高信噪比，但也会放大背景噪声。我可以使用imrotate来旋转内核吗？我的主要问题是旋转矩阵，为什么以及什么是正确的代码来实现它。

Answer 1

从每次旋转的解析表达式而不是使用imrotate构建过滤器的原因是，过滤器的范围不是圆形的，因此旋转会带来“新的”像素值，并将其他一些像素推出内核。此外，旋转按此处构建的核(沿一个方向平滑过渡，沿另一个维度的步长边缘)需要沿每个维度使用不同的插值方法，而imrotate无法做到这一点。由此产生的旋转内核总是错误的。

在显示您创建的内核和两个旋转版本时，可以很容易地看到这两个问题：

​

​

这种显示带来了一个额外的问题:内核没有以像素为中心，导致它将输出移动了半个像素。

另请注意，当减去平均值时，重要的是仅在滤波器的原始域上计算此平均值，并且用于将此域填充到矩形形状的任何零都保持为零(这些零不应变为负数)。

旋转的内核可以按如下方式构造：

m = max(ceil(3*s),(L-1)/2);
[x,y] = meshgrid(-m:m,-m:m); % non-rotated coordinate system, contains (0,0)
t = pi/6;                    % angle in radian
u = cos(t)*x - sin(t)*y;     % rotated coordinate system
v = sin(t)*x + cos(t)*y;     % rotated coordinate system
N = (abs(u) <= 3*s) & (abs(v) <= L/2);   % domain
k = exp(-u.^2/(2*s.^2));     % kernel
k = k - mean(k(N));
k(~N) = 0;                   % set kernel outside of domain to 0

这是上述示例中使用的三次旋转的结果(内核边缘周围的灰色对应于值0，黑色像素具有负值)：

​

​

另一个问题是，您将conv2与默认的'full'输出形状一起使用，在这里您应该使用'same'，以便过滤器的输出与输入匹配。

请注意，在计算每个滤波器响应时，计算最大值要容易得多，而不是计算所有滤波器响应，然后再计算最大值。以上所有代码都指向以下代码：

img = im2double(rgb2gray(img));

s = 1.5; %sigma
L = 7;
theta = 0:15:165; %different rotations

out = zeros(size(img));

m = max(ceil(3*s),(L-1)/2);
[x,y] = meshgrid(-m:m,-m:m); % non-rotated coordinate system, contains (0,0)
for t = theta
   t = t / 180 * pi;        % angle in radian
   u = cos(t)*x - sin(t)*y; % rotated coordinate system
   v = sin(t)*x + cos(t)*y; % rotated coordinate system
   N = (abs(u) <= 3*s) & (abs(v) <= L/2); % domain
   k = exp(-u.^2/(2*s.^2)); % kernel
   k = k - mean(k(N));
   k(~N) = 0;               % set kernel outside of domain to 0

   res = conv2(img,k,'same');
   out = max(out,res);
end

out = out/max(out(:)); % force output to be in [0,1] interval that MATLAB likes
imwrite(out,'so_result.png')

我得到以下输出：

​

​