2013-12-08 18:02:49 thnh169 阅读数 24318

       1.高通滤波器

       首先,对一副图像进行如下二维傅里叶变换。

我们将u=0和v=0带上式,我们可以得到如下式子。

根据上式,可以到F(0,0)的值是非常大的。这里,我们将F(0,0)称为直流分量,直流分量比其他的成分要大好几个数量级。所以,这也就是傅里叶谱为什么需要使用对数变换才能看清楚的原因。
       这里,对于高通滤波器而言,由于直流分量被衰减,所以,所得到的图像的动态范围是非常狭窄的,也就造成了图像偏灰。进一步而言,保持直流(DC)分量,对别的部分进行增幅,可以增强图像的细节。这样的滤波器称为锐化滤波器。这一小节主要介绍高通滤波器与锐化滤波器。

        1.1理想高通滤波器

        
        这里的D0是滤波器的阻带半径,而D(u,v)是点到滤波器中央的距离。理想高通的滤波器的振幅特性如下所示。

用这个滤波器对图像进行处理,可得到如下所示的结果。我们可以看到,与理想的低通滤波器一样,所得到的图像有很明显的振铃现象。结果图像从视觉上来看,有些偏暗,这是因为图像的直流分量被滤掉的原因。


       1.2巴特沃斯高通滤波器

       同样的,巴特沃斯高通滤波器也可以通过改变次数n,对过度特性进行调整。过大的n会造成振铃现象。


       1.3高斯高通滤波器

       高斯滤波器的过度特性很好,所以不会发生振铃现象。


       1.4 上述三种滤波器的实现Matlab代码

close all;
clear all;

%% ---------Ideal Highpass 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 = ones(P,Q);
H_2 = ones(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) = 0; end    
        if(D <= 160) H_2(x+(P/2)+1,y+(Q/2)+1) = 0; 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 Highpass 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 Highpass 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 Highpass 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_0/D)^2);   
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^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;

