原文地址：homography_estimation.pdf
1. 从三维坐标到二维坐标
在单应性变换下，我们可以将点从相机1的三维坐标到相机2的三维坐标的变换表示为：
\[\textbf{X}_2 = H\textbf{X}_1 \quad \textbf{X}_1, \textbf{X}_2 \in \mathbb{R}^3 \tag{1}
\]
在图像平面上，使用齐次坐标，我们有
\[\lambda_1 \textbf{x}_1 = \textbf{X}_1, \quad \lambda_2 \textbf{x}_2 = \textbf{X}_2 \tag{2}
\]
因此
\[\lambda_2 \textbf{x}_2 = H \lambda_1 \textbf{x}_1
\]
这意味着由于尺度歧义性，\(\textbf{x}_2\) 与 \(H\textbf{x}_1\) 在尺度上相等。注意，\(\textbf{x}_2 \sim H\textbf{x}_1\) 是图像平面上点之间的直接映射。如果已知某些点在图像中都位于同一平面上，则无需恢复和处理三维坐标即可直接对图像进行校正。
2. 单应性估计
为了估计单应矩阵 \(H\)，我们从方程 \(\textbf{x}_2 \sim H\textbf{x}_1\) 出发。在齐次坐标下，按元素展开可得到以下约束：
\[\begin{bmatrix} x_2 \\ y_2 \\ z_2 \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_1 \\ y_1 \\ z_1 \end{bmatrix} \quad \Leftrightarrow \quad \textbf{x}_2 = H\textbf{x}_1 \quad \tag{3}
\]
在非齐次坐标下（\(x_2' = x_2 / z_2\)且\(y_2' = y_2 / z_2\)），有
\[x_2' = \frac{H_{11}x_1 + H_{12}y_1 + H_{13}z_1}{H_{31}x_1 + H_{32}y_1 + H_{33}z_1} \quad \tag{4}
\]
\[y_2' = \frac{H_{21}x_1 + H_{22}y_1 + H_{23}z_1}{H_{31}x_1 + H_{32}y_1 + H_{33}z_1} \quad \tag{5}
\]
不失一般性，设 \(z_1 = 1\) 并重新整理可得：
\[x_2'(H_{31}x_1 + H_{32}y_1 + H_{33}) = H_{11}x_1 + H_{12}y_1 + H_{13} \quad \tag{6}
\]
\[y_2'(H_{31}x_1 + H_{32}y_1 + H_{33}) = H_{21}x_1 + H_{22}y_1 + H_{23} \quad \tag{7}
\]
我们希望求解 \(H\)。尽管这些非齐次方程涉及非线性坐标，但\(H\)的系数呈线性出现。重新整理方程(6)和(7)可得：
\[\textbf{a}_x^T \textbf{h} = 0 \quad \tag{8}
\]
\[\textbf{a}_y^T \textbf{h} = 0 \quad \tag{9}
\]
其中
\[\textbf{h} = (H_{11}, H_{12}, H_{13}, H_{21}, H_{22}, H_{23}, H_{31}, H_{32}, H_{33})^T \tag{10}
\]
\[\textbf{a}_x = (-x_1, -y_1, -1, 0, 0, 0, x_2'x_1, x_2'y_1, x_2')^T \tag{11}
\]
\[\textbf{a}_y = (0, 0, 0, -x_1, -y_1, -1, y_2'x_1, y_2'y_1, y_2')^T \tag{12}
\]
给定一组对应点，我们可以构建如下线性方程组：
\[A\textbf{h} = 0 \quad \tag{13}
\]
其中
\[A = \begin{pmatrix} \textbf{a}_{x1}^T \\ \textbf{a}_{y1}^T \\ \vdots \\ \textbf{a}_{xN}^T \\ \textbf{a}_{yN}^T \end{pmatrix} \quad \tag{14}
\]
方程(13)可通过齐次线性最小二乘法求解，具体内容将在下一节讨论。
3. 齐次线性最小二乘法
我们经常会遇到形如
\[A\textbf{x} = 0 \quad \tag{15}
\]
的问题，这类问题被称为齐次线性最小二乘问题。它与非齐次线性最小二乘问题
\[A\textbf{x} = b \tag{16}
\]
在形式上相似，后者通常使用伪逆或矩阵求逆的方法求解 \(\textbf{x}\)，但该方法不适用于方程(15)。相反，我们使用奇异值分解（SVD）来解决此问题。
从前一节的方程(13)出发，我们首先计算矩阵 \(A\) 的奇异值分解：
\[A = U\Sigma V^\top = \sum_{i=1}^9 \sigma_i \textbf{u}_i \textbf{v}_i^\top \tag{17}
\]
在Matlab中执行奇异值分解时，奇异值 \(\sigma_i\) 将按降序排列，因此 \(\sigma_9\) 是最小的奇异值。\(\sigma_9\) 的取值存在以下三种情况：
- 若单应性被精确确定，则 \(\sigma_9 = 0\)，此时存在一个能完美拟合所有点的单应矩阵。
- 若单应性为超定情况，则 \(\sigma_9 \geq 0\)，此时 \(\sigma_9\)表示“残差”或拟合优度。
- 我们不处理单应性为欠定的情况。
从奇异值分解结果中，我们选取与最小奇异值 \(\sigma_9\) 对应的“右奇异向量”（即矩阵 \(V\) 的某一列），该向量即为解 \(\textbf{h}\)，其包含了最佳拟合给定点的单应矩阵系数。我们将 \(\textbf{h}\) 重塑为矩阵 \(H\)，并构建方程 \(\textbf{x}_2 \sim H\textbf{x}_1\)。
4. 齐次线性最小二乘法的另一种推导
再次从方程 \(Ah = 0\) 出发，误差平方和可表示为：
\[\begin{align*}
f(\textbf{h}) &= \frac{1}{2}(A\textbf{h} - 0)^T(A\textbf{h} - 0) \tag{18}\\
f(\textbf{h}) &= \frac{1}{2}(A\textbf{h})^T(A\textbf{h}) \tag{19}\\
f(\textbf{h}) &= \frac{1}{2}\textbf{h}^T A^T A \textbf{h} \tag{20}\\
\end{align*}
\]
对 \(f\)关于\(h\) 求导并令结果为零，可得：
\[\begin{align*}
\frac{d}{d\textbf{h}}f = 0 &= \frac{1}{2}(A^T A + (A^T A)^T)\textbf{h} \tag{21} \\
0 &= A^T A \textbf{h} \tag{22}
\end{align*}
\]
从 \(A^T A\) 的特征分解可以看出，\(\textbf{h}\) 应等于 \(A^T A\) 中具有零特征值（或在存在噪声时最接近零的特征值）的特征向量。
这一结果与使用奇异值分解得到的结果一致，这可从以下事实轻松看出：
事实1 给定矩阵\(A\)的奇异值分解 \(A = U\Sigma V^T\)，矩阵 \(V\) 的列对应于 \(A^T A\) 的特征向量。
浙公网安备 33010602011771号