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?](http://stackoverflow.com/questions/8873462/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: 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
    

