【图像处理lec7】图像恢复、去噪

一、图像退化与恢复概述

二、图像退化中的噪声模型

1、使用 imnoise 函数添加噪声

（1）imnoise 函数的概述

（2）函数语法

（3）支持的噪声类型与具体语法

（4）噪声类型的详细说明

（5）示例与实践指导

2、生成具有特定分布的随机噪声

（1）随机噪声生成的理论基础

（2）随机变量的类型及其特性

i、分布类型及其特性

ii、各分布的特点与应用场景

（3）各噪声MATLAB代码示例与直方图展示

3、周期噪声

（1）周期噪声信号及其傅里叶变换

（2）MATLAB 代码生成周期噪声

4、噪声参数估计

（1）独立噪声对图像的影响

（2）周期噪声对图像的影响

（3）噪声参数估计

（4）独立噪声与周期噪声的对比

（5）噪声参数估计的重要性

三、空间域滤波处理噪声

1、噪声模型与空间域滤波

2、均值滤波器

（1）几种均值滤波器

（2）滤波示例

3、顺序统计滤波器

（1）几种顺序统计滤波器

（2）顺序统计滤波效果

（3）自适应中值滤波器

i、原理与算法

ii、自适应中值滤波效果

iii、椒盐噪声处理总结

4、综合总结

四、频域滤波去除周期噪声

1、带阻滤波器（Bandreject Filters）

2. 带阻滤波器的应用

3. 带通滤波器（Bandpass Filters）

4. 陷波滤波器（Optimum Notch Filtering）

（1）基本概念

（2）数学表示

（3）频率域滤波步骤

（4）应用示例

（5）最优陷波滤波器（不懂）

（6）结论

五、退化函数估计

1、图像退化函数的估计

2、大气湍流模型

3、运动模糊建模

（1）假设图像的模糊是由于线性匀速运动造成的

（2）傅里叶变换的推导

（3）运动模糊的传递函数

（4）不同的运动模型解析

i、匀速直线运动沿 x 方向

ii、匀速直线运动沿任意方向

5、综合分析

六、直接逆滤波恢复图像

1. 直接反滤波的基本理论

2. 图像退化与噪声问题

七、维纳滤波（Wiener Filtering）

1、Wiener 滤波的理论基础

2、Wiener 滤波的优点与应用

3、实际效果对比

4、MATLAB 实现

5、总结

八、约束最小二乘滤波

1、约束最小二乘滤波的理论基础

2、约束最小二乘滤波的数学推导

3、噪声水平的估计与约束参数

4、约束最小二乘滤波的实验效果

附录

1、随机噪声生成的理论基础

（1）理论背景与核心内容

（2）原理推导

（3）为什么逆变换法有效？

（4）定量证明：服从的分布

1. 问题描述

2. 基本思路

3. 证明过程

（5）定性分析

一、图像退化与恢复概述

退化与恢复的概念：

退化：通过退化函数H 和噪声 $\eta$ 模拟图像退化。
恢复：通过恢复滤波器，恢复被退化的图像。

图像退化公式：

$g(x,y) = H[f(x,y)] + \eta(x,y)$
$g(x,y) = h(x,y) * f(x,y) + \eta(x,y)$
$G(u,v) = H(u,v)F(u,v) + N(u,v)$

图像恢复：

改善图像，使其在某种预定义的标准下更好。
根据退化模式，选择不同的恢复滤波器

二、图像退化中的噪声模型

1、使用 `imnoise` 函数添加噪声

（1）`imnoise` 函数的概述

imnoise 是 MATLAB 中用于向图像添加噪声的函数，具有强大的功能和广泛的应用。它通过模拟现实环境中的噪声干扰（如传感器误差、传输噪声等），为图像处理算法的开发和测试提供数据支持。

（2）函数语法

$g = \text{imnoise}(f, \text{type}, \text{parameters})$

f：输入图像（可以是灰度图像或RGB图像）。
type：噪声类型，指定添加哪种噪声（例如高斯噪声、椒盐噪声）。
parameters：根据噪声类型不同，需要传入不同的参数来控制噪声的特性。
g：输出噪声图像。

（3）支持的噪声类型与具体语法

噪声类型	函数语法	参数含义	解释与特点
高斯噪声	`G = imnoise(f, 'gaussian', m, var)`	m：均值，var：方差。	模拟随机高斯噪声，广泛存在于传感器图像中。
局部方差噪声	`G = imnoise(f, 'localvar', V)`	V：局部方差分布矩阵。	适用于局部噪声强度变化的场景，噪声非均匀分布。
椒盐噪声	`G = imnoise(f, 'salt & pepper', d)`	d：噪声密度，范围 [0, 1]。	随机替换部分像素为最小值（黑色）或最大值（白色）。
斑点噪声	`G = imnoise(f, 'speckle', var)`	var：噪声方差。	乘性噪声，图像噪声与强度值成比例，常见于超声成像。
泊松噪声	`G = imnoise(f, 'poisson')`	无需参数，由图像像素值自动决定。	基于泊松分布，常见于光学成像（光子计数噪声）。

（4）噪声类型的详细说明

高斯噪声 (Gaussian Noise)
- 数学模型：高斯噪声服从正态分布： $P(x) = \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{(x - \mu)^2}{2\sigma^2}}$ 其中 μ 为均值， $\sigma^2$ 为方差。
- 特点：
  - 噪声均匀分布在整个图像中，适合模拟传感器噪声。
  - 图像整体出现细微随机抖动。
- 应用：
  - 模拟真实图像采集中的噪声干扰。
  - 测试去噪算法的效果（如均值滤波、维纳滤波）。
局部方差噪声 (Local Variance Noise)
- 特点：
  - 噪声强度在局部区域内变化，与图像的纹理和局部强度分布相关。
- 优势：
  - 适合模拟复杂噪声环境，如自然图像或光学系统中非均匀噪声。
椒盐噪声 (Salt & Pepper Noise)
- 数学描述：
  - 噪声以概率 d 随机替换像素为0（黑色）或1（白色）。
- 视觉效果：图像中出现随机的黑点和白点。
- 应用：
  - 模拟传输错误、像素损坏导致的噪声。
  - 测试中值滤波去噪算法的效果。
