高通巴特沃斯滤波器在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直接在幅度上做,你将得到一个除了中间以外零点的可视化.这是因为DC分量如此之大,相比之下其余频谱很小.

让我举个例子.这是摄影师图像,它是图像处理工具箱的一部分:

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

现在,让我们可视化光谱并确保DC分量位于图像的中心 - 您已经这样做了fftshift.将图像转换double为最佳数据精度也是一个好主意.此外,请确保您申请abs查找幅度:

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

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

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

现在,这是很多更好.现在我们将介绍您的过滤机制.您的巴特沃斯滤波器略有不正确.该meshgrid坐标是稍有不妥.在-1这在结束区间操作需要到外面去:

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

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

D(u,v)是在频域中Do距图像中心的距离,是截止距离,B而是控制在截止距离处所需增益的控制比例因子. n是过滤器的顺序. Do在你的情况下是d = 50.在实践中,B = sqrt(2) - 1使得在截止距离处Do,D(u,v) = 1 / sqrt(2) = 0.707这是在电子电路滤波器中最常见的3dB截止频率.有时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],然后uint8在乘以255后进行投射.这样,图像中的任何变化都不会显示为灰色,负面变化会显示为黑暗,正面变化可视化为白色....所以你会做这样的事情:

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

我们得到这个: