频域处理与频率滤波

Posted on

频域处理

二维DFT
1
2
3
4
5
6
7
8
9
10
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0403(a)(image).tif');
figure,imshow(f);
F=fft2(f)
figure,imshow(abs(F),[])
Fc=fftshift(F)%中心化
figure,imshow(abs(Fc),[])
figure,imshow(log(1+abs(Fc)),[])%对数放大
F=ifftshift(Fc)%反中心化
phi=angle(F);%相角
figure,imshow(phi,[])

原图及DFT：

image-20200422215336399image-20200422215357520

频谱居中及其对数放大，相角

image-20200422215420505image-20200422215434693image-20200422215449459

频域滤波
未填充
1
2
3
4
5
6
7
8
9
10
11
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0405(a)(square_original).tif');
figure,imshow(f);
[M,N]=size(f);
F=fft2(f);
figure,imshow(log(1+abs(fftshift(F)),[]);
sig=10;
H=lpfilter('gaussian',M,N,sig);
figure,imshow(abs(fftshift(H)),[]);%对数放大
G=H.*F;
g=ifft2(G);
figure,imshow(g,[])

原图及其傅里叶变换：

image-20200422221855557image-20200422222633292

高斯滤波器：

image-20200422222705552

滤波结果：左右边缘并没有计算到

image-20200422222719660

有填充

填充函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering. 
%   PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
%   computes the two-element size vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
%   vectors, computes the two-element size vector PQ.  The elements
%   of PQ are the smallest even integers greater than or equal to 
%   AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]). 
    
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $

if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 & ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.
   
   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
elseif nargin == 3
   m = max([AB CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
else 
   error('Wrong number of inputs.')
end
1
2
3
4
5
6
7
8
9
10
11
12
13
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0405(a)(square_original).tif');
figure,imshow(f);
[M,N]=size(f);
PQ=paddedsize(size(f));
Fp=fft2(f,PQ(1),PQ(2));
sig=10;
Hp=lpfilter('gaussian',PQ(1),PQ(2),2*sig);
figure,imshow(abs(fftshift(Hp)),[]);
Gp=Hp.*Fp
gp=ifft2(Gp);
figure,imshow(gp,[])
gpc=gp(1:M,1:N)
figure,imshow(gpc,[])

原图：

image-20200422224052640

高斯滤波器：

image-20200422224104340

image-20200422224119283

填充过后滤波：

image-20200422224126930

切割得到最终结果：此时左右边缘均模糊

image-20200422224139562

频域滤波封装成一个函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function g = dftfilt(f, H,classout)

[f,revertClass]=tofloat(f);

% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));

% Perform filtering. 
g = ifft2(H.*F);

% Crop to original size.
g = g(1:size(f,1), 1:size(f,2));

if nargin==2||strcmp(classout,'original')
    g=revertClass(g);
elseif strcmp(classout,'fltpoint')
        return
else
    error('undefined class for the out put')
end
从空域滤波器获得频域滤波器
1
2
3
4
5
6
figure
h=fspecial('sobel')%[1 0 -1;2 0 -2;1 0 -1]
H=freqz2(h);
H1=abs(H);
mesh(double(H1(1:1:64,1:1:64)))
axis tight

image-20200422231230346

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0409(a)(bld).tif');
figure,imshow(f);
f=tofloat(f);
F=fft2(f);
S=fftshift(log(1+abs(F)));%对数拉伸
figure,imshow(S,[])
h=fspecial('sobel')%[1 0 -1;2 0 -2;1 0 -1]
PQ=paddedsize(size(f));
H=freqz2(h,PQ(1),PQ(2));
H1=ifftshift(H);
figure,imshow(abs(H),[])
figure,imshow(abs(H1),[])
gs=imfilter(f,h)
figure,imshow(abs(gs)>0.2*abs(max(gs(:))),[])
gf=dftfilt(f,H1);
figure,imshow(abs(gf)>0.2*abs(max(gf(:))),[])

原图及其傅里叶变换：

image-20200422230834371image-20200422230849528

sobel滤波器及其中心化：

image-20200422230902229image-20200422230933140

采用空域滤波与频域滤波：两者一致

image-20200422231129912image-20200422231143826

在频域直接生成滤波器

建立网格数组：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function [U, V] = dftuv(M, N)
%DFTUV Computes meshgrid frequency matrices.
%   [U, V] = DFTUV(M, N) computes meshgrid frequency matrices U and
%   V.  U and V are useful for computing frequency-domain filter
%   functions that can be used with DFTFILT.  U and V are both
%   M-by-N.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.3 $  $Date: 2003/04/16 22:30:34 $

% 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);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
>> [U,V]=dftuv(5,5);
>> DSQ=U.^2+V.^2

DSQ =

     0     1     4     4     1
     1     2     5     5     2
     4     5     8     8     5
     4     5     8     8     5
     1     2     5     5     2
>> D=hypot(U,V)%开根号，运行速度比sqrt快

D =

         0    1.0000    2.0000    2.0000    1.0000
    1.0000    1.4142    2.2361    2.2361    1.4142
    2.0000    2.2361    2.8284    2.8284    2.2361
    2.0000    2.2361    2.8284    2.8284    2.2361
    1.0000    1.4142    2.2361    2.2361    1.4142
