ffriend - 1 month ago 5x
Python Question

# Procrustes Analysis with NumPy?

Is there something like Matlab's

function in NumPy/SciPy or related libraries?

For reference. Procrustes analysis aims to align 2 sets of points (in other words, 2 shapes) to minimize square distance between them by removing scale, translation and rotation warp components.

Example in Matlab:

``````X = [0 1; 2 3; 4 5; 6 7; 8 9];   % first shape
R = [1 2; 2 1];                  % rotation matrix
t = [3 5];                       % translation vector
Y = X * R + repmat(t, 5, 1);     % warped shape, no scale and no distortion
[d Z] = procrustes(X, Y);        % Z is Y aligned back to X
Z

Z =

0.0000    1.0000
2.0000    3.0000
4.0000    5.0000
6.0000    7.0000
8.0000    9.0000
``````

``````X = arange(10).reshape((5, 2))
R = array([[1, 2], [2, 1]])
t = array([3, 5])
Y = dot(X, R) + t
Z = ???
``````

Note: I'm only interested in aligned shape, since square error (variable
`d`
in Matlab code) is easily computed from 2 shapes.

I'm not aware of any pre-existing implementation in Python, but it's easy to take a look at the MATLAB code using `edit procrustes.m` and port it to Numpy:

``````def procrustes(X, Y, scaling=True, reflection='best'):
"""
A port of MATLAB's `procrustes` function to Numpy.

Procrustes analysis determines a linear transformation (translation,
reflection, orthogonal rotation and scaling) of the points in Y to best
conform them to the points in matrix X, using the sum of squared errors
as the goodness of fit criterion.

d, Z, [tform] = procrustes(X, Y)

Inputs:
------------
X, Y
matrices of target and input coordinates. they must have equal
numbers of  points (rows), but Y may have fewer dimensions
(columns) than X.

scaling
if False, the scaling component of the transformation is forced
to 1

reflection
if 'best' (default), the transformation solution may or may not
include a reflection component, depending on which fits the data
best. setting reflection to True or False forces a solution with
reflection or no reflection respectively.

Outputs
------------
d
the residual sum of squared errors, normalized according to a
measure of the scale of X, ((X - X.mean(0))**2).sum()

Z
the matrix of transformed Y-values

tform
a dict specifying the rotation, translation and scaling that
maps X --> Y

"""

n,m = X.shape
ny,my = Y.shape

muX = X.mean(0)
muY = Y.mean(0)

X0 = X - muX
Y0 = Y - muY

ssX = (X0**2.).sum()
ssY = (Y0**2.).sum()

# centred Frobenius norm
normX = np.sqrt(ssX)
normY = np.sqrt(ssY)

# scale to equal (unit) norm
X0 /= normX
Y0 /= normY

if my < m:
Y0 = np.concatenate((Y0, np.zeros(n, m-my)),0)

# optimum rotation matrix of Y
A = np.dot(X0.T, Y0)
U,s,Vt = np.linalg.svd(A,full_matrices=False)
V = Vt.T
T = np.dot(V, U.T)

if reflection is not 'best':

# does the current solution use a reflection?
have_reflection = np.linalg.det(T) < 0

# if that's not what was specified, force another reflection
if reflection != have_reflection:
V[:,-1] *= -1
s[-1] *= -1
T = np.dot(V, U.T)

traceTA = s.sum()

if scaling:

# optimum scaling of Y
b = traceTA * normX / normY

# standarised distance between X and b*Y*T + c
d = 1 - traceTA**2

# transformed coords
Z = normX*traceTA*np.dot(Y0, T) + muX

else:
b = 1
d = 1 + ssY/ssX - 2 * traceTA * normY / normX
Z = normY*np.dot(Y0, T) + muX

# transformation matrix
if my < m:
T = T[:my,:]
c = muX - b*np.dot(muY, T)

#transformation values
tform = {'rotation':T, 'scale':b, 'translation':c}

return d, Z, tform
``````
Source (Stackoverflow)