斑点噪声 (Speckle Noise)
- 特点：
  - 噪声与图像强度成正比，表现为斑点状纹理。
  - 乘性噪声模型： $g(x,y) = f(x,y) + f(x,y) \cdot \eta(x,y)$
  - 其中 $\eta(x,y)$ 是均值为0的高斯噪声。
- 应用：
  - 模拟超声图像、SAR雷达图像中的噪声干扰。
泊松噪声 (Poisson Noise)
- 数学背景：基于泊松分布： $P(k;\lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \; k = 0, 1, 2, \dots$ 其中 $\lambda$ 与图像的像素值成比例。
- 特点：
  - 噪声强度随像素值增加而增加。
- 应用：
  - 模拟光学成像中的光子噪声，例如天文图像。

（5）示例与实践指导

通过示例代码向图像添加噪声：

% 读取输入图像
f = imread('cameraman.tif');% 添加不同类型的噪声
G1 = imnoise(f, 'gaussian', 0, 0.01);       % 高斯噪声
G2 = imnoise(f, 'salt & pepper', 0.02);    % 椒盐噪声
G3 = imnoise(f, 'speckle', 0.05);          % 斑点噪声
G4 = imnoise(f, 'poisson');                % 泊松噪声% 显示噪声效果
figure;
subplot(2,2,1); imshow(G1); title('Gaussian Noise');
subplot(2,2,2); imshow(G2); title('Salt & Pepper Noise');
subplot(2,2,3); imshow(G3); title('Speckle Noise');
subplot(2,2,4); imshow(G4); title('Poisson Noise');

输出：

原图和带有不同噪声的图像直观对比，清晰展示各类噪声的特征。

2、生成具有特定分布的随机噪声

（1）随机噪声生成的理论基础

详细分析见附录1

随机噪声的生成：
如果 w 是一个区间 (0,1) 内均匀分布的随机变量，可以通过累积分布函数 (CDF) $F_Z$ 的逆函数 $F_Z^{-1}$ 来生成服从指定分布的随机变量 z。
公式如下：
$z = F_Z^{-1}(w)$
Rayleigh分布的示例：
累积分布函数 (CDF) 定义为：
$F_Z(z) = \begin{cases} 1 - e^{-\frac{(z - a)^2}{b}} & z \geq a \\ 0 & z < a \end{cases}$
由此，解出 z：
$z = a + \sqrt{-b \ln(1 - w)}$

分析：

理论支撑：
- 提供了从均匀分布生成指定分布随机数的通用方法，通过CDF逆变换实现。
- 公式清晰，Rayleigh分布作为示例有代表性，说明了如何从均匀分布生成Rayleigh噪声。
Rayleigh分布的具体实现：
- a：平移参数。
- b：控制噪声的扩展程度。

（2）随机变量的类型及其特性

i、分布类型及其特性

名称 (Name)	PDF (概率密度函数)	均值与方差 (Mean & Variance)	CDF (累积分布函数)	生成方法
Uniform (均匀分布)	$p_Z(z) = \frac{1}{b-a} \quad a \leq z \leq b$	$m = \frac{a+b}{2},$ $\, \sigma^2 = \frac{(b-a)^2}{12}$	$F_Z(z) = \frac{z-a}{b-a}$	MATLAB 函数：`rand`
Gaussian (高斯分布)	$p_Z(z) = \frac{1}{\sqrt{2\pi b}} e^{-\frac{(z-a)^2}{2b}}$	$m = a,$ $\, \sigma^2 = b$	$F_Z(z) = \int_{-\infty}^z p_Z(x)dx$	MATLAB 函数：`randn`
Salt & Pepper (椒盐噪声)	$p_Z(z) = \begin{cases} P_a & z = a \\ P_b & z = b \\ 0 & \text{otherwise} \end{cases}$	$m = aP_a + bP_b,$ $\sigma^2 = (a-m)^2P_a + (b-m)^2P_b$	$F_Z(z) =$ 分段定义	MATLAB 函数：`rand` + 逻辑条件
Lognormal (对数正态分布)	$p_Z(z) = \frac{1}{z\sqrt{2\pi b}} e^{-\frac{(\ln z - a)^2}{2b}}$	$m = e^{a + b/2},$ $\, \sigma^2 = e^{2a+b}(e^b - 1)$	$F_Z(z) = \int_{0}^z p_Z(x)dx$	通过指数变换 $z = e^{a + bN(0,1)}$
Rayleigh (瑞利分布)	$p_Z(z) = \frac{z-a}{b} e^{-\frac{(z-a)^2}{2b}} \, (z \geq a)$	$m = a + \sqrt{\frac{\pi b}{2}},$ $\, \sigma^2 = b(4-\pi)/4$	$F_Z(z) = 1 - e^{-\frac{(z-a)^2}{b}}$	$z = a + \sqrt{-b\ln(1-U)}$
Exponential (指数分布)	$p_Z(z) = \lambda e^{-\lambda z} \, (z \geq 0)$	$m = \frac{1}{\lambda},$ $\, \sigma^2 = \frac{1}{\lambda^2}$	$F_Z(z) = 1 - e^{-\lambda z}$	$z = -\frac{1}{\lambda} \ln(1-U)$
Erlang (厄兰分布)	$p_Z(z) = \frac{\lambda^k z^{k-1} e^{-\lambda z}}{(k-1)!}$	$m = \frac{k}{\lambda},$ $\, \sigma^2 = \frac{k}{\lambda^2}$	$F_Z(z) = 1 - \sum_{n=0}^{k-1} \frac{(\lambda z)^n}{n!} e^{-\lambda z}$	通过累加 k 个指数分布变量生成

ii、各分布的特点与应用场景

提供了6种概率分布的PDF图示：

高斯分布 (Gaussian)
瑞利分布 (Rayleigh)
伽马分布 (Gamma)
指数分布 (Exponential)
均匀分布 (Uniform)
椒盐分布 (Impulse)

均匀分布 (Uniform Distribution)
- 特点：随机变量在区间 [a,b] 内等概率分布。
- 应用：
  - 基础随机数生成。
  - 用作生成其他分布（如CDF逆变换法中的输入）。
- 生成：rand 在MATLAB中可直接生成。
高斯分布 (Gaussian Distribution)
- 特点：中心对称的钟形曲线，由均值 a 和方差 b 控制形态。
- 应用：
  - 噪声建模：传感器噪声、自然界随机噪声。
  - 数据分析与统计建模。
- 生成：randn 在MATLAB中可直接生成标准高斯分布。
椒盐噪声 (Salt & Pepper Noise)
- 特点：图像中随机出现“黑点 (0)” 和 “白点 (255)” 噪声。
- 应用：
  - 图像处理噪声模拟。
  - 用于测试中值滤波等去噪算法。
- 生成：结合 rand 和逻辑条件实现。
对数正态分布 (Lognormal Distribution)
- 特点：随机变量的对数服从正态分布。
- 应用：
  - 生物统计：描述生长现象、经济数据等。
  - 信号幅度建模。
- 生成：通过 $z = e^{a + bN(0,1)}$ 。
瑞利分布 (Rayleigh Distribution)
- 特点：右偏分布，常用于幅度建模。
- 应用：
  - 无线通信：信号衰减建模。
  - 噪声模拟：超声图像中的噪声。
- 生成：通过 $z = a + \sqrt{-b\ln(1-U)}$ 。
指数分布 (Exponential Distribution)
- 特点：单调递减分布，用于建模时间间隔。
- 应用：
  - 等待时间建模：系统故障时间、排队问题。
  - 噪声模拟。
- 生成：通过 $z = -\frac{1}{\lambda} \ln(1-U)$ 。
厄兰分布 (Erlang Distribution)
- 特点：由 k 个独立同分布的指数随机变量累加而成。
- 应用：
  - 排队论与服务时间建模。
  - 通信系统中的随机过程建模。

（3）各噪声MATLAB代码示例与直方图展示

图中展示了6种随机噪声的直方图：
a. 高斯噪声
b. 均匀噪声
c. 对数正态噪声
d. 瑞利噪声
e. 指数噪声
f. 厄兰噪声
每个直方图清晰呈现了不同噪声的分布特征。

3、周期噪声

（1）周期噪声信号及其傅里叶变换

周期噪声信号

其中：
- 定义公式：
- $A$ ：振幅，控制噪声的强度。
- $u_0, v_0$ ：频率参数，决定噪声在 $x$ 和 $y$ 方向的周期性。
- $B_x, B_y$ ：相位偏移。
- $M, N$ ：图像的大小。
周期噪声的傅里叶变换 (FT)
- 傅里叶变换结果 $R(u,v)$ ：
- 这里的 $\delta$ 表示傅里叶域中的脉冲，出现在 $(u + u_0, v + v_0)$ 和 $(u - u_0, v - v_0)$ 两个对称位置。

分析与解释

周期噪声的特点：
- 周期噪声在空间域表现为正弦波（周期性），在傅里叶域表现为两个对称的脉冲。
- 这些脉冲的位置和强度由噪声的频率 $u_0, v_0$ 和振幅 $A$ 决定。
傅里叶域表示的物理意义：
- 在 $(u + u_0, v + v_0)$ 和 $(u - u_0, v - v_0)$ 的脉冲对应于正弦噪声的频率分量。
- 相位偏移 $B_x$ 和 $B_y$ 通过指数项引入，影响噪声在时域的位置。

（2）MATLAB 代码生成周期噪声

代码示例
- 定义周期噪声的参数矩阵 C：
  - 这里的 C 定义了不同的噪声频率分量。
- 使用 imnoise3 函数生成噪声：
```
C=[0 64;0 128;32 32;64 0;128 0;-32 32];
[r,R,S]=imnoise3(512,512,C);
imshow(S,[]);
figure,imshow(r,[]);
```
结果展示
- 左侧图像：傅里叶域的噪声谱，表现为若干个点状脉冲。
- 右侧图像：空间域的周期噪声图案，呈现为规则的网格状或条纹状结构。

分析与解释

噪声生成的原理：
- 通过指定傅里叶域的频率位置 $(u_0, v_0)$ ，再进行逆变换，得到空间域的周期噪声。
- C 矩阵的行定义了不同的噪声频率分量，因此生成的噪声具有多重周期性。
MATLAB imnoise3 函数：
- 该函数通过傅里叶域的脉冲生成对应的周期噪声，是噪声仿真中的一种重要工具。
视觉结果：
- 空间域中的网格状图案是正弦波的叠加。
- 傅里叶域的脉冲清晰展示了噪声的频率分布。

4、噪声参数估计

（1）独立噪声对图像的影响

这些图展示了6种独立噪声对图像的污染效果，包括：
1. Gaussian 噪声（高斯噪声）
2. Rayleigh 噪声（瑞利噪声）
3. Gamma 噪声
4. Exponential 噪声（指数噪声）
5. Uniform 噪声（均匀噪声）
6. Salt & Pepper 噪声（椒盐噪声）
每种噪声图像下方显示了对应的直方图，反映了图像的灰度级分布。

分析与解释

高斯噪声 (Gaussian Noise)：
- 图像：随机噪声均匀分布于整个图像，呈现出“雪花”状噪声。
- 直方图：呈现典型的钟形曲线，符合高斯分布的特性。
- 应用：模拟传感器的随机噪声或自然界噪声。
瑞利噪声 (Rayleigh Noise)：
- 图像：噪声的随机性较明显，强度分布不均匀。
- 直方图：显示出右偏分布，符合Rayleigh分布的特征。
- 应用：常用于建模幅度噪声，如超声图像中的噪声。
Gamma噪声 (Gamma Noise)：
- 图像：噪声分布不均匀，有一定的随机性。
- 直方图：显示出偏态分布，具有较强的尾部特性。
- 应用：用于建模信号噪声或光照分布不均匀的情况。
指数噪声 (Exponential Noise)：
- 图像：噪声强度逐渐衰减。
- 直方图：显示出单调递减的指数形态。
- 应用：描述时间间隔、等待时间的噪声。
均匀噪声 (Uniform Noise)：
- 图像：噪声随机分布，均匀覆盖整个图像。
- 直方图：呈现平坦分布，表明灰度级分布一致。
- 应用：基础随机噪声，用于测试滤波器。
椒盐噪声 (Salt & Pepper Noise)：
- 图像：噪声表现为随机的黑点和白点。
- 直方图：直方图上出现极端峰值，分别对应最小值 (0) 和最大值 (255)。
- 应用：模拟传输错误或传感器故障导致的噪声。

（2）周期噪声对图像的影响

内容总结

图像被正弦噪声污染，表现为规律的条纹状干扰。
傅里叶频谱展示了周期噪声的频率成分，形成一对共轭脉冲。

分析与解释

空间域表现：
- 周期噪声在空间域中表现为规则的网格或条纹，这是由于正弦波的周期性。
傅里叶域特征：
- 在傅里叶域中，周期噪声对应于频谱图中的共轭脉冲，位置由噪声的频率确定。
实际应用：
- 周期噪声常见于扫描图像或电磁干扰，识别并去除周期噪声对于图像质量至关重要。

（3）噪声参数估计

如何估计噪声参数，主要包括：

傅里叶域方法：
- 检测周期噪声的频率分量，通过观察傅里叶谱的脉冲位置来确定噪声参数。
统计参数估计：
- 通过图像的灰度级直方图估计噪声的均值和方差。
实现步骤：
- 选择感兴趣区域 (ROI) 进行噪声分析（roipoly）。
- 使用 MATLAB 函数计算中心矩与均值。
公式与方法：
- 中心矩定义：
  - $n=0$ ：归一化常数。
  - $n=1$ ：均值。
  - $n=2$ ：方差。

图像与直方图对应关系

图像：展示噪声对图像的影响。
直方图：反映图像灰度级分布，帮助识别噪声类型。
- 高斯噪声 → 钟形分布。
- 均匀噪声 → 平坦分布。
- 指数噪声 → 右偏分布。
- 椒盐噪声 → 极端峰值。

（4）独立噪声与周期噪声的对比

类别	空间域表现	傅里叶域表现	常见去噪方法
独立噪声	随机分布的点状噪声，无明显规律。	频谱中无明显脉冲，呈现散射分布。	空间域滤波：均值滤波、中值滤波、维纳滤波等。
周期噪声	规则的条纹或网格状噪声。	频谱中出现对称的脉冲。	频域滤波：带阻滤波器、陷波滤波器等。

（5）噪声参数估计的重要性

准确建模：
- 通过估计噪声的均值、方差、分布形态等参数，可以准确建模噪声，为去噪提供理论基础。
去噪算法设计：
- 不同噪声对应不同的滤波器：
  - 高斯噪声 → 高斯滤波、维纳滤波。
  - 椒盐噪声 → 中值滤波。
  - 周期噪声 → 频域滤波。
实现步骤：
- 选择图像区域 (ROI) 进行分析。
- 计算灰度直方图并估计均值与方差。
- 比较直方图形状，确定噪声分布类型。

三、空间域滤波处理噪声

1、噪声模型与空间域滤波

内容总结

噪声模型：
- 给定图像模型：其中：
  - $f(x,y)$ ：原始图像。
  - $\eta(x,y)$ ：噪声（通常为随机噪声）。
- 在频域表示： $G(u,v) = F(u,v) + N(u,v)$
噪声滤波器：
- 空间噪声滤波器：直接在空间域处理噪声，例如均值滤波、中值滤波等。
- 自适应空间滤波器：根据图像局部统计特性自适应地滤除噪声。

2、均值滤波器

（1）几种均值滤波器

算术均值滤波器 (Arithmetic Mean Filter)：
- 基本的均值滤波，适用于去除高斯噪声。
几何均值滤波器 (Geometric Mean Filter)：
- 对乘积计算，适用于细微噪声。
调和均值滤波器 (Harmonic Mean Filter)：
- 对盐噪声效果好，但对椒噪声无效。
逆调和均值滤波器 (Contraharmonic Mean Filter)：
- $Q > 0$ ：去除椒噪声； $Q < 0$ ：去除盐噪声。

应用与局限性：

算术均值和几何均值适用于随机噪声（如高斯噪声）。
调和均值和逆调和均值适用于椒盐噪声，但无法同时处理两种噪声。

（2）滤波示例

示例1

a 图是原图
b 图是加入高斯噪声的图
c 图是3×3核尺寸的算数平均滤波
d 图是3×3核尺寸的几何平均滤波

结果分析：算术均值 和 几何均值 滤波后噪声显著降低，但图像细节模糊。

示例2

a 图是加入椒噪声
b 图是加入盐噪声
c 图是对 a 图进行逆调和均值滤波，核尺寸3×3，参数 Q=1.5
d 图是对 b 图进行逆调和均值滤波，核尺寸3×3，参数 Q=-1.5

结果分析：逆调和均值滤波 可以有效去除椒盐噪声，但要根据噪声类型调整 Q 值。如果Q值选择错误的话，滤波后效果会恶化，如下图：

第一个图是对上面的加入椒噪声的 a 图进行 Q=-1.5 的滤波
第二个图是对上面的加入盐噪声的 b 图进行 Q=1.5 的滤波

3、顺序统计滤波器

（1）几种顺序统计滤波器

中值滤波器 (Median Filter)：
- 定义： $\hat{f}(x,y) = \text{median}_{(s,t) \in S_{xy}} g(s,t)$
- 特点：对椒盐噪声效果最好，保留边缘细节。
最大和最小滤波器 (Max & Min Filters)：
- 最大滤波器：
  - 去除椒噪声。
- 最小滤波器：
  - 去除盐噪声。
中点滤波器 (Midpoint Filter)：
- 计算最大和最小值的平均： $\hat{f}(x,y) = \frac{1}{2} \left[ \max_{(s,t) \in S_{xy}} g(s,t) + \min_{(s,t) \in S_{xy}} g(s,t) \right]$
Alpha裁剪均值滤波器 (Alpha-Trimmed Mean Filter)：
- 去除部分极端值： $\hat{f}(x,y) = \frac{1}{mn - d} \sum_{(s,t) \in S_{xy}} g_r(s,t)$ 其中 d 表示被裁剪的元素数量。

（2）顺序统计滤波效果

示例1

中值滤波：
- (b)、(c)、(d) 展示了中值滤波连续多次应用的效果，随着迭代次数增加，椒盐噪声逐渐去除。

示例2

最大与最小滤波：

a图，最大滤波器去除了椒噪声，但会引入图像膨胀效果。
b图，最小滤波器去除了盐噪声，但会引入图像收缩效果。

（3）自适应中值滤波器

i、原理与算法

自适应中值滤波器的核心思想是：根据噪声的局部统计特性，自适应地调整滤波窗口大小，进一步提高去噪效果。

符号定义：
- $Z_{\min}$ ：邻域内的最小灰度值。
- $Z_{\max}$ ：邻域内的最大灰度值。
- $Z_{\text{med}}$ ：邻域内的中值。
- $Z_{xy}$ ：当前像素点的灰度值。
算法步骤：
- Level A：
  判断：
  - 如果满足，继续 Level B。
  - 否则增加窗口大小，直到最大窗口 $S_{\max}$ 。
- Level B：
  判断：
  - 如果满足，输出 $Z_{xy}$ 。
  - 否则输出 $Z_{\text{med}}$ 。

ii、自适应中值滤波效果

(a) 输入图像被椒盐噪声污染。
(b) 使用标准中值滤波器进行处理，部分噪声被去除。
(c) 使用自适应中值滤波器（窗口最大尺寸 $S_{\text{max}} = 7$ ），噪声被有效去除，图像细节得到较好保留。

优势：
- 能够处理椒盐噪声密度较高的图像。
- 自适应调整窗口大小，避免过度平滑图像细节。
局限性：
- 对于噪声密度极高的图像，可能需要更大窗口，导致计算量增加。

iii、椒盐噪声处理总结

低密度噪声：
- 使用中值滤波器即可达到良好的效果。
高密度噪声：
- 使用自适应中值滤波器，能够有效去噪并保持图像细节。
单一噪声：
- 使用逆调和均值滤波器，根据噪声类型设定 Q 值（Q>0：椒噪声；Q<0：盐噪声）。

4、综合总结

噪声与滤波器选择对比

噪声类型	最佳滤波器	优势	局限性
高斯噪声	算术均值、几何均值、调和均值滤波器	有效降低随机噪声	图像细节模糊
盐噪声	中值滤波、逆调和均值滤波 (Q < 0)	有效去除白色点噪声	逆调和均值需选择合适的 Q
椒噪声	中值滤波、逆调和均值滤波 (Q > 0)	有效去除黑色点噪声	逆调和均值无法同时处理盐椒噪声
椒盐噪声	中值滤波	同时处理盐噪声和椒噪声	多次迭代可能导致细节丢失

总结滤波器的使用策略

随机噪声（如高斯噪声）：使用均值滤波器。
椒盐噪声：使用中值滤波器，或逆调和均值滤波器（需要根据 QQQ 值调整）。
噪声参数估计：通过直方图或局部统计量确定噪声类型。

四、频域滤波去除周期噪声

1、带阻滤波器（Bandreject Filters）

带阻滤波器的目标是去除特定频带的噪声信号，同时尽量保留其他有效信号。
三种带阻滤波器：
1. 理想带阻滤波器：
  - 特征：具有严格的频带截止，但会产生振铃效应。
2. 巴特沃斯带阻滤波器（Butterworth）：
  - 特征：频带过渡较平滑，减少振铃效应。
3. 高斯带阻滤波器（Gaussian）：
  - 特征：具有平滑的频率响应，没有振铃效应。
图示（Figure 5.15）：
- 左侧：理想带阻滤波器（锐利的截断环）。
- 中间：巴特沃斯带阻滤波器（平滑的过渡）。
- 右侧：高斯带阻滤波器（平滑且柔和的响应）。

2. 带阻滤波器的应用

应用示例

噪声图像：图像被周期性正弦噪声污染（Figure 5.16(a)）。
频谱分析：在频域中观察到明显的周期性噪声频率尖峰（Figure 5.16(b)）。
带阻滤波器应用：
- 使用巴特沃斯带阻滤波器（Figure 5.16(c)）去除噪声频率成分。
- 滤波结果（Figure 5.16(d)）：噪声被去除，图像得到有效恢复。

3. 带通滤波器（Bandpass Filters）

带通滤波器是带阻滤波器的补集，其传递函数为： $H_{\text{bp}}(u,v) = 1 - H_{\text{br}}(u,v)$
应用：
- 滤出（Figure 5.16（a））的噪声成分以便分析（Figure 5.17）。
- 滤波后观察得到的噪声图案验证了噪声的频率分布。

4. 陷波滤波器（Optimum Notch Filtering）

（1）基本概念

陷波滤波器（Notch Filter）是一种特殊类型的滤波器，用于抑制特定频率或频率范围内的信号，而不影响其他频率成分。它主要用于去除图像中的周期性噪声，例如网格状或条纹噪声。

理想、巴特沃斯、高斯陷波滤波器（Figure 5.18）。

（2）数学表示

陷波滤波器的频率域函数 $H(u,v)$ 表达如下：

理想陷波滤波器（Ideal Notch Filter）
- 其中 $D_1(u,v)$ 和 $D_2(u,v)$ 是到噪声频率点的距离。
- $D_0$ 是一个阈值，定义了“陷波”的半径。
巴特沃斯陷波滤波器（Butterworth Notch Filter）
- n 是滤波器的阶数，决定了滤波器过渡的平滑程度。
- 阶数越高，滤波器越接近理想滤波器。
高斯陷波滤波器（Gaussian Notch Filter）
- 高斯滤波器具有无穷平滑的过渡特性。
距离函数
两个距离函数 $D_1(u,v)$ 和 $D_2(u,v)$ 分别表示到噪声频率点 $(u_0, v_0)$ 和 $(-u_0, -v_0)$ 的距离：
$D_1(u,v) = \left[ \left( u - \frac{M}{2} - u_0 \right)^2 + \left( v - \frac{N}{2} - v_0 \right)^2 \right]^{0.5}$ $D_2(u,v) = \left[ \left( u - \frac{M}{2} + u_0 \right)^2 + \left( v - \frac{N}{2} + v_0 \right)^2 \right]^{0.5}$
其中 M 和 N 是图像的尺寸， $(u_0, v_0)$ 是噪声频率位置。

（3）频率域滤波步骤

傅里叶变换
将图像 $g(x,y)$ 转换到频率域 $G(u,v)$ 。
设计陷波滤波器 $H(u,v)$
根据噪声频率位置和滤波器类型（理想、巴特沃斯或高斯）设计适当的陷波滤波器。
频率域滤波
将频率域图像与滤波器相乘：
$G_f(u,v) = H(u,v)G(u,v)$
反傅里叶变换
对滤波后的频率域图像进行逆傅里叶变换，得到处理后的图像 $g_f(x,y)$ 。

（4）应用示例

示例1：

图像被周期性噪声（例如传感器扫描线）污染（Figure 5.19(a)）。
频谱显示噪声对应的频率成分（Figure 5.19(b)）。
陷波滤波器应用去除噪声成分（Figure 5.19(c, d)）。
最终结果（Figure 5.19(e)）：去除了周期性噪声，图像得到恢复。

（5）最优陷波滤波器（不懂）

针对复杂噪声，最优陷波滤波器通过最小化局部方差的方式达到最佳效果：

目标函数
$\sigma^2 = \frac{1}{(2a+1)(2b+1)} \sum_{s=-a}^{a} \sum_{t=-b}^{b} \left[ g(x+s,y+t) - w(x+s,y+t)\eta(x+s,y+t) \right]^2$
求导并最小化
通过对 $w(x,y)$ 求导并设置为零，得到最优权重函数 $w(x,y)$ ：
$w(x,y) = \frac{\overline{g(x,y)\eta(x,y)} - \overline{g(x,y)}\overline{\eta(x,y)}}{\overline{\eta^2(x,y)} - \overline{\eta(x,y)}^2}$
该权重可以自适应调整，进一步提高噪声去除效果。

（6）结论

陷波滤波器通过选择性地抑制特定噪声频率，能有效去除周期性噪声。

理想陷波滤波器 适用于简单噪声，但易引入伪影。
巴特沃斯和高斯陷波滤波器 平滑过渡，适用于更复杂的噪声环境。
最优陷波滤波器 通过最小化局部方差，自适应地调整权重，达到最佳去噪效果。

图中的示例清楚地展示了陷波滤波器在去除周期性噪声中的高效性，尤其是在频域分析和滤波后的图像恢复方面。

五、退化函数估计

1、图像退化函数的估计

概述

退化模型：图像退化是图像质量下降的过程，通常由大气湍流、运动模糊等原因引起，这里的图像退化并不是噪声引起的。
观察到的图像： $G_s(u,v)$ 表示已观测到的图像在频域中的子图像。
退化函数估计：
- 从观测图像 $G_s(u,v)$ 和退化模型 $\hat{F}_s(u,v)$ 的比值推导出退化函数 $H_s(u,v)$ 。
- 如果已知观测图像与噪声较小，可用 $H_s(u,v) = \frac{G_s(u,v)}{\hat{F}_s(u,v)}$ 估计。

建模与实验方法

实验方法： $H(u,v) = \frac{G(u,v)}{A}$ ，其中 A 是观测数据的幅度。
基于物理模型的估计：
- 以 Hufnagel 提出的湍流模型为例，退化函数表示为： $H(u,v) = \exp \left( -k (u^2 + v^2)^{5/6} \right)$
其中 k 表示大气湍流程度。

2、大气湍流模型

图示分析：
- 图5.24：展示了一个光脉冲通过湍流后发生的扩散现象。
  - (a) 原始的光脉冲，几乎无扩散。
  - (b) 通过大气湍流后，脉冲显著模糊。
- 图5.25：展示了不同湍流条件下的地面图像：
  - (a) 无湍流，图像清晰。
  - (b) 严重湍流（k=0.0025），图像模糊。
  - (c) 中等湍流（k=0.001）。
  - (d) 轻微湍流（k=0.00025）。
分析：湍流程度越大，图像的高频成分被更多地衰减，导致细节丢失。

3、运动模糊建模

（1）假设图像的模糊是由于线性匀速运动造成的

原始图像在时间内发生了位移，定义模糊后的图像为：
- $x_0(t)$ 和 $y_0(t)$ 分别表示图像在 x 和 y 方向上的位移。
- T 是运动的持续时间。
- $f(x - x_0(t), y - y_0(t))$ 是在 $t$ 时刻，图像 $f$ 相对于 $(x, y)$ 的位移。

（2）傅里叶变换的推导

傅里叶变换将空间域表示的图像转换到频域，得到 $G(u, v)$ ：

$G(u, v) = \iint_{-\infty}^\infty g(x, y) \exp(-j 2\pi (ux + vy)) \, dx \, dy$

由时域平移性质：

$G(u, v) = F(u, v) \int_0^T \exp(-j 2\pi (u x_0(t) + v y_0(t))) \, dt$

（3）运动模糊的传递函数 $H(u, v)$

由上面推导可以看出，模糊的频域表示为：

$G(u, v) = F(u, v) \cdot H(u, v)$

其中：

$H(u, v) = \int_0^T \exp(-j 2\pi (u x_0(t) + v y_0(t))) \, dt$

（4）不同的运动模型解析

i、匀速直线运动沿 x 方向

当 $x_0(t) = \frac{a t}{T}$ 且 $y_0(t) = 0$ 时：

$H(u, v) = \int_0^T \exp\left( -j 2\pi u \frac{a t}{T} \right) \, dt$

积分结果为：

$H(u, v) = \frac{T}{\pi u a} \sin(\pi u a) \exp(-j \pi u a)$

这个结果表明，模糊函数 $H(u, v)$ 在频域上是一个正弦函数包络，造成了频域的衰减和周期性零点。

ii、匀速直线运动沿任意方向

当 $x_0(t) = \frac{a t}{T}, y_0(t) = \frac{b t}{T}$ 时：

$H(u, v) = \frac{T}{\pi (u a + v b)} \sin(\pi (u a + v b)) \exp(-j \pi (u a + v b))$

这里的 a 和 b 分别表示在 x 和 y 方向上的速度分量。
结果同样是一个正弦函数包络，但方向由 a 和 b 控制，模糊沿着运动方向扩展。

物理意义：

模糊现象：运动模糊会导致图像的高频成分被平均，形成方向性的模糊。
运动长度与角度：退化函数的幅度与 a,b 决定运动的长度和方向，进而影响模糊的程度和方向。

5、综合分析

退化函数估计：
- 通过观测图像或物理建模可估算退化函数 H(u,v)，进一步应用于图像复原。
大气湍流建模：
- 湍流的退化函数会衰减高频分量，导致图像模糊，退化程度随湍流参数 k 增大而加剧。
运动模糊建模：
- 线性运动导致图像高频成分沿运动方向平均化，形成方向性模糊。
MATLAB 实验：
- 利用 MATLAB 模拟运动模糊及噪声，可以验证理论模型的有效性。

六、直接逆滤波恢复图像

1. 直接反滤波的基本理论

原理介绍：

直接反滤波是对退化图像的一种简单恢复方法，它基于退化函数 $H(u,v)$ ，这里并不是去噪，是恢复其他类型的退化，如模糊、运动等。
在频域中，对于退化图像 $G(u,v)$ ，原图像 $F(u,v)$ 可以通过如下公式进行估计：

$\hat{F}(u,v) = \frac{G(u,v)}{H(u,v)}$

如果图像受到退化的影响（如模糊、运动等），直接反滤波尝试通过对退化函数的逆操作恢复原始图像。

2. 图像退化与噪声问题

实际情况下，由于噪声的存在，恢复公式会引入误差，特别是在 $H(u,v)$ 的某些值接近零或很小的情况下，反滤波会放大噪声，导致恢复图像质量下降。

噪声影响公式：

$\hat{F}(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)}$

其中 $N(u,v)$ 是噪声的频域表示。如果 $H(u,v)$ 很小，那么 $N(u,v)/H(u,v)$ 的贡献会显著放大。

解决方法：

限制滤波频率：如图中所示，通过在频域中限制滤波器的应用范围，避免高频噪声被放大。
滤波器半径限制：对滤波器 H 施加截止半径（如 40、70、85），可以改善恢复效果。

结果分析：

图示描述：

(a) 直接反滤波的结果：未对滤波器进行频率截止，图像中噪声被显著放大，图像模糊且噪声较大。
(b) 截止频率 40 的结果：一定程度改善，部分噪声被抑制，但图像仍较模糊。
(c) 截止频率 70 的结果：图像质量进一步改善，细节开始显现。
(d) 截止频率 85 的结果：恢复效果较好，但边缘和细节可能仍有轻微模糊。

七、维纳滤波（Wiener Filtering）

1、Wiener 滤波的理论基础

目标：通过 最小化统计误差函数 $e^2 = E\{(f - \hat{f})^2\}$ 来求得最优滤波器，恢复被退化的图像 f^\hat{f}f^。
统计误差函数：
$e^2 = E\{(f - \hat{f})^2\}$
其中 $E$ 为期望值操作， $f$ 是理想图像， $\hat{f}$ 是恢复的图像。
频域表示：维纳滤波器的解可以在频域中表示为：

