高通巴特沃斯滤波器在MATLAB中对图像进行滤波 [英] High Pass Butterworth Filter on images in MATLAB

问题描述

我需要在MATLAB中实现一个高通Butterworth滤波器，用于图像滤波。我已经实现了一个，但看起来它不起作用。这是我写的代码。谁能告诉我什么是错的？

I need to implement a high pass Butterworth filter in MATLAB for the purposes of image filtering. I have implemented one but it looks like it doesn't work. Here is the code I have written. Can anyone tell me what is wrong?

n=1;
d=50;
A=1.5;
im=imread('imagex.jpg');
h=size(im,1);
w=size(im,2);
[x y]=meshgrid(-floor(w/2):floor(w-1/2),-floor(h/2):floor(h-1/2));
hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));
image_2Dfilter=fftshift(fft2(im));
Image_butterworth=image_2Dfilter;
imshow(Image_butterworth);
ifftshow(Image_butterworth);

推荐答案

首先，没有这样的命令叫 ifftshow 。其次，你没有过滤任何东西。你正在做的就是可视化图像的光谱。

For one thing, there is no such command called ifftshow. Secondly, you aren't filtering anything. All you're doing is visualizing the spectrum of the image.

在可视化光谱方面，你现在如何做到这一点非常危险。首先，您可视化每个空间频率分量的系数，其本质上是复值。如果您希望以对我们大多数人有意义的方式可视化频谱，最好先看一下幅度或阶段。但是，由于这是Butterworth滤波器，因此最好将其应用于滤波器的幅度。

In terms of visualizing the spectrum, how you're doing it right now is very dangerous. For one thing, you are visualizing the coefficients at each spatial frequency component which is complex-valued in nature. If you want to visualize the spectrum in a way that makes sense to most of us, it's better to take a look at either the magnitude or phase. However, because this is a Butterworth filter, it's best to apply it to the magnitude of the filter.

您可以使用<$ c找到频谱的幅度$ c> abs 功能。即使你这样做，如果你直接在幅度上做了 imshow ，你将获得一个零点的可视化，除了以外的。这是因为直流分量太大而其余光谱相比较小。

You can find the magnitude of the spectrum by using the abs function. Even when you do that, if you did imshow directly on the magnitude, you will get a visualization that is zero everywhere except for the middle. This is because the DC component is so large and the rest of the spectrum is small in comparison.

让我举个例子。这是摄像师图像，它是图像处理工具箱的一部分：

Let me show you an example. This is the cameraman image that is part of the image processing toolbox:

im = imread('cameraman.tif');
figure;
imshow(im);

现在，让我们可视化光谱并确保DC分量位于图像的中心 - 您已经使用 fftshift 。 将图像转换为 double 也是一个好主意，以确保最佳的数据精度。此外，请确保您应用 abs 来查找幅度：

Now, let's visualize the spectrum and ensuring that the DC component is in the centre of the image - you already did this with fftshift. It's also a good idea to cast the image to double to ensure the best precision of data. In addition, make sure you apply abs to find the magnitude:

fftim = fftshift(fft2(double(im)));
mag = abs(fftim);
figure;
imshow(mag, []);

正如您所看到的，由于我提到的原因，它不是很有用。可视化图像光谱的更好方法通常是对光谱应用 log 变换。如果您想要取消平均值或删除平均值以使动态范围更适合显示，这也很有用。换句话说，您可以将1加到幅度上，然后将对数应用于幅度，以便更高的值可以逐渐减小。你使用哪个基数并不重要，所以我只使用由 log 命令封装的自然对数：

As you can see, it's not very useful due to the reason that I mentioned. A better way to visualize the spectrum of the image is usually to apply a log transformation to the spectrum. This is also useful if you want to de-mean or remove the mean so that the dynamic range fits better for display. In other words, you would add 1 to the magnitude, then apply a logarithm to the magnitude so that higher values can taper off. It doesn't matter which base you use, so I'll just use the natural logarithm which is encapsulated by the log command:

figure;
imshow(log(1 + mag), []);

现在更多更好。现在我们将介绍您的过滤机制。您的巴特沃斯滤波器略有不正确。坐标的 meshgrid 略有错误。在结束时间间隔的 -1 操作需要走出去：

Now that's much better. Now we'll get onto your filtering mechanism. Your Butterworth filter is slightly incorrect. The meshgrid of coordinates is slightly wrong. The -1 operation that's at the ending interval needs to go outside:

[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);

请记住，您正在定义关于图像中心的对称间隔，以及您最初的原因是什么不对。我还想提一下，这看起来像是一个高通滤波器，因此输出应该看起来像边缘检测。此外，Butterworth高通滤波器的定义不正确。频域中滤波器的正确定义是：

Remember, you are defining a symmetric interval about the centre of the image, and what you had originally wasn't correct. I'd also like to mention that this looks like a high-pass filter, so the output should look like an edge detection. In addition, the definition of the Butterworth high pass filter is incorrect. The correct definition of the filter in frequency domain is:

D（u，v）是频域中图像中心的距离， 是 B 时的截止距离是一个控制比例因子，用于控制在截止距离处所需的增益。 n 是过滤器的顺序。 你的情况是 d = 50 。在实践中， B = sqrt（2） - 1 ，以便在的截止距离， D（u，v）= 1 / sqrt（2）= 0.707 ，这是电子电路滤波器中最常见的3dB截止频率。有时为了简单起见，你会看到 B 设置为1，但通常将其设置为 B = sqrt（2） - 1 。

D(u,v) is the distance from the centre of the image in frequency domain, Do is the cutoff distance while B is a controlling scale factor controlling what the desired gain would be at the cutoff distance. n is the order of the filter. Do in your case is d = 50. In practice, B = sqrt(2) - 1 so that at the cutoff distance of Do, D(u,v) = 1 / sqrt(2) = 0.707, which is the 3 dB cutoff frequency mostly seen in electronics circuit filters. Sometimes you'll see B being set to 1 for simplicity, but it's common to set this to B = sqrt(2) - 1.

但是，您当前的代码没有进行任何过滤。要在频域中进行滤波，只需将图像的频谱与滤波器本身的频谱相乘即可。这相当于空间域中的卷积。完成后，您只需撤消对图像执行的 fftshift ，进行逆FFT，然后消除由于数值不精确而导致的任何虚构成分。另外，让我们转换为 uint8 以确保我们尊重原始图片类型。

However, your current code isn't doing any filtering. To filter in the frequency domain, you simply multiply the spectrum of the image with the spectrum of the filter itself. This is equivalent to convolution in the spatial domain. Once you do that, you simply undo the fftshift that was performed on the image, take the inverse FFT and then eliminate any imaginary components that are due to numerical imprecision. Also, let's cast to uint8 to make sure that we respect the original image type.

这可以像所以：

%// Your code with meshgrid fix
n=1;
d=50;
h=size(im,1);
w=size(im,2);
fftim = fftshift(fft2(double(im)));
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));

%%%%%%// New code
B = sqrt(2) - 1; %// Define B
D = sqrt(x.^2 + y.^2); %// Define distance to centre
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n)));
out_spec_centre = fftim .* hhp;

%// Uncentre spectrum
out_spec = ifftshift(out_spec_centre);

%// Inverse FFT, get real components, and cast
out = uint8(real(ifft2(out_spec)));

%// Show image
imshow(out);

如果您想查看过滤后的频谱，请执行以下操作：

If you want to see what the filtered spectrum looks like, just do this:

figure;
imshow(log(1 + abs(out_spec_centre)), []);

我们得到：

这是有道理的。您可以看到，在光谱的中间，与光谱的外边缘相比，它稍微更暗。这是因为使用高通巴特沃斯滤波器，你正在放大更高频率的项，并且它可视化为更高的强度。

This makes sense. You see that in the middle of the spectrum, it's slightly darker in comparison to the outer edges of the spectrum. That's because with the high-pass Butterworth filter, you are amplifying the higher frequency terms and it gets visualized to be a higher intensity.

现在， out 包含您的过滤图像，我们终于得到了这个：

Now, out contains your filtered image, and we finally get this:

看起来效果不错！但是，天真地将图像转换为 uint8 会将任何负值截断为0，并将任何正值大于255到255.由于这是边缘检测，因此您希望同时检测到负面和正面转换...所以一个好主意是规范化输出，使其范围从 [0,1] ，然后在乘以255后，使用 uint8 进行投射。这样，图像中的任何更改都不会显示为灰色，负面更改会显示为黑色，而正面更改会显示为白色。 ..所以你做这样的事情：

That looks like a fine result! However, naively casting the image to uint8 truncates any negative values to 0 and any positive values greater than 255 to 255. Because this is an edge detection, you want to detect both the negative and positive transitions... so a good idea would be to normalize the output so that it ranges from [0,1], and then cast with uint8 after you multiply by 255. This way, no changes in the image get visualized to gray, negative changes get visualized as dark and positive changes get visualized as white.... so you'd do something like this:

%// Your code with meshgrid fix
n=1;
d=50;
h=size(im,1);
w=size(im,2);
fftim = fftshift(fft2(double(im)));
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));

%%%%%%// New code
B = sqrt(2) - 1; %// Define B
D = sqrt(x.^2 + y.^2); %// Define distance to centre
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n)));
out_spec_centre = fftim .* hhp;

%// Uncentre spectrum
out_spec = ifftshift(out_spec_centre);

%// Inverse FFT, get real components
out = real(ifft2(out_spec));

%// Normalize and cast
out = (out - min(out(:))) / (max(out(:)) - min(out(:)));
out = uint8(255*out);

%// Show image
imshow(out);

我们得到这个：

这篇关于高通巴特沃斯滤波器在MATLAB中对图像进行滤波的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！