Series · Linear Algebra · Chapter 17

计算机视觉中的线性代数 -- 从像素到三维重建

计算机视觉几乎完全建立在线性代数之上:图像是矩阵,几何变换是矩阵乘法,相机成像是投影矩阵,三维重建是求解线性方程组。本章带你串起齐次坐标、单应矩阵、相机标定、对极几何、SfM 与 SLAM 的完整脉络。

计算机视觉的核心任务是让机器"看懂"图像。让人惊讶的是,整个学科几乎都建立在线性代数之上:图像本身就是矩阵,几何变换是矩阵乘法,相机成像是一个 $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 章)


1 图像:矩阵、张量、向量

从像素到矩阵

相机传感器是一格格"光桶"。每个桶在曝光期间收集到的光子数,归一化到固定区间后,就是图像在相应位置的像素值。一张 $H \times W$ 的灰度图像,就是一个矩阵:

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

故事到这里其实就讲完了。后面所有图像处理操作,都只是对这个矩阵做某种线性代数运算。

灰度图是矩阵,彩色图是 H x W x 3 的张量。

彩色图:三通道张量

彩色相机通过红绿蓝三种滤色片,每个像素记录三个亮度值,得到一个三维数组:

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

NumPy 用 (H, W, 3) 形状存储;PyTorch 习惯 (3, H, W)。把彩色图转成灰度时,通常按人眼亮度感知加权:

$$ 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

把图像看成高维向量

很多机器学习算法(PCA、SVM、全连接层)只接受向量输入。我们把 $H \times W$ 的图像按列(或按行)拉平:

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

一张 $640 \times 480$ 的图,就成为 30 多万维空间里的一个点,每个坐标轴对应一个像素位置。两张图片的相似度可以直接用余弦相似度度量:

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

图像检索、人脸识别都是这么做的。但维度太高也意味着 “邻近” 的几何意义不太可靠 – 这正是后面要用 PCA / SVD 把它压到一个真正有视觉意义的低维子空间的原因。


2 几何变换 = 矩阵乘法

为什么用矩阵

一个几何变换把每个像素坐标 $(x, y)$ 映射到新位置。用矩阵描述变换有两个无可替代的好处:

  • 组合就是相乘:先做 $\mathbf{A}$ 再做 $\mathbf{B}$,整体效果就是 $\mathbf{B}\mathbf{A}$;
  • 逆变换就是矩阵求逆:撤销操作不用重新推公式。

旋转、缩放、错切

绕原点逆时针旋转角度 $\theta$ 的矩阵是:

$$ \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

3 齐次坐标:把平移也变成线性

平移的尴尬

旋转、缩放、错切都让原点保持不动,所以是线性变换 $\mathbf{y} = \mathbf{A}\mathbf{x}$。但平移 $\mathbf{y} = \mathbf{x} + \mathbf{t}$ 把原点搬走了,不再线性 – 这意味着平移不能写成普通的矩阵乘法。我们希望所有几何变换都用同一种语言描述,怎么办?

升一维的小魔法

把每个 2D 点末尾加一个 $1$:

$$ \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}. $$

几何上,我们把 2D 平面"嵌入"到 3D 空间的 $z = 1$ 切片里。原本看似平移的操作,其实是 3D 里的一个错切;而错切是线性的。所有仿射变换都装进同一个 $3 \times 3$ 矩阵:

$$ \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)$ 标记着"这是个仿射变换"。

旋转、缩放、平移在齐次坐标下都是 3x3 矩阵。

绕任意点旋转

最常见的扩展:绕图像中心 $(c_x, c_y)$ 而不是绕原点旋转。三步搞定:

$$ \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
10
11
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))
# OpenCV 用的是 2x3 的前两行:
# warped = cv2.warpAffine(img, M[:2, :], (W, H))

4 透视变换与单应矩阵

仿射变换不够用

