利用MATLAB对图像进行高通巴特沃斯滤波

Question

为了实现图像滤波，我需要在MATLAB中实现一个高通巴特沃斯滤波器。我已经实现了一个，但它看起来不起作用。这是我写的代码。有人能告诉我是怎么回事吗？

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。其次，你没有过滤任何东西。你所做的就是把图像的光谱可视化。

在可视化光谱方面，你现在怎么做是非常危险的。首先，你正在可视化每个空间频率分量上的系数，这在本质上是complex-valued。如果你想以一种对我们大多数人来说都有意义的方式来可视化光谱，最好还是看看震级或者相位。但是，因为这是一个巴特沃斯过滤器，所以最好把它应用到过滤器的大小上。

您可以使用abs函数来找出光谱的大小。即使你这样做，如果你直接在震级上做了imshow，你将得到一个可视化，除了中间的之外，在任何地方都是零。这是因为直流分量是如此大，其余的频谱是很小的比较。

让我给你举个例子。这是作为图像处理工具箱的一部分的摄影师图像：

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

请记住，您正在定义一个关于图像中心的对称间隔，而最初的内容是不正确的。我还想提到，这看起来像一个高通滤波器，所以输出应该看起来像边缘检测。此外，巴特沃斯高通滤波器的定义是不正确的。滤波器在频域的正确定义是：

​

​

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]到[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);

我们得到了这个：

​

​

Answer 2

我认为你应该做一些不同的工作

n=1;
D0=50; % change the name for d0, d is usuaally the (u²+v²)⁽1/2)

A=1.5; % normally the amplitude is 1

im=imread('cameraman.jpg');

[M,N]=size(im); % is easy to get the h and w like this

% compute the 2d fourier transform in order to multiply

F=fft2(double(im));

% compute your filter and do the meshgrid for your matrix but it is M*n, and get only the real part

u=0:(M-1);
v=0:(N-1);

idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H =A * (1./(1 + (D0./D).^(2*n)));

% multiply element by element

G=H.*F;
g=real(ifft2(double(G)));
subplot(1,2,1); imshow(im); title('Input image');
subplot(1,2,2); imshow(g,[ ]); title('filtered image');