%% ---------Gaussian Highpass 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) = 1 - exp(-(D*D)/(2*D_0*D_0));   
        D_0 = 160;
        H_2(x+(P/2)+1,y+(Q/2)+1) = 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-------
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 Highpass (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 Highpass (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');

       1.5 锐化滤波器

       按照之前所说的,锐化滤波器是将傅里叶谱的直流分量保留,然后将其余的成分增幅。使用锐化滤波器,可以对图像的细节进行增强,使得细节凸显出来。锐化滤波器的表达式如下所示。

      其实上式的目的很明显,就是先将原图的傅里叶保留下来,然后叠加上高通滤波器的结果,所得到的图像就是锐化后的图像了。这里为了调整锐化程度,引入了两个变量可以调整直流分量的衰减程度,可以调整高频分量的增幅程度。

       下面是代码。
close all;
clear all;
clc;

%% ---------The High-Fre-Emphasis Filters (Fre. Domain)------------
f = imread('blurry_moon.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_HP = 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 = 80;
        H_HP(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0));   
     end
end

G_HP = H_HP .* F;

G_HFE = (1 + 1.1 * H_HP) .* F;

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

g_2 = real(ifft2(G_HFE));
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

g = histeq(g_2);

%g_1 = mat2gray(g_1);
%% -----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(g_1,[0 1]);
xlabel('c).Result image of High-pass Filter');

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

figure();
subplot(1,2,1);
imshow(g_2,[0 1]);
xlabel('e).Result image of High-Fre-Emphasis Filter');

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

       2.带阻滤波器

       同样的,带阻滤波器也有三种特性。高斯、巴特沃斯和理想,三种类型,其数学表达式如下所示。

       其带通滤波器可以使用上面的表格转化而得。

       带阻滤波器可以用于去除周期性噪声,为了体现带阻滤波器的特性,我们先对一幅图像增加很严重的噪声。


       在原图的傅里叶谱上添加了几个很明显的亮点。在对其做IDFT,可以看到,原图被严重的周期噪声污染了。此时,我们可以使用带阻滤波器,可以有很好的去噪效果。为了避免振铃现象,选择使用如下所示巴特沃斯带阻滤波器,所用滤波器的次数为2次。使用空间域的操作,要去除这种噪声基本是不可能的,这也是频域内的操作的优点。

  
close all;
clear all;
clc;
%% ---------------------Add Noise-------------------------
f = imread('left.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_NP = ones(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = 2;
        
        v_k = -200; u_k = 150;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = 200; u_k = 150;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = 0; u_k = 250;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = 250; u_k = 0;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        
        H_NP(x+(P/2)+1,y+(Q/2)+1) = 1 + 700*(1 - H_NP(x+(P/2)+1,y+(Q/2)+1));
     end
end

G_Noise = F .* H_NP;

g_noise = real(ifft2(G_Noise));
g_noise = g_noise(1:1:M,1:1:N);     

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


%% ---------Bondpass Filters (Fre. Domain)------------
H_1 = ones(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 = 250;
        W = 30;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+((D*W)/((D*D) - (D_0*D_0)))^6);   
     end
end

G_1 = H_1 .* G_Noise;

g_1 = real(ifft2(G_1));
g_1 = g_1(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);
    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,2);
imshow(log(1 + abs(G_Noise)),[ ]);
xlabel('c).Fourier spectrum of b');

subplot(1,2,1);
imshow(g_noise,[0 1]);
xlabel('b).Result of add noise');


figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('d).Butterworth Bandpass filter(D=355 W=40 n =2)');

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


figure();
subplot(1,2,2);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('e).Fourier spectrum of f');

subplot(1,2,1);
imshow(g_1,[0 1]);
xlabel('f).Result of denoise');

       3.陷波滤波器(Notch Filter)

       陷波滤波器也用于去除周期噪声,虽然带阻滤波器也能可以去除周期噪声,但是带阻滤波器对噪声以外的成分也有衰减。而陷波滤波器主要对,某个点进行衰减,对其余的成分不损失。使用下图将更容易理解。

       从空间域内看,图像存在着周期性噪声。其傅里叶频谱中,也能看到噪声的所在之处(这里我用红圈标注出来了,他们不是数据的一部分)。我们如果使用带阻滤波器的话,是非常麻烦的,也会使得图像有损失。陷波滤波器,能够直接对噪声处进行衰减,可以有很好的去噪效果。
       其表达式如下所示,陷波滤波器可以通过对高通滤波器的中心,进行位移就可以得到了。

这里,由于傅里叶的周期性,傅里叶频谱上不可能单独存在一个点的噪声,必定是关于远点对称的一个噪声对。这里的需要去除的噪声点,其求取的方式如下所示。
       针对于上图,我们设计如下滤波器,去进行去噪。
(图片下标错了,后续更新改过来!)
       所得到的结果,如下所示。噪声已经被去除了,画质得到了很大的改善。

                                                              (图片下标错了,后续更新改过来!)
close all;
clear all;
clc;

%% ---------Butterworth Notch filter (Fre. Domain)------------
f = imread('car_75DPI_Moire.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_NF = ones(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = 30;
        
        v_k = 59; u_k = 77;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = 59; u_k = 159;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = -54; u_k = 84;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        
        v_k = -54; u_k = 167;
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
        D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
     end
end

G_1 = H_NF .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(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);
    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_NF,[0 1]);
xlabel('c).Butterworth Notch filter(D=30 n=2)');

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

figure();
subplot(1,2,2);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('e).Fourier spectrum of d');

subplot(1,2,1);
imshow(g_1,[0 1]);
xlabel('d).Result image');
原文发于博客:http://blog.csdn.net/thnh169/ 

2015-04-23 20:13:23 chanxiaogang 阅读数 17245

1.图像频域滤波的基本概念
空间域和频域滤波的基础都是卷积定理。该定理可以写为
这里写图片描述
这里写图片描述
其中,符号*表示两个函数的卷积,双箭头两边的表达式组成了傅里叶变换对。第一个表达式表明两个空间函数的卷积可以通过计算两个傅里叶变换函数的乘积逆变换得到。在滤波问题上,我们更关注第一个表达式,它构成了整个频域滤波的基础,式中的乘积实际上就是两个二维矩阵F(u,v)和H(u,v)对应元素之间的乘积。
傅立叶变换相关函数
MATLAB提供了几个和傅里叶变换相关的函数。其说明如下:
fft2(F); 二维傅立叶变换
abs(F); 获得傅立叶频谱
fftshift(F); 将变换的原点移至频率矩形的中心
ifft2(F); 二维傅立叶反变换
real(ifft2(F)); 提取变换后的实部
imag(ifft2(F)); 提取变换后的虚部
2. 频域滤波的基本步骤
频域滤波通常应遵循以下步骤。
(1)计算原始图像f(x,y)的DFT,得到F(u,v)。
(2)将频谱F(u,v)的零频点移动到频谱图的中心位置。
(3)计算滤波器函数H(u,v)与F(u,v)的乘积G(u,v)。
(4)将频谱G(u,v)的零频点移回到频谱图的左上角位置。
(5)计算第(4)步计算结果的傅里叶反变换g(x,y)。
(6)取g(x,y)的实部作为最终滤波后的结果图像。
按照该步骤,在MATLAB中很容易编程实现频域滤波。滤波能否取得理想结果的关键取决于频域滤波函数H(u,v),常常称之为滤波器,或滤波器传递函数。因为它在滤波中抑制或滤除了频谱中某些频率的分量,而保留其他一些频率不受影响。
频域滤波的方框图表示
3 频域滤波器
1 理想低通滤波器(ILPF)
这里写图片描述

D(u,v)是点(u,v)距频率原点的距离,D0为截止频率。如果图像大小为M×N,其变换亦为M×N。中心化之后,矩形中心在(M/2,N/2)

这里写图片描述

说明:在半径为D0的圆内,所有频率没有衰减地通过滤波器,而在此半径的圆之外的所有频率完全被衰减掉。虽然这个滤波器不能使电子元件来模拟实现,但可以在计算机中用传递函数实现。由于理想低通滤波器的过度特性过于急峻,所以会产生了振铃现象。
理想低通滤波器三维透视图
2 n阶巴特沃斯低通滤波器(BLPF)
这里写图片描述

n表示的是巴特沃斯滤波器的次数。随着次数的增加,振铃现象会越来越明显。与理想低通滤波器不同的是,巴特沃斯低通滤波器传递函数并不是在D0处突然不连续。对于具有平滑传递函数的滤波器,我们通常要定义一个截止频率,在该点处H(u,v)会降低为其最大值的某个给定比例。例如当D(u,v)=D0时,H(u,v)=0.5。通常,BLPF的平滑效果好于ILPF。
巴特沃斯低通滤波器三维透视图
3 高斯低通滤波器(GLPF)
这里写图片描述
D0为截止频率,当D(u,v)=D0 时,滤波器下降到其最大值的0.607,无振铃现象。
高斯低通滤波器三维透视图
4 理想高通滤波器(IHPF)
这里写图片描述
D(u,v)是点(u,v)距频率原点的距离,D0为截止频率。如果图像大小为M×N,其变换亦为M×N。中心化之后,矩形中心在(M/2,N/2)
这里写图片描述

在半径为D0的圆外,所有频率没有衰减地通过滤波器,而在此半径的圆之内的所有频率完全被衰减掉。虽然这个滤波器不能使电子元件来模拟实现,但可以在计算机中用传递函数实现。由于理想高通滤波器的过度特性过于急峻,所以会产生了振铃现象。
![理想高通滤波器三维透视图](http://img.blog.csdn.net/20150423211457107)
5  n阶巴特沃斯高通滤波器(BHPF)

这里写图片描述
n表示的是巴特沃斯滤波器的次数。随着次数的增加,振铃现象会越来越明显。与理想高通滤波器不同的是,巴特沃斯高通滤波器传递函数并不是在D0处突然不连续。对于具有平滑传递函数的滤波器,我们通常要定义一个截止频率,在该点处H(u,v)会降低为其最大值的某个给定比例。例如当D(u,v)=D0时,H(u,v)=0.5。通常,BHPF的平滑效果好于IHPF。
巴特沃斯高通滤波器三维透视图
6 高斯高通滤波器(GHPF)
这里写图片描述
高斯高通滤波器三维透视图

程序代码如下:

function [U,V] = dftuv( M,N )
%DFTUV computes meshgrid frequency matrices
%set up range of variables
u=0:(M-1);
v=0:(N-1);
%compute the indices for use in meshgrid
idx=find(u > M/2);
u(idx)=u(idx)-M;
idy=find(v > N/2);
v(idy)=v(idy)-N;
%compute the meshgrid arrays
[V,U]=meshgrid(v,u);
End
function [H,D] = lpfilter( type,M,N,D0,n)
%lpfilter computes frequency domin lowpass filters
%H creates the transfer function of a lowpass filter

%use function dftuv to set up the meshgrid arrays needed for computing 
%the required distances
[U,V]=dftuv(M,N);
D=sqrt(U.^2+V.^2);
%begin filter computations
switch type
    case 'ideal'
        H=double(D<=D0);
    case 'btw'
        if nargin ==4
            n=1;
        end
        H=1./(1+(D./D0).^(2*n));
    case 'gaussian'
        H=exp(-(D.^2)./(2*(D0^2)));
    otherwise
        error('unknown filter type')  
end        
end
function H = hpfilter( type,M,N,D0,n )
%hpfilter computes frequency domin highpass filters
%H creates the transfer function of a highpass filter
%the transfer function Hhp of a highpass filter is 1-Hlp
%we can use function lpfilter to generate highpass filters.
if nargin==4
    n=1; 
end
%generate highpass filter.
hlp=lpfilter(type,M,N,D0,n);
H=1-hlp;

end
%图像频域高通滤波主程序
clear all;
I=imread('man.bmp'); %读入原图像
F=fft2(double(I));   %快速傅里叶变换
F=fftshift(F);       %图像零频点移动到频谱图中心

[M,N]=size(F);       %计算矩阵F的大小
%调用高通滤波器函数
H=hpfilter('ideal',M,N,40); 
%H=hpfilter('btw',M,N,40);
%H=hpfilter('gaussian',M,N,40);

H=fftshift(H);
G =F .* H;        %应用滤镜
%傅里叶反变换
g=ifft2(G);
%显示并比较结果
figure(1),imshow(I);  
title('显示原图像');
figure(2),imshow(log(1+abs(F)),[]);  
title('傅立叶变换后的频谱图');
figure(3),imshow(H,[]);
title('高通滤波器频域图像');
figure(4),imshow(log(1+abs(G)),[]);
title('高通滤波器处理后的频谱图');
figure(5),imshow(abs(real(g)),[]);
title('傅里叶反变换后的图像');

结论
在频谱中,低频主要对应图像在平滑区域的总体灰度级分布,而高频对应图像的细节部分,如边缘和噪声。

2017-11-16 10:38:43 hhyqhh 阅读数 1741

课后作业,实现“理想、巴特沃斯、高斯”高通低通滤波器。


代码基于Matlab实现。完整代码及处理结果见:GitHub

步骤

  1. 加载图像
  2. 中心化图像
  3. 傅里叶变换
  4. 与滤波器做运算(空域的卷积运算对应频域的乘法运算)
  5. 傅里叶反变换
  6. 裁剪图像

低通滤波器

理想低通滤波器

实现代码

close all;
clear all;
f = imread('cameraman.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);
H_3 = 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 <= 10)  H_1(x+(P/2)+1,y+(Q/2)+1) = 1; end    
        if(D <= 50)  H_2(x+(P/2)+1,y+(Q/2)+1) = 1; end
        if(D <= 150)  H_3(x+(P/2)+1,y+(Q/2)+1) = 1; end
     end
end

G_1 = H_1 .* F;
G_2 = H_2 .* F;
G_3 = H_3 .* 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);         

g_3 = real(ifft2(G_3));
g_3 = g_3(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);
        g_3(x,y) = g_3(x,y) * (-1)^(x+y);
    end
end

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

subplot(2,2,2);
imshow(g_1,[0 1]);
xlabel('b).Ideal_Lowpass/D0=10');

subplot(2,2,3);
imshow(g_2,[0 1]);
xlabel('c).Ideal_Lowpass/D0=50');

subplot(2,2,4);
imshow(g_3,[0 1]);
xlabel('d).Ideal_Lowpass/D0=150');

结果

可以看到D0的值越少,即允许用过的频率越低。高频信息损失越多,表现在图像上是边缘信息的丢失,看起来图像模糊。如图b。

巴特沃斯低通滤波器

这里只展示核心的滤波器实现代码。完整代码见GitHub

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 = 10;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^2);
        D_0 = 50;
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^2);
        D_0 = 150;
        H_3(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^2);    
     end
