【图像处理lec5、6】频域滤波与图像增强

一维傅里叶变换：
- 正变换 $F(u) = \int_{-\infty}^{\infty} f(x) e^{-j2\pi ux} dx$
- 逆变换 $f(x) = \int_{-\infty}^{\infty} F(u) e^{j2\pi ux} du$
二维傅里叶变换：
- 用于处理图像，公式为： $F(u,v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y) e^{-j2\pi(ux + vy)} dx dy$
- 逆变换： $f(x,y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u,v) e^{j2\pi(ux + vy)} du dv$
解析：
- 图像在频域表示为一系列频率分量，方便进行滤波和增强。

从图中可以看出，方波的傅里叶变换是Sa函数，且时域越 “胖” ，频域就越 “瘦” 。

3、极坐标表示与功率谱

傅里叶变换在极坐标下可表示为： $F(u) = |F(u)| e^{j\phi(u)}$
成分解析：
- $|F(u)|= \sqrt{R^2(u) + I^2(u)}$ ：幅度或频谱，表示各频率分量的强度。
- $\phi(u)= \tan^{-1} \left( \frac{I(u)}{R(u)} \right)$ ：相位谱，表示各分量的相位关系。
功率谱：
- $P(u) = |F(u)|^2 = R^2(u) + I^2(u)$
- 反映频率分量的能量分布。

4、二维离散傅里叶变换 (DFT)

（1）公式

离散傅里叶变换（DFT）：
- 一维DFT：将连续变换离散化，用于数字信号处理。 $F(u) = \frac{1}{M} \sum_{x=0}^{M-1} f(x) e^{-j2\pi ux/M}$
- 二维DFT：扩展到图像。 $F(u,v) = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{-j2\pi (ux/M + vy/N)}$
频谱表示：
- 幅度 $|F(u)| = \sqrt{R^2(u) + I^2(u)}$
- 相位角 $|F(u)| = \sqrt{R^2(u) + I^2(u)}$
- 功率谱 $P(u) = |F(u)|^2$

（2）离散傅里叶变换 (DFT) 的频谱展示

一维

显示了一个傅里叶频谱的特性，分为：
1. (a) 原始傅里叶频谱：展示了对称的频谱，两侧的半周期相对称。0频率在0处和M-1处。
2. (b) 居中频谱：通过将输入信号乘以 $(-1)^{x+y}$ ，将频谱的零频成分移动到中心。（这里的中心是指 $[0,2\pi ]$ 的中心）
理解：
- DFT 的频谱是周期性的，即频谱重复出现。
- 中心频谱通过平移将零频移到中心，有助于观察频谱的结构。

二维

显示了一个傅里叶频谱的特性，分为：
1. (a) 原始傅里叶频谱：展示了对称的频谱，两侧的半周期相对称。(0频在原点和M-1处，即在正方形的边缘)
2. (b) 居中频谱：通过将输入信号乘以 $(-1)^{x+y}$ ，将频谱的零频成分移动正方形的到中心。

需要频谱中心化的原因主要是因为变换后的原点不在正方形的中心，而在正方形的左上角；当然如果变换时在原点周围截取一个正方形的话，原点就直接在正方形中心了，就不用频谱中心化了。

理解：
- DFT 的频谱是周期性的，即频谱重复出现。
- 中心频谱通过平移将零频移到中心，有助于观察频谱的结构。

（3）常见二维的傅里叶变换对

