线性代数（十七）：计算机视觉中的线性代数——从像素到三维重建

计算机视觉就是教机器学会“看”。有意思的是，这个领域几乎所有问题最后都能归结为线性代数：图像可以看作矩阵，几何变换是矩阵相乘，相机成像用一个 $3 \times 4$ 的投影矩阵表示，两视图几何关系简化为方程 $\mathbf{x}_2^\top \mathbf{F}\, \mathbf{x}_1 = 0$ ，三维重建则是一个稀疏线性最小二乘问题。一旦从这个角度理解 CV，原本看似杂乱无章的算法其实只是反复应用了少数几个线性代数的核心思想。

你会学到什么
图像如何表示为矩阵、张量和高维向量
旋转、缩放、错切和平移如何通过矩阵乘法实现，齐次坐标统一了这些操作
为什么透视变换需要单应矩阵 $\mathbf{H}$ ，以及如何从点对应关系中求解它
针孔相机模型与投影矩阵 $\mathbf{P} = \mathbf{K}[\mathbf{R}\,|\,\mathbf{t}]$
对极几何、基础矩阵 / 本质矩阵、三角测量
基于 SVD 的图像压缩与 Eckart-Young 定理
卷积的本质是矩阵乘法； Harris 角点检测和光流问题归结为 $2 \times 2$ 线性系统

前置知识
线性变换（第 3 章）、特征分解（第 6 章）、 SVD （第 9 章）

图像：矩阵、张量与向量#

从像素到矩阵#

\mathbf{I} \in \mathbb{R}^{H \times W}, \qquad I_{ij} \in [0, 1].

到这里，故事已经讲清楚了。后续所有对图像的操作，都可以看作是对这个矩阵进行线性代数变换。

彩色图像：三通道张量#

\mathcal{I} \in \mathbb{R}^{H \times W \times 3}.

Y = 0.299\,R + 0.587\,G + 0.114\,B.

1
2
3
4
import numpy as np, cv2
img = cv2.imread('photo.jpg')           # (H, W, 3)，注意 BGR 顺序
B, G, R = cv2.split(img.astype(float) / 255)
gray = 0.299 * R + 0.587 * G + 0.114 * B

图像作为高维向量#

\mathrm{vec}(\mathbf{I}) \in \mathbb{R}^{HW}.

\cos\theta = \frac{\mathbf{a}^\top \mathbf{b}}{\lVert\mathbf{a}\rVert\,\lVert\mathbf{b}\rVert}.

这种方法广泛应用于图像检索和人脸识别。但维度太高会导致几何邻近关系失去意义，这正是我们需要用 PCA 或 SVD 将其压缩到低维子空间的原因。

几何变换 = 矩阵乘法#

为什么用矩阵#

几何变换会把每个像素坐标 $$(x, y)$$ 映射到一个新位置。用矩阵表示这种变换有两个巨大的优势：

组合就是相乘。先应用 $\mathbf{A}$ ，再应用 $\mathbf{B}$ ，结果直接是 $\mathbf{B}\mathbf{A}$ 。
逆变换就是求逆矩阵。撤销操作时，完全不用重新推导公式。

旋转、缩放、错切#

\mathbf{R}(\theta) = \begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{bmatrix}.

旋转矩阵有三个特别让人省心的性质：

它是正交矩阵（ $\mathbf{R}^\top\mathbf{R} = \mathbf{I}$ ），所以逆矩阵就是转置， $\mathbf{R}^{-1} = \mathbf{R}^\top$ ；
行列式 $\det \mathbf{R} = 1$ ，这意味着它既保面积又保朝向；
旋转角可以直接相加， $\mathbf{R}(\alpha)\mathbf{R}(\beta) = \mathbf{R}(\alpha + \beta)$ 。

\mathbf{S} = \begin{bmatrix}s_x & 0 \\ 0 & s_y\end{bmatrix}, \qquad \mathbf{H}_k = \begin{bmatrix}1 & k \\ 0 & 1\end{bmatrix}.

顺序很重要#

变换组合是从右往左读的： $\mathbf{T} = \mathbf{R}\,\mathbf{S}$ 意思是“先缩放，再旋转”。如果调换顺序，结果通常不一样。唯一例外是等比缩放，它和旋转可以交换。

1
2
3
4
5
6
def rotation_matrix(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c, -s], [s, c]])

S = np.diag([2.0, 2.0])
T = rotation_matrix(np.pi / 4) @ S        # 先 S 再 R

