二维匹配滤波器

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 对齐时匹配的滤波器响应最大（横截面分布是高斯分布，而不是容器本身）。

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

我的 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 保留高斯形状。那么使用Rotation matrix有什么好处呢？

输入图像是：

输出：

匹配过滤有助于提高 SNR，但也会放大背景噪声。 我可以使用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')

我得到以下输出：