PLUS-WAVE's WORLD

DLT算法求解单应性矩阵

 
  • 2024年 03月13日
  • 0 阅读

DLT算法求解单应性矩阵

原理：

单应性矩阵描述了两个图像之间的投影变换关系，即从一张图到另一张图的变换。

下面是DLT算法的基本原理：

  1. 构建投影方程： 对于两个图像中的对应点 $(x, y, 1)$ 和 $(u, v, 1)$ ，投影关系可以用齐次坐标表示为 $c \begin{bmatrix} u \ v \ 1 \end{bmatrix} = H \begin{bmatrix} x \ y \ 1 \end{bmatrix}$ （ $c$ 只是一个常数，由于是齐次坐标，所以不影响 ）。这里的 $H$（ $3 \times 3$ 矩阵）是我们要求解的单应性矩阵
\[H =\begin{bmatrix} h1 & h2 & h3 \\ h4 & h5 & h6 \\ h7 & h8 & h9 \end{bmatrix}\]
  1. 构建矩阵 $A$： 将上面的投影方程展开成 $Ah = 0$ 的形式，其中 $A$ 是一个 $2n \times 9$ 的矩阵， $h$ 是包含矩阵 $H$ 所有元素的列向量
\[A = \begin{bmatrix} -x_1 & -y_1 & -1 & 0 & 0 & 0 & u_1x_1 & u_1y_1 & u_1 \\ 0 & 0 & 0 & -x_1 & -y_1 & -1 & v_1x_1 & v_1y_1 & v_1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ -x_n & -y_n & -1 & 0 & 0 & 0 & u_nx_n & u_ny_n & u_n \\ 0 & 0 & 0 & -x_n & -y_n & -1 & v_nx_n & v_ny_n & v_n \\ \end{bmatrix}\] \[h=\begin{bmatrix}h1&h2&h3&h4&h5&h6&h7&h8&h9\end{bmatrix}\]

  1. 奇异值分解（SVD）： 对矩阵 $A$ 进行奇异值分解，得到 $A = U \Sigma V^T$。取 $V^T$ 的最后一列作为 $h$ 的估计

    方程的最小二乘解有一个既定的结论，即对 $A$ 进行SVD分解，得到的 ==$V^T$ 的最后一行== 即是 $h$ 的解，对 $h$ 做 reshape 得到 $H$ 。

实现：

根据你提供的信息，DLT（Direct Linear Transform）算法用于通过最小二乘法来估计单应性矩阵 $H$，以拟合两组特征点之间的关系。下面是DLT算法的具体步骤：

  1. 构建矩阵 $A$： 对于每一对特征点 $(x, y, 1)$ 和 $(u, v, 1)$，构建一个对应的矩阵 $A_i$。将所有这些矩阵==堆叠==成一个大矩阵 $A$

     def get_Ai(xi_vector, xi_prime_vector):
     assert (xi_vector.shape == (3,) and xi_prime_vector.shape == (3,))
	
     Ai = np.zeros((2, 9))
     Ai[0] = np.array(
         [-xi_vector[0], -xi_vector[1], -1, 0, 0, 0, xi_vector[0] * xi_prime_vector[0],
          xi_vector[1] * xi_prime_vector[0], xi_prime_vector[0]])
     Ai[1] = np.array(
         [0, 0, 0, -xi_vector[0], -xi_vector[1], -1, xi_vector[0] * xi_prime_vector[1],
          xi_vector[1] * xi_prime_vector[1], xi_prime_vector[1]])
	
     assert (Ai.shape == (2, 9))
     return Ai

    \(A_i = \begin{bmatrix} -x & -y & -1 & 0 & 0 & 0 & ux & vy & uv \\ 0 & 0 & 0 & -x & -y & -1 & ux & vy & uv \\ \end{bmatrix}\)

     def get_A(points_source, points_target):
     N = points_source.shape[0]
     A = np.zeros((2 * N, 9))
     for i in range(len(points_target)):
         Ai = get_Ai(points_source[i], points_target[i])
         A[2 * i:2 * i + 2] = Ai
	
     assert (A.shape == (2 * N, 9))
     return A

  2. SVD分解： 对矩阵 $A$ 进行奇异值分解（SVD），得到 $A = U \Sigma V^T$ 。取 $V^T$ 矩阵的最后一行作为矩阵 $h$

  3. Reshape： 将向量 $h$ reshape 为 $3 \times 3$ 的单应性矩阵 $H$

     def get_homography(points_source, points_target):
     A = get_A(points_source, points_target)
     _, _, Vt = np.linalg.svd(A)
     Homo = Vt[-1].reshape((3, 3))
     assert (Homo.shape == (3, 3))
     return Homo
搜索