仿射变换保持平行线平行 – 但现实世界里,平行的铁轨在远方明明会汇聚。同样地,斜着拍一张平面物体(如一张纸),照片上的矩形会变成不规则四边形。要把这种"近大远小"翻译成数学,需要更一般的投影变换

$3 \times 3$ 的单应矩阵

单应矩阵 (homography) 是一个一般的 $3 \times 3$ 可逆矩阵 $\mathbf{H}$,作用在齐次坐标上:

$$ \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'$,把平行线弯成了汇聚的样子。

由于矩阵整体只在尺度意义下确定,单应矩阵有 $9 - 1 = 8$ 个自由度。每对对应点 $(x_i, y_i) \leftrightarrow (u_i, v_i)$ 提供两个约束方程,因此理论上四对点就够了。实际中我们用更多噪声点,把所有约束堆成

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

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

从一张倾斜拍摄的平面照片,通过单应矩阵 H 矫正回正视图。

1
2
3
4
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)
# corrected = cv2.warpPerspective(img, H, (W, H))

它出现在哪些场景: 全景拼接(相机只旋转时,每对相邻图像之间就是一个 $\mathbf{H}$);文档扫描 App(把斜拍的纸面纠正成正视图);增强现实(把虚拟物体贴到检测出来的平面标志上);自动驾驶里把车前视图变成俯视 BEV 图。


5 针孔相机与投影矩阵

针孔几何

最简单的相机模型把镜头想象成一个无穷小的针孔位于原点,成像平面位于距离 $f$ 处(焦距)。相机坐标系下的 3D 点 $(X, Y, Z)$ 投影到像平面上:

$$ 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$ 是传感器的倾斜量(现代相机基本为零)。

外参与完整投影方程

世界坐标系一般不等于相机坐标系,二者之间是一个刚体运动:先旋转 $\mathbf{R}$,再平移 $\mathbf{t}$。把所有部分串起来:

$$ \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 个平移分量)。

针孔相机投影:3D 点经 P = K[R|t] 变成 2D 像素。

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 一行调用就能跑完整套流程。


6 两视图:对极几何

对极约束

同一个 3D 点经过两个相机分别拍到了 $\mathbf{x}_1$ 和 $\mathbf{x}_2$。这两个像点之间不是无关的 – 3D 点、两个相机光心、两个像点这五个点必然共面(这个面叫对极平面)。把"共面"这个事翻译成代数,就是

$$ \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}} = \mathbf{K}^{-1}\mathbf{x}$,相同的约束就变成了本质矩阵

$$ \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)

7 三维重建:三角测量、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 取最小奇异值对应的右奇异向量,再做齐次除法就得到 3D 坐标。

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

Structure from Motion

SfM 把三角测量推广到多视图。常见的增量式流程:

  1. 特征检测与匹配:所有图像两两之间用 SIFT / ORB 匹配;
  2. 初始化:选两张视差适中的图,用 5 点法 + RANSAC 估计 $\mathbf{E}$,恢复相对位姿,三角测量出初始点云;
  3. 逐张加入:每加入一张新图,用 PnP 求解新相机位姿,再三角测量出能新看到的 3D 点;
  4. 全局精化:调用 Bundle Adjustment 同时优化所有变量。

Bundle Adjustment

多视图几何里的旗舰优化问题:联合优化所有相机参数 $\{\mathbf{P}_i\}$ 和所有 3D 点 $\{\mathbf{X}_j\}$,使总重投影误差最小:

$$ \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 补就能利用这种稀疏性,让笔记本电脑也跑得动。


8 SLAM:实时跑的线性代数

SLAM 问题

机器人在未知环境中移动,一边读传感器数据,一边同时估计自己的轨迹和环境地图,这就是 SLAM(Simultaneous Localization And Mapping)。现代 SLAM 本质上就是在线版的 Bundle Adjustment,里面有两块特别有意思的线性代数。

用李群表示位姿

3D 刚体位姿生活在矩阵群

$$ \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)

位姿图优化

