image - 在MATLAB中对图像进行高通巴特沃斯滤波器

Question

为了图像过滤的目的，我需要在MATLAB中实现一个高通Butterworth滤波器。我已经实现了一个，但看起来似乎行不通。这是我编写的代码。谁能告诉我这是怎么回事？

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);

Answer 1

一方面，没有这样的命令ifftshow。其次，您没有过滤任何内容。您要做的只是可视化图像的光谱。

就可视化频谱而言，您现在的工作方式非常危险。一方面，您正在可视化每个空间频率分量的系数，这些系数本质上是复数值。如果您想以对我们大多数人有意义的方式可视化频谱，则最好查看幅度或相位。但是，由于这是一个Butterworth滤波器，所以最好将其应用于滤波器的幅度。

您可以使用abs函数找到频谱的幅度。即使执行此操作，如果直接在幅度上执行imshow，则在的所有地方(中间的除外)也将得到可视化的零。这是因为相比之下，直流分量很大，而频谱的其余部分很小。

让我给你看一个例子。这是摄影师图像，它是图像处理工具箱的一部分:

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

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

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

如您所见，由于我提到的原因，它不是很有用。可视化图像光谱的更好方法通常是将 log 变换应用于光谱。如果您想降低平均数或去除均值，以使动态范围更适合显示，这也很有用。换句话说，您可以将幅度加1，然后对幅度应用对数，以便可以逐渐减小较高的值。使用哪种基础都没有关系，所以我只使用log命令封装的自然对数:

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

现在好很多。现在，我们将介绍您的过滤机制。您的Butterworth筛选器略有错误。坐标的meshgrid略有错误。结束时间间隔处的-1操作需要移出外部:

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

请记住，您正在定义一个围绕图像中心的对称间隔，而最初的设置是不正确的。我还要提及的是，这看起来像是高通滤波器，因此输出应该看起来像是边缘检测。此外，巴特沃思高通滤波器的定义不正确。频域中滤波器的正确定义是:


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

但是，您当前的代码未执行任何过滤。要在频域中进行滤波，只需将图像的频谱与滤波器本身的频谱相乘即可。这等效于空间域中的卷积。完成此操作后，您只需撤消在图像上执行的fftshift，进行逆FFT，然后消除由于数值不精确而引起的任何虚部。另外，让我们强制转换为uint8以确保我们尊重原始图像类型。

可以这样做:

%// 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);

如果您想查看过滤后的光谱是什么样子，请执行以下操作:

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

我们得到:



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

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



看起来效果不错!但是，将图像天真地转换为uint8会将任何负值都截断为0，并将任何大于255到255的正值都截断。由于这是边缘检测，因此您希望同时检测负向和正向过渡...因此，一个好主意是到标准化输出，使其在[0,1]范围内，然后乘以255后再用uint8进行转换。这样，图像中的任何变化都不会显示为灰色，负的变化会显示为暗，正的变化会显示为白色....因此您将执行以下操作:

%// 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);

我们得到这个: