MATLAB 中图像的高通巴特沃斯滤波器
High Pass Butterworth Filter on images in MATLAB
我需要在 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);
一方面，没有名为 ifftshow 的命令。其次，你没有过滤任何东西。您所做的只是可视化图像的光谱。
就频谱可视化而言，您现在的做法非常危险。一方面，您正在可视化每个空间频率分量的系数，这些分量本质上是 complex-valued。如果您想以对我们大多数人都有意义的方式可视化频谱，最好看一下 magnitude 或 phase .但是，因为这是一个巴特沃斯滤波器，所以最好将它应用到滤波器的幅度上。
您可以使用 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 变换。如果您想要 de-mean 或删除均值以使动态范围更适合显示，这也很有用。换句话说，您可以将幅度加 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);
请记住，您定义的是关于图像中心的对称间隔，而您最初的定义是不正确的。我还想提一下，这看起来像 high-pass 过滤器 ，因此输出应该看起来像边缘检测。另外，巴特沃斯高通滤波器的定义不正确。滤波器在频域的正确定义是：
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)), []);
我们得到：
这是有道理的。你会看到在光谱的中间，与光谱的外边缘相比，它稍微暗一些。那是因为使用 high-pass Butterworth 滤波器，您正在放大更高频率的项，并且它被可视化为更高的强度。
现在，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);
我们得到这个：
我觉得你的工作应该有点不同
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');
我需要在 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);
一方面，没有名为 ifftshow 的命令。其次，你没有过滤任何东西。您所做的只是可视化图像的光谱。
就频谱可视化而言，您现在的做法非常危险。一方面，您正在可视化每个空间频率分量的系数，这些分量本质上是 complex-valued。如果您想以对我们大多数人都有意义的方式可视化频谱，最好看一下 magnitude 或 phase .但是，因为这是一个巴特沃斯滤波器，所以最好将它应用到滤波器的幅度上。
您可以使用 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 变换。如果您想要 de-mean 或删除均值以使动态范围更适合显示，这也很有用。换句话说，您可以将幅度加 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);
请记住，您定义的是关于图像中心的对称间隔，而您最初的定义是不正确的。我还想提一下，这看起来像 high-pass 过滤器 ，因此输出应该看起来像边缘检测。另外，巴特沃斯高通滤波器的定义不正确。滤波器在频域的正确定义是：
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)), []);
我们得到：
这是有道理的。你会看到在光谱的中间，与光谱的外边缘相比，它稍微暗一些。那是因为使用 high-pass Butterworth 滤波器，您正在放大更高频率的项，并且它被可视化为更高的强度。
现在，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);
我们得到这个：
我觉得你的工作应该有点不同
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');