5.2.1 图像退化的物理过程#

一张图像从拍摄到保存，可能受到：

1. 光学模糊：焦外、衍射、景深
2. 运动模糊：曝光期间相机/被摄物移动
3. 大气扰动：湍流、雾、水
4. 传感器噪声：热噪声、暗电流、放大噪声
5. 压缩伪影：JPEG 块效应
6. 传输噪声：椒盐、条纹

5.2.2 通用模型#

定义 5.1（线性退化 + 加性噪声）：

其中：

• $f$：原图（未知，想要恢复）
• $h$：点扩散函数 (Point Spread Function, PSF)，表征退化系统
• $\eta$：加性噪声
• $g$：观测到的图像（已知）

频域形式：

其中 $H$ 称为光学传递函数 (Optical Transfer Function, OTF)。

5.2.3 PSF 的意义#

PSF = 系统对一个"点光源"的响应。

一个理想点源（δ 冲激）经过退化系统，会被"扩散"成一个小斑点 $h(x, y)$，这就是 PSF。

整幅图像可看作无数点源的叠加（线性），因此：

PSF 的典型形态：

5.2.4 复原问题的核心困难#

已知 $g, h$（甚至不知道 $h$），要求 $f$：

困难 1：$h$ 在某些频率接近零（比如运动模糊在运动方向的高频），除法会放大噪声到无穷。

困难 2：噪声 $\eta$ 未知，只能估计其统计特性。

困难 3：病态 (Ill-posed) ——小扰动 $\Delta g$ 可能导致巨大的 $\Delta f$。

复原算法的所有技巧，都在对抗这三个困难。

5.3.1 常见噪声的 PDF#

5.3.2 关键统计量#

5.3.3 从图像估计噪声参数#

找一块均匀区域（比如天空、墙面），计算其像素值直方图：

1. 直方图对称钟形 → 高斯，拟合 $\mu, \sigma$
2. 直方图偏态 → 瑞利 / Gamma
3. 直方图两端有尖峰 → 椒盐

鲁棒估计：用中位绝对偏差 (MAD)：

MAD 对异常值鲁棒（比标准差好）。

5.3.4 周期噪声#

与上面"随机噪声"不同，周期噪声有确定的空间频率（如扫描仪机械振动、电源耦合），只能在频域识别（亮斑对），空域滤波几乎无效 → 陷波滤波（第 4 章）。

5.4.1 均值类滤波#

5.4.2 排序统计类（非线性）#

5.4.3 自适应局部降噪（经典技巧）#

5.4.4 自适应中值滤波（强鲁棒）#

1
Level A:
2
  A1 = z_med - z_min
3
  A2 = z_med - z_max
4
  如果 A1 > 0 且 A2 < 0 （即 z_med 不是极值）
5
      → 进入 Level B
6
  否则
7
      → 扩大窗口 +2；若超过最大尺寸 → 输出 z_med
8

9
Level B:
10
  B1 = z_xy - z_min
11
  B2 = z_xy - z_max
12
  如果 B1 > 0 且 B2 < 0 （中心不是极值）
13
      → 输出 z_xy（保留）
14
  否则
15
      → 输出 z_med（被替换）

5.5.1 最优陷波滤波（Gonzalez）#

简单陷波会留下残影。最优陷波：

1. 提取噪声分量：$\hat{\eta} = \mathcal{F}^{-1}\{G \cdot (1 - H)\}$
2. 从 $g$ 减去 $w \cdot \hat{\eta}$：

3. $w$ 通过最小化邻域方差求得（使 $\hat{f}$ 局部最"干净"）：

（上划线表示邻域平均）

5.6.1 逆滤波（最朴素方案）#

由 $G = H F + N$：

致命问题：

• 当 $|H(u, v)|$ 小（高频位置很常见）时，$N / H$ 被放大到无穷
• 原图信号 $F$ 被彻底淹没

1
直接逆滤波的结果：看起来像一团彩色雪花。

改进 1：低通截断

只在 $|H| > T$ 处做除法，其余置 0（或低通外等价）：

改进 2：Tikhonov 正则化

$\alpha$ 越大越稳定（偏向平滑），越小越贴近逆滤波。实际上就是维纳滤波的简化。

5.6.2 维纳滤波（Wiener，MMSE 最优）#

目标：最小化均方误差 $E[\|f - \hat{f}\|^2]$。

推导（简要）： 假设 $\hat{F} = W \cdot G$。对 $W$ 求导使 MSE 最小，得：

公式：

• $S_\eta$：噪声功率谱
• $S_f$：原信号功率谱
• $K = S_\eta / S_f$：噪信比（工程上常设常数）

两种极端：

• $K = 0$（无噪声）：$W = 1/H$ → 退化为逆滤波
• $H = 0$（某频率被完全抑制）：$W = 0$ → 不放大噪声

1
def wiener(g, H, K=0.01):
2
    G = np.fft.fft2(g)
3
    H_mag2 = np.abs(H)**2
4
    W = np.conj(H) / (H_mag2 + K)
5
    return np.real(np.fft.ifft2(G * W))

5.6.3 约束最小二乘（CLS）#

维纳滤波需要知道两个功率谱。若无法估计，改用 CLS。

目标：最小化"平滑度"，约束残差等于噪声能量：

$C$ 常取拉普拉斯算子（梯度平方度量平滑度）。

拉格朗日法求得：

其中 $P$ 是拉普拉斯的 FT，$\gamma$ 是拉格朗日系数。

与 Wiener 对比：

• Wiener：分母加 $S_\eta / S_f$（依赖功率谱）
• CLS：分母加 $\gamma |P|^2$（只需一个系数，工程更实用）

迭代调整 $\gamma$：

1. 初始 $\gamma$
2. 算 $\hat{F}, \hat{f}$
3. 计算残差 $r = \|g - H \hat{f}\|^2$
4. 如果 $r$ 太小 → 增大 $\gamma$；太大 → 减小
5. 直到 $r \approx \|\eta\|^2$

5.6.4 Lucy-Richardson 迭代（极大似然）#

假设噪声是泊松分布（光子计数场景），用 EM 算法导出迭代：

特点：

• 保证非负
• 收敛慢，但质量高
• 经典应用：天文望远镜（修复哈勃太空望远镜早期图像）、荧光显微镜

1
from skimage.restoration import richardson_lucy
2
deconv = richardson_lucy(blurry, psf, num_iter=30)

5.6.5 盲反卷积（Blind Deconvolution）#

现实中 $h$ 往往未知。

盲反卷积同时估计 $f$ 和 $h$。这是超病态问题，需要额外先验：

• 自然图像的梯度分布：稀疏（重尾）
• PSF 的形状先验：非负、有限支集、归一化

经典方法：

• Fergus et al. 2006（贝叶斯 + 变分推断）
• Krishnan 2009（归一化稀疏正则）
• Levin et al. 2011（MAP 问题的理论分析）

深度学习完全碾压传统方法：

• DeblurGAN、MPRNet、NAFNet 等
• 端到端学习 $g \to \hat{f}$，无需显式 $h$

5.8.1 畸变模型（Brown-Conrady）#

设归一化坐标 $(x, y)$，畸变后 $(x_d, y_d)$：

其中 $r = \sqrt{x^2 + y^2}$。

• $(k_1, k_2, k_3)$：径向
• $(p_1, p_2)$：切向

5.8.2 相机标定（张正友法简述）#

用棋盘格拍多张不同角度照片：

1. 自动检测角点 → 每角点给出 (像素坐标, 世界坐标) 对
2. 估计每张图的单应 $H_i$
3. 由多个 $H_i$ 约束解内参 $K$ 和畸变系数
4. 非线性细化

1
ret, K, dist, rvecs, tvecs = cv2.calibrateCamera(obj_pts, img_pts, gray.shape[::-1], None, None)
2
undistorted = cv2.undistort(img, K, dist)

5.8.3 几何校正的一般流程#

1
畸变图 → 逆向映射 → 无畸变坐标 → 插值取值 → 校正图

关键：逆向映射（每个目标像素反查源）。

5.9.1 Radon 变换#

定义 5.2：对给定角度 $\theta$ 和距离 $\rho$，沿直线 $x\cos\theta + y\sin\theta = \rho$ 的线积分：

物理意义：X 光穿过物体，射线强度衰减 $\propto \int f$。每个投影就是一个 $\theta$ 下的"影子"。

正弦图 (Sinogram)：$R(\rho, \theta)$ 作为 $\rho, \theta$ 的函数图——每个物体点对应正弦图中一条正弦曲线。

5.9.2 中心切片定理（Fourier Slice Theorem）#

定理 5.1（核心）：$R(\rho, \theta)$ 关于 $\rho$ 的 1D 傅里叶变换，等于原函数 $f(x, y)$ 的 2D 傅里叶变换 $F(u, v)$ 沿过原点、方向 $\theta$ 的一维切片。

证明思路：

几何直觉：

1
2D Fourier plane:           Sinogram:
2
    ▲ v                           R(ρ, θ)
3
    │
4
  切片 ─→  = 1D FFT of  ─→ projection at θ
5
    │     /
6
    │   /
7
    │ / θ
8
────┼────→ u

推论：只要我们收集足够多角度的投影，就能"填满"整个 2D 频谱，再 IFT 即得原图。

5.9.3 直接反投影（失败方案）#

