数字图像处理(5):图像复原

Last updated on November 20, 2025 pm

这是SJTU-CS3324《数字图像处理》课程的知识点整理系列。本文整理部分为“第 5 章:图像复原”。

图像复原的目标与基本概念

  • 图像复原的目标: 提升图像质量(可以是主观的或客观的)
  • 图像增强与图像复原的区别:
    • 图像增强 (Image Enhancement):
      • 主观性 (Subjective): 这是一个主观的过程,目的是让处理后的图像在视觉上更适合于某个特定应用
      • 无先验模型: 通常不考虑图像是如何“变坏”的,没有一个具体的退化模型
    • 图像复原 (Image Restoration):
      • 客观性 (Objective): 这是一个客观的过程,试图将退化(degraded)的图像恢复到其原始状态
      • 基于先验模型: 图像复原试图对退化过程进行建模,需要事先知道(或估计出)图像退化的原因,例如模糊的类型、噪声的统计特性等
      • 逆过程: 复原可以看作是退化过程的逆向操作 (reverse of the degradation)

图像退化/复原过程的模型 (Model of the Image Degradation/Restoration Process)

模型结构

  • 退化过程 (DEGRADATION):
    1. 原始的、清晰的图像 f(x,y)f(x,y) 首先经过一个退化函数 HH(Degradation function)的处理
    • 这个函数 H 模拟了造成图像质量下降的各种因素,最典型的是模糊
    1. 然后,一个与图像内容无关的噪声项 η(x,y)\eta(x,y)(Noise)以加性的方式被引入
    2. 最终得到我们观察到的、质量下降的图像 g(x,y)g(x,y)
  • 复原过程 (RESTORATION):
    • 将退化图像 g(x,y)g(x,y) 输入到一个或多个复原滤波器(Restoration filter(s))中
    • 滤波器的目标是产生一个对原始图像的最佳估计 f^(x,y)\hat{f}(x,y)
    • 理想情况下,我们希望 f^(x,y)\hat{f}(x,y) 尽可能地接近 f(x,y)f(x,y)

模型的数学表示

  • 通用形式:

g(x,y)=H[f(x,y)]+η(x,y)g(x,y) = H[f(x,y)] + \eta(x,y)

  • 线性、位置不变 (LPI) 系统下的简化: 如果退化函数 HH 是一个线性、位置不变 (Linear, Position-Invariant, LPI) 的系统,模型可以被大大简化
    • 空间域 (Spatial Domain): 当 HH 是 LPI 系统时,退化过程在空间域表现为卷积

      g(x,y)=h(x,y)f(x,y)+η(x,y)g(x,y) = h(x,y) * f(x,y) + \eta(x,y)

    • 频率域 (Frequency Domain): 根据卷积定理,空间域的卷积等价于频率域的乘法

      G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v) = H(u,v)F(u,v) + N(u,v)

噪声模型与噪声去除 (Noise Model & Noise Reduction)

这一部分聚焦于退化模型中的噪声项 η(x,y)\eta(x,y),探讨它的来源、数学模型以及如何去除它。

噪声的来源与分类

  • 主要来源:
    • 图像采集 (Image Acquisition): 图像传感器的性能和工作环境是噪声的主要来源,例如电路噪声、传感器因光照不足或温度过高产生的噪声
    • 图像传输 (Image Transmission): 图像在信道中传输时,可能会受到电磁干扰等,从而引入噪声
    • 随机噪声: 当退化只有噪声时(HH 为单位算子),主要使用空间滤波 (Spatial Filtering) 来处理
    • 周期噪声: 主要使用频率域滤波 (Frequency Domain Filtering) 来处理

随机噪声的数学模型 (Noise Models for Random Noise)

  • 基本假设:
    • 噪声在空间上是独立的(周期噪声除外),并且与图像内容本身也无关
    • 噪声被看作是随机变量,其特性可以通过其概率密度函数 (Probability Density Function, PDF) 来描述

  • 常见的噪声 PDF 模型:
    • 高斯噪声 (Gaussian Noise):
      • 成因:来源于电子电路噪声、传感器在光照不足或高温下的噪声
      • PDF: 呈经典的钟形曲线
    • 瑞利噪声 (Rayleigh Noise):
      • 成因: 主要用于描述距离图像 (range imaging) 中的噪声
      • PDF: 是向右偏斜的
    • 爱尔兰(伽马)噪声 (Erlang (Gamma) Noise) & 指数噪声 (Exponential Noise):
      • 成因: 常用于描述激光成像中的噪声
      • PDF: 都是向右偏斜的,指数噪声是伽马噪声的一个特例
    • 均匀噪声 (Uniform Noise):
      • 成因: 来源于量化过程中的误差。
      • PDF: 在某个范围内,每个灰度值出现的概率是均等的
    • 脉冲噪声(椒盐噪声) (Impulse (Salt-and-Pepper) Noise):
      • 成因: 由快速的瞬时干扰引起
      • PDF: 只有两个值(通常是 0 和 255),表现为图像中的纯黑点(“胡椒”)和纯白点(“盐”)