冲激函数与常数：
- 冲激函数 $\delta(x,y)$ 的傅里叶变换是常数 1： $\mathfrak{F}[\delta(x,y)] \Leftrightarrow 1$
高斯函数的傅里叶变换：
- 空间域中的高斯函数在频域中仍为高斯函数，且宽度与标准差 $\sigma$ 有反比关系： $A2\pi\sigma^2 \exp(-2\pi^2\sigma^2(x^2 + y^2)) \Leftrightarrow A \exp\left(-\frac{u^2 + v^2}{2\sigma^2}\right)$
正弦和余弦的变换：
- 空间正弦/余弦在频域中表现为两个对称的冲激函数： $\cos(2\pi u_0 x + 2\pi v_0 y) \Leftrightarrow \frac{1}{2} \left[ \delta(u+u_0, v+v_0) + \delta(u-u_0, v-v_0) \right]$
- 正弦会引入相位变化： $\sin(2\pi u_0 x + 2\pi v_0 y) \Leftrightarrow \frac{1}{2} j \left[ \delta(u+u_0, v+v_0) - \delta(u-u_0, v-v_0) \right]$
总结：
- 常见傅里叶对提供了在空间域与频域间进行信号变换的基础规则。
- 对称性质和冲激函数的变换是频域分析的关键。

（4）二维傅里叶变换示例

图像与傅里叶谱：
- 左图：空间域中一个 $20 \times 40$ 的白色矩形位于 $512 \times 512$ 的黑色背景上。
- 右图：傅里叶变换后，显示出明显的中心对称 高频分量 和 低频分量。
特征解释：
- 低频分量：傅里叶变换结果中，图像的能量主要集中在中心区域。这表示图像的全局形状由低频分量表示。
- 高频分量：矩形的边缘和细节表现为远离中心的高频分量。
- 对称性：傅里叶谱具有中心对称性质（负频率和正频率对称）。

5、二维离散傅里叶变换（DFT）的性质

（1）时间移位和频率移位性质

时间移位
若 $f(x - x_0, y - y_0)$ 是图像 $f(x, y)$ 的移位版本，则对应的频域表示为：
$\mathcal{F}[f(x - x_0, y - y_0)] = F(u, v)e^{-j2\pi \left( \frac{ux_0}{M} + \frac{vy_0}{N} \right)}$
说明时域的平移会引入频域中的相位变化，但幅度保持不变。
频域平移
若时域信号 $f(x, y)$ 乘以 $e^{-j2\pi \left( \frac{u_0 x}{M} + \frac{v_0 y}{N} \right)}$ ，则其傅里叶变换在频域平移为 $F(u - u_0, v - v_0)$ ：
$\mathcal{F} \left[ f(x, y) e^{-j2\pi \left( \frac{u_0 x}{M} + \frac{v_0 y}{N} \right)} \right] = F(u - u_0, v - v_0)$
中心化处理
若时域信号 $f(x, y)$ 乘以 $(-1)^{x+y}$ ，对应的频域 $F(u, v)$ 平移到中心 $\left( u - \frac{M}{2}, v - \frac{N}{2} \right)$ ：
$\mathcal{F} \left[ f(x, y)(-1)^{x+y} \right] = F \left( u - \frac{M}{2}, v - \frac{N}{2} \right)$

（2）平均值与对称性性质

平均值
图像的平均值可以通过 $F(0,0)$ 得到，表示频域中 $F(0,0)$ 的值等于时域中 $f(x,y)$ 的平均值：
$F(0,0) = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)$
共轭对称性
对于实值图像，频域的性质：
$F(u, v) = F^*(-u, -v)$
幅度满足对称性：
$|F(u, v)| = |F(-u, -v)|$

（3）可分离性性质

二维DFT可以分解为两次一维DFT：

步骤：
1. 对每一行图像 $f(x, y)$ 进行一维 DFT，得到 $F(u, y)$ 。
2. 对结果 $F(u, y)$ 的每一列进行一维 DFT，得到 $F(u, v)$ 。
这种可分离性使得二维 DFT的计算效率更高，可以使用FFT（快速傅里叶变换）实现。

（4）旋转性质

旋转性质表明，当时域图像 $f(r, \theta + \theta_0)$ 旋转一定角度 $\theta_0$ 时，对应的频域 $F(\omega, \phi + \theta_0)$ 也发生相同角度的旋转。
直观效果：时域图像的直线旋转，频域中的条纹图案也随之旋转。

（5）周期性和线性性质