end

高斯低通滤波器

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 = 10;
        H_1(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));
        D_0 = 50;
        H_2(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));
        D_0 = 150;
        H_3(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));
     end
end

高通滤波器

高通滤波器和低通滤波器在形式上相反,即之前允许通过的频率截止,之间截止的频率允许通过。公式上表示为高通滤波器=1-低通滤波器的关系。

理想高通滤波器

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 >= 10)  H_1(x+(P/2)+1,y+(Q/2)+1) = 1; end    
        if(D >= 50)  H_2(x+(P/2)+1,y+(Q/2)+1) = 1; end
        if(D >= 150)  H_3(x+(P/2)+1,y+(Q/2)+1) = 1; end
     end
end

巴特沃斯高通滤波器

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 = 10;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^2);
        D_0 = 50;
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^2);
        D_0 = 150;
        H_3(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^2);       
     end
end

高斯高通滤波器

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 = 10;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1-exp(-(D*D)/(2*D_0*D_0));
        D_0 = 50;
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1-exp(-(D*D)/(2*D_0*D_0));
        D_0 = 150;
        H_3(x+(P/2)+1,y+(Q/2)+1) = 1-exp(-(D*D)/(2*D_0*D_0));
     end
end

2018-09-23 20:19:06 iiiil7 阅读数 4990

