Calculating rotation and translation of axes given original and transformed points

Question

How can we derive the angle by which the origin is rotated and axes are shifted given the list of points before the transformation was done and after the transformation is done.

example:

Before transformation - [(-3173.24, 1503.76), (-3173.24, 1599.98), (-2921.24, 1941.7), (-2777.24, 1941.7), (-2969.24, 1905.7), (-2969.24, 1941.7), (-3017.24, 1941.7), (-3065.24, 1905.7), (-3161.24, 2013.6), (-3179.24, 2049.6), (-2759.24, 1905.7), (-3017.24, 1803.81), (-3113.24, 1803.81), (-3161.24, 1803.81), (-3179.24, 1839.81), (-2759.24, 1803.81), (-2777.24, 1839.81), (-2789.24, 1623.98), (-2789.24, 1527.76), (-2760.86, 1737.92)]

After transformation - [(52.12, 146.39), (52.12, 242.61), (592.12, 584.33), (448.12, 584.33), (256.12, 548.33), (640.12, 584.33), (688.12, 584.33), (160.12, 548.33), (64.12, 656.22), (46.12, 692.22), (466.12, 548.33), (208.12, 446.44), (112.12, 446.44), (64.12, 446.44), (46.12, 482.44), (466.12, 446.44), (448.12, 482.44), (436.12, 266.61), (436.12, 170.39), (464.5, 380.54)]

How can we generate the correlation between these points so that given an older point(x, y) we can find the derived point (x', y')

Origin details are not known. We are using a Cartesian coordinate system having x and y axes.

Answer 1

New Edit.

import numpy as np
import matplotlib.pyplot as plt

x = np.array([(-3173.24, 1503.76), (-3173.24, 1599.98), (-2921.24, 1941.7), (-2777.24, 1941.7), (-2969.24, 1905.7), 
              (-2969.24, 1941.7), (-3017.24, 1941.7), (-3065.24, 1905.7), (-3161.24, 2013.6), (-3179.24, 2049.6), 
              (-2759.24, 1905.7), (-3017.24, 1803.81), (-3113.24, 1803.81), (-3161.24, 1803.81), (-3179.24, 1839.81), 
              (-2759.24, 1803.81), (-2777.24, 1839.81), (-2789.24, 1623.98), (-2789.24, 1527.76), (-2760.86, 1737.92)])
y = np.array([(52.12, 146.39), (52.12, 242.61), (592.12, 584.33), (448.12, 584.33), (256.12, 548.33), (640.12, 584.33), 
              (688.12, 584.33), (160.12, 548.33), (64.12, 656.22), (46.12, 692.22), (466.12, 548.33), (208.12, 446.44), 
              (112.12, 446.44), (64.12, 446.44), (46.12, 482.44), (466.12, 446.44), (448.12, 482.44), (436.12, 266.61), 
              (436.12, 170.39), (464.5, 380.54)])

def centroid(s):
    # Centroid of a set s
    c = np.sum(s, axis=0) / len(s)
    return c[np.newaxis,:]

def rotation_translation(x, y):
    # Centroids of each set
    xG = centroid(x)
    yG = centroid(y)
    # Shape of each set
    Mx = (x-xG).T.dot(x-xG)
    My = (y-yG).T.dot(y-yG)
    # Eigendecomposition
    ex, Ux = np.linalg.eig(Mx)
    ey, Uy = np.linalg.eig(My)
    # rotation matrix
    U = Uy.dot(Ux.T)
    # translation vector
    t = yG - xG.dot(U.T)
    return U, t, xG, yG

def matrix_vector(x, y):
    # Centroids of each set
    xG = centroid(x)
    yG = centroid(y)
    # Shape of each set
    Mx = (x-xG).T.dot(x-xG)
    My = (y-yG).T.dot(y-yG)
    # Eigendecomposition: eigeinvalues, rotation matrices
    ex, Ux = np.linalg.eig(Mx)
    ey, Uy = np.linalg.eig(My)
    ex = np.sqrt(ex)
    ey = np.sqrt(ey)
    ex = ex[:, np.newaxis]
    ey = ey[:, np.newaxis]
    # linear matrix, rotation-scaling-rotation
    U = (ey/ex) * Ux.T
    U = Uy.dot(U)
    # trasnlation vector
    t = yG - xG.dot(U.T)
    return U, t, xG, yG