周期性
$f(x, y)$ 和 $F(u, v)$ 都是周期的，具有如下性质：
$f(x + M, y) = f(x, y), \quad F(u + M, v) = F(u, v)$
线性性质
傅里叶变换是线性的，线性组合在时域对应于频域的线性组合：
$\mathcal{F}[af(x, y) + bg(x, y)] = a\mathcal{F}[f(x, y)] + b\mathcal{F}[g(x, y)]$

（6）微分性质

在时域进行微分操作，频域中对应乘以 $j2\pi u$ 的幂次项：
$\mathcal{F} \left[ \frac{\partial^n f(x, y)}{\partial x^n} \right] = (j2\pi u)^n \mathcal{F}[f(x, y)]$
此性质说明傅里叶变换在边缘检测和高通滤波等微分运算中有重要应用。
拉普拉斯算子
对于二阶微分拉普拉斯算子 $\nabla^2 f(x, y)$ ，频域中变换为：
$\mathcal{F}(\nabla^2 f(x, y)) = -4\pi^2 (u^2 + v^2) F(u, v)$
拉普拉斯变换放大高频分量，适用于图像锐化。

（7）卷积性质

空间域卷积 $f(x,y) * g(x,y)$ 在频域中变成 频域乘法 $F(u,v)G(u,v)$ 。
数学公式： $F(u,v)G(u,v)$
结论：在时域/空间域进行卷积等同于在频域进行简单乘法运算。

（8）相关性质

相关运算的频域表示中涉及共轭运算。
自相关： $f(x,y) \circ f(x,y)$ 的傅里叶变换等于 功率谱 $|F(u,v)|^2$ 。
数学公式： $\mathfrak{F} [f(x,y) \circ f(x,y)] = |F(u,v)|^2$
结论：傅里叶变换中相关运算可直接通过频谱幅度计算。

（9）相似性性质

空间域缩放 $f(ax, by)$ 会导致频域中缩放反向变化并乘以比例系数 $\frac{1}{|ab|}$ .
数学公式： $\mathfrak{F}[f(ax,by)] = \frac{1}{|ab|} F\left(\frac{u}{a}, \frac{v}{b}\right)$
结论：缩放信号的空间坐标会影响频率分布的位置和密度。

6、二维（DFT）在 MATLAB 中的实现

（1）通过 `fft2` 函数计算二维图像的 DFT

通过 MATLAB 中的 fft2 函数计算二维图像的 DFT。
公式： F = fft2(f)，其中 f 是输入图像，F 是其 DFT 结果。
说明：DFT 返回的结果是一个大小与输入图像相同的复数矩阵，数据起始点位于图像的左上角。

代码步骤：

f = imread(...);       %读取图像
F = fft2(f);           %计算 DFT
S = abs(F);            %获取幅度谱
imshow(S, []);         %显示结果

结果：直接计算的 DFT 图像低频分量位于四个角落，高频分量位于中间，如下图：

（2）DFT 结果中心化

使用 fftshift 对 DFT 结果进行中心化处理：
- 函数：Fc = fftshift(F)
- 作用：将 DFT 结果的低频分量移动到图像中心，高频分量移到边缘。
- 说明：
  - 使用 fftshift 相当于将输入图像乘以 $(-1)^{x+y}$ ，实现低频分量的中心化。
  - 原始 fft2 结果不适合直接观察，fftshift 更直观显示频谱特征。
- 代码：

Fc=fftshift(F);
imshow(abs(Fc),[])

- 结果：中心化后的图像显示了十字状的低频能量集中在中心，如下图：

（3）对数变换增强频谱可视化

操作：S2 = log(1 + abs(Fc))
作用：通过取对数操作，将幅度谱数据进行压缩，使得低幅度和高幅度的值都能更清晰地展示出来。
代码：

S2=log(1+abs(Fc)); 
imshow(S2,[])

结果：对数变换后的频谱图显示出更清晰的结构，如下图：

（4）逆变换与重建

逆 DFT 及注意事项：
- 函数：ifft2(F) 计算 DFT 的逆变换。
- 语法： f = ifft2(F)。
- 实际输出：由于计算精度问题，逆变换结果通常包含非常小的虚数部分，因此需要使用 real 函数提取实部：
```
f = real(ifft2(F))
```
- 结果：输入图像可以通过逆变换重建，同时去除微小误差。