周期噪声 (Periodic Noise)

  • 成因: 主要来源于图像采集过程中电力或机电设备的干扰
  • 特性:
    • 空间依赖性 (Spatial Dependence): 与随机噪声不同,周期噪声在空间上不是独立的,而是呈现周期性重复的模式
    • 频域表现: 周期噪声在频率域表现得非常明显,通常是在频谱图上出现成对的、孤立的亮点(脉冲)

噪声参数的估计 (Estimation of Noise Parameters)

  • 目的: 为了选择合适的滤波器并设置其参数,我们需要估计出图像中噪声的统计特性,如均值、方差等
  • 估计方法:
    • 周期噪声: 直接观察图像的频谱图 (Frequency Spectrum),找到其中的亮点,即可确定其频率和方向
    • 随机噪声:
      • 情况1:成像系统可用: 拍摄一张“平坦”环境的图像,理论上该图像所有像素值应相同,实际的波动就完全由噪声引起,直接分析这张噪声图像的直方图和统计量即可
      • 情况2:只有带噪图像: 在带噪图像中,选取一小块内容相对恒定的区域(即图像本身的灰度变化很小),然后分析这个小块的直方图
        • 这个小块的直方图形状可以近似看作是噪声的 PDF,通过计算这个小块的均值和方差,就可以估计出噪声的参数

仅有加性噪声时的复原 - 空间滤波

当退化模型简化为 g(x,y)=f(x,y)+η(x,y)g(x,y) = f(x,y) + \eta(x,y) 时,可以使用各种空间滤波器来抑制噪声 η\eta,这些技术与上一章“图像增强”中的平滑滤波非常相似。

  • 均值滤波器 (Mean Filters): 算术均值滤波器、几何均值滤波器、谐波均值滤波器、逆谐波均值滤波器
  • 顺序统计滤波器 (Order-statistics Filters): 中值滤波器、最大/最小值滤波器、中点滤波器、阿尔法截尾均值滤波器
  • 自适应滤波器 (Adaptive Filters): 滤波器的行为会根据邻域的局部统计特性进行自适应调整

周期噪声的去除 - 频率域滤波

由于周期噪声在频谱图上表现为孤立的亮点,因此可以在频率域精确地将其滤除。

  • 带阻滤波器 (Bandreject Filters): 拒绝特定频率范围内的信号通过
    • 如理想、巴特沃斯和高斯三种类型的带阻滤波器,它们可以在以原点为中心的环形区域内抑制频率分量
  • 陷波滤波器 (Notch Filters): 带阻滤波器的推广,它可以在频谱图的任意位置拒绝(或通过)一个预定义邻域内的频率

线性、位置不变的退化 (Linear, Position-Invariant Degradations)

  • 线性 (Linear): H[af1(x,y)+bf2(x,y)]=aH[f1(x,y)]+bH[f2(x,y)]H[a*f_1(x,y) + b*f_2(x,y)] = a*H[f_1(x,y)] + b*H[f_2(x,y)]
  • 位置不变 (Position-Invariant): 退化效果与图像内容的位置无关,在 A 处造成的模糊和在 B 处造成的模糊是完全一样的
    • 数学上讲,如果对输入的平移会导致输出产生完全相同的平移,那么该系统就是位置不变的

      H[f(xα,yβ)]=g(xα,yβ)H[f(x-\alpha, y-\beta)] = g(x-\alpha, y-\beta)

      其中 (α,β)(\alpha, \beta) 是任意的空间位移
  • 将退化过程表达为卷积: 当 H 是LPI系统时,退化过程在空间域表现为卷积

    g(x,y)=h(x,y)f(x,y)+η(x,y)g(x,y) = h(x,y) * f(x,y) + \eta(x,y)

    其中 h(x,y)h(x,y) 是退化函数 HH 的脉冲响应,也称为点扩散函数 (Point Spread Function, PSF),它描述了系统对一个单位点光源的响应,即一个理想的点光源在经过退化系统后会变成什么样子

连续图像的复原 (Restoration for Continuous Images)

复原问题的核心与困难

  • 核心问题: 给定观察到的退化图像 g(x,y)g(x,y),以及退化模型 g(x,y)=h(x,y)f(x,y)+η(x,y)g(x,y) = h(x,y) * f(x,y) + \eta(x,y),并且我们对退化函数 h(x,y)h(x,y) 和噪声 η(x,y)\eta(x,y) 的统计特性有所了解,我们的目标是找到原始图像 f(x,y)f(x,y) 的一个最优估计 f^(x,y)\hat{f}(x,y)
  • 病态问题 (Ill-conditioned Problem): 图像复原在数学上是一个典型的病态问题,这意味着解可能不存在、不唯一

