空间域恢复与边缘保持滤波

June 16, 2026 物理世界开发计算机视觉, 图像处理, 反卷积, 双边滤波计算机视觉系列 4334 字 9 分钟阅读

上一篇介绍了频率域的图像处理方法——通过傅里叶变换将图像转换到频域，然后进行滤波、复原等操作。但频率域方法有时不太直观，尤其是我们更习惯直接操作像素。

空间域方法直接在像素空间中处理图像，不需要进行变换。本章介绍空间域的图像恢复、边缘保持滤波和锐化技术。

Lucy-Richardson 迭代反卷积

反卷积（deconvolution）的目标是从退化图像中恢复原始图像。退化过程通常建模为卷积：退化图像 = 原始图像卷积模糊核 + 噪声。

直接反卷积（逆滤波）在噪声存在时极不稳定——噪声会被高频部分的除法放大。Lucy-Richardson 算法提供了一种更稳健的迭代方法。

最大似然原理

Lucy-Richardson 算法基于最大似然估计，假设图像噪声服从泊松分布。泊松噪声的特点是方差等于均值——这在低光成像、X 射线成像等场景中很常见。

迭代公式

Lucy-Richardson 的迭代更新公式是：

$$ f_{k+1}(x, y) = f_k(x, y) \cdot \left( \frac{g(x, y)}{h * f_k(x, y)} * h(-x, -y) \right) $$

其中：

$f_k$ 是第 $k$ 次迭代的估计图像
$g$ 是观测到的退化图像
$h$ 是模糊核（point spread function）
$*$ 表示卷积
$h(-x, -y)$ 是模糊核的翻转

迭代行为

直觉上，这个公式在做两件事：

计算当前估计与观测的比值 $g / (h * f_k)$ —— 如果当前估计偏暗，比值 > 1，需要增强
用翻转的模糊核对比值进行反向扩散，修正像素值

迭代次数的选择很关键：

迭代太少：图像仍然模糊，恢复不足
迭代太多：细节恢复充分，但噪声被放大，出现过拟合

早期停止（early stopping）是常用策略——在噪声过度放大前停止迭代。

图1 - Lucy-Richardson 迭代过程：

mermaid
flowchart TD
    A["初始估计 f_0<br/>通常为均匀分布或退化图像"] --> B["计算卷积 h * f_k"]
    B --> C["计算比值 g / (h * f_k)"]
    C --> D["反向扩散: 比值 * h(-x, -y)"]
    D --> E["更新: f_{k+1} = f_k × 反向扩散"]
    E --> F{迭代次数<br/>足够?}
    F -->|否| B
    F -->|是| G["输出恢复图像 f_k"]

    classDef start fill:#FF9800,color:#fff
    classDef proc fill:#2196F3,color:#fff
    classDef check fill:#9C27B0,color:#fff
    classDef end fill:#4CAF50,color:#fff

    class A start
    class B,C,D,E proc
    class F check
    class G end

实际使用时，模糊核 $h$ 可以通过实验测量（如拍摄点光源），也可以用模糊边缘分析等方法估计。

全变分正则化

反卷积问题本质上是不适定问题——即使知道模糊核，从退化图像恢复原始图像也没有唯一解。正则化（regularization）通过引入先验约束，让解更稳定。

为什么选择 L1 范数

最直观的正则化是能量最小化：让恢复图像的总变分（total variation, TV）最小。总变分衡量图像的梯度总和：

$$ TV(f) = \sum_{x,y} \sqrt{(\nabla_x f)^2 + (\nabla_y f)^2} $$

这里用 L1 范数（梯度的绝对值）而不是 L2 范数（梯度的平方），是因为：

L2 范数会过度惩罚强梯度——边缘的梯度会被平滑掉
L1 范数对大梯度更宽容——允许边缘保持锐利

直觉上：L2 范数希望梯度处处都小（平滑），而 L1 范数允许梯度在某些地方很大（边缘），只要其他地方梯度足够小（平坦区域）。

优化目标

带 TV 正则化的图像恢复目标函数：

$$ \min_f | g - h * f |_2^2 + \lambda \cdot TV(f) $$

第一项是数据保真度（fidelity term）——恢复图像必须与观测图像一致。第二项是正则化项（regularization term）——恢复图像的总变分要小。$\lambda$ 是平衡参数：$\lambda$ 大则更平滑，$\lambda$ 小则更贴近观测。

这个优化问题通常用梯度下降或分裂 Bregman等迭代算法求解。实际使用时可以调用 scikit-image 的 restoration.denoise_tv_chambolle 或 restoration.unsupervised_wiener 等函数。

TV 正则化的局限

TV 正则化虽然能保持边缘，但有两个常见问题：

阶梯效应（staircase effect）——平坦区域的恢复会出现阶梯状伪影，因为 L1 范数倾向于产生分段常数区域
纹理丢失——细微纹理被过度平滑

更高级的正则化方法（如非局部均值、稀疏表示）可以缓解这些问题，但计算复杂度更高。

双边滤波

边缘保持滤波（edge-preserving filtering）的核心目标是：平滑噪声的同时不模糊边缘。传统高斯滤波无法做到这一点——它对所有像素一视同仁，边缘和平坦区域都会被平滑。

双边滤波（bilateral filter）通过引入像素值相似度作为额外权重，实现了边缘保持。

两维权重

双边滤波的每个像素有两个权重：

空间权重 $G_s$：由像素间的几何距离决定，距离越近权重越大
颜色权重 $G_r$：由像素值的差异决定，值越相似权重越大

完整公式：

$$ f(x, y) = \frac{1}{W} \sum_{(i,j) \in \Omega} I(i, j) \cdot G_s(|(x,y)-(i,j)|) \cdot G_r(|I(x,y) - I(i,j)|) $$

其中 $W$ 是归一化因子：

$$ W = \sum_{(i,j) \in \Omega} G_s(|(x,y)-(i,j)|) \cdot G_r(|I(x,y) - I(i,j)|) $$

$G_s$ 和 $G_r$ 通常都是高斯函数：

$$ G_s(d) = \exp\left(-\frac{d^2}{2\sigma_s^2}\right), \quad G_r(v) = \exp\left(-\frac{v^2}{2\sigma_r^2}\right) $$

$\sigma_s$ 控制空间邻域大小，$\sigma_r$ 控制颜色相似度的阈值。

为什么能保持边缘

双边滤波的魔法在于颜色权重 $G_r$：

在平坦区域，邻域像素值相近，$G_r \approx 1$，双边滤波退化为高斯滤波——平滑噪声。

在边缘处，邻域像素值差异大，$G_r \approx 0$，跨边缘的权重被抑制——边缘不被模糊。

图2 - 双边滤波的权重机制：

mermaid
flowchart TD
    CENTER["中心像素 P(x,y)<br/>I = 128"]
    P1["P1: I=125, 距离1.0<br/>颜色差3 → 权重0.93"]
    P2["P2: I=120, 距离1.4<br/>颜色差8 → 权重0.77"]
    P3["P3: I=30, 距离1.0<br/>颜色差98 → 权重≈0"]
    CENTER --> P1
    CENTER --> P2
    CENTER --> P3

    classDef center fill:#FF9800,color:#fff
    classDef neighbor fill:#2196F3,color:#fff
    classDef edge fill:#f44336,color:#fff
    class CENTER center
    class P1,P2 neighbor
    class P3 edge

上图显示，P3 虽然空间距离很近（$G_s \approx 0.95$），但因为颜色差异太大（跨越边缘），颜色权重 $G_r \approx 0$，最终权重接近 0——边缘被保留。

实际应用

双边滤波常用于：

高动态范围（HDR）成像：将多曝光图像融合后再用双边滤波调整局部对比度
卡通化效果：强双边滤波后量化颜色，产生卡通风格
图像抠图：用双边滤波优化 alpha 遮罩

OpenCV 提供 cv2.bilateralFilter() 函数：

python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import cv2
import numpy as np

# 读取图像
img = cv2.imread('input.jpg')

# 双边滤波
# d: 邻域直径（必须为奇数）
# sigmaColor: 颜色空间的标准差
# sigmaSpace: 坐标空间的标准差
filtered = cv2.bilateralFilter(img, d=9, sigmaColor=75, sigmaSpace=75)

双边滤波的缺点是计算量大——每个像素需要遍历邻域内所有像素，时间复杂度 $O(N \cdot r^2)$，其中 $N$ 是像素数，$r$ 是邻域半径。实际使用时，$d$ 通常不超过 15，否则处理速度会变慢。

双边滤波 vs 高斯滤波

对比一下双边滤波和高斯滤波在边缘处的行为：

高斯滤波：不管边缘是否存在，所有邻域像素都参与平滑 → 边缘模糊
双边滤波：跨边缘像素的权重被抑制 → 边缘保留

可以用实验验证：用 cv2.GaussianBlur() 和 cv2.bilateralFilter() 分别处理同一张图像，对比边缘区域——高斯滤波会让边缘变宽，双边滤波保持边缘锐利。

引导滤波

引导滤波（guided filter）是一种线性时间 $O(N)$ 的边缘保持滤波方法，由何恺明等人在 2010 年提出。相比双边滤波，它有两个优势：

速度快：线性时间复杂度，可快速处理大图像
无梯度反转：不会出现双边滤波的梯度反转 artifacts（梯度反转 artifact 指滤波后某些区域的梯度方向与原图相反）

局部线性模型

引导滤波的核心假设是：输出图像是引导图像的局部线性变换。对于每个局部窗口 $\omega_k$，有：

$$ q_i = a_k I_i + b_k, \quad \forall i \in \omega_k $$

其中：

$q_i$ 是输出图像的第 $i$ 个像素
$I_i$ 是引导图像的第 $i$ 个像素（引导图像可以是输入图像本身，也可以是另一幅图像）
$a_k$ 和 $b_k$ 是窗口 $\omega_k$ 内的线性系数（同一个窗口内所有像素共享 $a_k, b_k$）

如果 $a_k$ 小，说明输出图像在窗口内变化平缓（平坦区域）；如果 $a_k$ 大，说明输出图像跟随引导图像变化（边缘区域）。

系数求解

$a_k$ 和 $b_k$ 通过最小二乘拟合得到。目标是最小化：

$$ \min_{a_k, b_k} \sum_{i \in \omega_k} \left[ (a_k I_i + b_k - p_i)^2 + \epsilon a_k^2 \right] $$

其中 $p_i$ 是输入图像的第 $i$ 个像素。$\epsilon a_k^2$ 是正则化项，防止 $a_k$ 过大（过拟合噪声）。

解为：

$$ a_k = \frac{\frac{1}{|\omega|} \sum_{i \in \omega_k} I_i p_i - \mu_k \bar{p}_k}{\sigma_k^2 + \epsilon} $$

$$ b_k = \bar{p}_k - a_k \mu_k $$

其中：

$\mu_k$ 是引导图像在窗口 $\omega_k$ 内的均值
$\sigma_k^2$ 是引导图像在窗口 $\omega_k$ 内的方差
$\bar{p}_k$ 是输入图像在窗口 $\omega_k$ 内的均值
$|\omega|$ 是窗口内的像素数

直觉上：

方差 $\sigma_k^2$ 大（边缘区域），$a_k$ 偏向较大值 → 输出跟随引导图像
方差 $\sigma_k^2$ 小（平坦区域），$a_k$ 偏向 0 → 输出被平滑

实际应用

OpenCV 提供 cv2.ximgproc.guidedFilter() 函数（需要 opencv-contrib-python）：

python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
import cv2
import numpy as np

# 读取图像
img = cv2.imread('input.jpg')
guide = img  # 引导图像可以是原图或其他图像

# 引导滤波
filtered = cv2.ximgproc.guidedFilter(
    guide=guide,
    src=img,
    radius=8,      # 窗口半径
    eps=0.01       # 正则化参数
)

引导滤波常用于：

抠图细化：用原图作为引导图像，细化粗糙的 alpha 遮罩
去雾：用暗通道图作为引导，恢复清晰图像
HDR 压缩：用原图引导，调整局部对比度

双边滤波 vs 引导滤波

特性	双边滤波	引导滤波
时间复杂度	$O(N \cdot r^2)$	$O(N)$
边缘保持	是	是
梯度反转	可能出现	不会出现
可引导性	仅自身	可用其他图像引导

对于大图像或实时应用，引导滤波更合适。如果需要精细控制边缘保持程度，双边滤波的参数调节更直观。

拉普拉斯锐化

锐化（sharpening）与滤波相反——滤波是为了平滑，锐化是为了增强边缘对比度。拉普拉斯锐化是最基础的锐化方法。

二阶微分算子

拉普拉斯算子是一种二阶微分算子，定义为：

$$ \nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} $$

二阶导数对边缘（梯度变化剧烈处）响应强烈，对平坦区域响应接近 0。因此拉普拉斯响应可以看作是边缘强度的度量。

离散化与卷积核

在离散图像中，二阶导数用差分近似：

$$ \frac{\partial^2 f}{\partial x^2} \approx f(x+1, y) - 2f(x, y) + f(x-1, y) $$

$$ \frac{\partial^2 f}{\partial y^2} \approx f(x, y+1) - 2f(x, y) + f(x, y-1) $$

相加后得到4-邻域拉普拉斯核：

1
2
3
[ 0  -1   0]
[-1   4  -1]
[ 0  -1   0]

如果考虑对角线方向（8-邻域），核为：

1
2
3
[-1  -1  -1]
[-1   8  -1]
[-1  -1  -1]

锐化公式

将拉普拉斯响应加回原图，可以增强边缘对比度：

$$ g(x, y) = f(x, y) + c \cdot \nabla^2 f(x, y) $$

$c$ 是锐化强度：

$c > 0$：增强边缘（当核的中心系数为正，如 4 或 8）
$c < 0$：反向操作（当核的中心系数为负，如 -4 或 -8）

实践中，为了保持图像亮度不变，通常用中心系数为正的核（如 4-邻域核中心改为 4）：

1
2
3
[ 0  -1   0]
[-1   5  -1]
[ 0  -1   0]

这个核等价于 $g = f + \nabla^2 f$（因为 $5 = 4 + 1$）。

实际应用

OpenCV 用 cv2.filter2D() 应用自定义卷积核：

python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
import cv2
import numpy as np

img = cv2.imread('input.jpg')

# 4-邻域拉普拉斯锐化核
laplacian_kernel = np.array([
    [ 0, -1,  0],
    [-1,  5, -1],
    [ 0, -1,  0]
], dtype=np.float32)

# 8-邻域拉普拉斯锐化核
laplacian_kernel_8 = np.array([
    [-1, -1, -1],
    [-1,  9, -1],
    [-1, -1, -1]
], dtype=np.float32)

# 应用锐化
sharpened = cv2.filter2D(img, -1, laplacian_kernel)

锐化的局限

拉普拉斯锐化有几个问题：

噪声放大：锐化也会放大噪声，因为噪声在高频区域
过冲与下冲：过度锐化会在边缘两侧产生过冲（overshoot）和下冲（undershoot）——边缘周围出现亮环和暗环
不自然感：过度锐化的图像看起来不自然，像是过度调色

实际使用时，锐化前最好先降噪（如高斯滤波），锐化强度不要太大。

总结

空间域方法直接操作像素，更符合直觉，但也有计算量大、难以分析全局特性等局限。

方法	用途	核心思想	局限
Lucy-Richardson	图像去模糊	迭代最大似然估计	迭代次数难选择，噪声放大
TV 正则化	图像去噪/去模糊	最小化总变分（L1 梯度）	阶梯效应，纹理丢失
双边滤波	边缘保持降噪	空间权重 + 颜色权重	计算量大
引导滤波	边缘保持滤波/抠图	局部线性模型	参数调节敏感
拉普拉斯锐化	边缘增强	二阶导数检测边缘	噪声放大，过冲下冲