§1.2.9

SVD 与最小二乘求解 Homography？DLT 算法与归一化的必要性？

核心概念

单应性（Homography）是射影几何中的一个概念，描述了两个平面之间的一种射影变换关系。在计算机视觉中，它特指一个 3x3 的矩阵 H，用于将在一个平面上的点 p 映射到另一个图像中对应平面上的点 p'，即 p' = H * p。直接线性变换（Direct Linear Transform, DLT）是一种通过建立线性方程组 Ah=0 来求解单应性矩阵 H 的算法。该算法的核心思想是将非线性的映射问题转化为一个齐次线性方程组，然后通过奇异值分解（Singular Value Decomposition, SVD）在最小二乘意义下求解。

原理与推导

1. Homography 问题的数学表述

假设在第一个平面（图像）上有一个点 p = [x, y, 1]^T，在第二个平面（图像）上有其对应的点 p' = [x', y', 1]^T。两者都使用齐次坐标表示。单应性变换关系为：

s \begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \mathbf{H} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}

其中 s 是一个非零的尺度因子。我们的目标是求解这个 3x3 的矩阵 H。H 有 9 个元素，但由于尺度等价性（乘以任意非零常数 k 的 kH 仍然表示同一个变换），它实际上只有 8 个自由度。

2. DLT 算法：构建线性方程组 `Ah=0`

为了求解 H，我们需要将其转化为一个标准的线性方程组。一种优雅的方式是利用向量叉乘的性质。因为 p' 和 Hp 是平行的向量（仅相差一个尺度因子 s），它们的叉乘为零向量：

\mathbf{p'} \times (\mathbf{H} \mathbf{p}) = \mathbf{0}

令 p' = [x', y', w']^T（这里 w'=1）和 Hp = [r_1, r_2, r_3]^T，其中 r_i 是 H 的第 i 行与 p 的点积。展开叉乘得到：

\begin{bmatrix} y' r_3 - w' r_2 \\ w' r_1 - x' r_3 \\ x' r_2 - y' r_1 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}

这三个方程中只有两个是线性独立的。我们取前两个（或任意两个）来构建方程。令 h_i^T 为 H 的第 i 行（作为行向量），h 为 H 矩阵按行展开的 9x1 向量 h = [h_1^T, h_2^T, h_3^T]^T。将 r_1 = h_1^T p, r_2 = h_2^T p, r_3 = h_3^T p 和 w'=1 代入前两个方程：

y'(h_3^T p) - 1(h_2^T p) = 0
1(h_1^T p) - x'(h_3^T p) = 0

将 p^T = [x, y, 1] 代入并整理，可以得到关于 h 的两个线性方程：

0^T h_1 - p^T h_2 + y' p^T h_3 = 0
p^T h_1 + 0^T h_2 - x' p^T h_3 = 0

写成矩阵形式 A_i h = 0，其中 A_i 是一个 2x9 的矩阵：

\mathbf{A}_i = \begin{bmatrix} \mathbf{0}^T & -x & -y & -1 & y'x & y'y & y' \\ x & y & 1 & \mathbf{0}^T & -x'x & -x'y & -x' \end{bmatrix}