逆滤波 (Inverse Filtering)

逆滤波是最直接、最直观的复原方法。

  • 思想: 从频率域的退化模型 G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v) = H(u,v)F(u,v) + N(u,v) 出发
  • 推导:
    • 首先做一个理想化的假设:不存在噪声,即 N(u,v)=0N(u,v) = 0,模型简化为

      G(u,v)=H(u,v)F(u,v)G(u,v) = H(u,v)F(u,v)

    • 直接通过代数运算求解 F(u,v)F(u,v):

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

    • 最后,通过傅里叶逆变换得到复原后的图像估计:

      f^(x,y)=F1[G(u,v)H(u,v)]\hat{f}(x,y) = F^{-1}\left[\frac{G(u,v)}{H(u,v)}\right]

  • 过程: 这个过程可以看作是将退化图像 G(u,v)G(u,v) 通过一个传递函数为 1/H(u,v)1/H(u,v) 的“逆”滤波器
  • 缺陷:
    • 对噪声极其敏感: 现实中噪声 N(u,v)N(u,v) 总是存在的,除法操作 G/HG/H 会变成 (HF+N)/H=F+N/H(HF+N)/H = F + N/H,而退化函数 H(u,v)H(u,v) 在高频区域的值会非常小,这会导致 N/HN/H 项变得非常大,极大地放大了噪声
    • H(u,v)H(u,v) 为零的问题: 如果 H(u,v)H(u,v) 在某些频率上恰好为零,那么除法将无法进行,信息在这些频率上已经丢失

维纳滤波 (Wiener Filtering)

维纳滤波是一种更先进、更实用的方法,它克服了逆滤波对噪声敏感的缺点。

  • 思想: 维纳滤波不再追求对 f(x,y)f(x,y) 的精确还原,而是在存在噪声的情况下,寻找一个最优的估计 f^(x,y)\hat{f}(x,y)
  • 优化目标: 其优化准则是最小化均方误差 (Minimum Mean Square Error, MMSE),即使得估计图像 f^\hat{f} 与原始图像 ff 之间的差的平方的期望值最小:

Minimizee2=E[f(x,y)f^(x,y)2]\text{Minimize}\quad e^2 = \mathbb{E}\left[ \left|f(x,y) - \hat{f}(x,y)\right|^2 \right]

  • 维纳滤波器的表达式: 通过变分法推导,可以得到维纳滤波器在频率域的表达式:

    F^(u,v)=[1H(u,v)H(u,v)2H(u,v)2+Sn(u,v)Sf(u,v)]G(u,v)\hat{F}(u, v)=\left[\frac{1}{H(u, v)} \frac{|H(u, v)|^2}{|H(u, v)|^2+\frac{S_n(u, v)}{S_f(u, v)}}\right] G(u, v)

    其中各项的定义为:
    • H(u,v)H(u,v): 退化函数(的傅里叶变换)
    • H(u,v)2=H(u,v)H(u,v)|H(u,v)|^2 = H^*(u,v)H(u,v): 退化函数的功率谱,HH^*HH 的复共轭
    • Sn(u,v)=N(u,v)2S_n(u,v) = |N(u,v)|^2: 噪声的功率谱
    • Sf(u,v)=F(u,v)2S_f(u,v) = |F(u,v)|^2: 原始图像的功率谱
  • 特点与优势:
    • 自适应性: 维纳滤波器是自适应的。方括号中的项 […] 是一个调节因子:
      • 当信噪比 Sf/SnS_f/S_n 很高时,Sn/SfS_n/S_f 接近于0,整个表达式趋近于逆滤波 G/HG/H,此时以去模糊为主
      • 当信噪比很低时,Sn/SfS_n/S_f 很大,分母变得很大,导致整个表达式的值趋近于0,从而抑制了噪声的放大
    • 需要先验知识: 维纳滤波需要我们事先估计出噪声的功率谱 SnS_n 和原始图像的功率谱 SfS_f,这在实际应用中往往是困难的
    • 与逆滤波的关系: 如果噪声为零 Sn(u,v)=0S_n(u,v)=0,维纳滤波退化为逆滤波

离散图像的复原 (Restoration for Discrete Images)

离散退化模型

  • 思想: 将连续域的积分和函数,用离散域的求和与矩阵来表示
  • 模型表示:
    • 空间域: 连续的卷积运算变为离散卷积,这个过程可以用矩阵和向量的乘法来表示,即通过堆叠 (Stacking) 将图像和卷积核矩阵化

      g(m,n)=f(m,n)h(m,n)g=[h]fg(m,n) = f(m,n) * h(m,n) \quad \to \quad g = [h]f

    • 频率域: 离散傅里叶变换下,卷积定理依然成立

      G(u,v)=F(u,v)H(u,v)G=[H]FG(u,v) = F(u,v) · H(u,v) \quad \to \quad G = [H]F