齐次坐标：让平移变成线性#

平移的麻烦#

旋转、缩放和错切都不会动原点，所以它们是线性变换，公式是 $\mathbf{y} = \mathbf{A}\mathbf{x}$ 。但平移 $\mathbf{y} = \mathbf{x} + \mathbf{t}$ 不一样，它把原点移到别的地方，破坏了线性性质。我希望所有几何变换都能用矩阵乘法统一描述。

升维的小技巧#

\begin{bmatrix} x \\ y \end{bmatrix} \;\longrightarrow\; \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}.

\begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & t_x \\ 0 & 1 & t_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}.

\mathbf{M} = \begin{bmatrix} a_{11} & a_{12} & t_x \\ a_{21} & a_{22} & t_y \\ 0 & 0 & 1 \end{bmatrix}.

左上角的 $2 \times 2$ 块负责线性部分，右边一列是平移，最后一行 $$(0, 0, 1)$$ 标记着这是个仿射变换。

绕任意点旋转#

\mathbf{M} = \mathbf{T}(c_x, c_y)\,\mathbf{R}(\theta)\,\mathbf{T}(-c_x, -c_y).

1
2
3
4
5
6
7
8
9
def affine_matrix(rotation_deg, scale, translation):
    theta = np.radians(rotation_deg)
    c, s = np.cos(theta), np.sin(theta)
    sx, sy = scale; tx, ty = translation
    return np.array([[sx * c, -sy * s, tx],
                     [sx * s,  sy * c, ty],
                     [0,       0,      1]])

M = affine_matrix(30, (1.5, 1.5), (100, 50))

透视变换与单应矩阵#

为什么仿射变换不够用#

仿射变换能保持平行线不变，但现实世界不是这样。比如铁轨在远处会汇聚到地平线。拍一张平面物体的照片时，如果角度倾斜，照片里的平行线也会汇聚。要解决这个问题，需要一种更灵活的变换。

$3 \times 3$ 单应矩阵#

\begin{bmatrix} x' \\ y' \\ w' \end{bmatrix} = \mathbf{H} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}, \qquad u = x'/w', \quad v = y'/w'.

它和仿射矩阵的关键区别在于最后一行。仿射矩阵的最后一行固定为 $$(0, 0, 1)$$ ，而单应矩阵允许非零值。这些非零值引入了透视除法 $$w'$$ ，让平行线看起来汇聚在一起。

\mathbf{A}\,\mathbf{h} = \mathbf{0}

然后用 SVD 求解：取最小奇异值对应的右奇异向量作为 $\mathbf{h}$ 。这就是直接线性变换 (DLT) 算法，也是本章所有投影几何估计算法的原型。

1
2
3
src = np.float32([[100,100],[200,100],[200,200],[100,200]])
dst = np.float32([[120, 90],[220,110],[210,220],[ 90,210]])
H, _ = cv2.findHomography(src, dst, cv2.RANSAC, 5.0)

应用场景
全景拼接：相机只旋转时，相邻图像之间可以用一个 $\mathbf{H}$ 表示。
文档扫描：把斜拍的纸面矫正成正面视角。
增强现实：在检测到的平面标志上放置虚拟物体。
自动驾驶：将车前视角转换为俯视 BEV 图像。

针孔相机与投影矩阵#

针孔几何#

u = f\,\frac{X}{Z}, \qquad v = f\,\frac{Y}{Z}.

这里除以深度 $$Z$$ 是透视效果的来源：远处的物体在图像上显得更小。

内参：从物理单位到像素#

\mathbf{K} = \begin{bmatrix} f_x & s & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix},

其中 $$f_x, f_y$$ 是以像素为单位的焦距， $$(c_x, c_y)$$ 是主点的像素坐标， $$s$$ 是传感器的倾斜量（现代相机基本为零）。

外参与完整投影#

\lambda \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \mathbf{K}\,[\mathbf{R}\,|\,\mathbf{t}] \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix} = \mathbf{P}\,\mathbf{X}_w,

其中 $\mathbf{P} \in \mathbb{R}^{3 \times 4}$ 是投影矩阵， $\lambda$ 吸收了透视除法引入的比例因子。投影矩阵包含 11 个参数： 5 个内参加 6 个位姿参数（3 个旋转角和 3 个平移分量）。

