在计算机视觉与图像处理的浩瀚领域中,我们通常面临两类截然不同的任务:
本书专注于底层视觉。在这个领域,我们不仅是把图像当作像素矩阵,更是将其视为物理信号的离散采样。我们的目标是通过数学建模,逆转物理退化过程。
如果你想理解如何把模糊的监控画面变清晰,或者想知道 MRI(核磁共振)是如何从一堆毫无意义的频率数据中重建出大脑图像的,那么欢迎来到变分与优化的世界。在这里,我们将物理定律翻译成目标函数,将图像质量的追求转化为数学约束。
本章学习目标:
几乎所有的图像复原问题,都可以被浓缩进一个极其简洁的线性方程组中:
\[y = A(x) + n\]这个公式虽然简单,但它是连接现实世界与数学世界的桥梁。让我们逐一解剖其中的变量。
根据算子 $A$ 的不同构造,我们定义了不同的图像处理任务。理解 $A$ 的结构是建模的第一步。
| 任务 (Task) | 算子 $A$ 的性质 | 物理意义 | 矩阵形态 (直观) |
|---|---|---|---|
| 去噪 (Denoising) | 单位矩阵 $I$ | 理想成像,仅受噪声干扰 | 对角线为1 |
| 去模糊 (Deblurring) | 卷积矩阵 $H$ | 镜头失焦、物体运动、大气湍流 | 托普利兹 (Toeplitz) 或循环矩阵 |
| 超分辨率 (Super-Resolution) | 降采样 $D$ $\times$ 模糊 $H$ | 传感器像素密度不足 | 扁矩阵 (Fat Matrix),行数少 |
| 修复/补全 (Inpainting) | 掩膜矩阵 $M$ | 遮挡、坏点、文字水印覆盖 | 从 $I$ 中抽掉若干行 |
| 压缩感知 (Compressed Sensing) | 随机投影 $\Phi$ | 单像素相机、加速成像 | 随机高斯/伯努利矩阵 |
| 断层重建 (CT/MRI) | 变换域采样 $\mathcal{F}$ | 测量的是投影或频率 | 傅里叶/Radon 变换矩阵 |
Rule of Thumb (工程法则) 在实际代码中,永远不要显式生成矩阵 $A$! 对于一张 $1024 \times 1024$ 的图像,$N \approx 10^6$。矩阵 $A$ 的大小是 $10^6 \times 10^6$,存储它需要数 TB 的内存。 做法:我们实现两个函数(或类方法):
A(x):模拟前向过程。At(y):模拟 $A^T$(伴随算子/转置矩阵)的过程。 所有的优化算法只需要这两个函数即可运行。
[ 真实世界 x ] [ 算子 A 的作用 ] [ 观测 y ]
+---------+ (Convolution) +---------+
| Clear | Shift & Sum | Blurry |
| Image | --------> Kernel K --------> | Image |
+---------+ +---------+
| ^
| |
v |
(Vectorized) (Add Noise)
[ x_1 ] |
[ x_2 ] x [ A (Huge Matrix) ] = [ y_clean ] + n
[ ... ] [ (Implicitly structured) ]
[ x_N ]
既然 $y = Ax$,那为什么不直接计算 $x = A^{-1}y$ 呢?这引出了逆问题的核心困难——病态性。 法国数学家 Hadamard 定义了“适定问题”(Well-posed problem)必须满足三个条件:
图像逆问题通常违反这三条中的至少一条,尤其是第三条。
假设 $A$ 可逆,但它通常是一个“平滑”算子(如模糊)。这意味着 $A$ 的奇异值(Singular Values)衰减非常快,包含很多接近于 0 的微小值。 根据 SVD 分析,$A^{-1}$ 的奇异值是 $A$ 的倒数。 \(\hat{x} = A^{-1}(y_{true} + n) = x_{true} + \underbrace{A^{-1}n}_{\text{Explosion!}}\) 如果 $A$ 有一个特征值 $\sigma_i = 10^{-6}$,那么对应的噪声分量将被放大 $10^6$ 倍!恢复出来的图像将被巨大的噪声淹没,完全不可用。
在超分辨率或修复问题中,$A$ 是不可逆的“扁”矩阵(欠定系统,Under-determined)。这意味着有无穷多个 $x$ 都能满足 $Ax = y$。 例如,填补缺失的像素:我们可以填黑色,也可以填白色,方程都成立。但只有一种填法看起来像“真实的图像”。
为了从无数可能的解中通过“选美”挑出最合理的一个,或者为了压制噪声的爆炸,我们需要引入正则化(Regularization)。
我们利用贝叶斯定理: \(P(x|y) = \frac{P(y|x) P(x)}{P(y)}\) 我们想要最大化后验概率(MAP): \(\hat{x}_{MAP} = \arg\max_x P(y|x) P(x)\) 对两边取负对数(Negative Log-Likelihood),问题转化为最小化能量函数: \(\min_x \underbrace{-\log P(y|x)}_{\text{数据保真项 } D(Ax,y)} + \underbrace{-\log P(x)}_{\text{正则项 } \lambda R(x)}\)
这就构成了本书通用的变分模型(Variational Model):
\[\min_x \quad D(Ax, y) + \lambda R(x)\]其中 $\lambda$ 是平衡参数,控制了我们多大程度上相信观测数据,多大程度上相信先验知识。
这一项由噪声的物理统计特性决定。
如果 $n \sim \mathcal{N}(0, \sigma^2I)$,则 $P(y|x) \propto \exp(-|Ax-y|_2^2 / 2\sigma^2)$。 \(D(Ax, y) = \frac{1}{2} \|Ax - y\|_2^2\)
如果噪声服从拉普拉斯分布(尾巴比高斯长),或者包含椒盐噪声。 \(D(Ax, y) = \|Ax - y\|_1 = \sum_i |(Ax)_i - y_i|\)
常见于低光照成像、天文摄影、显微成像(光子计数过程)。 \(D(Ax, y) = \langle Ax, 1 \rangle - \langle y, \log(Ax) \rangle\)
这一项决定了我们恢复出的图像长什么样。这是过去几十年图像处理研究的核心。
\(R(x) = \| \Gamma x \|_2^2\) 通常 $\Gamma$ 是梯度算子。
\(R(x) = \| \nabla x \|_1 = \sum_{i,j} \sqrt{(\partial_x x)_{ij}^2 + (\partial_y x)_{ij}^2}\)
\(R(x) = \| \Psi x \|_1\) 其中 $\Psi$ 是某种变换(小波 Wavelet、DCT、Curvelet)。
\(R(X) = \| X \|_* = \sum \sigma_i(X)\) 针对矩阵形式的图像或视频数据。
一旦模型确定(即选定了 $D$ 和 $R$),剩下的就是纯粹的数学计算问题:求最小值。本书后续章节将详细介绍以下几类算法:
[ 优化方法族谱 ]
1. 光滑无约束 (Smooth Unconstrained)
|-- 梯度下降 (Gradient Descent) -> 也就是现在的深度学习训练基础
|-- 牛顿法 / 拟牛顿法 (L-BFGS) -> 收敛快,但内存消耗大
|-- 共轭梯度 (Conjugate Gradient)-> 求解线性方程组 (Ax=b) 的神器
2. 近端方法 (Proximal Methods) -> 处理 R(x) 不可导的情况 (如 L1, TV)
|-- ISTA (Iterative Shrinkage-Thresholding)
|-- FISTA (Fast ISTA) -> 加速版,图像处理标配
3. 算子分裂 / 增广拉格朗日 (Splitting / ALM) -> 处理复杂约束
|-- ADMM (交替方向乘子法) -> 极其灵活,万能胶水
|-- Split Bregman -> 专门针对 L1/TV 问题的高效解法
|-- Primal-Dual (Chambolle-Pock) -> 数学优美,处理鞍点问题
4. 非凸优化 (Non-convex) -> 更好的先验,更难的求解
|-- IRLS (迭代重加权最小二乘)
|-- MM (Majorization-Minimization)
选择指南:
Q1.1 (算子的理解): 对于“去模糊”任务,如果模糊核是 $3 \times 3$ 的均匀模糊(每个元素都是 $1/9$),对于 $100 \times 100$ 的图像。
Q1.2 (L2 vs L1): 考虑一维信号恢复问题。观测值 $y=3$,模型 $x$。 情况 A: $D(x) = (x-3)^2$。 情况 B: $D(x) = |x-3|$。 如果现在增加一个离群观测 $y_{outlier} = 100$,两个观测一起考虑。 A 的目标: $\min_x (x-3)^2 + (x-100)^2$ B 的目标: $\min_x |x-3| + |x-100|$ 请分别求出最优解,并说明哪个对离群值更鲁棒。
Q1.3 (正则化参数): 在超分辨率任务中,如果 $\lambda \to \infty$,恢复出来的图像 $x$ 会变成什么样?(假设 $R(x) = |\nabla x|_2^2$)。
Q1.4 (贝叶斯与 MAP): 证明:如果先验 $P(x) \propto e^{-\lambda |x|1}$(拉普拉斯分布),似然 $P(y|x) \propto e^{-\frac{1}{2\sigma^2}|Ax-y|_2^2}$(高斯分布)。推导出的 MAP 估计等价于 Lasso 问题($\ell_2$ 数据项 + $\ell_1$ 正则)。并给出 $\lambda{opt}$ 与噪声方差 $\sigma^2$ 以及拉普拉斯参数 $b$ 的关系。
Q1.5 (逆犯罪 Inverse Crime): 在学术论文的仿真实验中,如果你用模型 $y = A_{sim}x + n$ 生成数据,然后用完全相同的 $A_{sim}$ 去进行恢复 $\hat{x} = \arg\min |A_{sim}x - y|^2 + R(x)$,你的结果通常会好得出奇。 但在真实世界中,真实的退化过程 $A_{real}$ 永远无法被完美的数学模型 $A_{sim}$ 描述(例如离散化误差、未建模的光学畸变)。 这种现象被称为“逆犯罪”。如何设计实验避免这个问题?
img[row, col]。edgetaper (Matlab) 或在变分模型中引入边界处理。