MATLAB频域图像处理

为什么要在频率域研究图像增强

可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务，在频率域中变得非常普通滤波在频率域更为直观，它可以解释空间域滤波的某些性质

给出一个问题，寻找某个滤波器解决该问题，频率域处理对于试验、迅速而全面地控制滤波器参数是一个理想工具

一旦找到一个特殊应用的滤波器，通常在空间域采用硬件实现它

傅里叶变换的频率分量和图像空间特征之间的联系

当从变换的原点移开时，低频对应着图像的慢变化分量，如图像的平滑部分

进一步离开原点时，较高的频率对应图像中变化越来越快的灰度级，如边缘或噪声等尖锐部分

二维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 horizontalf = 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 imageg_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滤波器。滤波器在频域中的描述如下：