图像处理高通频域滤波器

2013-12-08 18:02:49 thnh169 阅读数 26194
  • H.264/AVC的背景与基本概念

    通过本系列课程,观众有望理解视频压缩编码技术的整体发展,精通H.264视频编码技术的框架与细节,为进一步研究H.265/HEVC编码标准、音视频流媒体、视频直播点播等技术奠定坚实的基础。本系列课程适合多媒体方向的...

    2777课时 0分钟 5067人学习 殷汶杰
    免费试看

       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/ 

2020-04-11 16:28:10 qq_44132878 阅读数 71
  • H.264/AVC的背景与基本概念

    通过本系列课程,观众有望理解视频压缩编码技术的整体发展,精通H.264视频编码技术的框架与细节,为进一步研究H.265/HEVC编码标准、音视频流媒体、视频直播点播等技术奠定坚实的基础。本系列课程适合多媒体方向的...

    2777课时 0分钟 5067人学习 殷汶杰
    免费试看

上一期:医学图像处理——彩色图像的显示

“ 在上上期的讲解中,小白已经跟大家讲解了如何构造频域低通滤波器,今天我们就接着讲讲它的“同胞兄第”——高通滤波器。”

 

      就像低通滤波会使图像变得模糊那样,高通滤波通过削弱傅里叶变换的低频而保持高频相对不变,会使得图像变得更加清晰(锐化)在这一期中,我们将考虑几种高通滤波的方法中最简单的一种—— 基本高通滤波器。

     给定一个低通滤波器的传递函数Hlp(u,v),通过使用如下的简单关系,我们可以获得相应高通滤波器的传递函数:

 

代码步骤具体如下:

  1. 按照低通滤波器中的理论构造lpfilter.

function [H,D]=lpfilter(type,M,N,D0,n)

u=-M/2:(M/2-1);v=-N/2:(N/2-1);

idx =find(u > M/2) ;

u(idx) = u(idx)-M;

idy =find(v > N/2) ;

u(idy) = u(idy)-N;

[V,U]=meshgrid(v,u);%画网格准备计算U.^2+V.^2,即公式中的D(u,v)的平方

D = sqrt(V.^2 + U.^2) ;

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



2.构造高通滤波器。

function H = hpfilter (type, M,N, DO,n)

if nargin==4

    n =1;

end

    Hlp = lpfilter(type, M,N, DO,n) ;

    H=1-Hlp;

      可以看到我们构造的函数中可产生三种滤波器,分别是理想、巴氏、高斯高通滤波器。

     最后我们可以看到构造出的理想高通滤波器的频域形态。

H = hpfilter('ideal',500,500,50);

mesh (H(1:10:500,1:10:500)) ;

axis([0 50 0 50 0 1])

axis off

grid off

 

 

实战演练

 

下面我们用高斯高通滤波器来对一幅图进行锐化。

f=imread("Fig0313(a).tif");

PQ = paddedsize (size(f)) ;

DO = 0.05* PQ(1) ;

H = hpfilter('ideal', PQ(1), PQ(2), DO) ;

F=fft2(f,size(H,1),size(H,2));

G=H.*fftshift(F);

g=real(ifft2(fftshift(G)));

g=g(1:size(f,1),1:size(f,2));

figure(1);

subplot(3,2,1);

imshow(f,[]);

title("原图");subplot(3,2,2);

imshow(log(1+abs(fftshift(F))),[]);

title("原图频谱");

subplot(3,2,3);

imshow(H,[]);

title("滤波器频谱");subplot(3,2,4);

imshow(g,[]);

title("滤波后图");

subplot(3,2,5);

imshow(log(1+abs(G)),[]);

title("滤波后图的频谱");

结果如下:

  可以看到滤波后的图片的确被锐化了。

更多干货请关注公众微信号:医电小白的进阶之路

2019-12-14 21:47:59 qq_43571150 阅读数 484
  • H.264/AVC的背景与基本概念

    通过本系列课程,观众有望理解视频压缩编码技术的整体发展,精通H.264视频编码技术的框架与细节,为进一步研究H.265/HEVC编码标准、音视频流媒体、视频直播点播等技术奠定坚实的基础。本系列课程适合多媒体方向的...

    2777课时 0分钟 5067人学习 殷汶杰
    免费试看

Matlab-频域增强实验-彩色图像的频域滤波器

代码链接:https://download.csdn.net/download/qq_43571150/12033263

问题
自行选择一种频域的高通滤波器对彩色图像进行滤波操作, 取3组不同的参数进行实验,根据实验效果进行参数的比较分析。

结果图像👇
在这里插入图片描述
巴特沃兹高通滤波代码👇

clc;
clear all;
close all;
J=imread('05.jpg');
if size(J, 3)==3
    J = rgb2gray(J);
end 
subplot(3,3,1);imshow(uint8(J));xlabel('input');

J=double(J);
f=fft2(J);      %采用傅里叶变换
g=fftshift(f);   %数据矩阵平衡
[M,N]=size(f);
n1=floor(M/2);
n2=floor(N/2);
n=2;

d1=5;%截止频率
for i=1:M        %进行巴特沃兹高通滤波和巴特沃兹高通加强滤波
    for j=1:N
        d=sqrt((i-n1)^2+(j-n2)^2);
        if d==0
            h1=0;
            h2=0.5;
        else
            h1=1/(1+(d1/d)^(2*n));
            h2=1/(1+(d1/d)^(2*n))+0.5;
        end
        gg1(i,j)=h1*g(i,j);
        gg2(i,j)=h2*g(i,j);
    end
end
gg1=ifftshift(gg1);
gg1=uint8(real(ifft2(gg1))); 
subplot(3,3,2);imshow(gg1);   %显示巴特沃兹高通滤波
xlabel('巴特沃兹高通滤波 5');
imwrite(gg1,'05 巴特沃兹高通滤波 5.jpg');
gg2=ifftshift(gg2); 
gg2=uint8(real(ifft2(gg2))); 
subplot(3,3,3);imshow(gg2);   %显示巴特沃兹高通加强滤波
xlabel('巴特沃兹高通加强滤波 5');
imwrite(gg2,'05 巴特沃兹高通加强滤波 5.jpg');


d3=50;%截止频率
for i=1:M        %进行巴特沃兹高通滤波和巴特沃兹高通加强滤波
    for j=1:N
        d=sqrt((i-n1)^2+(j-n2)^2);
        if d==0
            h5=0;
            h6=0.5;
        else
            h5=1/(1+(d3/d)^(2*n));
            h6=1/(1+(d3/d)^(2*n))+0.5;
        end
        gg5(i,j)=h5*g(i,j);
        gg6(i,j)=h6*g(i,j);
    end
end
gg5=ifftshift(gg5);
gg5=uint8(real(ifft2(gg5))); 
subplot(3,3,5);imshow(gg5);   %显示巴特沃兹高通滤波
xlabel('巴特沃兹高通滤波 50');
imwrite(gg5,'05 巴特沃兹高通滤波 50.jpg');
gg6=ifftshift(gg6); 
gg6=uint8(real(ifft2(gg6))); 
subplot(3,3,6);imshow(gg6);   %显示巴特沃兹高通加强滤波
xlabel('巴特沃兹高通加强滤波 50');
imwrite(gg6,'05 巴特沃兹高通加强滤波 50.jpg');



d5=250;%截止频率
for i=1:M        %进行巴特沃兹高通滤波和巴特沃兹高通加强滤波
    for j=1:N
        d=sqrt((i-n1)^2+(j-n2)^2);
        if d==0
            h7=0;
            h8=0.5;
        else
            h9=1/(1+(d5/d)^(2*n));
            h10=1/(1+(d5/d)^(2*n))+0.5;
        end
        gg9(i,j)=h9*g(i,j);
        gg10(i,j)=h10*g(i,j);
    end