对每个 $(\rho, \theta)$，把 $R$ 的值均匀涂回对应直线上：

结果：图像模糊，有"星芒伪影"（因为重复叠加低频）。不可用。

5.9.4 滤波反投影（FBP）#

核心思想：反投影前，先对每个投影做1D 滤波（Ram-Lak 斜坡滤波）。

算法：

1
1. 对每个角度 θ:
2
   1a. 对 R(·, θ) 做 1D FFT → R̂(ω, θ)
3
   1b. 乘以 ramp 滤波器 |ω|   (+ 可选窗函数如 Hamming)
4
   1c. IFFT → R̃(ρ, θ)  (滤波后的投影)
5
2. 反投影:
6
   f(x, y) = ∫ R̃(x cosθ + y sinθ, θ) dθ

为什么乘 $|\omega|$？ 将直接反投影对应的低通响应 $1/|\omega|$ 校正为平坦响应，防止低频过度叠加。

FBP 是 CT 医疗设备的标准算法（40 年主力）。

1
from skimage.transform import radon, iradon
2
sinogram = radon(img, theta=np.arange(180))
3
recon = iradon(sinogram, theta=np.arange(180), filter_name='ramp')

5.9.5 代数重建（ART）与迭代重建#

将投影过程写成矩阵方程：

$A$ 是"投影矩阵"，每行对应一条射线。

ART 迭代（Kaczmarz）：

对每条射线 $i$ 逐一更新，遍历所有射线为一轮。

SIRT（Simultaneous）：一次用所有射线，更稳定但慢。

优势：

• 不需要 $\theta$ 密集采样
• 可加先验（非负、平滑、稀疏）
• 低剂量 CT、稀疏视图重建的主力

5.9.6 现代方向：深度学习重建#

• FBPConvNet：FBP + U-Net 去伪影
• LEARN / MoDL：把迭代算法展开为网络
• NeRF-based CT：神经辐射场重建 3D 体积
• 效果显著优于 FBP，尤其在低剂量下

✅ 必须掌握#

• 退化模型 $g = h * f + \eta$ 及频域形式
• 常见噪声 PDF 与识别方法
• 逆滤波为何失败，如何改进
• Wiener 滤波公式推导与 $K$ 的意义
• 运动模糊的 PSF 空域与频域形式
• 中心切片定理
• FBP 流程（尤其为何要 ramp 滤波）

💡 高频面试题#

Q1. 逆滤波为何不可用？如何改进？

答：退化模型 $G = HF + N$。逆滤波 $\hat{F} = G/H = F + N/H$。在 $H$ 接近 0 的频率，$N/H$ 被无限放大。 改进：

• 低通截断（只在 $|H| > T$ 处做除法）
• 加正则项（Tikhonov）
• 维纳滤波：分母加 $S_N/S_f$ 防爆

Q2. 推导维纳滤波公式。

答：设 $\hat{F} = W G$，目标 $\min E[|F - \hat{F}|^2]$。对 $W$ 求导置零，得

Q3. 如何估计运动模糊的 PSF？

答：

• 先验法：若已知曝光时长与相机运动（IMU），直接画 PSF
• 频域法：观察 $G$ 频谱中的黑条纹（sinc 零点），沿运动方向均匀分布。 从零点间距 $1/a$ 反推运动长度 $a$，条纹方向给出运动角度。

Q4. 中心切片定理是什么？

答：投影数据 $R(\rho, \theta)$ 关于 $\rho$ 做 1D FFT 得到的结果，等于原图 $f$ 的 2D FFT 沿方向 $\theta$ 过原点的切片。意义：多角度投影可"填满"2D 频谱。

Q5. FBP 里为什么要乘 ramp 滤波器？

答：直接反投影相当于在频域每个频率累加 $1/|\omega|$ 次（角度积分），低频被过度强调。乘 $|\omega|$ 校正后响应平坦。

Q6. 自适应中值滤波的优势？

答：椒盐噪声比例 > 30% 时，标准中值也会被污染（窗口内多数噪声）。自适应中值先扩大窗口直到中值不是极值（即不受椒盐干扰），再判断中心像素是否要替换。对极端椒盐强鲁棒。

Q7. 椒盐 + 高斯混合噪声用什么？

答：先去椒盐再去高斯。椒盐是脉冲，若先做均值/高斯会将脉冲扩散成圆斑；先用自适应中值去椒盐（信号几乎不受影响），再用 NLM/BM3D 去高斯。也可以直接 Alpha-trimmed。

Q8. Lucy-Richardson 和 Wiener 的区别？

答：

• Wiener 在频域线性最优（MMSE），假设高斯噪声
• Lucy-Richardson 迭代 ML 估计，假设泊松噪声（光子计数），保证非负；在天文、显微更适合