低通图像滤波器问题

Question

我需要在图像上应用一个简单的低通滤波器，但这需要在频域中完成。然而最终的结果是一个非常嘈杂的图像，这表明我的过程在某处不正确。

我的思路一直是

1) 将图像居中 (Testpatt) 并应用 2D FFT

tpa1 = size(Testpatt, 1); tpa2 = size(Testpatt, 2);
dTestpatt = double(Testpatt);


for i = 1:tpa1
    for j = 1:tpa2

        dTestpatt(i,j) = (dTestpatt(i,j))*(-1)^(i + j);
    end
end

fftdTestpatt = fft2(dTestpatt);

2) 生成低通滤波器，并与傅里叶变换后的图像相乘（注意：低通滤波器需要是半径Do内为1的圆）

lowpfilter = zeros(tpa1, tpa2);
Do = 120;

for i = 1:tpa1
    for j = 1:tpa2


       if sqrt((i - tpa1/2)^2 + (j - tpa2/2)^2) <= Do
           lowpfilter(i,j) = 1;
       end

   end
end


filteredmultip = lowpfilter*fftdTestpatt;

3) 取逆二维 FFT 的实部并恢复居中。

filteredimagecent = real(ifft2(filteredmultip));


for i = 1:tpa1
    for j = 1:tpa2

        filteredimage(i,j) = (filteredimagecent(i,j))*(-1)^(i + j);
    end
end

任何帮助或 cmets 将不胜感激。

Answer 1

我很惊讶为什么这没有给你一个错误，但你正在执行 矩阵乘法，而不是频域中的元素乘法。回想一下，频域中的滤波需要两个变换信号的元素相乘，以在空间域中执行等效的卷积运算。您只需要更改乘法语句，使其在元素方面等效：

filteredmultip = lowpfilter.*fftdTestpatt;

此外，请确保在显示之前将图像转换回与原始图像相同的数据类型。例如，如果您的图像是uint8，您将需要使用im2uint8 来转换它。这对于显示目的和将图像写入文件很重要。如果您将其保留为double，就像您在代码中所做的那样，则显示图像将大部分显示为白色，因为显示双精度类型的图像假定值的范围来自[0,1]。

作为旁注，我怀疑您没有收到错误的原因是因为您的图像和过滤器是方形的，因此矩阵乘法的尺寸是有效的。如果您决定对非方形图像执行此操作，则肯定会出错。

警告 - 使用理想的低通滤波器

您正在实施的是使用理想低通滤波器进行滤波，因此当您使用低通滤波器时会出现振铃效应。原因是因为这直接来自信号处理理论。它需要空间域中的大量（或相当无限......）数量的正弦曲线来实现频域中的硬边缘。请记住，FFT 是根据正弦曲线对信号进行分解。当您使用此低通滤波器过滤您的图像时，这被可视化为重建图像中的振铃，因为硬边缘需要大量的正弦曲线总和（因此是振铃）来创建它们。

为了向您展示这些效果，让我们以示例图像进行演示。我将使用作为图像处理工具箱一部分的摄影师图像。如果我使用Do = 40 的半径并运行您的代码（已更正），这是原始图像，这是我过滤后得到的：

如您所见，这非常糟糕。振铃来自频域中滤波器的硬边缘。当您远离中心时，人们倾向于使用更逐渐减小的幅度。 高斯模糊就是一个很好的例子。相反，您要做的是定义一个标准偏差，然后创建一个与您的图像大小相同的掩码，以表示 2D 高斯。

回想一下，二维高斯可以表示为：

^{来源：Wikipedia}

您只需遍历所有像素坐标并计算上述等式。我没有乘以前面的比例因子，因为你真的不需要这样做。因此，您可以使用此代码生成一个高斯掩码，这相当于您使用两个 for 循环所拥有的，但它的矢量化程度更高。我定义了一个 2D 坐标网格，该网格与您的图像大小相同，然后针对图像中的每个位置运行方程。我们当然需要将内核居中，因此您必须在图像的中间偏移，就像您在之前的低通滤波代码中所做的那样。我还决定使用您的 Do 变量作为标准差。

Do = 40;
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = exp(-((X - (tpa2/2)).^2 + (Y - (tpa1/2)).^2) / (2*Do*Do));

所以我们现在得到这个作为过滤器响应：

...当我过滤时，我得到：

与原图对比：

如您所见，模糊效果要好得多。没有响铃。因此，当您过滤图像时，请确保避免过滤器中出现硬边缘，因为这会在空间域中产生振铃伪影。

更多建议

我还有一些建议可以让你快速运行这段代码。

避免使用循环来居中图像

正如您从信号处理理论中已经知道的那样，如果您将图像中的像素值预乘以 (-1)^(i+j)，其中 (i,j) 是您要过滤的像素的空间位置，这会使您的图像以频率为中心领域。我实际上会避免使用循环来执行此操作，并首先进行 FFT 然后 居中图像。您可以使用名为fftshift 的函数为您执行居中。首先对你的图像进行 FFT，然后调用这个函数：

fftdTestpatt = fft2(dTestpatt);
fftdTestpatt = fftshift(fftdTestpatt); % Centre the FFT

避免使用循环来生成过滤器

我也会避免使用循环来生成您的过滤器。具体来说，对于具有理想过滤器的代码，将您的代码替换为遵循与高斯过滤器相同的逻辑的代码。我们还可以删除sqrt 操作，并确保您与半径的平方进行比较以避免任何不必要的计算：

[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = (X - tpa1/2).^2 + (Y - tpa2/2).^2 <= Do^2;
lowpfilter = double(lowpfilter);

请特别注意上述代码中的元素级幂运算 (.^)。最后一条语句很重要，因为我们需要将结果转换回double，因为首先生成过滤器实际上会给你一个logical 类型作为结果。我们需要 double 类型来确保 FFT 的精度得到尊重。

过滤后避免循环使图像不居中

完成过滤后，再次乘以 (-1)^(i+j) 以使图像不居中。使用 FFT 过滤后，使用相关的ifftshift 函数使图像不居中：

filteredmultip = ifftshift(filteredmultip); % Uncentre first
filteredimage = real(ifft2(filteredmultip));