end
gg9=ifftshift(gg9);
gg9=uint8(real(ifft2(gg9))); 
subplot(3,3,8);imshow(gg9);   %显示巴特沃兹高通滤波
xlabel('巴特沃兹高通滤波 250');
imwrite(gg9,'05 巴特沃兹高通滤波 250.jpg');
gg10=ifftshift(gg10); 
gg10=uint8(real(ifft2(gg10))); 
subplot(3,3,9);imshow(gg10);   %显示巴特沃兹高通加强滤波
xlabel('巴特沃兹高通加强滤波 250');
imwrite(gg10,'05 巴特沃兹高通加强滤波 250.jpg');

2017-11-16 10:38:43 hhyqhh 阅读数 2035
  • H.264/AVC的背景与基本概念

    通过本系列课程,观众有望理解视频压缩编码技术的整体发展,精通H.264视频编码技术的框架与细节,为进一步研究H.265/HEVC编码标准、音视频流媒体、视频直播点播等技术奠定坚实的基础。本系列课程适合多媒体方向的...

    2777课时 0分钟 5067人学习 殷汶杰
    免费试看

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


代码基于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

2014-02-23 21:37:24 u012627502 阅读数 19414
  • H.264/AVC的背景与基本概念

    通过本系列课程,观众有望理解视频压缩编码技术的整体发展,精通H.264视频编码技术的框架与细节,为进一步研究H.265/HEVC编码标准、音视频流媒体、视频直播点播等技术奠定坚实的基础。本系列课程适合多媒体方向的...

    2777课时 0分钟 5067人学习 殷汶杰
    免费试看

1、主要步骤:空域(傅里叶变换、卷积)>>>频域(与转移函数相乘、处理、傅里叶反变换)>>>空域

2、常用频域增强方法:

巴特沃斯滤波器:阶为n,截断频率为D0的转移函数为:

(1)低通滤波:

低通巴特沃斯滤波器在高低频率间的过渡比较光滑,所以得到的输出图其“振铃”现象不明显。

频域低通滤波器能消除虚假轮廓。

(2)高通滤波:就是利用滤波器的频率特性,让高频的通过,低频的无法通过,就好比在频率域设置阈值,频率域每一个频率分量有一个“幅值”,滤波器就好比在不同的频率分量给这个幅值乘以不同的增益,高通就好比高频部分增益为1,低频部分增益为0,当然这是理想高通。高斯高通滤波器就是频域每一个频率分量的增益的连接而成的曲线是一个高斯曲线

高通巴特沃斯滤波器.

G( x, y ) = g ( x, y ) + c * f( x, y )    c = 0.5,0<=c<=1

(3)带通和带阻滤波

带阻滤波器:阻止一定频率范围内的信号通过而允许其他频率范围内的信号通过。

带通滤波器:允许一定频率范围内的信号通过而阻止其它频率范围内的信号通过。带通和带阻互补。


低通、高通、带通、带阻等线性滤波器可以较好地消除线性叠加在图像上的加性噪声。但实际中,噪声和图像也常常是以非线性的方式结合。例如光源照明成像的情况,其中光的入射和景物的反射是以相乘的形式对成像做出贡献的,这样成像中的噪声与景物也是相乘的关系。


同态滤波:

是一种在频域中同时将图像亮度范围进行压缩和将图像对比度进行增强的方法。

f(x, y) = i(x, y) r(x, y)(一幅图像可以看做是  照度分量和反射分量的乘积)


ln f(x, y) = ln i(x, y) + ln r(x, y)(为了分开处理,两边分别取对数)

cvLog(ImgLog, ImgLog);


F(u, v) = I(u, v) + R(u, v)  (傅里叶变换)

cvDFT(Fourier, Fourier, CV_DXT_FORWARD);


H(u, v)F(u, v) = H(u, v)I(u, v) + H(u, v)R(u, v) (H(u, v)同态滤波函数,可分别作用在照度分量上和反射分量上)

cvMul(ImageRe,matH,ImageRe);
cvMul(ImageIm,matH,ImageIm);


hf (x,y) = hi (x,y)+ hr (x,y) (反变换)

