[数字图像处理]频域滤波(1)--基础与低通滤波器

时间  2020-07-25 标签 数字 图像处理 频 域 滤波 1 基础 低 通 滤波器

       以前的博文主要介绍了空间域内的滤波器，本文主要从频域的角度进行分析。主要使用傅里叶变换，将空间域的图像转换到频域内，在频域内进行数字图像处理。这部分的内容及其重要，频域内的处理能够解决空间域内没法完成的图像加强。本文首先从数学角度，对图像的频域内的性质进行分析，而后在着重介绍滤波器在频域内的性质。函数

        1.傅里叶变换与频域

        在以前的文中，咱们已经进行过一些基本的图像处理。好比，使用低通滤波能够将图像模糊，也有些许降噪的做用。这些都是在空间域内进行的滤波处理，这个处理主要是依靠卷积来进行计算的。首先，从连续的一维卷积入手，以下所示。spa


       将上式进行傅里叶变换，能够获得以下结果。.net


        从这个式子，咱们能够获得一个重要的结论。也就是，函数与卷积的傅里叶变换所获得的结果，是函数与的傅里叶变换与的乘积。再将其总结得简单易懂一些，有以下结论。3d


        在将其扩展到二维的形况下，假设尺寸为MxN的图像，以下关系是成立的。code


       其实到这，基本的原理就明了的。咱们所看到的图像，均为空间域内的表现形式，咱们没法辨识出频域内的图像。要进行频域内的滤波器处理，首先就须要进行傅里叶变换，而后直接进行滤波处理，最后再用反傅里叶变换倒回到空间域内。blog

       到此，已经能够开始空间域内的滤波处理了。可是，还有一点须要注意的地方。使用某个一维信号来举例子，一维信号的傅里叶变换是以2π为周期的函数。因此，咱们经常使用的范围[-π,π]来表示这个信号的傅里叶变换，以下所示。get


        这样作的好处是，靠近0的成分就是低频，靠近-π与π的成分就表示高频。而对于图像而言，在Matlab中，咱们使用fft2()这个函数来求取图像的傅里叶变换。博客

g = fft2(f);
        上面这个代码求取的，实际上是 范围[0,π]内的傅里叶变换。为了方便理解，下图画出了本行代码所求取的图像的傅里叶变换的范围（右）和与其等效的一维傅里叶变换的范围（左）。


       很显然，这并非但愿的范围，下面这个代码能够求取[0,2π]内的傅里叶变换。数学

P = 2*M;
Q = 2*N;
F = fft2(f,P,Q);
       下图画出了本行代码所求取的图像的傅里叶变换的范围（右）和与其等效的一维傅里叶变换的范围 （左） 。所获得的图像F(u,v)的尺寸为PxQ。


         咱们须要对其移动一下，以下图所示，咱们须要的是粉色范围的区域。io


下面，从数学上分析一下，如何得到这个部分的频谱。对于傅里叶变换，有以下性质。


这个特性称为平移特性，粉色部分的频谱，将带入上式，咱们能够获得以下式子。


为次，咱们已经获得了粉色范围的频谱。越靠近傅里叶频谱图像中间的成分，表明了低频成分。其Matlab代码以下所示。

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

        代码所获得的结果，以下图所示。


        接下来，咱们总结一下频域滤波的步骤：

        ①：先将图像作频域内的水平移动，而后求原图像f(x,y)的DFT，获得其图像的傅里叶谱F(u,v)。


        ②：与频域滤波器作乘积，

        ③：求取G(u,v)的IDFT，而后再将图像作频域内的水平移动（移动回去），其结果可能存在寄生的虚数，此时忽略便可。


        ④：这里使用ifft2函数进行IDFT变换，获得的图像的尺寸为PxQ。切取左上角的MxN的图像，就能获得结果了。

        2.低通滤波器

        2.1理想的低通滤波器

       其中，D0表示通带的半径。D(u,v)的计算方式也就是两点间的距离，很简单就能获得。

       使用低通滤波器所获得的结果以下所示。低通滤波器滤除了高频成分，因此使得图像模糊。因为理想低通滤波器的过分特性过于急峻，因此会产生了振铃现象。


         

        2.2巴特沃斯低通滤波器

       一样的，D0表示通带的半径，n表示的是巴特沃斯滤波器的次数。随着次数的增长，振铃现象会愈来愈明显。

   

       2.3高斯低通滤波器
       D0表示通带的半径。高斯滤波器的过分特性很是平坦，所以是不会产生振铃现象的。

       3.实现代码

close all;
clear all;

%% ---------Ideal Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        if(D <= 60)  H_1(x+(P/2)+1,y+(Q/2)+1) = 1; end    
        if(D <= 160)  H_2(x+(P/2)+1,y+(Q/2)+1) = 1; end
     end
end

G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end

%% -----show-------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c).Ideal Lowpass filter(D=60)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Ideal Lowpass filter(D=160)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');
close all;
clear all;

%% ---------Butterworth Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        D_0 = 100;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^2);   
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^6);
     end
end


G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end

%% -----show-------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c)Butterworth Lowpass (D_{0}=100,n=1)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Butterworth Lowpass (D_{0}=100,n=3)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');
close all;
clear all;
clc;
%% ---------Gaussian Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        D_0 = 60;
        H_1(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));   
        D_0 = 160;
        H_2(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));
     end
end

G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end


%% -----show-------
close all;

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c)Gaussian Lowpass (D_{0}=60)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Gaussian Lowpass (D_{0}=160)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');

      原文发于博客：http://blog.csdn.net/thnh169/

联系我们 沪ICP备2021010478号-7