一,实验原理

设置低通滤波截止频率为Dl0=50,即:

高通滤波器截止频率为Dh0=30,即:

二,实验过程和结果

clc;
clear;
close all;

I=imread('test5.gif','gif');            %自己设置路径
subplot(2,2,1);
imshow(I);
title('test5原图');
f=double(I);     % chage into double as MATLAB doesn’t suppor calculation
  % of image in unsigned int type
  subplot(2,2,2);
  G=fft2(f);
imshow(log(abs(G)),[]);
title('test5幅度谱图');
g=fft2(f);       % fourier transform
g=fftshift(g);  % zero-frequency area centralized
[M,N]=size(g);
dl0=50;            %cutoff frequency
m=fix(M/2); n=fix(N/2);
for i=1:M
       for j=1:N
           d=sqrt((i-m)^2+(j-n)^2);
           if(d<=dl0)
           h=1;
           else h=0;
                  end
           result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
Gl1=ifft2(result);
Gl2=uint8(real(Gl1));
subplot(2,2,3);
imshow(Gl2) ;
title('test5低通滤波图');

dh0=30;
for i=1:M
       for j=1:N
           d=sqrt((i-m)^2+(j-n)^2);
           if(d>=dh0)
           h=1;
           else h=0;
                  end
           result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
Gh1=ifft2(result);
Gh2=uint8(real(Gh1));
subplot(2,2,4);
imshow(Gh2) ;
title('test5高通滤波图');

三,实验总结

低通滤波器由于滤掉了高频成分,高频成分含有大量边缘信息,所以造成了一定程度的图像模糊。高通滤波器滤掉了低频成分,保留了高频成分,即保留了边界信息,所以显示出原图像的边界。

2016-05-25 16:48:40 u014593748 阅读数 633

图像的低频与图像中变化缓慢的灰度分量有关,图像的高频与图像中变化快速的灰度分量有关。频率域相乘等同于空间域卷积。

使用频率域滤波器平滑图像

频率域平滑可以通过对高频的衰减来达到,即低通滤波。

使用频率域滤波器锐化图像

频率域的锐化可以通过低频的衰减来达到,即高通滤波。

选择性滤波

带通滤波、带阻滤波。

参考书目:数字图像处理(中译第三版)[冈萨雷斯]

没有更多推荐了,返回首页