1
2
3
4
5
6
7
def project(X_world, K, R, t):
    Xc = X_world @ R.T + t                 # 世界坐标 -> 相机坐标
    x  = Xc[:, 0] / Xc[:, 2]                # 透视除法
    y  = Xc[:, 1] / Xc[:, 2]
    u  = K[0, 0] * x + K[0, 2]
    v  = K[1, 1] * y + K[1, 2]
    return np.stack([u, v], axis=1)

张正友标定法#

用一组已知尺寸的棋盘格图像估计 $\mathbf{K}$ 和镜头畸变系数，这就是相机标定。张正友在 1999 年提出的算法分为三步：每张棋盘图估计一个单应矩阵，利用旋转矩阵列向量正交的约束提取关于 $\mathbf{K}^{-\top}\mathbf{K}^{-1}$ 的线性方程，最后用非线性优化细化结果。 OpenCV 的 calibrateCamera 函数一行代码就能完成整个流程。

两视图：对极几何#

对极约束#

\mathbf{x}_2^\top\,\mathbf{F}\,\mathbf{x}_1 = 0,

其中 $\mathbf{F}$ 是 $3 \times 3$ 的基础矩阵。它只和两台相机的相对位姿有关，与场景无关。如果知道第一张图像中的一个点 $\mathbf{x}_1$ ，那么它在第二张图像中的对应点一定落在一条直线上。这条直线叫对极线，计算公式是 $\boldsymbol{\ell}_2 = \mathbf{F}\,\mathbf{x}_1$ 。这样一来，立体匹配的搜索范围就从二维降到了一维。

本质矩阵 vs 基础矩阵#

\hat{\mathbf{x}}_2^\top\,\mathbf{E}\,\hat{\mathbf{x}}_1 = 0, \qquad \mathbf{F} = \mathbf{K}_2^{-\top}\,\mathbf{E}\,\mathbf{K}_1^{-1}.

本质矩阵有特殊形式 $\mathbf{E} = [\mathbf{t}]_\times \mathbf{R}$ ，其中 $[\mathbf{t}]_\times$ 是平移向量的反对称矩阵。对 $\mathbf{E}$ 做 SVD 分解，可以恢复出两台相机的相对位姿 $(\mathbf{R}, \mathbf{t})$ 。不过， $\mathbf{t}$ 的绝对长度永远无法确定，这是单目视觉固有的尺度不确定性。

1
2
3
F, _ = cv2.findFundamentalMat(pts1, pts2, cv2.FM_RANSAC, 3.0, 0.99)
E, _ = cv2.findEssentialMat(pts1, pts2, K, cv2.RANSAC, 0.999, 1.0)
_, R, t, _ = cv2.recoverPose(E, pts1, pts2, K)

三维重建：三角测量、 SfM 和 Bundle Adjustment#

三角测量是一个线性系统#

已知两个投影矩阵 $\mathbf{P}_1, \mathbf{P}_2$ ，以及一对对应点 $\mathbf{x}_1 \leftrightarrow \mathbf{x}_2$ 。每幅图像的方程 $\lambda_i \mathbf{x}_i = \mathbf{P}_i \mathbf{X}$ 消去未知深度 $\lambda_i$ 后，会为 3D 点 $\mathbf{X}$ 提供两个线性约束。把这些约束堆叠起来，就得到一个 $4 \times 4$ 的线性系统 $\mathbf{A}\mathbf{X} = \mathbf{0}$ 。用 SVD 求解时， $\mathbf{X}$ 是最小奇异值对应的右奇异向量，经过齐次除法后就能恢复出 3D 坐标。

1
2
X4 = cv2.triangulatePoints(P1, P2, pts1.T, pts2.T)
X3 = (X4[:3] / X4[3]).T

运动恢复结构#

SfM 把三角测量扩展到多视图场景。标准的增量式流程如下：

特征检测与匹配：在所有重叠图像对之间用 SIFT 或 ORB 检测和匹配特征点。
初始化：选两张视差适中的图像，用 5 点法加 RANSAC 估计本质矩阵 $\mathbf{E}$ ，恢复相对位姿，再通过三角测量生成初始点云。
逐步扩展：每次加入一张新图像，用 PnP（Perspective-n-Point）求解新相机的位姿，然后三角测量出新可见的 3D 点。
全局优化：调用 Bundle Adjustment 对所有变量进行联合优化。

光束法平差#

