使用 Python 从 3D 中的六个点确定齐次仿射变换矩阵

Posted 2023-02-21

【中文标题】使用 Python 从 3D 中的六个点确定齐次仿射变换矩阵

【英文标题】：Determining a homogeneous affine transformation matrix from six points in 3D using Python

【发布时间】：2014-12-18 11:58:31

【问题描述】：

我得到了三个点的位置：

p1 = [1.0, 1.0, 1.0]
p2 = [1.0, 2.0, 1.0]
p3 = [1.0, 1.0, 2.0]

及其转化的对应物：

p1_prime = [2.414213562373094,  5.732050807568877, 0.7320508075688767]
p2_prime = [2.7677669529663684, 6.665063509461097, 0.6650635094610956]
p3_prime = [2.7677669529663675, 5.665063509461096, 1.6650635094610962]

仿射变换矩阵的形式

trans_mat = np.array([[…, …, …, …],
                      […, …, …, …],
                      […, …, …, …],
                      […, …, …, …]])

这样与

import numpy as np

def transform_pt(point, trans_mat):
    a  = np.array([point[0], point[1], point[2], 1])
    ap = np.dot(a, trans_mat)[:3]
    return [ap[0], ap[1], ap[2]]

你会得到：

transform_pt(p1, trans_mat) == p1_prime
transform_pt(p2, trans_mat) == p2_prime
transform_pt(p3, trans_mat) == p3_prime

假设变换是齐次的（仅由旋转和平移组成），我如何确定这个变换矩阵？

通过 CAD 程序，我知道矩阵是：

trans_mat = np.array([[0.866025403784, -0.353553390593, -0.353553390593, 0],
                      [0.353553390593,  0.933012701892, -0.066987298108, 0],
                      [0.353553390593, -0.066987298108,  0.933012701892, 0],
                      [0.841081377402,  5.219578794378,  0.219578794378, 1]])

我想知道如何找到它。

【参考方案1】：

仅六个点不足以唯一确定仿射变换。但是，根据您之前在一个问题中提出的问题（在它被删除之前不久）以及your comment，您似乎不仅在寻找仿射变换，而且在寻找齐次仿射变换 em>。

This answer by robjohn 提供了问题的解决方案。虽然它解决了一个更普遍的问题，但可以在答案的最底部找到 6 分的解决方案。我将在此处以对程序员更友好的格式对其进行转录：

import numpy as np

def recover_homogenous_affine_transformation(p, p_prime):
    '''
    Find the unique homogeneous affine transformation that
    maps a set of 3 points to another set of 3 points in 3D
    space:

        p_prime == np.dot(p, R) + t

    where `R` is an unknown rotation matrix, `t` is an unknown
    translation vector, and `p` and `p_prime` are the original
    and transformed set of points stored as row vectors:

        p       = np.array((p1,       p2,       p3))
        p_prime = np.array((p1_prime, p2_prime, p3_prime))

    The result of this function is an augmented 4-by-4
    matrix `A` that represents this affine transformation:

        np.column_stack((p_prime, (1, 1, 1))) == \
            np.dot(np.column_stack((p, (1, 1, 1))), A)

    Source: https://math.stackexchange.com/a/222170 (robjohn)
    '''

    # construct intermediate matrix
    Q       = p[1:]       - p[0]
    Q_prime = p_prime[1:] - p_prime[0]

    # calculate rotation matrix
    R = np.dot(np.linalg.inv(np.row_stack((Q, np.cross(*Q)))),
               np.row_stack((Q_prime, np.cross(*Q_prime))))

    # calculate translation vector
    t = p_prime[0] - np.dot(p[0], R)

    # calculate affine transformation matrix
    return np.column_stack((np.row_stack((R, t)),
                            (0, 0, 0, 1)))

对于您的样本输入，这将恢复与您从 CAD 程序中获得的完全相同的矩阵：

