为什么要在频率域研究图像增强
可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务,在频率域中变得非常普通滤波在频率域更为直观,它可以解释空间域滤波的某些性质
给出一个问题,寻找某个滤波器解决该问题,频率域处理对于试验、迅速而全面地控制滤波器参数是一个理想工具
一旦找到一个特殊应用的滤波器,通常在空间域采用硬件实现它
傅里叶变换的频率分量和图像空间特征之间的联系
当从变换的原点移开时,低频对应着图像的慢变化分量,如图像的平滑部分
进一步离开原点时,较高的频率对应图像中变化越来越快的灰度级,如边缘或噪声等尖锐部分
二维DFT(二维离散傅里叶变换)的数学定义
二维DFT考虑M列N行的图像f(x,y)。在我们的符号中,x是水平坐标,x=0,1,2,…,M-1,y是垂直坐标,y=0,1,…,N-1。图像f(x,y)的二维离散傅立叶变换(2-D DFT)是一个矩阵,其项定义为:
其中u和v是整数频率指数,u=0,1,…,M-1和v=0,1,…,N-1。
二维DFT F(u,v)被认为是图像F(x,y)的频率表示。
公式中的归一化频率,u/M称为归一化水平频率。标准化水平频率有时简单地表示为u。
v/N称为标准化垂直频率。归一化垂直频率有时简单地表示为v。
区分这些表示:标准化频率的值介于[0,1)或[-0.5,0.5]之间,而频率指数具有整数值。
逆二维DFT
从其频率表示F(u,v)出发,应用逆二维DFT可以重建原始图像,定义为
注意方程中的正向变换和方程中的逆变换之间指数项的符号差。
分析二维DFT
一般来说,变换后的图像F(u,v)具有复值:,其中R(u,v)是实部,I(u,v)是虚部。
F(u,v)的幅度称为输入图像的幅度谱。这个幅度谱可以揭示输入图像的显著频率成分。
MATLAB中的二维DFT
a)在MATLAB中,可以使用fft2函数找到图像的二维DFT:
F = fft2(f);
(6.6)如果输入图像f具有N行M列,则输出F也将是具有N行M列的矩阵。 F中的每个条目都称为输入图像f的频率分量。
在解释函数fft2的输出时必须小心。 这是因为在MATLAB中,矩阵索引从1开始,而在数学定义(6.1)中,u和v从0开始。因此,如果F是函数fft2的MATLAB输出,则
F(1,1)是对应于v = 0和u = 0的频率分量。也就是说,F(1,1)是输入图像的DC频率分量。
F(3,8)是对应于v = 2和u = 7的频率分量,或者垂直归一化频率为2 / N,水平归一化频率为7 / M。 使用相同的规则来解释F的其他条目。
b)明确指定DFT大小
fft2函数具有另一种语法,可让我们指定转换后图像的大小:
F = fft2(f,H,W);
在这种情况下,输出矩阵F将具有H行和W列。 在评估DFT之前,输入图像f用零填充以创建具有H行和W列的图像。 高×宽称为DFT尺寸。 为了加快DFT计算(使用快速2D傅立叶变换),通常将H和W的值选择为2的幂。
c)分析二维DFT下面显示了如何在MATLAB中将二维DFT应用于图像。 执行程序并研究其步骤。
注意:在MATLAB中,可以分别使用函数abs和angle计算复数的幅度和相位。
由于幅度谱具有很大的动态范围,因此通常使用强度调整技术(对数和缩放)来可视化幅度谱。
% Compute 2D DFT of an image
f = double(imread('Lena.bmp')); % read an image
[H, W ] = size(f); % image size
F = fft2(f); % compute 2-D DFT
F_mag = abs(F); % compute the magnitude spectrum
% Visualize the magnitude spectrum
u = 0:1/W:1-1/W; % Horizontal normalized frequencies
v = 0:1/H:1-1/H; % Vertical normalized frequencies
imagesc(v, u, log10(1 + F_mag)); % Magnitude spectrum versus normalized frequencies
colorbar, colormap gray, axis image
title('Magnitude Spectrum');
xlabel('Horizontal normalized frequency u');
ylabel('Vertical normalized frequency v');
d)在MATLAB中计算逆二维DFT
,逆DFT可用ifft2函数计算:
f1 = real(ifft2(F));(6.8)
理论上,通过ifft2重建的图像与方程(6.6)中的原始图像F相同。
然而,由于数值舍入误差,函数ifft2可能产生复值矩阵。
因此,函数real()仅用于提取实值部分并丢弃虚部。
e) 移动上面中的DC组件,DC组件显示在图的左上角。为了可视化的目的,有时会移动频谱,使直流分量位于图形的中心。
这可以使用MATLAB fftshift函数:Fc=fftshift(F)来实现;fftshift函数的效果如图6.1所示。
从数学上讲,用fftshift得到的矩阵Fc(u,v)是以下图像的二维DFT:
中心变化
Fc = fftshift(F); % Shift DC component
% Horizontal normalized frequencies
u = [0:1/H:1-1/H] - 1/2; % Vertical normalized frequencies
v = [0:1/W:1-1/W] - 1/2;
Fc_mag = abs(Fc); % Magnitude spectrum
% Magnitude spectrum versus normalized frequencies
imagesc(v, u, log10(1 + Fc_mag));
colorbar, colormap gray, axis image, title('Magnitude spectrum');
xlabel('Horizontal normalized frequency u');
练习计算并显示三幅图像的幅度谱:hstripes.bmp, vstripes.bmp,和dstripes.bmp. 解释每个图像与其幅度谱之间的关系。
计算通过对hstripes.bmp. 重建图像与原始图像的中误差是多少hstripes.bmp?
使用二维DFT
a)空间域卷积的线性图像滤波等同于频域中的乘法:
- 假设f和h是两个图像,g是通过f和h的卷积获得的图像。
- 设F、H、G分别为图像F、H、G的DFT。
- 然后,通过将两个矩阵F和H逐项相乘得到G:(6.9)
该特性用于在频域中实现线性图像滤波。
计算图像f和卷积核h之间的卷积输出有四个主要步骤
第一步:计算图像f的二维DFT。
第二步:计算卷积核h的二维DFT。
将第1步和第2步中的输入和输出相乘。
第四步:对G进行二维DFT逆运算,得到输出图像G。
为了避免混叠,必须适当选择DFT大小,
- 假设图像f有Hf行和Wf列。
- 假设卷积核h有Hh行和Wh列。
- 然后,我们必须使用至少(Hf+Hh–1)行和(Wf+Wh–1)列的DFT大小。
因此,必须使用方程式(6.7)中所示的fft2的扩展语法
function g = imfilter_dft(f, h)
[Hf, Wf] = size(f); % input image size
[Hh, Wh] = size(h); % convolution kernel size
P = Hf + Hh - 1; % DFT size vertical
Q = Wf + Wh - 1; % DFT size horizontal
f = double(f);
F = fft2(f, P, Q); % DFT of input image
H = fft2(h, P, Q); % DFT of convolution kernel
G = F .* H; % DFT of output image
g_inverse = real(ifft2(G)); % Inverse DFT
% Extract the output image of size HfxWf from g_inverse
offset_y = floor((Hh + 1)/2);
offset_x = floor((Wh + 1)/2);
g = g_inverse(offset_y:offset_y + Hf - 1, ...
offset_x:offset_x + Wf - 1);
g = uint8(g); % Convert to unsigned 8-bit image
c) 利用imfilter_dft.m函数,编写MATLAB程序计算并显示在图像上应用5×5平均滤波器时的输出图像lena.bmp.
d) 比较获得的输出图像(b)和在实验室5中使用函数imfilter.m获得的卷积输出。
e) 利用MATLAB函数tic和toc,测量imfilter_dft和imfilter对图像进行滤波所需的时间lena.bmp. 阅读有关如何使用tic和toc的在线文档。
f) 比较平均滤波器尺寸为99×99像素时imfilter_dft和imfilter的处理时间。
频域通用图像滤波器
a)执行MATLAB命令open('视频脚本.slx)观看频域图像滤波的演示;此演示使用MATLAB计算机视觉工具箱实现。单击“开始模拟”按钮。
b) 当滤波器在频域内给定时,用二维DFT进行滤波更为合适。
在这里,我们学习几种常见的图像滤波器,它们的规格在频域中给出:
- 低通图像滤波器衰减图像的高频分量。也就是说,它将平滑图像。
- 高通图像滤波器衰减低频分量。也就是说,它将锐化图像中的边缘。
- 带通图像滤波器只保留属于频带的频率分量。
- 带阻图像滤波器拒绝属于频带的频率分量
c)6.4中的MATLAB程序生成一个具有截止频率的理想低通滤波器,并将其应用于噪声图像。它使用一个名为dftuv的函数;该函数包含在文件夹“Lab-06”中。
函数dftuv生成标准化频率的二维网格,其范围为,因此可以计算滤波器规范中的表达式。函数dftuv的语法是[v,u]=dftuv(H,W),其中H×W是DFT大小(H=行数,W=列数)。
执行清单6.4中的代码并解释其步骤。注意滤波器的频率响应是如何在矩阵Hd中创建的
E。频域中的图像滤波
a)修改清单6.4中的代码以使用理想的高通滤波器(见表6.1)。找到合适的截止频率并显示输出图像。
输入图像是lena.bmp
b) 在图像Lena-gaussian上应用三阶Butterworth低通滤波器计算输出图像-噪声.bmp. 确定合适的截止频率。
c) 编写了一个MATLAB程序,在频域内实现了一个高boost滤波器。滤波器在频域中的描述如下: