单应性估计

原文地址: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\) 是图像平面上点之间的直接映射。如果已知某些点在图像[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\) 的特征向量。


  1. 对于通过纯旋转关联的相机,每个场景点可视为位于无穷远处的平面上。 ↩︎

posted @ 2025-06-04 23:07  GShang  阅读(95)  评论(0)    收藏  举报