其中：
- $H(u,v)$ ：退化函数
- $S_\eta(u,v)$ ：噪声功率谱
- $S_f(u,v)$ ：原始图像的功率谱
- $G(u,v)$ ：观测到的退化图像的频域表示
- $H^*(u,v)$ ： $H(u,v)$ 的复共轭

2、Wiener 滤波的优点与应用

去噪能力强：通过平衡噪声功率谱和图像功率谱，维纳滤波器能有效减少噪声影响。
优化性能：通过最小化误差函数，使恢复图像接近真实图像。
适应退化：不仅适用于退化图像，还能处理加性噪声。

3、实际效果对比

示例1：

恢复效果对比：

(a)：使用 全逆滤波 恢复图像，导致噪声被放大，效果差。
(b)：限制 H 在某个半径范围外的使用，减少了高频噪声放大，但恢复效果仍不理想。
(c)：维纳滤波 恢复效果明显优于逆滤波，噪声被有效抑制，图像质量更高。

示例2：

不同条件下的恢复效果：

第一行：图像受运动模糊和噪声影响严重。
第二行：逆滤波虽然恢复了一部分细节，但放大了噪声。
第三行：维纳滤波结果更为平滑，噪声被有效减少，图像可读性提高。

4、MATLAB 实现

MATLAB 命令：

