低通图像滤波器的问题
Issue with lowpass image filter
我需要在图像上应用一个简单的低通滤波器，但这需要在频域中完成。然而最终的结果是一个非常嘈杂的图像，表明我的过程在某处不正确。
我的思路是
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
如有任何帮助或意见，我们将不胜感激。
我很惊讶为什么这没有给你一个错误，但你在频域中执行 矩阵乘法 ，而不是 element-wise 乘法。回想一下，频域中的滤波需要两个 t运行 变形信号的 element-wise 相乘，以在空间域中执行等效的卷积运算。您只需更改乘法语句，使其等效于 element-wise：
filteredmultip = lowpfilter.*fftdTestpatt;
此外，请确保在显示之前将图像转换回与原始图像相同的数据类型。例如，如果您的图像是 uint8，您将需要使用 im2uint8 来转换它。这主要对于显示目的和将图像写入文件很重要。如果您像在代码中所做的那样将其保留为 double，则显示图像将被可视化为大部分白色，因为显示 double 类型的图像假定值的 运行ge 来自 [0,1].
附带说明一下，我怀疑您没有收到错误的原因是因为您的图像和滤镜是正方形的，因此矩阵乘法的维度是有效的。如果您决定对 non-square 图像执行此操作，您肯定会遇到错误。
警告 - 使用理想的低通滤波器
您正在实施的是使用理想 低通滤波器进行过滤，因此您将在低通滤波器时看到振铃效应。原因是因为这直接来自信号处理理论。它需要空间域中的大量（或者更确切地说是无限的...）正弦曲线来实现频域中的硬边缘。请记住，FFT 是根据正弦波分解信号。当您使用此低通滤波器过滤图像时，这在重建图像中被可视化为振铃，因为硬边缘需要大量正弦波求和（因此是振铃）才能创建它们。
为了向您展示这些效果，让我们以示例图像进行演示。我将使用属于图像处理工具箱一部分的摄影师图像。如果我使用 Do = 40 的半径和 运行 你的代码（更正），这是原始图像，这是我过滤后得到的：
如您所见，这很糟糕。振铃来自频域中滤波器的硬边缘。当您远离中心时，人们倾向于使用幅度逐渐减小的方式。一个很好的例子是 高斯模糊 。相反，您要做的是定义一个标准差，然后创建一个与代表 2D 高斯图像的图像大小相同的掩码。
回想一下，二维高斯可以表示为：
来源：Wikipedia
您只需遍历所有像素坐标并计算上述等式。我没有乘以前面的比例因子，因为你真的不需要这样做。因此，您可以使用此代码生成一个高斯掩码，这与您的两个 for 循环所具有的等效，但它的矢量化程度更高。我定义了一个二维坐标网格，其跨度与您的图像大小相同，然后 运行 通过图像中每个位置的方程式。我们当然需要 centre 内核，因此您必须像使用之前的低通滤波代码那样在图像中间进行偏移。我还决定使用您的 Do 变量作为标准偏差。
Do = 40;
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = exp(-((X - (tpa2/2)).^2 + (Y - (tpa1/2)).^2) / (2*Do*Do));
所以我们现在得到这个作为过滤器响应：
...当我过滤时，我得到：
与原图比较：
如您所见，模糊效果好多了。没有铃声。因此，当您过滤图像时，请确保避免过滤器中的硬边，因为这会在空间域中产生振铃效应。
更多建议
我还有一些建议可以让您更快地完成此代码运行。
避免使用循环使图像居中
正如您从信号处理理论中已经知道的那样，如果您 pre-multiply 图像中的像素值 (-1)^(i+j) 其中 (i,j) 是您要过滤的像素的空间位置，这会使您的图像在频域中居中。我实际上会避免使用循环来执行此操作并首先使用 FFT 然后 将图像居中。您可以使用名为 fftshift 的函数来为您执行居中操作。首先对图像进行 FFT，然后在以下之后调用此函数：
fftdTestpatt = fft2(dTestpatt);
fftdTestpatt = fftshift(fftdTestpatt); % Centre the FFT
避免使用循环生成过滤器
我也会避免使用循环来生成您的过滤器。具体来说，对于具有理想滤波器的代码，请将代码替换为遵循与高斯滤波器相同逻辑的代码。我们还可以删除 sqrt 操作并确保您正在与 avo 的半径的平方进行比较d 任何不必要的计算：
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = (X - tpa1/2).^2 + (Y - tpa2/2).^2 <= Do^2;
lowpfilter = double(lowpfilter);
特别注意上面代码中的element-wise幂运算(.^)。最后一条语句很重要，因为我们需要将结果转换回 double，因为首先生成过滤器实际上会给你一个 logical 类型作为结果。我们需要 double 类型来确保 FFT 的精度得到尊重。
避免在过滤后循环使图像不居中
完成过滤后，您再次乘以 (-1)^(i+j) 以使图像不居中。在使用 FFT 过滤后，使用相关的 ifftshift 函数使图像不居中：
filteredmultip = ifftshift(filteredmultip); % Uncentre first
filteredimage = real(ifft2(filteredmultip));
我需要在图像上应用一个简单的低通滤波器，但这需要在频域中完成。然而最终的结果是一个非常嘈杂的图像，表明我的过程在某处不正确。
我的思路是
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
如有任何帮助或意见，我们将不胜感激。
我很惊讶为什么这没有给你一个错误，但你在频域中执行 矩阵乘法 ，而不是 element-wise 乘法。回想一下，频域中的滤波需要两个 t运行 变形信号的 element-wise 相乘，以在空间域中执行等效的卷积运算。您只需更改乘法语句，使其等效于 element-wise：
filteredmultip = lowpfilter.*fftdTestpatt;
此外，请确保在显示之前将图像转换回与原始图像相同的数据类型。例如，如果您的图像是 uint8，您将需要使用 im2uint8 来转换它。这主要对于显示目的和将图像写入文件很重要。如果您像在代码中所做的那样将其保留为 double，则显示图像将被可视化为大部分白色，因为显示 double 类型的图像假定值的 运行ge 来自 [0,1].
附带说明一下，我怀疑您没有收到错误的原因是因为您的图像和滤镜是正方形的，因此矩阵乘法的维度是有效的。如果您决定对 non-square 图像执行此操作，您肯定会遇到错误。
警告 - 使用理想的低通滤波器
您正在实施的是使用理想 低通滤波器进行过滤，因此您将在低通滤波器时看到振铃效应。原因是因为这直接来自信号处理理论。它需要空间域中的大量（或者更确切地说是无限的...）正弦曲线来实现频域中的硬边缘。请记住，FFT 是根据正弦波分解信号。当您使用此低通滤波器过滤图像时，这在重建图像中被可视化为振铃，因为硬边缘需要大量正弦波求和（因此是振铃）才能创建它们。
为了向您展示这些效果，让我们以示例图像进行演示。我将使用属于图像处理工具箱一部分的摄影师图像。如果我使用 Do = 40 的半径和 运行 你的代码（更正），这是原始图像，这是我过滤后得到的：
如您所见，这很糟糕。振铃来自频域中滤波器的硬边缘。当您远离中心时，人们倾向于使用幅度逐渐减小的方式。一个很好的例子是 高斯模糊 。相反，您要做的是定义一个标准差，然后创建一个与代表 2D 高斯图像的图像大小相同的掩码。
回想一下，二维高斯可以表示为：
来源：Wikipedia
您只需遍历所有像素坐标并计算上述等式。我没有乘以前面的比例因子，因为你真的不需要这样做。因此，您可以使用此代码生成一个高斯掩码，这与您的两个 for 循环所具有的等效，但它的矢量化程度更高。我定义了一个二维坐标网格，其跨度与您的图像大小相同，然后 运行 通过图像中每个位置的方程式。我们当然需要 centre 内核，因此您必须像使用之前的低通滤波代码那样在图像中间进行偏移。我还决定使用您的 Do 变量作为标准偏差。
Do = 40;
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = exp(-((X - (tpa2/2)).^2 + (Y - (tpa1/2)).^2) / (2*Do*Do));
所以我们现在得到这个作为过滤器响应：
...当我过滤时，我得到：
与原图比较：
如您所见，模糊效果好多了。没有铃声。因此，当您过滤图像时，请确保避免过滤器中的硬边，因为这会在空间域中产生振铃效应。
更多建议
我还有一些建议可以让您更快地完成此代码运行。
避免使用循环使图像居中
正如您从信号处理理论中已经知道的那样，如果您 pre-multiply 图像中的像素值 (-1)^(i+j) 其中 (i,j) 是您要过滤的像素的空间位置，这会使您的图像在频域中居中。我实际上会避免使用循环来执行此操作并首先使用 FFT 然后 将图像居中。您可以使用名为 fftshift 的函数来为您执行居中操作。首先对图像进行 FFT，然后在以下之后调用此函数：
fftdTestpatt = fft2(dTestpatt);
fftdTestpatt = fftshift(fftdTestpatt); % Centre the FFT
避免使用循环生成过滤器
我也会避免使用循环来生成您的过滤器。具体来说，对于具有理想滤波器的代码，请将代码替换为遵循与高斯滤波器相同逻辑的代码。我们还可以删除 sqrt 操作并确保您正在与 avo 的半径的平方进行比较d 任何不必要的计算：
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = (X - tpa1/2).^2 + (Y - tpa2/2).^2 <= Do^2;
lowpfilter = double(lowpfilter);
特别注意上面代码中的element-wise幂运算(.^)。最后一条语句很重要，因为我们需要将结果转换回 double，因为首先生成过滤器实际上会给你一个 logical 类型作为结果。我们需要 double 类型来确保 FFT 的精度得到尊重。
避免在过滤后循环使图像不居中
完成过滤后，您再次乘以 (-1)^(i+j) 以使图像不居中。在使用 FFT 过滤后，使用相关的 ifftshift 函数使图像不居中：
filteredmultip = ifftshift(filteredmultip); % Uncentre first
filteredimage = real(ifft2(filteredmultip));