每个点对 (p, p') 提供两个方程。由于 H 有 8 个自由度，我们至少需要 4 个不共线的点对来唯一确定 H（2N >= 8，N=4）。当有 N 个点对时，我们将 N 个 A_i 矩阵堆叠起来，得到一个 2N x 9 的大矩阵 A：

\mathbf{A} = \begin{bmatrix} \mathbf{A}_1 \\ \mathbf{A}_2 \\ \vdots \\ \mathbf{A}_N \end{bmatrix}

我们的问题就变成了求解齐次线性方程组：

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

3. SVD 求解最小二乘问题

在实际情况中，由于噪声的存在，Ah=0 通常没有精确解。因此，我们寻求一个最小二乘解，即最小化 ||Ah||^2，同时为了避免 h=0 的平凡解，我们增加一个约束 ||h||=1。

\min_{\mathbf{h}} ||\mathbf{A}\mathbf{h}||^2 \quad \text{subject to} \quad ||\mathbf{h}|| = 1

这个问题可以通过 SVD 完美解决。对矩阵 A 进行奇异值分解：

\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T

其中 V 是一个 9x9 的正交矩阵，其列是 A^T A 的特征向量。||Ah||^2 = h^T A^T A h。解 h 就是 V 中对应最小奇异值的列向量。在 Σ 中，奇异值是按降序排列的，所以最小的奇异值是 σ_9，对应的右奇异向量是 V 的最后一列。因此，h 的解就是 V 的最后一列。

算法复杂度:

构建 A 矩阵: O(N)
对 2N x 9 的矩阵 A 进行 SVD: O(N * 9^2) = O(N) (因为 9 是常数)。
总体时间复杂度为 O(N)。

4. 归一化 (Normalization) 的必要性

问题: DLT 算法最小化的是代数误差 ||Ah||^2，这并没有明确的几何意义。当输入点的坐标值很大或分布不均时（例如，图像左上角的点 (10, 20) 和右下角的点 (1900, 1070)），A 矩阵的行之间数值差异会非常大。例如，A_i 中包含 x, y 这样的项，也包含 x'x, x'y 这样的二次项。这会导致 A 矩阵的条件数 (condition number) 非常大，成为一个病态矩阵 (ill-conditioned matrix)。

后果:

数值不稳定性: 对微小的输入噪声非常敏感，SVD 求解的结果可能非常不稳定和不准确。
次优解: 最小化代数误差不等于最小化几何误差（即，在图像上的重投影误差）。

解决方案: Hartley 和 Zisserman 在其经典著作中提出了归一化 DLT 算法。

中心化: 对两组点分别进行平移，使它们的质心（centroid）移动到原点 (0,0)。
缩放: 对平移后的点进行缩放，使得它们到原点的平均距离为 sqrt(2)。（对于 2D 点 (x,y)，到原点的距离是 sqrt(x^2+y^2)，选择 sqrt(2) 是因为这使得一个单位正方形 [-1,1]x[-1,1] 的顶点的平均距离是 sqrt(2)）。

归一化步骤:

对第一组点 {p_i}，计算变换矩阵 T，使得归一化后的点 p_i_norm = T * p_i。
对第二组点 {p'_i}，计算变换矩阵 T'，使得归一化后的点 p'_i_norm = T' * p'_i。
使用归一化后的点对 (p_i_norm, p'_i_norm) 运行标准 DLT 算法，求解得到归一化后的单应性矩阵 H_norm。
将 H_norm 反归一化，得到最终的单应性矩阵 H: p'_i_norm = H_norm * p_i_norm (T' * p'_i) = H_norm * (T * p_i) p'_i = (T')^{-1} * H_norm * T * p_i 所以，H = (T')^{-1} * H_norm * T。

几何直观解释: 归一化将所有点都变换到一个统一的、"表现良好"的坐标系下（大致在 [-1, 1] 范围内）。在这个空间里，不同项的数值大小接近，A 矩阵的条件数得到极大改善。此时，代数误差 ||A_norm * h_norm||^2 成为几何误差的一个更好的近似，因此求解结果在几何上更加精确和稳定。

代码实现

下面提供了使用 NumPy 实现的 DLT 和归一化 DLT 算法。

python

 1import numpy as np
 2import cv2
 3
 4def solve_homography_dlt(src_pts, dst_pts):
 5    """
 6    使用 DLT 算法求解单应性矩阵 H (p_dst = H * p_src)。
 7    :param src_pts: 源点, shape (N, 2)
 8    :param dst_pts: 目标点, shape (N, 2)
 9    :return: H, 3x3 的单应性矩阵
10    """
11    if src_pts.shape[0] < 4:
12        raise ValueError("求解单应性矩阵至少需要4个点对")
13
14    num_points = src_pts.shape[0]
15    A = np.zeros((2 * num_points, 9))
16
17    # 构建 A 矩阵
18    # 为什么这样做：每个点对提供两个线性独立的方程，用于约束 H 的9个未知数
19    for i in range(num_points):
20        x, y = src_pts[i]
21        xp, yp = dst_pts[i]
22        
23        # 第一个方程: 0^T h1 - p^T h2 + y'p^T h3 = 0
24        A[2 * i] = [0, 0, 0, -x, -y, -1, yp * x, yp * y, yp]
25        # 第二个方程: p^T h1 + 0^T h2 - x'p^T h3 = 0
26        A[2 * i + 1] = [x, y, 1, 0, 0, 0, -xp * x, -xp * y, -xp]
27
28    # 使用 SVD 求解 Ah = 0
29    # 为什么这样做：SVD 是求解 ||Ah||^2 最小化问题的标准方法，解是对应最小奇异值的右奇异向量
30    _, _, V_T = np.linalg.svd(A)
31    
32    # h 是 V_T 的最后一行 (即 V 的最后一列)
33    h = V_T[-1]
34    
35    # 将 1x9 的向量 h 变回 3x3 的矩阵 H
36    H = h.reshape((3, 3))
37    
38    return H
39
40def solve_homography_normalized_dlt(src_pts, dst_pts):
41    """
42    使用归一化的 DLT 算法求解单应性矩阵 H。
43    :param src_pts: 源点, shape (N, 2)
44    :param dst_pts: 目标点, shape (N, 2)
45    :return: H, 3x3 的单应性矩阵
46    """
47    if src_pts.shape[0] < 4:
48        raise ValueError("求解单应性矩阵至少需要4个点对")
49
50    # --- 归一化步骤 ---
51    # 为什么这样做：改善 A 矩阵的条件数，提高数值稳定性和求解精度
52    
53    # 1. 计算源点的归一化矩阵 T_src
54    src_centroid = np.mean(src_pts, axis=0)
55    # 计算点到质心的平均距离
56    src_dist = np.mean(np.linalg.norm(src_pts - src_centroid, axis=1))
57    # 缩放因子，使得平均距离为 sqrt(2)
58    scale_src = np.sqrt(2) / src_dist
59    
60    T_src = np.array([
61        [scale_src, 0, -scale_src * src_centroid[0]],
62        [0, scale_src, -scale_src * src_centroid[1]],
63        [0, 0, 1]
64    ])
65
66    # 2. 计算目标点的归一化矩阵 T_dst
67    dst_centroid = np.mean(dst_pts, axis=0)
68    dst_dist = np.mean(np.linalg.norm(dst_pts - dst_centroid, axis=1))
69    scale_dst = np.sqrt(2) / dst_dist
70    
71    T_dst = np.array([
72        [scale_dst, 0, -scale_dst * dst_centroid[0]],
73        [0, scale_dst, -scale_dst * dst_centroid[1]],
74        [0, 0, 1]
75    ])
76
77    # 3. 将点转换为齐次坐标并进行归一化
78    src_pts_h = np.hstack((src_pts, np.ones((src_pts.shape[0], 1))))
79    dst_pts_h = np.hstack((dst_pts, np.ones((dst_pts.shape[0], 1))))
80    
81    norm_src_pts_h = (T_src @ src_pts_h.T).T
82    norm_dst_pts_h = (T_dst @ dst_pts_h.T).T
83    
84    # 转换回非齐次坐标以传入 DLT 函数
85    norm_src_pts = norm_src_pts_h[:, :2] / norm_src_pts_h[:, 2:]
86    norm_dst_pts = norm_dst_pts_h[:, :2] / norm_dst_pts_h[:, 2:]
87
88    # 4. 对归一化后的点使用 DLT 算法求解 H_norm
89    H_norm = solve_homography_dlt(norm_src_pts, norm_dst_pts)
90
91    # 5. 反归一化得到最终的 H
92    # 为什么这样做：将变换从归一化坐标系转换回原始坐标系
93    # H_norm * T_src * p_src = T_dst * p_dst  =>  p_dst = (T_dst)^-1 * H_norm * T_src * p_src
94    H = np.linalg.inv(T_dst) @ H_norm @ T_src
95    
96    # 保证 H_33 为 1 (可选，但常用)
97    H /= H[2, 2]
98    
99    return H
100
101if __name__ == '__main__':
102    # --- 测试代码 ---
103    # 1. 定义一组源点 (例如一个矩形的四个角点)
104    src_pts = np.array([
105        [100, 100],
106        [500, 100],
107        [500, 400],
108        [100, 400]
109    ], dtype=np.float32)
110
111    # 2. 定义一个真实的 H_gt，并生成目标点
112    # 这是一个轻微的透视变换
113    H_gt = np.array([
114        [1.1, 0.2, -50],
115        [-0.1, 1.2, 20],
116        [0.0005, 0.0002, 1]
117    ])
118    
119    src_pts_h = np.hstack((src_pts, np.ones((4, 1))))
120    dst_pts_h = (H_gt @ src_pts_h.T).T
121    dst_pts = dst_pts_h[:, :2] / dst_pts_h[:, 2:]
122
123    # 3. 添加一点噪声模拟真实情况
124    noise = np.random.normal(0, 1.5, dst_pts.shape)
125    noisy_dst_pts = dst_pts + noise
126
127    # 4. 使用不同方法求解 H
128    # OpenCV 的实现 (作为参考基准，它内部使用了 RANSAC 和归一化)
129    H_cv, _ = cv2.findHomography(src_pts, noisy_dst_pts)
130    
131    # 基础 DLT
132    H_dlt = solve_homography_dlt(src_pts, noisy_dst_pts)
133    H_dlt /= H_dlt[2, 2] # 归一化以便比较
134
135    # 归一化 DLT
136    H_norm_dlt = solve_homography_normalized_dlt(src_pts, noisy_dst_pts)
137
138    print("--- Ground Truth H ---")
139    print(H_gt / H_gt[2, 2])
140    print("\n--- Basic DLT Result ---")
141    print(H_dlt)
142    print("\n--- Normalized DLT Result ---")
143    print(H_norm_dlt)
144    print("\n--- OpenCV findHomography Result ---")
145    print(H_cv)
146
147    # 5. 评估误差 (重投影误差)
148    def reprojection_error(H, src, dst):
149        src_h = np.hstack((src, np.ones((src.shape[0], 1))))
150        pred_dst_h = (H @ src_h.T).T
151        pred_dst = pred_dst_h[:, :2] / pred_dst_h[:, 2:]
152        return np.mean(np.linalg.norm(pred_dst - dst, axis=1))
153
154    err_dlt = reprojection_error(H_dlt, src_pts, noisy_dst_pts)
155    err_norm_dlt = reprojection_error(H_norm_dlt, src_pts, noisy_dst_pts)
156    err_cv = reprojection_error(H_cv, src_pts, noisy_dst_pts)
157
158    print(f"\nReprojection Error (Basic DLT): {err_dlt:.4f} pixels")
159    print(f"Reprojection Error (Normalized DLT): {err_norm_dlt:.4f} pixels")
160    print(f"Reprojection Error (OpenCV): {err_cv:.4f} pixels")
161    # 通常，归一化 DLT 的误差会更小，结果更接近 Ground Truth

工程实践

RANSAC (Random Sample Consensus): DLT 对离群点 (outliers) 非常敏感，一个错误的点对就能让结果面目全非。在实际应用中，DLT 几乎总是与 RANSAC 结合使用。RANSAC 的流程是：
- 随机采样一个最小子集（对 Homography 是 4 个点对）。
- 使用这个子集通过 DLT 估计一个模型 H_candidate。
- 计算所有点对中，有多少是 "内点" (inliers)，即它们的重投影误差小于一个阈值 t。
- 重复以上步骤 K 次，选择拥有最多内点的 H_candidate 作为最佳模型。
- （可选）使用所有内点重新估计一次 H，得到更精确的结果。
超参数选择:
- RANSAC 阈值 t: 通常设置为 1-5 个像素。它取决于特征点检测器的精度。SIFT/ORB 等特征的匹配精度通常在像素级。
- RANSAC 迭代次数 K: 可以通过公式 K = log(1-p) / log(1 - w^n) 估算，其中 p 是期望的成功率（如 0.99），w 是内点比例，n 是采样点数（这里是 4）。实践中，通常直接设为一个较大的数，如 500 或 1000。
非线性优化: DLT/归一化 DLT 提供了很好的初始解。为了达到最高精度，可以在此基础上进行非线性优化，直接最小化几何误差（重投影误差）。例如，使用 Levenberg-Marquardt (LM) 算法，将 DLT 的解作为初始值，以 H 的 8 个参数为变量，最小化 sum(||p'_i - proj(H, p_i)||^2)。OpenCV 的 findHomography 在 RANSAC 之后就执行了这一步。
部署与性能:
- Homography 的计算（包括 RANSAC 和 DLT）通常在 CPU 上完成，对于实时应用（如 AR、SLAM）的几百个点对，这个计算速度足够快。
- 计算出的 3x3 矩阵 H 非常轻量。在推理阶段，它主要用于图像的透视变换（warping），这个操作在 GPU 上有高效的实现，如 cv2.warpPerspective 或 PyTorch 的 F.grid_sample。

常见误区与边界情况

误区: "点越多，DLT 结果越好"
- 纠正: 不一定。如果新增的点是离群点，结果会迅速恶化。DLT 没有鲁棒性。正确的说法是："内点越多，（在 RANSAC 框架下）结果越好"。
误区: "DLT 最小化的是图像上的像素距离"
- 纠正: 不是。DLT 最小化的是代数误差 ||Ah||^2，这是一个数学上方便但缺乏直观几何意义的量。归一化 DLT 的结果在几何上更优，因为它使得代数误差更接近几何误差。真正的几何误差最小化需要非线性优化。
边界情况: 共线点
- 如果用于估计 H 的 4 个源点中有 3 个或 4 个是共线的，A 矩阵会变得秩亏（rank deficient），导致无法得到唯一解。几何上，你无法从一条线（或其一部分）的变换确定整个平面的变换。cv2.findHomography 在这种情况下会报错或返回一个不正确的 H。
边界情况: 变换的物理有效性
- DLT 求解出的 H 只是一个数学矩阵，它不保证物理上的可实现性。例如，它可能导致一个 "翻转" 的图像（H 的行列式为负）。在某些应用中（如姿态估计），需要对 H 进行分解，检查其物理合理性。
面试追问:
- "如果不能用 SVD，怎么解 Ah=0？": 可以求解 A^T A 的特征值和特征向量。A^T A 是一个 9x9 的对称半正定矩阵，它的最小特征值对应的特征向量就是 h 的解。SVD 在数值上更稳定，是首选方法。
- "Homography 的自由度是多少？为什么？": 8 个。因为 3x3 矩阵有 9 个元素，但变换对全局尺度是不变的（p' = sHp），所以有一个自由度是冗余的。通常通过约束 h_33=1 或 ||h||=1 来消除这个歧义。
- "什么场景下两张图之间存在精确的 Homography 关系？": 1. 场景是一个平面。 2. 相机围绕其光心进行了纯旋转（无平移）。这是全景图拼接（Panorama Stitching）的基础。