二、频域滤波与图像增强

1、频域滤波基本概念

在频域进行图像增强的关键步骤：
1. 傅里叶变换：将图像从时域转换到频域。
2. 频域处理：应用滤波器（低通、高通、带通等）。
3. 逆傅里叶变换：将处理后的图像转换回时域。
频域滤波的类型：
- 低通滤波：保留低频分量，平滑图像，去除高频噪声。
- 高通滤波：保留高频分量，突出图像边缘和细节。
- 带通滤波：仅保留某一频段的分量，适用于特定频率的增强。

2、频域滤波原理

卷积定理：
线性滤波的基础是卷积定理，即： $f(x,y) * h(x,y) \iff F(u,v)H(u,v)$ 这意味着在频域中，空间域的卷积操作可以简化为乘法操作。
滤波思想：
在频域滤波中，通过设计一个滤波器传递函数 $H(u,v)$ 并与图像的频谱 $F(u,v)$ 相乘来修改频率成分。
低通滤波器示例（图中图4.4）：
- 图4.4(a) 显示了低通滤波器的理想中心格式。
- 图4.4(b) 显示了DFT实现中使用的频域滤波器。
- 如果原图DFT后进行了中心化，则使用（a）图进行滤波；如果没有进行中心化，则使用（b）图进行滤波。

3、卷积与包裹误差

包裹误差（Wraparound Error）：
由于DFT对输入图像进行周期性延拓，非零部分的周期性延拓会相互干扰，形成误差。
解决方法：
为避免这种干扰，需对输入函数进行零填充（Padding）。
低通滤波示例：
图4.5显示了使用低通滤波器的结果：
- 左侧图：未进行零填充的滤波结果。
- 右侧图：进行了零填充的滤波结果，避免了边界干扰。这里并没有绝对避免，0填充相当于都是黑色填充，对黑色区域避免了干扰，对白色区域并没有避免干扰，反而还加大了干扰。

图4.6：
- 左侧：未进行零填充的图像频谱，显示了周期性延拓的干扰。
- 右侧：对图像进行零填充后，频谱没有边界干扰，便于滤波操作。
图4.7：
对于零填充后的图像，经过ifft2变换得到一个完整的、无边界干扰的滤波结果。（这里有疑问？）

4、低通滤波的MATLAB实现

通过MATLAB的具体代码实现频域低通滤波：

f = imread('Fig0405(a)(square original).tif');
[m, n] = size(f);        % 获取图像尺寸
F = fft2(f);             % 计算2D FFT
H = lpfilter('gaussian', m, n, 10); % 生成高斯低通滤波器
G = H .* F;              % 滤波器与频谱相乘
g = real(ifft2(G));      % 逆FFT还原图像，取实部
imshow(g, []);

结果：

展示了MATLAB滤波后得到的图像，说明在未零填充的情况下，图像的频域滤波效果。

5、DFT滤波基本步骤

DFT滤波的步骤：
1. 零填充：使用paddedsizes函数计算合适的填充尺寸。
```
PQ = paddedsize(size(f));
```
2. 计算傅里叶变换：
```
F = fft2(f, PQ(1), PQ(2));
```
3. 生成滤波器（如低通滤波器 H）。
4. 与频谱相乘：
```
G = H .* F;
```
5. 计算逆FFT：
```
g = real(ifft2(G));
```
6. 裁剪图像：
```
g = g(1:size(f,1), 1:size(f,2));
```
图4.8：展示了频域滤波的完整流程图，包括预处理、傅里叶变换、滤波、逆变换和后处理。

6、频域滤波器相较于空间域滤波器的优势

效率：频域滤波将复杂的空间域卷积简化为频域乘法，计算速度更快。
解释：
- 空间域卷积计算复杂度较高，尤其是大尺寸图像。而频域通过傅里叶变换（FFT）将图像和滤波器转到频域，再通过简单的乘法完成滤波，效率显著提高。

