matlab - 低通图像滤镜问题

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）取反2D FFT的实部，恢复对中。

filteredimagecent = real(ifft2(filteredmultip));


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

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

如有任何帮助或意见，将不胜感激。

Answer 1

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

filteredmultip = lowpfilter.*fftdTestpatt;

另外，请确保在显示前将图像转换回与原始图像相同的数据类型例如，如果您的图像是uint8，则需要使用im2uint8来转换它这对于显示和将图像写入文件非常重要如果您像在代码中所做的那样将其保留为double，则显示图像将显示为主要为白色，因为显示double类型的图像假定值的范围是从[0,1]开始的。
顺便说一下，我怀疑你没有得到错误的原因是因为你的图像和滤波器是正方形的，所以矩阵乘法的维数是有效的如果你决定在非正方形图像上这样做，你肯定会得到一个错误。
警告-使用理想低通滤波器
你要实现的是用一个理想的低通滤波器进行滤波，所以当你使用低通滤波器时，你会看到振铃效应原因是这直接来自信号处理理论它需要在空间域中有大量（或者更确切地说是无限）的正弦波才能在频域中实现硬边记住FFT是信号的正弦分解当你使用这个低通滤波器来过滤你的图像时，这在重建图像中被视为振铃，因为硬边需要一个大的正弦和（因此振铃的事实）来创建它们。
为了向您展示这些效果，让我们演示一个示例图像我将使用摄像图像，这是图像处理工具箱的一部分如果我使用半径Do = 40并运行代码（已更正），这是原始图像，这是过滤后得到的结果：


如你所见，那很糟糕铃声来自你的滤波器在频域的硬边当你远离中心时，人们倾向于使用一个更为渐进的幅度递减高斯模糊就是一个很好的例子你要做的是定义一个标准差，然后创建一个与你的图像大小相同的遮罩，它代表一个二维高斯。
回想一下，二维高斯可以表示为：
来源：Wikipedia
你只需循环所有的像素坐标并计算上面的方程我没有乘以前面的比例因子，因为你真的不需要这么做因此，您可以使用此代码生成一个高斯掩码，它相当于您的两个for循环所具有的功能，但这样做会更加矢量化我定义了一个二维坐标网格，它的跨度与图像的大小相同，然后运行图像中每个位置的方程式当然，我们需要将内核居中，这样您就必须像处理以前的低通滤波代码那样，偏移图像的中间我还决定用你的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));