Affine 3D transformation in Python

Question

I'm programming a function in Python in Autodesk Maya (using PyMel for Maya)

I have three 3D points; p0, p1, p2.

Then they make a rigid transformation, so after the transformation (an affine transformation) I have their new positions; q0, q1, q2.

I also have a fourth point before the transformation; p3. I want to calculate its position after the same transformation; q4.

So I need to calculate the transformation matrix, and then apply it to p4. I don't know how to do either. List = an array of objects

import pymel.core as pm
import pymel.core.datatypes as dt

p0 = dt.Vector(pm.getAttr(list[0]+".tx"), pm.getAttr(list[0]+".ty"), pm.getAttr(list[0]+".tz"))
p1 = dt.Vector(pm.getAttr(list[1]+".tx"), pm.getAttr(list[1]+".ty"), pm.getAttr(list[1]+".tz"))
p2 = dt.Vector(pm.getAttr(list[2]+".tx"), pm.getAttr(list[2]+".ty"), pm.getAttr(list[2]+".tz")
p3 = dt.Vector(pm.getAttr(list[3]+".tx"), pm.getAttr(list[3]+".ty"), pm.getAttr(list[3]+".tz"))

The 3D points are read from animated objects in the Maya scene. So at another frame, I run this code to get

q0 = dt.Vector(pm.getAttr(list[0]+".tx"), pm.getAttr(list[0]+".ty"), pm.getAttr(list[0]+".tz"))
q1 = dt.Vector(pm.getAttr(list[1]+".tx"), pm.getAttr(list[1]+".ty"), pm.getAttr(list[1]+".tz"))
q2 = dt.Vector(pm.getAttr(list[2]+".tx"), pm.getAttr(list[2]+".ty"), pm.getAttr(list[2]+".tz"))
#q3 = TransformationMatrix between (p0,p1,p2) and (q0,q1,q2), applied to p3

I tried to calculate with vectors, but I ended up with errors due to divisions by zero... So I figured that a transformation matrix should solve it without problems.

I've got a deadline not far ahead and I REALLY need help solving this! PLEASE HELP!

Edit: how to perform coordinates affine transformation using python?

I need this function "solve_affine", but it should take only 3 points from each set instead of 4. And I can't use numpy...

Answer 1

Here's a solution using numpy and scipy. scipy is mostly used to generate random rotations, except for scipy.linalg.norm which is easy to code oneself. The main things used from numpy are cross product and matrix multiplication, which are also easy to code oneself.

The basic idea is this: given three non-collinear points x1,x2,x3, it's possible to find an orthogonal triple of vectors (axes) v1,v2,v3, with v1 in the direction of x2-x1, v2 in the plane spanned by (x2-x1) and (x3-x1), and v3 completing the triple.

The points y1,y2,y3 are rotated and translated relative to x1,x2,x3. The axes w1,w2,w3 generated from y1,y2,y3 are rotated (i.e., no translation) from v1,v2,v3. These two sets of triples are each orthogonal, so it's easy to find the rotation from them: R = W * transpose(V)

Once we have the rotation, finding the translation is simple: y1 = R*x + t, so t = y1 - R*x. It might be a better to use a least-squares solver and combine all three points to get an estimate of t.

import numpy
import scipy.linalg


def rand_rot():
    """Return a random rotation

    Return a random orthogonal matrix with determinant 1"""
    q, _ = scipy.linalg.qr(numpy.random.randn(3, 3))
    if scipy.linalg.det(q) < 0:
        # does this ever happen?
        print "got a negative det"
        q[:, 0] = -q[:, 0]
    return q


def rand_noncollinear():
    """Return 3 random non-collinear vectors"""
    while True:
        b = numpy.random.randn(3, 3)
        sigma = scipy.linalg.svdvals(b)
        if sigma[2]/sigma[0] > 0.1:
            # "very" non-collinear
            break
        # "nearly" collinear; try again

    return b[:, 0], b[:, 1], b[:, 2]


def normalize(a):
    """Return argument normalized"""
    return a/scipy.linalg.norm(a)


def colstack(a1, a2, a3):
    """Stack three vectors as columns"""
    return numpy.hstack((a1[:, numpy.newaxis],
                         a2[:, numpy.newaxis],
                         a3[:, numpy.newaxis]))


def get_axes(a1, a2, a3):
    """Generate orthogonal axes from three non-collinear points"""
    # I tried to do this with QR, but something didn't work
    b1 = normalize(a2-a1)
    b2 = normalize(a3-a1)
    b3 = normalize(numpy.cross(b1, b2))
    b4 = normalize(numpy.cross(b3, b1))
    return b1, b4, b3

# random rotation and translation
r = rand_rot()
t = numpy.random.randn(3)

# three non-collinear points
x1, x2, x3 = rand_noncollinear()
# some other point
x4 = numpy.random.randn(3)

# the images of the above in the transformation.
# y4 is for checking only -- won't be used to estimate r or t
y1, y2, y3, y4 = [numpy.dot(r, x) + t
                  for x in x1, x2, x3, x4]


v1, v2, v3 = get_axes(x1, x2, x3)
w1, w2, w3 = get_axes(y1, y2, y3)

V = colstack(v1, v2, v3)
W = colstack(w1, w2, w3)

# W = R V, so R = W * inverse(V); but V orthogonal, so inverse(V) is
# transpose(V):
rfound = numpy.dot(W, V.T)

# y1 = R x1 + t, so...
tfound = y1-numpy.dot(r, x1)

# get error on images of x2 and x3, just in case

y2err = scipy.linalg.norm(numpy.dot(rfound, x2) + tfound - y2)
y3err = scipy.linalg.norm(numpy.dot(rfound, x3) + tfound - y3)

# and check error image of x4 -- getting an estimate of y4 is the
# point of all of this
y4err = scipy.linalg.norm(numpy.dot(rfound, x4) + tfound - y4)

print "y2 error: ", y2err
print "y3 error: ", y3err
print "y4 error: ", y4err

Answer 2

Both the description and your code are confusing. Description is a bit vague while the code examples are missing important bits and pieces. So here is how I understand the question:

Knowing three points in two spaces how to construct a transform from space A to space B?

Image 1: How to form a transformation between 2 spaces.

The answer depends on the type of transform the spaces have. You see three points always form a planar span. This means that you can know what the rotation, transform, and uniform scale of the new space is. You can also know the shear on the plane, as well as nonuniform scale. However, you can not know what the shear or nonuniform scale would be in the plane normal direction.

Therefore to make sense the question mutates into how to rotate and translate two spaces to match? this is pretty easy to do Translation part is directly:

trans = q0 - p0

That leaves you with rotation which has been explained in several posts:

• python + maya: Rotate Y axis to be along vector
• How to convert three dimensional vector to an Euler rotation in software like Maya using python

You can also calculate a scaling factor after this.

Answer 3

I've figured it out

p0p1 = p1-p0
p0p2 = p2-p0
p0p3 = p3-p0

q0q1 = q1-q0
q0q2 = q2-q0
q0q3 = q3-q0

before = dt.Matrix([p0.x, p0.y, p0.z, 0],[p1.x, p1.y, p1.z, 0],[p2.x, p2.y, p2.z, 0], [0,0,0,1]);
after = dt.Matrix([q0.x, q0.y, q0.z, 0],[q1.x, q1.y, q1.z, 0],[q2.x, q2.y, q2.z, 0], [0,0,0,1]);

normal = p0p1.cross(p0p2).normal()
dist = p0p3.dot(normal)
q3 = p3 - dist*normal

transformMatrix = before.inverse()*after
solve = dt.Matrix(q3.x, q3.y, q3.z, 1)*transformMatrix

q3 =  dt.Vector(solve[0][0], solve[0][1], solve[0][2])

newNormal = q0q1.cross(q0q2).normal()
q3 = q3 + newNormal*dist

pm.move(list[3], q3, r=False)

The transformation matrix only worked for points that are within the plane p0p1p2. So I solved it by transforming the projected point of p3, then move it out from the plane by the same distance.

If you have a solution that only involves a matrix, feel free to share, it may still help me! :)