Procrustes Analysis with NumPy?

Question

Is there something like Matlab's procrustes 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

Same task in NumPy:

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.

Answer 1

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 != '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

Answer 2

There is a Scipy function for it: scipy.spatial.procrustes

I'm just posting its example here:

>>> import numpy as np
>>> from scipy.spatial import procrustes

>>> a = np.array([[1, 3], [1, 2], [1, 1], [2, 1]], 'd')
>>> b = np.array([[4, -2], [4, -4], [4, -6], [2, -6]], 'd')
>>> mtx1, mtx2, disparity = procrustes(a, b)
>>> round(disparity)
0.0

Answer 3

You can have both Ordinary Procrustes Analysis and Generalized Procrustes Analysis in python with something like this:

import numpy as np
        
def opa(a, b):
    aT = a.mean(0)
    bT = b.mean(0)
    A = a - aT 
    B = b - bT
    aS = np.sum(A * A)**.5
    bS = np.sum(B * B)**.5
    A /= aS
    B /= bS
    U, _, V = np.linalg.svd(np.dot(B.T, A))
    aR = np.dot(U, V)
    if np.linalg.det(aR) < 0:
        V[1] *= -1
        aR = np.dot(U, V)
    aS = aS / bS
    aT-= (bT.dot(aR) * aS)
    aD = (np.sum((A - B.dot(aR))**2) / len(a))**.5
    return aR, aS, aT, aD 
        
def gpa(v, n=-1):
    if n < 0:
        p = avg(v)
    else:
        p = v[n]
    l = len(v)
    r, s, t, d = np.ndarray((4, l), object)
    for i in range(l):
        r[i], s[i], t[i], d[i] = opa(p, v[i]) 
    return r, s, t, d

def avg(v):
    v_= np.copy(v)
    l = len(v_) 
    R, S, T = [list(np.zeros(l)) for _ in range(3)]
    for i, j in np.ndindex(l, l):
        r, s, t, _ = opa(v_[i], v_[j]) 
        R[j] += np.arccos(min(1, max(-1, np.trace(r[:1])))) * np.sign(r[1][0]) 
        S[j] += s 
        T[j] += t 
    for i in range(l):
        a = R[i] / l
        r = [np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]
        v_[i] = v_[i].dot(r) * (S[i] / l) + (T[i] / l) 
    return v_.mean(0)

For testing purposes, the output of each algorithm can be visualized as follows:

import matplotlib.pyplot as p; p.rcParams['toolbar'] = 'None';

def plt(o, e, b):
    p.figure(figsize=(10, 10), dpi=72, facecolor='w').add_axes([0.05, 0.05, 0.9, 0.9], aspect='equal')
    p.plot(0, 0, marker='x', mew=1, ms=10, c='g', zorder=2, clip_on=False)
    p.gcf().canvas.set_window_title('%f' % e)
    x = np.ravel(o[0].T[0])
    y = np.ravel(o[0].T[1])
    p.xlim(min(x), max(x)) 
    p.ylim(min(y), max(y))
    a = []
    for i, j in np.ndindex(len(o), 2):
        a.append(o[i].T[j])    
    O = p.plot(*a, marker='x', mew=1, ms=10, lw=.25, c='b', zorder=0, clip_on=False)
    O[0].set(c='r', zorder=1)
    if not b:
        O[2].set_color('b')
        O[2].set_alpha(0.4)
    p.axis('off')     
    p.show()

# Fly wings example (Klingenberg, 2015 | https://en.wikipedia.org/wiki/Procrustes_analysis)
arr1 = np.array([[588.0, 443.0], [178.0, 443.0], [56.0, 436.0], [50.0, 376.0], [129.0, 360.0], [15.0, 342.0], [92.0, 293.0], [79.0, 269.0], [276.0, 295.0], [281.0, 331.0], [785.0, 260.0], [754.0, 174.0], [405.0, 233.0], [386.0, 167.0], [466.0, 59.0]])
arr2 = np.array([[477.0, 557.0], [130.129, 374.307], [52.0, 334.0], [67.662, 306.953], [111.916, 323.0], [55.119, 275.854], [107.935, 277.723], [101.899, 259.73], [175.0, 329.0], [171.0, 345.0], [589.0, 527.0], [591.0, 468.0], [299.0, 363.0], [306.0, 317.0], [406.0, 288.0]])

def opa_out(a):
    r, s, t, d = opa(a[0], a[1])
    a[1] = a[1].dot(r) * s + t
    return a, d, False
plt(*opa_out([arr1, arr2, np.matrix.copy(arr2)]))

def gpa_out(a):
    g = gpa(a, -1) 
    D = [avg(a)]
    for i in range(len(a)):
        D.append(a[i].dot(g[0][i]) * g[1][i] + g[2][i])
    return D, sum(g[3])/len(a), True 
plt(*gpa_out([arr1, arr2]))

Answer 4

Probably you want to try this package with various flavors of different Procrustes methods, https://github.com/theochem/procrustes.