deconvwnr：实现维纳滤波。
fspecial('motion', len, theta)：生成运动模糊的点扩散函数 (PSF)。
fft2 和 ifft2：对图像进行傅里叶变换和逆变换。

代码示例：

PSF = fspecial('motion', 7, 45); % 运动模糊 PSF
g = imread('Fig0507(d).tif'); % 读取退化图像
fr1 = deconvwnr(g, PSF); % 基本维纳滤波Sn = abs(fft2(noise)).^2; % 噪声功率谱
nA = sum(Sn(:)) / prod(size(noise)); % 噪声均方值
Sf = abs(fft2(f)).^2; % 原始图像功率谱
fA = sum(Sf(:)) / prod(size(f)); % 原始图像均方值
R = nA / fA; % 噪声与信号的比例
fr2 = deconvwnr(g, PSF, R); % 加噪声功率谱的维纳滤波

功能说明：
- deconvwnr 可接受不同参数（噪声功率、图像自相关）以优化恢复效果。
- edgetaper 函数用于减少边缘效应。

5、总结

Wiener 滤波优点：
- 能同时处理退化和噪声影响。
- 自适应于噪声和信号的比例，最小化恢复误差。
- 对比逆滤波，维纳滤波能有效控制噪声放大，恢复效果更好。
实际应用：
- 运动模糊恢复：运动模糊在特定方向的退化，通过 PSF 和 Wiener 滤波可进行有效恢复。
- 噪声抑制：适用于加性高斯噪声等不同类型噪声。
限制：
- 需要事先估计噪声功率谱和退化函数 $H(u,v)$ 。
- 当图像中噪声水平未知时，恢复效果可能受限。