现代 SLAM 框架(g2o、GTSAM、Ceres)把问题建成因子图:节点是位姿和路标,边是带协方差 $\boldsymbol{\Omega}_k$ 的相对测量。优化目标:

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

用 Gauss-Newton 迭代求解,每一步都要解一个稀疏线性方程组:

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

其中 $\mathbf{H} = \mathbf{J}^\top \boldsymbol{\Omega}\,\mathbf{J}$ 继承了图的稀疏结构。在 $\mathbf{H}$ 上跑稀疏 Cholesky 分解,几乎是所有现代 SLAM 系统的内层循环。


9 图像滤波 = 矩阵乘法

三种经典 3x3 卷积核(边缘、模糊、锐化)的输入输出对比。

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

二维卷积 $\mathbf{G} = \mathbf{I} \ast \mathbf{K}$,把图像拉成向量后写成矩阵乘法 $\mathbf{g} = \mathbf{T}\mathbf{i}$,这里 $\mathbf{T}$ 是从卷积核构造出来的双重块 Toeplitz 矩阵。我们当然不会把 $\mathbf{T}$ 真的存下来 – 那要 $(HW)^2$ 个元素 – 但它的存在让我们可以用线性代数的语言分析卷积(事实上 $\mathbf{T}$ 的特征值正好是卷积核的离散傅里叶变换)。

三个该背下来的卷积核

均值核(低通,模糊)。在邻域内平均;权重之和为 1,所以不改变整体亮度。

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

拉普拉斯核(带通,边缘)。离散二阶导;权重之和为 0,平坦区直接归零,只有突变处才被点亮。

$$ \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}]. $$

设计滤波器,本质上就是设计它的频率响应:低通去噪、高通提边、带通抓纹理。


10 两个看似不同、其实是同一个 $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$ 描述了窗口内图像沿两个主方向的变化幅度:

特征值局部结构
都很小平坦区
一大一小边缘
都很大角点

Harris 巧妙地避开了显式求特征值,引入响应函数

$$ 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

光流:同一个矩阵又出现了

光流估计的是相邻两帧之间每个像素的位移 $(u, v)$。假设亮度守恒,把第一帧到第二帧的灰度变化做一阶泰勒展开,得到亮度恒定方程

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

一个标量方程,两个未知数 – 局部欠定(这就是著名的光圈问题)。Lucas-Kanade 的解法:假设 $(u, v)$ 在小窗口内是常数,对窗口内每个像素列一个方程,得到一个超定线性系统,用最小二乘求解:

$$ \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 它的光流估计可不可信。

光流:从亮度恒定方程局部最小二乘解出来的位移向量场。


11 SVD 图像压缩

任何 $H \times W$ 灰度图都可以做奇异值分解:

$$ \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)、图像去噪都源于同一个思想。

SVD 在三个不同秩下的低秩重建及奇异值谱。


12 现代深度视觉里的线性代数

im2col:把卷积写成矩阵乘法

GPU 在稠密矩阵乘法上跑得飞快。im2col 技巧把每次卷积重写成一次巨大的 gemm:把每个感受野展开成列向量,把所有卷积核展开成行向量,做一次矩阵乘法,再 reshape 回特征图。机制上看,CNN 就是一连串矩阵乘法 + 逐元素非线性。

自注意力也是矩阵乘法

Transformer 的自注意力模块就是

$$ \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$。这个变换可以在部署时与前一层卷积合并,从计算图里彻底消失。


练习题

基础题

  1. 写出对灰度图 $\mathbf{I}$ 做下列操作的矩阵表达式:(a) 上下翻转;(b) 左右翻转;(c) 亮度 $+50$;(d) 对比度 $\times 1.5$。
  2. 证明二维旋转矩阵:(a) $\det \mathbf{R} = 1$;(b) $\mathbf{R}^{-1} = \mathbf{R}^\top$;(c) $\mathbf{R}(\alpha)\mathbf{R}(\beta) = \mathbf{R}(\alpha + \beta)$。
  3. 用 $3 \times 3$ 齐次矩阵分别表示并合并:先平移 $(2, 3)$,再绕原点旋转 $45^\circ$,再缩放 $2$ 倍。

