Determining a homogeneous affine transformation matrix from six points in 3D using Python

Question

I am given the locations of three points:

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

and their transformed counterparts:

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

The affine transformation matrix is of the form

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

such that with

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]]

you would get:

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

Assuming the transformation is homogeneous (consists of only rotations and translations), how can I determine this transformation matrix?

From a CAD program, I know the matrix is:

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]])

I'd like to know how this can be found.

Answer 1

Six points alone is not enough to uniquely determine the affine transformation. However, based on what you had asked in a question earlier (shortly before it was deleted) as well as your comment, it would seem that you are not merely looking for an affine transformation, but a homogeneous affine transformation.

This answer by robjohn provides the solution to the problem. Although it solves a more general problem with many points, the solution for 6 points can be found at the very bottom of the answer. I shall transcribe it here in a more programmer-friendly format:

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)))

For your sample inputs, this recovers the exact same matrix as what you had obtained from the CAD program:

>>> 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.        ]])

Answer 2

Finding a transformation is like solving any system of equations with unknown. First you have to write down the equation, which in your case means that you must know what transformation you are looking for. E.g. a rigid translation takes three parameters (x, y, and z) so you must have at least three parameters. General rotation takes another three parameters, which give you six unknowns. Scaling give you another three parameters for a total of 9 parameters. Since you state only three points, that yield nine unknows, this is the transformation that you are looking for. This means that you are ignoring any shearing and projection. You should always know the type of transformation that you are looking for.

Once you know the type of transformation you should write down the matrix equation, and then solve for the unknowns. This can be done with a linear algerbra library through a matrix multiplication, e.g. by numpy.

Answer 3

This problem is called point-to-point registration or point-set registration.

For a rigid transform, ie ignoring shearing and scaling, I like this tutorial. I involved finding the centroids and applying Singular Value Decomposition.

Note that for your particular case, with exactly three points, you could find a closed form solution.

OpenCV is good to help with this.

Oh, check out Finding translation and scale on two sets of points to get least square error in their distance?

Answer 4

It is possible to determine transformation matrix if original data (p1,p2,p3 in your case) and transformed data (p1_prime,p2_prime,p3_prime) are given as shown below:

>>> 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]])

In your case since there is some unknown data transformed ap[3] values for all the three points, the transformation matrix cannot be obtained. It can only be obtained if these three values are known.