【六】频率域图像增强

1 频率域滤波

从空间域变换到频率域，傅里叶变换可以做到转换过程不丢失任何信息。

2 傅里叶变换

傅里叶级数：满足狄利赫里的正弦函数都可以用正弦函数和余弦函数构成的无穷级数。

傅里叶变换：傅里叶变换在这里不做介绍，它的本质就是基的变换，数字图像处理中主要关注二维离散函数的傅里叶变换(DFT)。

快速傅里叶变换：利用傅里叶变换中变量的周期性和对称性，可以进行简化运算，便于实现。（FFT）

2.1 频谱增强

Matlab实现

I1 = imread('bacteria.bmp');

fcoef = fft2(double(I1));%做fft变换
spectrum = fftshift(fcoef);%将零点移到中心
temp = log(1+abs(spectrum));%对幅值做对数变换以压缩动态范围

subplot(1,4,1),imshow(I1),title('source');
subplot(1,4,2),imshow(temp,[]),title('FFT1');

I2 = imread('babyNew.bmp');

fcoef2 = fft2(double(I2));%做fft变换
spectrum2 = fftshift(fcoef2);%将零点移到中心
temp2 = log(1+abs(spectrum2));%对幅值做对数变换以压缩动态范围

subplot(1,4,3),imshow(I2),title('source');
subplot(1,4,4),imshow(temp2,[]),title('FFT2');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

第一张图中对应FFT图的高频分量多 (如水平竖直的光)，第二张图对应FFT图的高频分量少

2.2 频谱和相位谱

Matlab实现

A = imread('beauty.jpg');
B = imread('cat.jpg');

% 求傅立叶变换
Af = fft2(double(A));
Bf = fft2(double(B));

% 分别求幅度谱和相位谱
AfA = abs(Af);
AfB = angle(Af);

BfA = abs(Bf);
BfB = angle(Bf);

% 交换相位谱并重建复数矩阵
AfR = AfA .* cos(BfB) + AfA .* sin(BfB) .* i;
BfR = BfA .* cos(AfB) + BfA .* sin(AfB) .* i;

% 傅立叶反变换
AR = abs(ifft2(AfR));
BR = abs(ifft2(BfR));

% 显示图像
subplot(2,2,1),imshow(A);
title('美女原图像');

subplot(2,2,2),imshow(B);
title('猫的原图像');

subplot(2,2,3),imshow(AR, []);
title('美女的幅度谱和猫的相位谱组合');

subplot(2,2,4),imshow(BR, []);
title('猫的幅度谱和美女的相位谱组合');
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

可以看出相位谱决定的是图像结构，幅度谱决定图像灰度分布。

3 频域滤波

这个要结合后面的滤镜才能起作用。
Matlab实现
频域滤波函数imfreqfilt.m

function out = imfreqfilt(I, ff)
% imfreqfilt函数			对灰度图像进行频域滤波
% 参数I				输入的空域图像
% 参数ff				应用的与原图像等大的频域滤镜

if (ndims(I)==3) && (size(I,3)==3)   % RGB图像
    I = rgb2gray(I);
end

if (size(I) ~= size(ff))
    msg1 = sprintf('%s: 滤镜与原图像不等大，检查输入', mfilename);
    msg2 = sprintf('%s: 滤波操作已经取消', mfilename);
    eid = sprintf('Images:%s:ImageSizeNotEqual',mfilename);
    error(eid,'%s %s',msg1,msg2);
end

% 快速傅立叶变换
f = fft2(double(I));

% 移动原点
s = fftshift(f);

% 应用滤镜及反变换
out = s .* ff; %对应元素相乘实现频域滤波
out = ifftshift(out);
out = ifft2(out);

% 求模值
out = abs(out);

% 归一化以便显示
out = out/max(out(:));
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

4 频域低通滤波器

频谱中，低频对应图像的整体灰度分布，高频对应细节。平滑是通过衰减高频部分实现。
Matlab实现
imidealflpf.m

function out = imidealflpf(I, freq)
% imidealflpf函数			构造理想的频域低通滤波器
% I参数				输入的灰度图像
% freq参数				低通滤波器的截止频率
% 返回值：out – 指定的理想低通滤波器
[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        if (sqrt(((i-M/2)^2+(j-N/2)^2))>freq)
            out(i,j)=0;
        end
    end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14

4.1 理想低通滤波器

I = imread('baby_noise.bmp'); %读入原图像

% 生成滤镜
ff = imidealflpf(I, 20);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,1);
imshow(I);
title('Source');

% 计算FFT并显示
temp = fft2(double(I));
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,1);
imshow(temp, []);
title('Source');

figure (1);
subplot(2,2,2);
imshow(out);
title('Ideal LPF, freq=20');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,2);
imshow(temp, []);
title(' Ideal LPF, freq=20');

% 生成滤镜
ff = imidealflpf(I, 40);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,3);
imshow(out);
title('Ideal LPF, freq=40');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,3);
imshow(temp, []);
title(' Ideal LPF, freq=40');

% 生成滤镜
ff = imidealflpf(I, 60);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,4);
imshow(out);
title('Ideal LPF, freq=60');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,4);
imshow(temp, []);
title(' Ideal LPF, freq=60');

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

4.2 高斯低通滤波器

Matlab实现
imgaussflpf.m