八、约束最小二乘滤波

1、约束最小二乘滤波的理论基础

主要内容：
- 问题：恢复退化图像需要了解退化函数 H，这是常规方法的难点。
- Wiener滤波 的缺点：需要知道未退化图像和噪声的功率谱，实际应用中可能难以获得。
- 约束最小二乘滤波 的优势：只需知道噪声的均值和方差，而无需详细的功率谱信息。
- 模型公式：
  - $g(x,y)$ ：观测到的退化图像
  - $h(x,y)$ ：退化函数（点扩散函数，PSF）
  - $f(x,y)$ ：原始图像
  - $\eta(x,y)$ ：噪声
关键要点：
- 知道退化函数 H 是关键，但约束最小二乘方法可以通过平滑度约束在未知 HHH 的情况下实现图像复原。
- 优势：只需噪声均值和方差，降低了参数依赖性。

2、约束最小二乘滤波的数学推导

公式和推导：
- 目标是最小化平滑度约束下的代价函数 C：
  - 其中 $\nabla^2$ 是图像的二阶导数（拉普拉斯算子），用于衡量图像的平滑度。
  - 约束条件：
    - $H \hat{f}$ ：滤波后的估计图像
    - $\eta$ ：噪声水平
- 频域解：
  在频域中，最小二乘解为：
  - $H^*(u,v)$ ：退化函数的共轭
  - $\gamma$ ：平滑度约束参数
  - $P(u,v)$ ：拉普拉斯算子的频域表示，用于平滑度约束。
拉普拉斯模板：
- 拉普拉斯模板在空间域中体现了平滑约束。

3、噪声水平的估计与约束参数

残差 $r$ ：
- 残差的平方和 $\| r \|^2$ 表明了重建图像与退化图像的差异。
噪声方差的估计：
- 其中 $\sigma_\eta^2$ ：噪声方差， $m_\eta$ ：噪声均值。
噪声方差计算公式：
$\sigma_\eta^2 = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} [\eta(x,y) - m_\eta]^2$ $m_\eta = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} \eta(x,y)$
关键点：
- 通过估计噪声水平，可以调整平滑度约束参数 $\gamma$ ，以获得最优结果。

4、约束最小二乘滤波的实验效果

图5.29 (a)由运动模糊及加性噪声污染的图像，(b)逆滤波的结果，(c)维纳滤波的结果。(d)~(f)同样的序列,但是噪声幅度的方差小了一个数量级,,(g)-(i)同样的序列，但是噪声方差比(a)小5个数景级,注意,在(h)中,去模糊图像透过噪声“宽帘”清晰可见。

图5.30 约東最小二乘方滤波的结果，用图5.29(c)，(f)和(i)中的结果分别与维纳滤波(a),(b),(c)相比较

直接逆滤波 (Inverse Filtering)：
- 图5.29(b)：通过直接逆滤波恢复图像。
- 特点：直接恢复，但对噪声极其敏感，容易放大高频噪声，导致恢复质量差。
- 图中表现：图像呈现出严重的噪声和模糊。
维纳滤波 (Wiener Filtering)：
- 图5.29(c, f, i) ：维纳滤波在频域中平衡退化函数和噪声功率谱，得到估计图像。
- 特点：通过引入噪声功率与信号功率比，降低了对噪声的敏感度，恢复效果较好。
- 图中表现：相较于直接逆滤波，噪声有所抑制，但高频细节仍然可能丢失。
约束最小二乘滤波 (Constrained Least Squares Filtering)：
- 图5.30(a, b，c)：约束最小二乘滤波的结果。
- 特点：基于图像平滑度约束（如拉普拉斯约束），通过调节参数 $\gamma$ 使得噪声影响最小化。
- 图中表现：恢复图像具有更好的平滑性与锐度，相较于维纳滤波（图5.29(c, f, i) ），细节保留更好，噪声较少。