7、如何实现频率滤波器

MATLAB中将空间域滤波器转换为频域滤波器的完整过程：

f = imread('Fig0409(a)(bld).tif');   %读取图像：
F = fft2(f);                         % 2D快速傅里叶变换
h = fspecial('sobel');               % 通过 fspecial 创建空间域滤波器（例如 Sobel）
H = freqz2(h);                       % 使用 freqz2 将其转换到频域
H1 = ifftshift(H);                   %使用 ifftshift 将滤波器中心化
mesh(abs(H1));                       %可视化滤波器，使用 mesh 可视化滤波器的频率响应

可视化滤波器响应：
- 使用 mesh 显示频域滤波器的三维响应图，直观展示滤波器的作用范围。

8、频域滤波器的应用与可视化

如何通过傅里叶变换查看图像的频率分布，并验证滤波器效果。

步骤：
- 对图像进行傅里叶变换，使用对数变换增强显示：
```
f=imread(’Fig0409(a)(bld).tif’);
F=fft2(f);
S=fftshift(log(1+abs(F)));
S=gscale(S);
imshow(S)
```
- 观察频谱中心（低频分量）与边缘（高频分量）。

9、典型的频域低通滤波器

（1）理想低通滤波器 (ILPF)

定义：理想低通滤波器 $H(u,v)$ 的传递函数如下：
$H(u,v) = \begin{cases} 1 & \text{if } D(u,v) \leq D_0 \\ 0 & \text{if } D(u,v) > D_0 \end{cases}$
其中 $D(u,v)$ 是频率点 $(u,v)$ 到频域中心 $(M/2, N/2)$ 的距离：
$D(u,v) = \left[ \left( u - \frac{M}{2} \right)^2 + \left( v - \frac{N}{2} \right)^2 \right]^{\frac{1}{2}}$
解释：
- $D_0$ ：截止频率，控制滤波器允许的频率范围。
- 其在频域中表现为一个“硬边界”圆盘，在 $D_0$ 范围内保持 1，外部为 0。
- 图中展示了 ILPF 的三维透视图、平面图（圆形滤波区域）和沿径向的剖面线。
特性：
- 理想但不实际，因为其会产生明显的振铃效应，或说年轮效应（频域硬截止会导致时域的边缘振荡）。

理想低通滤波器的时域是Sa函数，旁瓣面积大，容易出现年轮效应。

（2）巴特沃斯低通滤波器 (BLPF)

定义：巴特沃斯低通滤波器 $H(u,v)$ 具有平滑的频率衰减特性：
- $D_0$ ：截止频率，决定滤波器的半径大小。
- $n$ ：滤波器的阶数，控制频率响应的平滑程度。

图示：
- 透视图显示 $n$ 的影响，随着 $n$ 增加，滤波器从平滑过渡到接近理想低通滤波器。
- $n$ 小：滤波器的响应较平滑，频域衰减缓慢。
- $n$ 大：响应更加接近理想低通，但边缘会产生轻微振铃。
特性：
- 巴特沃斯滤波器避免了 ILPF 的硬边界，提供了较平滑的频率衰减。
- 高阶 $n$ 会引入振铃现象，但相对于 ILPF 较轻。
阶数的影响
- 随着滤波器阶数 $n$ 的增加，滤波器的截止区域变得更清晰，滤波器的过渡带更窄。
- 但高阶滤波器会在频域边缘产生振荡，时域表现为边缘的振铃效应。
- 低阶滤波器：平滑但过渡带较宽，适用于对边缘细节要求不高的应用。
- 高阶滤波器：频域截止锐利，但易引起时域的振铃。

阶数越高，时域越接近Sa函数，旁瓣越大，卷积时越容易出现年轮效应。

（3）高斯低通滤波器

滤波器的定义

高斯低通滤波器的传递函数 $H(u,v)$ 由以下公式定义：

$H(u,v) = e^{-\frac{D^2(u,v)}{2D_0^2}}$

其中：