\min_{\{\mathbf{P}_i\},\,\{\mathbf{X}_j\}}\;\sum_{(i,j)\in\mathcal{V}} \rho\!\left(\lVert \mathbf{x}_{ij} - \pi(\mathbf{P}_i, \mathbf{X}_j)\rVert^2\right),

其中 $\pi$ 是投影函数， $\rho$ 是 Huber 鲁棒核函数，用来处理外点。这是一个超大规模的非线性最小二乘问题——现代重建中变量数可能达到百万级。但雅可比矩阵是块稀疏的：每个观测只涉及一个相机和一个 3D 点。 Levenberg-Marquardt 算法结合 Schur 补利用了这种稀疏性，让普通笔记本也能高效求解。

SLAM：实时中的线性代数#

SLAM 问题#

机器人在未知环境中移动，一边读取传感器数据，一边同时估计自己的轨迹和环境地图。这就是 SLAM（Simultaneous Localization And Mapping）。现代 SLAM 系统本质上是在线版的 Bundle Adjustment，其中包含两个非常有趣的线性代数部分。

用李群表示旋转#

\mathbf{T} = \begin{bmatrix}\mathbf{R} & \mathbf{t} \\ \mathbf{0}^\top & 1\end{bmatrix} \in SE(3)

其中 $\mathbf{R} \in SO(3)$ 。直接把 $\mathbf{R}$ 当成 9 个自由变量优化很麻烦，因为必须时刻保证 $\mathbf{R}^\top\mathbf{R} = \mathbf{I}$ 。李代数 $\mathfrak{se}(3) \cong \mathbb{R}^6$ 是一个无约束的向量空间，指数映射 $\mathbf{T} = \exp(\boldsymbol{\xi}^\wedge)$ 在两者之间切换。我们在 $\mathbb{R}^6$ 中做梯度下降，再映射回 $$SE(3)$$ —— 干净利落，没有约束。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def skew(v):
    return np.array([[0, -v[2], v[1]],
                     [v[2], 0, -v[0]],
                     [-v[1], v[0], 0]])

def so3_exp(phi):                    # 罗德里格斯公式
    theta = np.linalg.norm(phi)
    if theta < 1e-10:
        return np.eye(3)
    a = phi / theta; A = skew(a)
    return np.eye(3) + np.sin(theta) * A + (1 - np.cos(theta)) * (A @ A)

位姿图优化#

\min\;\sum_k \lVert \mathbf{e}_k\rVert^2_{\boldsymbol{\Omega}_k}.

\mathbf{H}\,\Delta\mathbf{x} = -\mathbf{b},

其中 $\mathbf{H} = \mathbf{J}^\top \boldsymbol{\Omega}\,\mathbf{J}$ 继承了图的稀疏结构。对 $\mathbf{H}$ 做稀疏 Cholesky 分解，几乎是所有现代 SLAM 系统的核心步骤。

图像滤波 = 矩阵乘法#

卷积是超大、稀疏、有结构的线性映射#

二维卷积 $\mathbf{G} = \mathbf{I} \ast \mathbf{K}$ ，把图像拉成向量后可以写成矩阵乘法 $\mathbf{g} = \mathbf{T}\mathbf{i}$ 。这里的 $\mathbf{T}$ 是一个由卷积核生成的双重块 Toeplitz 矩阵。我们不会真的去存储 $\mathbf{T}$ ，因为它需要 $$(HW)^2$$ 个元素。但它的存在让我们可以用线性代数工具分析卷积。比如， $\mathbf{T}$ 的特征值正好是卷积核的离散傅里叶变换。

必须记住的三个卷积核#

\mathbf{K}_{\text{blur}} = \tfrac{1}{9}\begin{bmatrix}1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1\end{bmatrix}

\mathbf{K}_{\text{edge}} = \begin{bmatrix}0 & -1 & 0\\ -1 & 4 & -1\\ 0 & -1 & 0\end{bmatrix}

\mathbf{K}_{\text{sharp}} = \begin{bmatrix}0 & -1 & 0\\ -1 & 5 & -1\\ 0 & -1 & 0\end{bmatrix} = \mathbf{I} + \mathbf{K}_{\text{edge}}.

卷积定理#

\mathcal{F}[\mathbf{I} \ast \mathbf{K}] = \mathcal{F}[\mathbf{I}] \cdot \mathcal{F}[\mathbf{K}].

设计滤波器其实就是设计它的频率响应：低通用来去噪，高通用来提取边缘，带通用来捕捉纹理。

两个看似不同、本质却是同一个 $2 \times 2$ 特征值问题的算法#