cvDFT(Fourier, Fourier, CV_DXT_INV_SCALE);
cvSplit(Fourier,ImageRe, ImageIm, 0, 0);

g(x,y)=exp|hf (x,y)|=exp|hi (x,y)|exp|hr (x,y)|(前面去对数,在转变回来)

 cvExp(ImageRe, ImageRe);
 cvExp(ImageIm, ImageIm);

 cvMul(ImageRe, ImageIm, ImgDst);

最后归一化,显示。


同态滤波消噪过程:先利用非线性的对数变换将乘性噪声转化为加性噪声,然后就可用线性滤波器进行消除,最后再进行非线性的指数反变换以获得原始的无噪声图像。


同态滤波处理的流程如下:

S(x,y)------>Log---->FFT---->高通滤波---->IFFT---->Exp---->T(x,y)

其中S(x,y)表示原始图像;T( x,y)表示处理后的图像;Log 代表对数运算;FFT 代表傅立叶变换;IFFT 代表傅立叶逆变换;Exp 代表指数运算。


因为一般照度分量在空间变化较缓慢,而反射分量在不同物体交界处会急剧变化,所以图像对数的傅里叶变换后的低频部分主要对应照度分量,而高频部分对应反射分量。

这样我们可以设计一个对傅里叶变换结果的高频和低频分量影响不同的滤波函数H(u,v)。


高斯型高通滤波器修改形式:

式中,常数c 被用来控制滤波器函数斜面的锐化,它在之间。

OpenCV中,构造高斯型高通滤波器

CvMat *matH = cvCreateMat( height, width, CV_32FC1);
HPF( matH );

void HPF(CvMat *matH, float D0 = 10,float rH = 2.0,float rL = 0.3,float c = 0)
{
    if (D0 < 0){
        qDebug()<<"ERROR!";
        return ;
    }

    for (int u = 0; u < matH->rows; ++u)
        for (int v =0 ;v < matH->cols; ++v){
            float  D = (u * u )+(v * v);
            *( (float*)CV_MAT_ELEM_PTR(*matH, u, v) ) =
                    (rH - rL)*( 1 - exp(-(c * D * D)/(D0*D0)) ) + rL;//高斯型高斯滤波器
        }
}

》》》下文摘自:http://course.cug.edu.cn/rs/COURSE/6-3-4-a.HTM《《《

同态滤波增强是把频率过虑和灰度变换结合起来的一种处理方法。它是把图像的照明反射模型作为频域处理的基础,利用压缩灰度范围和增强对比度来改善图像的一种处理技术。它在密度域中运用相当成功。
一幅图像f(x,y)可以看成由两个分量组合而成,即
f(x,y)=i(x,y).r(x,y)
i(x,y)为照明分量(入射分量),是入射到景物上的光强度;
r(x,y)为反射分量,是受到景物反射的光强度。
具体步骤如下:
(1)先对上式的两边同时取对数,即
Inf(x,y)=Ini(x,y)+Inr(x,y) 
(2)将上式两边取傅立叶变换,得
F(u,v)=I(u,v)+R(u,v)
(3)用一个频域函数H(u,v)处理F(u,v),可得到
H(u,v)F(u,v)=H(u,v)I(u,v)+H(u,v)R(u,v)
(4)逆傅立叶变换到空间域得
Hff(x,y)=hi(x,y)+hr(x,y)
可见增强后得图像是由对应照度分量与反射分量得两部分叠加而成。
(5)再将上式两边取指数,得

g(x,y)=exp|hff(x,y)|=exp|hi(x,y)|exp|hr(x,y)|

这里,称作同态滤波函数,它可以分别作用于照度分量和反射分量上。
                                           

    一幅图像得照明分量通常用慢变化来表征,而反射分量则倾向于急剧变换。所以图像取对数后得傅立叶变换的低频部分主要对应照度分量,而高频部分主要对应反射分量。适当的选择滤波器函数将会对傅立叶变换中的低频部分和高频部分产生不同的响应。处理结果会使像元灰度的动态范围或图像对比度得到增强。