def Transform(v, mtrx_vctr):
    return v.dot(mtrx_vctr[0].T) + mtrx_vctr[1] 


U, t, xG, yG = rotation_translation(x, y)
z_rot = Transform(x, (U, t))

U, t, xG, yG = matrix_vector(x, y)
z_aff = Transform(x, (U, t))


plt.plot(x[:,0], x[:,1], 'go')
plt.plot(y[:,0], y[:,1], 'bo')
plt.plot(z_rot[:,0], z_rot[:,1], 'ro')
plt.plot(z_aff[:,0], z_aff[:,1], 'yo')
plt.show()

Initial Answer.

These two triangles are not congruent, i.e. there is no translation and rotation that can transform x into y. There is an affine transformation.

import numpy as np
import matplotlib.pyplot as plt

def Affine_map(before, after):
    
    b0 = before - before[0,:][np.newaxis,:]
    b0 = b0[1:3,:]
    
    a0 = after - after[0,:][np.newaxis,:]
    a0 = a0[1:3,:]
    
    A = np.linalg.inv(b0).dot(a0) 
    t = after[0,:] - before[0,:].dot(A)
    return A, t
    

x = np.array([[-4648, 2144], [-4311, 2155], [-3589, 357]])
y = np.array([[591, 157], [1073, 168], [1596, -223]])

A, t = Affine_map(x, y)

# map is applied as y = x.dot(A) + t

print('')
print('Matrix of affine triansformation:')
print(A)
print('')
print('Vector of translation:')
print(t)
print('')
print('Check that y = Ax + t ')
print( y - x.dot(A) - t)
print('')

plt.plot(x[:,0], x[:,1], 'ro')
plt.plot(y[:,0], y[:,1], 'bo')
plt.show()

Answer 2

You are looking for a 2x2 rotation matrix M and a 2x1 translation vector T such as:

(x',y') = M.(x,y) + B 

for all of your points.

Let's say that M is the following Matrix:

( a b )
( c d )

And T the following vector:

(Tx, Ty)

Each one of your points before transformation P = (Px, Py) and after transformation P' = (P'x, P'y) gives you 2 equations:

P'x = a.Px + b.Py + Tx
P'y = c.Px + d.Py + Ty

• With a,b,c,d,Tx,Ty being the unknown elements.

• Gathering all the equations given by all your points, you can find those unknown elements using optimization algorithms like Levenberg Marquart or Gradient Descent ...

• Note that you have 6 unknowns, so you need at least 3 points which will give you 6 equations to solve this problem. (And if you feel like pen and paper, you can solve this by hand :P)

Once you solved the system of equation, you can get the rotation angle theta as the matrix M is the rotation matrix defined by:

( cos(theta) -sin(theta) )
( sin(theta)  cos(theta) )

• More details on rotation matrix and the relation between angle theta and matrix coefficients here: https://en.wikipedia.org/wiki/Rotation_matrix

Note 1: be aware that in the general case, "Rotation THEN Translation" is not equivalent to "Translation THEN Rotation", simply because if you rotate before, you won't translate in the same direction ;-)

• The equations I mentioned above correspond to "Rotation THEN Translation".

Note 2: my answer suppose that the transformation is really just rotation and translation but if you are not 100% sure that it is the case, the transformation can be an affine transformation (P' = M.P + T with M not being necessarily a rotation matrix) or even non-linear transformation (and in this case, it can't be represented with a matrix and I would use a 2D model like 2D polygons or stuff like that:

f(x,y) = a.x² + b.y² + c.xy + d.x + e.y + f and try a numerical method to find a,b,c,d,e,f) which are less straight forward.