>>> recover_homogenous_affine_transformation(
        np.array(((1.0,1.0,1.0),
                  (1.0,2.0,1.0),
                  (1.0,1.0,2.0))),
        np.array(((2.4142135623730940, 5.732050807568877, 0.7320508075688767),
                  (2.7677669529663684, 6.665063509461097, 0.6650635094610956),
                  (2.7677669529663675, 5.665063509461096, 1.6650635094610962))))
array([[ 0.8660254 , -0.35355339, -0.35355339,  0.        ],
       [ 0.35355339,  0.9330127 , -0.0669873 ,  0.        ],
       [ 0.35355339, -0.0669873 ,  0.9330127 ,  0.        ],
       [ 0.84108138,  5.21957879,  0.21957879,  1.        ]])

【参考方案2】：

寻找变换就像求解任何未知的方程组。首先，您必须写下方程式，这在您的情况下意味着您必须知道您正在寻找什么转换。例如。刚性平移需要三个参数（x、y 和 z），因此您必须至少具有三个参数。一般旋转需要另外三个参数，这给你六个未知数。缩放为您提供另外三个参数，总共 9 个参数。由于你只陈述了三个点，产生了九个不知道，这就是你正在寻找的转变。这意味着您忽略了任何剪切和投影。您应该始终知道您正在寻找的转换类型。

一旦你知道了变换的类型，你应该写下矩阵方程，然后求解未知数。这可以通过矩阵乘法使用线性代数库来完成，例如通过 numpy。

【参考方案3】：

可以确定变换矩阵如果原始数据（p1,p2,p3 in 你的情况）和变换数据（p1_prime,p2_prime,p3_prime )如下所示：

>>> p   # original data
array([[ 1.,  1.,  1.],
       [ 1.,  2.,  1.],
       [ 1.,  1.,  2.]])
>>> p_prime  # transformed data
array([[ 2.41421356,  5.73205081,  0.73205081],
       [ 2.76776695,  6.66506351,  0.66506351],
       [ 2.76776695,  5.66506351,  1.66506351]])
# Get transformation matrix
>>> trans = np.dot(np.linalg.inv(p),p_prime)
>>> trans  # transformation matrix
array([[ 1.70710678,  4.86602541, -0.13397459],
       [ 0.35355339,  0.9330127 , -0.0669873 ],
       [ 0.35355339, -0.0669873 ,  0.9330127 ]])
# obtain transformed data from original data and transformation matrix
>>> np.dot(a, trans)  
array([[ 2.41421356,  5.73205081,  0.73205081],
       [ 2.76776695,  6.66506351,  0.66506351],
       [ 2.76776695,  5.66506351,  1.66506351]])

在您的情况下，由于所有三个点都有一些未知数据转换ap[3] 值，因此无法获得转换矩阵。只有知道这三个值才能得到。

【讨论】：

CAD 程序提供的变换矩阵的第 4 维可能是由于使用了所谓的齐次坐标。它们通常用于计算机图形学中，因此可以将平移计算为矩阵乘法，从而与旋转变换相结合。要获得齐次坐标，将 1 作为第 4 维附加到所有点，本文中的计算应该是 OP 的矩阵中的结果。

在数学上转换为齐次坐标，将 3d 空间中的仿射变换 (y = Ax + b) 转换为线性变换 (y = Ax) 4d 空间。

@BHAT...如果该点不是最初的 6 个点之一怎么办？在这种情况下，这个解决方案似乎失效了。我们可以假设没有剪切或缩放……只有变换和旋转。

@Jenny_Winters，我没明白你的意思。 如果该点不是最初的 6 个点之一，您想问什么？

【参考方案4】：

此问题称为点对点注册或point-set registration。

对于刚性变换，即忽略剪切和缩放，我喜欢这个tutorial。我参与了寻找质心和应用奇异值分解。

请注意，对于您的特定情况，只需三点，您就可以找到一个封闭形式的解决方案。

OpenCV 可以很好地帮助解决这个问题。

哦，看看Finding translation and scale on two sets of points to get least square error in their distance?