DLT 算法计算变换矩阵及源码

彰阙有常

已于 2023-05-10 17:19:27 修改

阅读量864

点赞数

分类专栏：图像配准文章标签：矩阵算法线性代数计算机视觉

于 2022-05-11 11:57:30 首次发布

本文链接：https://blog.csdn.net/qq_38844437/article/details/124683679

版权

图像配准专栏收录该内容

1 篇文章 0 订阅

订阅专栏

假设原始图像上特征点坐标为 $x_i=(u_i,v_i,1)$ ，目标图像上的坐标为 $x_i'=(u_i',v_i',1)$
假设 $H$ 为两者之间的变换关系。
$Hx_i= \begin{bmatrix} h_1&h_2&h_3\\h_4&h_5&h_6\\h_7&h_8&h_9 \end{bmatrix} \begin{bmatrix} u_i\\v_i\\1 \end{bmatrix} =\begin{bmatrix} h^{1T}x_i\\h^{2T}x_i\\h^{3T}x_i\end{bmatrix}= \begin{bmatrix}x_i^Th^1\\x_i^Th^2\\x_i^Th^3\end{bmatrix}= \begin{bmatrix}u_i'\\v_i'\\1 \end{bmatrix}$

$h_i$ 代表 H 矩阵第 i 行的列向量表示。

可以通过 $x_i'\cdot x_i$ 得到三个线性方程：

$\begin{cases} \frac{x_i^Th^3}{x_i^Th_1}=\frac{1}{u_i'}\\ \frac{x_i^Th^3}{x_i^Th_2}=\frac{1}{v_i'}\\ \frac{x_i^Th^2}{x_i^Th_1}=\frac{v_i'}{u_i'} \end{cases} \rightarrow \begin{cases} v_i'x_ih^3-x_i^Th_2=0\\ u_i'x_i^Th_3-x_i^Th_1=0\\ u_i'x_i^Th_2-v_i'x_i^Th_1=0 \end{cases} \rightarrow \begin{bmatrix} 0_{3\times1}^T&-x_i^T&v_i'x_i^T\\x_i^T&0_{3\times1}^T&-u_i'x_i^T\\-v_i'x_i^T&u_i'x_i^T&0_{3\times1}^T \end{bmatrix} \begin{bmatrix} h^1\\h^2\\h^3 \end{bmatrix} =0$

即有 $A_i^{(3)}h=0$ 对于图像片的每一个点。 $A_i$ 是一个 $3\times9$ 矩阵，因为 $A$ 的最后一行是由另外两行确定，所以只需要计算线性方程 $\hat A_i \hat h=0$ 。这里的 $\hat A_i$ 由 $A_i^{(3)}$ 的前两行组成。A 的矩阵表示为：
在这里插入图片描述
在估计H时，可以使用更多的匹配点信息来减少估计误差。选取前两行组成一个独立的系数矩阵Ai，将所有的匹配点都考虑在内，形成一个2N × 9的系数矩阵a。用最小二乘法，h的解可以表示为

求解上述公式存在两种方式：

奇异值求解
通过SVD 求解 A 的最小右奇异向量作为 h的估计结果。
转化为非齐次线性最小二乘形式。
假设 h9=1，这时等式可以转换为

如果包含n 个匹配点

第一种方式源码如下：

# Ametrix to solv DLT
    def matrix_generate(sample_n, cf1, cf2):
        """
        AH=b
        生成A
        没对样本能生成一行
        [[0,0,0,-x,-y,1,xy',yy',y'],[x,y,1,0,0,0,-xx',-yx',-x']]
        :param sample_n:
        :param cf1: 点1
        :param cf2:
        :return: A [2N,9]
        """
        A = np.zeros([sample_n * 2, 9], dtype=np.float)
        for k in range(sample_n):
            # 每对点对应量
            A[2 * k, 0] = cf1[k, 0]
            A[2 * k, 1] = cf1[k, 1]
            A[2 * k, 2] = 1
            A[2 * k, 6] = (-cf2[k, 0]) * cf1[k, 0]
            A[2 * k, 7] = (-cf2[k, 0]) * cf1[k, 1]
            A[2 * k, 8] = (-cf2[k, 0])

            A[2 * k + 1, 3] = cf1[k, 0]
            A[2 * k + 1, 4] = cf1[k, 1]
            A[2 * k + 1, 5] = 1
            A[2 * k + 1, 6] = (-cf2[k, 1]) * cf1[k, 0]
            A[2 * k + 1, 7] = (-cf2[k, 1]) * cf1[k, 1]
            A[2 * k + 1, 8] = (-cf2[k, 1])
        return A

def compute_H(src_points,dst_points):
		A = self.matrix_generate(sample_n, cf1, cf2)

        # SVD 奇异值求解
        # Singular Value Decomposition
        W, U, V = cv2.SVDecomp(A)

        # get global-homography
        h = V[-1, :]
        h = h.reshape((3, 3))
        # 反归一化用与还原
        h = h / h[2, 2]
        return h