低通图像滤波器问题

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)生成低通滤波器并将其与傅里叶变换图像相乘(注意:低通滤波器需要在半径范围内为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

我很惊讶为什么这不会给你一个错误,但你正在执行矩阵乘法,而不是频域中的逐元素乘法.回想一下,频域中的滤波需要两个变换信号的逐元素乘法,以在空间域中执行等效的covnvolution操作.您只需更改乘法语句,使其与元素相当:

filteredmultip = lowpfilter.*fftdTestpatt;

此外,请确保在显示之前将图像转换回与原始图像相同的数据类型.uint8例如,如果您的图像是,您将要用im2uint8它来转换它.这主要用于显示目的和将图像写入文件.如果你像double在代码中那样离开它,显示图像将显示为大部分白色,因为显示double类型的图像假定值的范围来自[0,1].

作为旁注,我怀疑你没有得到错误的原因是因为你的图像和滤镜是方形的,所以矩阵乘法的尺寸是有效的.如果您决定在非方形图像上执行此操作,您肯定会收到错误.

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

您正在实现的是使用理想的低通滤波器进行滤波,因此当您低通滤波器时,您将看到振铃效应.之所以如此,是因为这直接来自信号处理理论.它需要空间域中的大量(或相当无限......)数量的正弦曲线来实现频域中的硬边缘.请记住,FFT是根据正弦曲线分解信号.当您使用此低通滤波器来过滤图像时,这可视化为重建图像中的振铃,因为硬边需要大量的正弦曲线(因此振铃的事实)来创建它们.

为了向您展示这些效果,让我们演示一个示例图像.我将使用摄影师图像作为图像处理工具箱的一部分.如果我使用半径Do = 40并运行你的代码(更正),这是原始图像,这是我在过滤后得到的:

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

回想一下2D高斯可以表示为:

^{资料来源:维基百科}

您只需循环遍历所有像素坐标并计算上述等式.我没有乘以前面的缩放因子,因为你真的不需要这样做.因此,您可以使用此代码生成高斯蒙板,这相当于您对两个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)图像.使用ifftshiftFFT过滤后,使用相关功能可以不影响图像:

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