>> fftshift(DSQ)%每个点到获得到中点的距离平方

ans =

     8     5     4     5     8
     5     2     1     2     5
     4     1     0     1     4
     5     2     1     2     5
     8     5     4     5     8

用生成的高斯滤波器滤波：

1
2
3
4
5
6
7
8
9
10
11
12
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
figure,imshow(f);
PQ=paddedsize(size(f));
[U,V]=dftuv(PQ(1),PQ(2));
D=hypot(U,V);
D0=0.05*PQ(2);
F=fft2(f,PQ(1),PQ(2));
H=exp(-(D.^2)/(2*(D0^2)));
g=dftfilt(f,H);
figure,imshow(fftshift(H));
figure,imshow(log(1+abs(fftshift(F))),[])
figure,imshow(g)

原图及其傅里叶变换：

image-20200422233840072image-20200422233947493

滤波器与滤波结果：

image-20200422233854576image-20200422234058343

封装函数：产生几种

低通滤波器
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function H = lpfilter(type, M, N, D0, n)
%LPFILTER Computes frequency domain lowpass filters.
%   H = LPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a lowpass filter, H, of the specified TYPE and size (M-by-N). To
%   view the filter as an image or mesh plot, it should be centered
%   using H = fftshift(H). 
%
%   Valid values for TYPE, D0, and n are:
%
%   'ideal'    Ideal lowpass filter with cutoff frequency D0. n need
%              not be supplied.  D0 must be positive.
%
%   'btw'      Butterworth lowpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0.  D0 must be
%              positive.
%
%   'gaussian' Gaussian lowpass filter with cutoff (standard
%              deviation) D0.  n need not be supplied.  D0 must be
%              positive. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.8 $  $Date: 2004/11/04 22:33:16 $

% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances. 
[U, V] = dftuv(M, N);

% Compute the distances D(U, V).
D = hypot(U,V);

% Begin filter computations.
switch type
case 'ideal'
   H = single(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
高通滤波器：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
function H = hpfilter(type, M, N, D0, n)
%HPFILTER Computes frequency domain highpass filters.
%   H = HPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a highpass filter, H, of the specified TYPE and size (M-by-N).
%   Valid values for TYPE, D0, and n are: 
%
%   'ideal'    Ideal highpass filter with cutoff frequency D0.  n
%              need not be supplied. D0 must be positive.
%
%   'btw'      Butterworth highpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0. D0 must be
%              positive.
%
%   'gaussian' Gaussian highpass filter with cutoff (standard
%              deviation) D0. n need not be supplied. D0 must be
%              positive.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/08/25 14:28:22 $

% The transfer function Hhp of a highpass filter is 1 - Hlp, 
% where Hlp is the transfer function of the corresponding lowpass 
% filter.  Thus, we can use function lpfilter to generate highpass 
% filters.

if nargin == 4
   n = 1; % Default value of n.
end

% Generate highpass filter.
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;

image-20200423002127952image-20200423002153603image-20200423002201364

高通滤波：

1
2
3
4
5
6
7
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
figure,imshow(f);
PQ=paddedsize(size(f));
D0=0.05*PQ(1);
H=hpfilter('gaussian',PQ(1),PQ(2),D0);
g=dftfilt(f,H);
figure,imshow(g)

image-20200423002422824image-20200423002436507

高频强调滤波：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
f=imread('C:\Users\HP\Desktop\ImageDatabase\dipum_images_ch04\Fig0419(a)(chestXray_original).tif');
figure,imshow(f);
PQ=paddedsize(size(f));
D0=0.05*PQ(1);
HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);
H=0.5+0.2*HBW;
gbw=dftfilt(f,HBW,'fltpoint');
gbw=gscale(gbw);
figure,imshow(gbw)
ghf=dftfilt(f,H,'fltpoint');
ghf=gscale(ghf);
figure,imshow(ghf)
ghe=histeq(ghf,256);
figure,imshow(ghe)

image-20200423003714304image-20200423003730349image-20200423003743152image-20200423003816656

选择性滤波

陷波带阻与带通滤波器：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
function H=bandfilter(type,band,M,N,D0,W,n)
    [U,V]=dftuv(M,N);
    D=hypot(U,V);
    if nargin<7
        n=1;
    end
    switch lower(type)
        case 'ideal'
            H=idealReject(D,D0,W);
        case 'btw'
            H=btwReject(D,D0,W,n);
        case 'gaussian'
            H=gaussReject(D,D0,W);
        otherwise
            error('Unknown filter type.')
    end
    if strcmp(band,'pass')
        H=1-H;
    end     
end
%--------------------------------------%
function  H=idealReject(D,D0,W)
    RI=D<D0-(W/2);
    RO=D>D0+(W/2);
    H=tofloat(RO|RI);
end
%--------------------------------------%
function H=btwReject(D,D0,W,n)
    H=1./(1+(((D*W)./(D.^2-D0.^2)).^2*n));
end
%--------------------------------------%
function H=gaussReject(D,D0,W)
    H=1-exp(-((D.^2-D0.^2)./(D.*W+eps)).^2);
end

‘ideal’

image-20200423010805885image-20200423010746577

‘btw’

image-20200423010949531image-20200423011006651

‘gaussian’

image-20200423011121141image-20200423011053633