how to perform coordinates affine transformation using python? part 2

Question

I have same problem as described here: how to perform coordinates affine transformation using python?

I was trying to use method described but some reason I will get error messages. Changes I made to code was to replace primary system and secondary system points. I created secondary coordinate points by using different origo. In real case for which I am studying this topic will have some errors when measuring the coordinates.

primary_system1 = (40.0, 1160.0, 0.0)
primary_system2 = (40.0, 40.0, 0.0)
primary_system3 = (260.0, 40.0, 0.0)
primary_system4 = (260.0, 1160.0, 0.0)

secondary_system1 = (610.0, 560.0, 0.0) 
secondary_system2 = (610.0,-560.0, 0.0) 
secondary_system3 = (390.0, -560.0, 0.0)
secondary_system4 = (390.0, 560.0, 0.0)

Error I get from when executing is following.

*Traceback (most recent call last):
  File "affine_try.py", line 57, in <module>
    secondary_system3, secondary_system4 )
  File "affine_try.py", line 22, in solve_affine
    A2 = y * x.I
  File "/usr/lib/python2.7/dist-packages/numpy/matrixlib/defmatrix.py", line 850, in getI
    return asmatrix(func(self))
  File "/usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py", line 445, in inv
    return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
  File "/usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py", line 328, in solve
    raise LinAlgError, 'Singular matrix'
numpy.linalg.linalg.LinAlgError: Singular matrix*

What might be the problem ?

Answer 1

The problem is that your matrix is singular, meaning it's not invertible. Since you're trying to take the inverse of it, that's a problem. The thread that you linked to is a basic solution to your problem, but it's not really the best solution. Rather than just inverting the matrix, what you actually want to do is solve a least-squares minimization problem to find the optimal affine transform matrix for your possibly noisy data. Here's how you would do that:

import numpy as np

primary = np.array([[40., 1160., 0.],
                    [40., 40., 0.],
                    [260., 40., 0.],
                    [260., 1160., 0.]])

secondary = np.array([[610., 560., 0.],
                      [610., -560., 0.],
                      [390., -560., 0.],
                      [390., 560., 0.]])

# Pad the data with ones, so that our transformation can do translations too
n = primary.shape[0]
pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])
unpad = lambda x: x[:,:-1]
X = pad(primary)
Y = pad(secondary)

# Solve the least squares problem X * A = Y
# to find our transformation matrix A
A, res, rank, s = np.linalg.lstsq(X, Y)

transform = lambda x: unpad(np.dot(pad(x), A))

print "Target:"
print secondary
print "Result:"
print transform(primary)
print "Max error:", np.abs(secondary - transform(primary)).max()

The reason that your original matrix was singular is that your third coordinate is always zero, so there's no way to tell what the transform on that coordinate should be (zero times anything gives zero, so any value would work).

Printing the value of A tells you the transformation that least-squares has found:

A[np.abs(A) < 1e-10] = 0  # set really small values to zero
print A

results in

[[  -1.    0.    0.    0.]
 [   0.    1.    0.    0.]
 [   0.    0.    0.    0.]
 [ 650. -600.    0.    1.]]

which is equivalent to x2 = -x1 + 650, y2 = y1 - 600, z2 = 0 where x1, y1, z1 are the coordinates in your original system and x2, y2, z2 are the coordinates in your new system. As you can see, least-squares just set all the terms related to the third dimension to zero, since your system is really two-dimensional.