function out = imgaussflpf(I, sigma)
% imgaussflpf函数     		构造频域高斯低通滤波器
% I参数				输入的灰度图像
% sigma参数			高斯函数的Sigma参数

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);
    end
end
1
2
3
4
5
6
7
8
9
10
11
12

高斯低通滤波器实现

I = imread('baby_noise.bmp'); %读入原图像

% 生成滤镜
ff = imgaussflpf (I, 20);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,1);
imshow(I);
title('Source');

% 计算FFT并显示
temp = fft2(double(I));
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,1);
imshow(temp, []);
title('Source');

figure (1);
subplot(2,2,2);
imshow(out);
title('Gauss LPF, sigma=20');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,2);
imshow(temp, []);
title(' Gauss LPF, sigma=20');

% 生成滤镜
ff = imgaussflpf (I, 40);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,3);
imshow(out);
title('Gauss LPF, sigma =40');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,3);
imshow(temp, []);
title(' Gauss LPF, sigma =40');

% 生成滤镜
ff = imgaussflpf (I, 60);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,4);
imshow(out);
title('Gauss LPF, sigma =60');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,4);
imshow(temp, []);
title(' Gauss LPF, sigma =60');
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

5 频率域高斯滤波器

图像锐化就是通过衰减图像中低频成分来实现的。
Matlab实现
imgaussfhpf.m函数

function out = imgaussfhpf(I, sigma)
% imgaussfhpf函数			构造频域高斯高通滤波器
% I参数				输入的灰度图像
% sigma参数			高斯函数的Sigma参数

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = 1 - exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);
    end
end
1
2
3
4
5
6
7
8
9
10
11
12

5.1 高斯高通滤波实现

I = imread('coins.png');

% 生成滤镜
ff = imgaussfhpf (I, 20);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,1);
imshow(I);
title('Source');

% 计算FFT并显示
temp = fft2(double(I));
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,1);
imshow(temp, []);
title('Source');

figure (1);
subplot(2,2,2);
imshow(out);
title('Gauss HPF, sigma=20');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,2);
imshow(temp, []);
title(' Gauss HPF, sigma=20');

% 生成滤镜
ff = imgaussfhpf (I, 40);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,3);
imshow(out);
title('Gauss HPF, sigma =40');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,3);
imshow(temp, []);
title(' Gauss HPF, sigma =40');

% 生成滤镜
ff = imgaussfhpf (I, 60);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,4);
imshow(out);
title('Gauss HPF, sigma =60');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,4);
imshow(temp, []);
title(' Gauss HPF, sigma =60');
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

5.2 拉普拉斯滤波实现

Matlab实现
imlapf.m函数

function out = imlapf(I)
% imlapff函数			构造频域拉普拉斯滤波器
% I参数				输入的灰度图像

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = -((i-M/2)^2+(j-N/2)^2);
    end
end
1
2
3
4
5
6
7
8
9
10
11

拉普拉斯滤波实现

I = imread('bacteria.bmp');

ff = imlapf (I);
out = imfreqfilt(I, ff);

figure (1);
subplot(1,2,1);
imshow(I);
title('Source');

temp = fft2(double(I));
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(1,2,1);
imshow(temp, []);
title('Source');

figure (1);
subplot(1,2,2);
imshow(out);
title('Laplace Filter');

temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(1,2,2);
imshow(temp, []);
title('Laplace Filter');
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

6 实例——频域滤波消除周期噪声

6.1 带阻滤波

高斯带阻滤波器
Matlab实现
imgaussbrf.m函数

function out = imgaussfbrf(I, freq, width)
% imidealflpf函数 		构造频域高斯带阻滤波器
% I参数				输入的灰度图像
% freq参数				阻带中心频率
% width参数			阻带宽度

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = 1-exp(-0.5*((((i-M/2)^2+(j-N/2)^2)-freq^2)/(sqrt(i.^2+j.^2)*width))^2);
    end
end
1
2
3
4
5
6
7
8
9
10
11
12
13

带阻滤波实现

O = imread('pout.tif');
[M,N] = size(O);
I = O;
%这一部分是在添加噪声
for i=1:M
    for j=1:N
        I(i,j) = I(i,j)+20*sin(20*i)+20*sin(20*j);
    end
end
figure;
subplot(1,2,1),imshow(O);
title('原图像');
subplot(1,2,2),imshow(I);
title('添加噪声后的图像');

%获得频谱
i_f = fft2(I);
i_f = fftshift(i_f);
i_f = abs(i_f);
i_f = log(1+i_f);

o_f = fft2(O);
o_f = fftshift(o_f);
o_f = abs(o_f);
o_f = log(1+o_f);

figure;
subplot(1,4,1),imshow(i_f,[]);
title('原图像的频谱');
subplot(1,4,2),imshow(o_f,[]);
title('添加噪声后图像的频谱');


%带阻滤波
ff = imgaussfbrf(I,50,5);%构造高斯带阻滤波器
subplot(1,4,3),imshow(ff,[]);
title('高斯带阻滤波器');
out = imfreqfilt(I,ff);%带阻滤波
subplot(1,4,4),imshow(out,[]);

figure;
subplot(1,2,1),imshow(I);
title('高斯带阻滤波前');

subplot(1,2,2),imshow(out);
title('高斯带阻滤波后');
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
46