逆滤波 (Inverse Filtering)

  • 思想: 与连续域完全相同,即在频率域直接做除法,但在离散情况下,我们通常从最小化代价函数 J[f^]J[\hat{f}] 的角度来推导
    • 假设无噪声,退化图像 gg 与复原图像 f^\hat{f} 经过同样退化 [h][h] 后的差 g[h]f^g - [h]\hat{f} 的范数(能量)应最小

      MinimizeJ[f^]=g[h]f^2\text{Minimize}\quad J[\hat{f}] = \left\Vert g - [h]\hat{f}\right\Vert²

    • f^\hat{f} 求导并令其为零,可以解出:

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

  • 问题: 对噪声极其敏感

维纳滤波 (Wiener Filtering)

  • 思想与优化目标: 与连续域相同,目标是最小化均方误差

e2=E[ff^2]e^2 = \mathbb{E}\left[ \left\Vert f - \hat{f} \right\Vert^2 \right]

  • 离散域下的公式: 最终得到的频率域滤波器表达式与连续情况完全一致

F^(u,v)=[1H(u,v)H(u,v)2H(u,v)2+Sn(u,v)Sf(u,v)]G(u,v)\hat{F}(u, v)=\left[\frac{1}{H(u, v)} \frac{|H(u, v)|^2}{|H(u, v)|^2+\frac{S_n(u, v)}{S_f(u, v)}}\right] G(u, v)

  • 挑战: 需要预先知道或估计噪声功率谱 SnS_n 和原始图像功率谱 SfS_f

约束最小二乘滤波 (Constrained Least Squares Filtering)

  • 思想: 维纳滤波需要知道噪声和图像的谱,这在实践中很难,我们能否找到一种只需要少量先验知识的复原方法?约束最小二乘法的思想是,我们希望找到一个解 $\hat{f},它在满足某个约束条件的同时,使某个目标函数最小化
  • 优化目标与约束:
    • 目标函数: 我们希望复原后的图像 f^\hat{f} 具有一定的平滑性,一个衡量图像平滑度的常用指标是其拉普拉斯算子的能量 2f^2\Vert \nabla^2 \hat{f}\Vert^2。因此,我们的目标是最小化 Qf^2\Vert Q \hat{f} \Vert^2,其中 QQ 是拉普拉斯算子
    • 约束条件: 复原后的图像 f^\hat{f} 经过退化 hh 后,与原始退化图像 gg 之间的残差 ghf^g - h\hat{f} 应该约等于噪声 nn。因此,我们约束这个残差的能量等于我们估计出的噪声能量 n2\Vert n \Vert^2

      ghf^2=n2\left\Vert g - h \hat{f} \right\Vert^2 = \Vert n \Vert^2

  • 求解与公式: 通过拉格朗日乘子法求解这个约束优化问题,可以得到频率域下的解

    F^(u,v)=[H(u,v)H(u,v)2+1αQ(u,v)2]G(u,v)\hat{F}(u, v)=\left[\frac{H^*(u, v)}{|H(u, v)|^2+\frac{1}{\alpha}|Q(u, v)|^2}\right] G(u, v)

    其中 α\alpha 是拉格朗日乘子,需要迭代调整以满足约束条件,Q(u,v)Q(u,v) 是拉普拉斯算子的傅里叶变换

图像重建 (Image Reconstruction)

基本概念

  • 核心问题: 图像重建的目标是从物体的“投影” (Projections) 数据来重建物体内部的“切片” (Slice) 图像
  • CT 简介: CT是一种革命性的医学成像技术,它从多个不同角度发射 X 射线穿过身体,并测量穿透后的射线强度,然后通过计算机算法来重建身体内部的横截面图像
  • CT 工作流程:
    1. 一个 X 射线源和一个探测器阵列环绕着被检测物体(病人)
    2. X 射线源发射一束窄射线,穿过物体,被另一端的探测器接收
    3. X 射线源和探测器同步旋转一个小的角度
    4. 重复步骤2和3,直到完成至少 180 度的扫描,从而获得一系列在不同角度下的投影数据
    5. 使用断层扫描算法(Tomography algorithms)处理这些投影数据,重建出代表物体内部一个“切片”的图像
    6. 通过在垂直于扫描环的方向上移动物体,可以获得一系列的切片图像,最终构成三维的内部结构

参考资料

本文参考上海交通大学电子工程系《数字图像处理》课程 CS3324 闵雄阔老师的 PPT 课件整理。


数字图像处理(5):图像复原
https://cny123222.github.io/2025/11/20/数字图像处理-5-:图像复原/
Author
Nuoyan Chen
Posted on
November 20, 2025
Licensed under