附录

1、随机噪声生成的理论基础

（1）理论背景与核心内容

这一部分基于 累积分布函数（CDF）逆变换法，描述了如何通过均匀分布随机变量 $w \in (0, 1)$ 生成服从指定分布的随机变量 z 。具体原理包括：

累积分布函数 $F_Z(z)$ ：
定义 $F_Z(z)$ 是随机变量 Z 的CDF，即：
$F_Z(z) = P(Z \leq z)$
它表示 Z 小于等于某个值 z 的概率。
逆变换法原理：
如果 $w \sim U(0,1)$ （均匀分布），则可以通过逆CDF $F_Z^{-1}$ 将 w 映射为服从指定分布的随机变量 z：
$z = F_Z^{-1}(w)$

（2）原理推导

均匀分布的特性：
均匀分布 $U(0,1)$ 的定义是：在区间 $(0, 1)$ 上，每个 w 的出现概率是等同的。换句话说，w 是一个在 $(0, 1)$ 区间上随机取值的变量。
CDF的性质：
- $F_Z(z)$ 是单调递增的，且取值范围为 $[0,1]$ 。
- 若我们设定 $w = F_Z(z)$ ，则 w 也会落在区间 (0,1) 上。
逆变换的关键：
- 已知 $w \sim U(0,1)$ ，我们需要找到满足 $F_Z(z) = w$ 的 z 。
- 由于 $F_Z(z)$ 是单调递增的，可以通过其逆函数 $F_Z^{-1}$ 求解 z：
这样 z 就服从 $F_Z(z)$ 对应的分布。
数学直观：
- 将均匀分布的随机变量 w 通过CDF $F_Z$ 映射到 [0,1] 区间上；
- 然后通过逆CDF $F_Z^{-1}$ 再映射回原始的随机变量 z。
- 这种映射保证了 z 的分布与 $F_Z$ 一致。

（3）为什么逆变换法有效？

保持分布一致性：
- 均匀分布随机变量 $w \sim U(0,1)$ 被逆CDF $F_Z^{-1}$ 转换后，其分布自然服从指定的 $F_Z$ 。
单调性保证唯一性：
- 因为 $F_Z(z)$ 是单调递增函数，逆函数 $F_Z^{-1}$ 也是单调的，因此对于任意 $w \in (0,1)$ ，对应的 z 是唯一的。
概率守恒：
- 通过逆CDF变换，w 的概率密度在 [0,1] 区间上的均匀性被“映射”到 z 上，保证了 z 的分布符合 $F_Z$ 。