Harris 角点：结构张量#

\mathbf{M} = \sum_{(x, y) \in W} w(x, y) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix}.

这是一个 $2 \times 2$ 的对称半正定矩阵。它的特征值 $\lambda_1 \ge \lambda_2 \ge 0$ 描述了窗口内图像沿两个主方向的变化情况：

特征值	含义
都很小	平坦区域
一大一小	边缘
都很大	角点

R = \det(\mathbf{M}) - k\,\mathrm{trace}(\mathbf{M})^2 = \lambda_1\lambda_2 - k(\lambda_1 + \lambda_2)^2,

只有当两个特征值都较大时， $$R$$ 才会显著增大，从而准确刻画角点。

1
2
3
4
5
6
7
8
def harris(img, k=0.04, win=2, ksize=3):
    Ix = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=ksize)
    Iy = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=ksize)
    Sxx = cv2.GaussianBlur(Ix * Ix, (2 * win + 1,) * 2, 0)
    Syy = cv2.GaussianBlur(Iy * Iy, (2 * win + 1,) * 2, 0)
    Sxy = cv2.GaussianBlur(Ix * Iy, (2 * win + 1,) * 2, 0)
    detM, trM = Sxx * Syy - Sxy ** 2, Sxx + Syy
    return detM - k * trM ** 2

光流：同一个矩阵再次登场#

I_x\,u + I_y\,v + I_t = 0.

\underbrace{\begin{bmatrix} \sum I_x^2 & \sum I_x I_y \\ \sum I_x I_y & \sum I_y^2 \end{bmatrix}}_{\text{正是结构张量 } \mathbf{M}} \begin{bmatrix} u \\ v \end{bmatrix} = -\begin{bmatrix} \sum I_x I_t \\ \sum I_y I_t \end{bmatrix}.

这个方程组是否能解得好，取决于 $\mathbf{M}$ 的条件数——换句话说，取决于当前区域是否为角点。同样的特征值分析告诉 Harris 哪里有角点，也告诉 Lucas-Kanade 它的光流估计是否可靠。

SVD 图像压缩#

\mathbf{I} = \sum_{i=1}^{r} \sigma_i\,\mathbf{u}_i\,\mathbf{v}_i^\top, \qquad \sigma_1 \ge \sigma_2 \ge \cdots \ge 0.

Eckart-Young 定理（第 9 章）指出，取前 $$k$$ 项的截断结果 $\mathbf{I}_k$ 是在 Frobenius 范数和谱范数意义下的最佳秩 $$k$$ 近似。自然图像的奇异值衰减得很快，大部分视觉信息集中在前几十个分量中。因此，即使 $$k$$ 很小，也能生成可辨认的图像。

存储需求从 $$HW$$ 减少到 $$k(H + W + 1)$$ 个数，压缩比大约是 $$k(H + W) / (HW)$$ 。 SVD 在自然图像上的效率远不如 JPEG （JPEG 利用了 DCT 系数的感知冗余），但它思路清晰，是许多视觉任务中低秩近似的金标准：比如人脸识别（Eigenfaces）、背景建模（RPCA）和图像去噪。

现代深度视觉中的线性代数#

卷积通过 im2col 转化为矩阵乘法#

GPU 在处理稠密矩阵乘法时性能极强。im2col 方法把卷积操作转化为一次大规模的 gemm：将每个感受野区域提取为一列，把所有卷积核展开成行，然后一次性完成矩阵乘法，最后再将结果重塑为特征图。从本质上讲， CNN 就是一系列矩阵乘法和逐元素非线性变换的交替操作。

自注意力机制的核心是矩阵乘法#

\mathrm{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \mathrm{softmax}\!\left(\frac{\mathbf{Q}\mathbf{K}^\top}{\sqrt{d_k}}\right)\mathbf{V},

其中 $\mathbf{Q}, \mathbf{K}, \mathbf{V}$ 是对序列进行线性投影的结果。 Vision Transformer (ViT) 把图像分割成小块作为 token 输入，整个前向传播过程就是一堆矩阵乘法堆叠起来的操作。

推理时 BatchNorm 是仿射变换#

推理阶段， BatchNorm 使用存储的运行统计量和学到的 $(\gamma, \beta)$ 参数，执行逐元素仿射变换 $y = \gamma\,\hat{x} + \beta$ 。在部署时，这个操作可以合并到前一层卷积中，从而完全从计算图中移除。