进阶题

  1. 为什么单应矩阵理论上需要 4 对点?写出 DLT 算法对应的 $\mathbf{A}\mathbf{h} = \mathbf{0}$ 系统。
  2. 证明本质矩阵可以分解为 $\mathbf{E} = [\mathbf{t}]_\times \mathbf{R}$。利用 SVD $\mathbf{E} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\top$ 恢复 $\mathbf{R}$ 与 $\mathbf{t}$(差一个符号)。
  3. 给定两个投影矩阵和一对像点,推导三角测量的线性方程组 $\mathbf{A}\mathbf{X} = \mathbf{0}$。

编程题

  1. 从零实现仿射变换:手工构造 $3 \times 3$ 矩阵,用反向映射 + 双线性插值完成图像变换。
  2. 实现八点法估计 $\mathbf{F}$,包含坐标归一化和 SVD 强制 rank-2 约束。
  3. 实现完整的 Harris 角点检测,含非极大值抑制。

应用题

  1. 设计一个全景拼接系统。为什么相机只做旋转时单应矩阵就够用?多张图累积下来的误差怎么处理?曝光不一致怎么处理?
  2. 实现一个简单的单目视觉里程计:提特征 → 帧间匹配 → 估 $\mathbf{E}$ → 恢复相对位姿 → 累加轨迹。讨论尺度漂移问题,以及如何检测跟踪失败。

本章小结

计算机视觉与线性代数密不可分。

表示:灰度图是矩阵,彩色图是三通道张量,拉平后是 $\mathbb{R}^{HW}$ 中的向量。

几何:旋转、缩放、错切是线性变换;齐次坐标把平移也提升成线性操作;透视则需要单应矩阵 – 它的"非平凡"之处就在最后一行。

相机与三维:针孔相机是 $3 \times 4$ 投影矩阵 $\mathbf{P} = \mathbf{K}[\mathbf{R}\,|\,\mathbf{t}]$;两视图几何浓缩在基础矩阵 $\mathbf{F}$ 里;三角测量、SfM、Bundle Adjustment 都是(稀疏)线性最小二乘。

状态估计:SLAM 在 $SE(3)$ 流形上做优化,李代数让旋转优化无约束化;位姿图优化的内层循环是稀疏 Cholesky。

滤波与特征:卷积是巨大的块 Toeplitz 乘法。Harris 与光流共享同一个 $2 \times 2$ 结构张量的特征值分析。

压缩与学习:SVD 给出图像的最优低秩近似;CNN、Transformer 都是矩阵乘法的塔;BatchNorm 在推理时退化为一个仿射映射。

把这套东西内化以后,再读 CV 的论文,你会觉得它们不再是一堆毫不相关的算法,而是同一套线性代数工具不断被重新组合的结果。


参考资料

  • Hartley, R., & Zisserman, A. Multiple View Geometry in Computer Vision (第 2 版). Cambridge University Press, 2004.
  • Szeliski, R. Computer Vision: Algorithms and Applications (第 2 版). Springer, 2022.
  • Forsyth, D. A., & Ponce, J. Computer Vision: A Modern Approach (第 2 版). Pearson, 2011.
  • Strang, G. Introduction to Linear Algebra (第 5 版). Wellesley-Cambridge Press, 2016.
  • Barfoot, T. D. State Estimation for Robotics. Cambridge University Press, 2017.
  • Zhang, Z. “A flexible new technique for camera calibration.” IEEE TPAMI 22(11), 2000.
  • Lucas, B., & Kanade, T. “An iterative image registration technique with an application to stereo vision.” IJCAI, 1981.
  • Harris, C., & Stephens, M. “A combined corner and edge detector.” Alvey Vision Conference, 1988.

系列导航

Liked this piece?

Follow on GitHub for the next one — usually one a week.

GitHub