$D(u,v)$ ：频域中点 $(u,v)$ 到中心点的欧几里得距离。
$D_0$ ：截止频率（控制滤波器的扩展半径）。

高斯低通滤波器的特性
- 平滑的频率衰减：
  - 由于高斯函数的性质，GLPF 的衰减是连续的，没有明显的截止点。
  - 这种平滑的衰减有效地减少了时域中的“振铃效应”。
- 无振铃效应：
  - 理想低通滤波器（ILPF）和高阶巴特沃斯滤波器在时域中会产生振铃现象，而高斯滤波器由于其平滑的响应，可以避免这种问题。
- 参数 $D_0$ 的控制：
  - $D_0$ 越大，滤波器允许通过的频率范围越广，滤波效果越弱。
  - $D_0$ 越小，滤波器的截止频率越低，平滑效果越强。

10、典型的频域高通滤波器（锐化滤波器）

（1）高通滤波器的基本概念与推导

数学定义：高通滤波器是通过低通滤波器的补集构建而来，公式为： $H_{hp}(u,v) = 1 - H_{lp}(u,v)$ 其中 $H_{lp}(u,v)$ 是低通滤波器。
物理意义：高通滤波器用于保留图像中的高频分量，抑制低频成分（如均匀背景），实现图像锐化。

（2）理想高通滤波器（Ideal Highpass Filter, IHPF）：

数学公式： $H(u,v) = \begin{cases} 0 & \text{if } D(u,v) \leq D_0 \\ 1 & \text{if } D(u,v) > D_0 \end{cases}$
特性：简单直接，但引起严重的振铃效应（Gibbs现象），如图 4.24。

（3）巴特沃斯高通滤波器（Butterworth Highpass Filter, BHPF）：

数学公式： $H(u,v) = \frac{1}{1 + [D_0 / D(u,v)]^{2n}}$
特性：平滑过渡，减少振铃效应，阶数 $n$ 控制滤波器的锐度，如图 4.25。

（4）高斯高通滤波器（Gaussian Highpass Filter, GHPF）：

数学公式： $H(u,v) = 1 - e^{-D^2(u,v)/2D_0^2}$
特性：无振铃效应，平滑度最好，如图 4.26。

（5）三种高通滤波器比较

二维和一维的时域图：

二维和一维的频域图

（6）高频强调滤波器（锐化）

高频强调滤波的定义

高频强调滤波是一种结合低频分量与高频分量的滤波方法，其目的是在保留低频分量的同时，强调图像中的高频细节，公式如下：

$H_{hfe}(u,v) = a + b H_{hp}(u,v)$

$a$ 是低频分量的权重系数。
$b$ 是高频分量的权重系数。
$H_{hp}(u,v)$ 是高通滤波器（如巴特沃斯高通滤波器）。

图像处理步骤

读取原始图像：
```
f = imread('Fig0419(a)(chestXray_original).tif');
```
原始图像为胸部 X 光片，如图 4.19(a)。
图像填充：
```
PQ = paddedsize(size(f));
```
为了避免频域卷积中的“环绕误差（Wraparound Error）”，对图像进行零填充操作。
构建高通滤波器：
```
D0 = 0.05 * PQ(1);
HBW = hpfilter('btw', PQ(1), PQ(2), D0, 2);
```
- 使用巴特沃斯高通滤波器（'btw'）：
  - 截止频率 $D_0 = 0.05 \times$ 图像尺寸。
  - 滤波器阶数 $n = 2$ 。
高频强调滤波：
```
H = 0.5 + 2 * HBW;
gbf = dftfilt(f, H);
```
- 构建最终滤波器，结合低频与高频分量：
  - $a = 0.5$ ，保留低频分量。
  - $b = 2$ ，放大高频分量。
- dftfilt：对图像进行频域滤波，输出滤波后的结果 $gbf$ 。
后处理：
```
ghf = gscale(gbf);
ghe = histeq(ghf, 256);
imshow(ghe);
```
- gscale：将图像归一化到可视范围。
- histeq：执行直方图均衡，增强图像对比度。
- imshow：显示最终图像。