Procrustes Superimposition - Again!

If you copy code, copy it right.

Posted by Falk Mielke on 7th Sep 2020

Procrustes Superimposition - Again!

If you copy code, copy it right.

Falk Mielke, Functional Morphology, Universiteit Antwerpen

You can download this notebook here.


It is hard to say whether any procedure has been more central to my work than Procrustes Superimposition (the only candidates are Fourier Series, MCMC sampling, and Principal Component Analysis). Nevertheless, I keep coming back to old Procrustes, as I did today. I have so far achieved two publications which introduce tweaks of the general procedure, and I'm currently working on a more subtle implementation which will use Procrustes as part of my own inverse dynamics workflow.

One reason might be the seemingly trivial goal of Procrustes Superimposition. When we see two forms, we compare them. Automatically. We compare size, rotation, position and shape. Even kids do that, all the time. Procrustes Superimposition is the quantitative way to compare, and it is surprisingly simple.

I got primed on the matter when preparing my master's thesis in John Nyakatura's lab, working all through the excellent explanations by Norm MacLeod in his PalaeoMaths 101 column. The title was "Who is Procrustes and What Has He Done With My Data?". Then, for application, it started as so often with code copied from a stack overflow response which someone got by inspiration from matlab (wouldn't matlab be better if it were all open source?) and which I've adapted for my purpose.

Initially, I understood close to nothing.

But over time, I (wrongly) adjusted the code, found mistakes, and understood more. Only today, by coincidence, I fully understood the last enigmatic component (the "reflection" issue with SVD).

My personal anecdote should serve as a motivation to students to go out there and use (and thereby gradually understand) the uncountable piles of code they find online. I'll also take the occasion to briefly explain my Procrustes toolbox in the short blog post below.

Enjoy!

# The only library required herein is numpy.
import numpy as NP

# However, plotting can be nice.
import matplotlib.pyplot as PLT

Data Properties

All data are points will be organized in an \(n\times m\) matrix, with n being the number of points, and m the dimension (always 3 with the code below). Mind your dimensions and signs! Some of the matrix operations further below are sensitive and can be found in inverse notation.

The following functions extract features of the data which are relevant.

def Centroid(pointset):
    # the centroid of the data.
    return NP.mean(pointset, axis = 0)
def RCS(pointset):
    # root centroid size, i.e. the scale of the data, defined as the root sum of 
    # euclidean distances of all points from their centroid.
    return NP.sqrt(NP.sum(NP.power(pointset-Centroid(pointset), 2), axis = 0 ))

Transformation Operations

These are rather trivial operations that change attributes of a data set. While doing that, they store the exact modifications in a dictionary (see below).

def Center(pointset, trafo_history = None):
    # move all points so that the centroid is at the origin.
    centroid = Centroid(pointset)

    if trafo_history is not None:
        MakeHistory(trafo_history, translate = -centroid, log = 'center')

    return pointset - centroid
def UnitScale(pointset, trafo_history = None):
    # Scale the points so that RCS is one.
    rcs = RCS(pointset)

    if trafo_history is not None:        
        MakeHistory(trafo_history, scale = 1/rcs, log = 'unitscale')

    return pointset / rcs
def Rotate(pointset, rotation_matrix, trafo_history = None):
    # rotate a set of points
    # Note: usually (but not generally) it makes sense to center the points.
    centroid = Centroid(pointset)

    if trafo_history is not None:
        MakeHistory(trafo_history, rotate = rotation_matrix, log = 'rotate')

    return NP.dot(pointset.copy(), rotation_matrix )

Finally, it is useful to have a combined function that "standardizes" a form.


def Standardize(pointset, trafo_history = None):
    # Standardize a point set, i.e. center and scale to unity.
    pointset = Center(pointset.copy(), trafo_history)
    pointset = UnitScale(pointset, trafo_history)

    return pointset

Making History

While modifying objects, we would like to keep track of the things that happen. Obviously, this facilitates bug fixing.

But what's more important, it allows to exapt a series of transformations to different point sets.

# all starts with an empty history.
def GetEmptyHistory():
    # as history, I use dicts with defined keywords.
    return {key: [] for key in ['rotate', 'scale', 'translate', 'log']}
def MakeHistory(hist, scale = 1., translate = NP.zeros((3,)), rotate = NP.eye(3), log = ''):
    # This convenience function is used to appended actions to a history dict.  
    hist['scale'].append(scale)
    hist['translate'].append(translate)
    hist['rotate'].append(rotate)
    hist['log'].append(log)
def InvertHistory(trafo_history):
    # this inverts the course of actions to exactly undo a series of transformations.
    # (not fully tested)
    inverted_history = GetEmptyHistory()
    for shift, scale, rotation_matrix, log in zip(trafo_history['translate'], trafo_history['scale'], trafo_history['rotate'], trafo_history['log']):
        MakeHistory(inverted_history, scale = 1/scale, translate = -shift, rotate = rotation_matrix.T, log = f'{log} (inverted)')

    return inverted_history

And here is the function that allows to apply a superimposition to a separate point set:

def ApplyTransformations(points, trafo_history):
    # This will apply a previously stored series of transformations to data.
    points = points.copy()
    for shift, scale, rotation_matrix in zip(trafo_history['translate'], trafo_history['scale'], trafo_history['rotate']):
        points = scale * Rotate(points, rotation_matrix) + shift

    return points

Hopefully few conceptual challenges so far. But maybe the procedures, when broken down into simple functions, just "loose their magic" quicker.

Toy Data

I will first go through the procedure step by step, before pasting a single function which you can copy and use elsewhere.

Procrustes Superimposition will remove differences between two objects in position, scale, and rotation, thereby leaving them to differ only in shape.

Historically, Procrustes specialists work on triangles. I prefer smileys.

data = {\
"face1": NP.array([ \
          [14.041144, 41.9160903] \
        , [23.635776, 42.1840953] \
        , [18.129748, 43.5327753] \
        , [17.430506, 46.9459693] \
        , [11.763092, 50.8475533] \
        , [12.2746, 51.6814123] \
        , [13.395186, 52.1052683] \
        , [14.999929, 52.2049983] \
        , [16.963911, 52.0664843] \
        , [19.162211, 51.7756033] \
        , [21.469911, 51.4182353] \
        , [23.762091, 51.0802583] \
        , [25.913831, 50.8475513] \
        ]) \
, "face2": NP.array([ \
          [91.405685, 13.157473] \
        , [101.000317, 12.467501] \
        , [96.732169, 14.142255] \
        , [97.082454, 17.007084] \
        , [89.127633, 19.904976] \
        , [91.119536, 20.980855] \
        , [93.025648, 21.712223] \
        , [94.856997, 22.117637] \
        , [96.624615, 22.215657] \
        , [98.339531, 22.024843] \
        , [100.012776, 21.563752] \
        , [101.655379, 20.850943] \
        , [103.278372, 19.904975] \
        ])
}

data = {k: NP.concatenate([v, NP.linspace(0, (float(k[-1])-1)*0.2, v.shape[0]).reshape([-1,1])], axis = 1) \
       for k, v in data.items()}
data['face1'] *= 1.5
data['face2'] -= 40
data['face1'][:, 1] -= 80

data
{'face1': array([[ 21.061716  , -17.12586455,   0.        ],
        [ 35.453664  , -16.72385705,   0.        ],
        [ 27.194622  , -14.70083705,   0.        ],
        [ 26.145759  ,  -9.58104605,   0.        ],
        [ 17.644638  ,  -3.72867005,   0.        ],
        [ 18.4119    ,  -2.47788155,   0.        ],
        [ 20.092779  ,  -1.84209755,   0.        ],
        [ 22.4998935 ,  -1.69250255,   0.        ],
        [ 25.4458665 ,  -1.90027355,   0.        ],
        [ 28.7433165 ,  -2.33659505,   0.        ],
        [ 32.2048665 ,  -2.87264705,   0.        ],
        [ 35.6431365 ,  -3.37961255,   0.        ],
        [ 38.8707465 ,  -3.72867305,   0.        ]]),
 'face2': array([[ 51.405685  , -26.842527  , -40.        ],
        [ 61.000317  , -27.532499  , -39.98333333],
        [ 56.732169  , -25.857745  , -39.96666667],
        [ 57.082454  , -22.992916  , -39.95      ],
        [ 49.127633  , -20.095024  , -39.93333333],
        [ 51.119536  , -19.019145  , -39.91666667],
        [ 53.025648  , -18.287777  , -39.9       ],
        [ 54.856997  , -17.882363  , -39.88333333],
        [ 56.624615  , -17.784343  , -39.86666667],
        [ 58.339531  , -17.975157  , -39.85      ],
        [ 60.012776  , -18.436248  , -39.83333333],
        [ 61.655379  , -19.149057  , -39.81666667],
        [ 63.278372  , -20.095025  , -39.8       ]])}

The first four points are "eyes" and "nose", the other are mouthpoints. Note that, ideally, points must be biologically homologous, and I will resist my urge to rant about the overuse of "semilandmarks" (here: the ones on the mouth line).

Our goal is to superimpose the second face onto the first.

Instead of calling them "face1" and "face2", or second and first, I prefer calling the reference one "reference". Deep copies with .copy() are used to make sure we don't accidentally change our raw data.

referenceform = data['face1'].copy()
focalform = data['face2'].copy()
def PlotForms(form1, form2):
    for form, color in zip([form1, form2], [(0.2, 0.2, 0.2), (0.2, 0.2, 0.8)]):
        PLT.scatter(form[:2,0], -form[:2,1], s = 320, edgecolor = color, marker = 'o', facecolor = 'none')
        PLT.plot(form[2:4,0], -form[2:4,1], color = color, lw = 2, ls = '-')
        PLT.plot(form[4:,0], -form[4:,1], color = color, lw = 1, ls = '-')

    PLT.gca().set_aspect('equal')
    PLT.show();
PlotForms(referenceform, focalform)

png

You can intuitively see that these "faces" differ in position, size, and in the way their mouths curl. It is harder to guess whether they differ in rotation.

And it is impossible to quantify "by eye" how different they are.

Procrustes Superimposition

Position and Scaling

This is trivial, because we can "center" and "unit scale" objects with the little helper functions above. The combination is sometimes called "standardization" (although a proper standardization should also involve a default rotation).

We will store the raw centroid and RCS for later.

ref_centroid = Centroid(referenceform)
ref_rcs = RCS(referenceform)

Then get both forms standardized.

referenceform = Standardize(referenceform)
focalform = Standardize(focalform)

PlotForms(referenceform, focalform)

png

Voila! Still, rotation might differ (see eyes and noses)...

Rotation

Finding the matchmaking rotation is a bit more tricky. I found this blog post showing the quickest solution and further reading in the Kwon archives (where else?).

Start by calculating a cross dispersion matrix, which is in fact the correlation matrix of the two point sets.

cross_dispersion = NP.dot(referenceform.T, focalform)

Next, perform singular value decomposition of that matrix. What is singular value decomposition? Well, I would argue that SVD is maybe the best reason to ever visit a lecture in linear algebra. It's relevant and has many applications (e.g. Principal Component Analysis). Spend some time to read it up, try it on your favorite matrices, and digest what you need for application.

U,singular_values,V = NP.linalg.svd(cross_dispersion, full_matrices = False)

One can get the rotation matrix R as follows:

rotation_matrix = NP.dot(V.T, U.T)

But note: the SVD can occasionaly yield a "reflected" outcome. I guess this happens if the forms are more similar if one of them is mirrored.

This can be intended in some applications. For me, it lead to nonsense results because mirroring made no sense; I had to debug for two days, which ultimately made me start this text.

You can identify such "reflection" cases by testing if the determinant of the rotation matrix is negative. Then, adjustment works by two sign flips.

# rectify the mirrored outcomes
if NP.linalg.det(rotation_matrix) < 0:
    rotation_matrix = NP.dot(NP.dot( V.T, NP.eye(3) * NP.array([1., 1., -1.]) ), U.T)
    singular_values[-1] *= -1

The singular values can be of further relevance: their magnitude (sum) might be different from \(1\), which must be regarded for scaling. And they are used to calculate procrustes distance.

singular_value_trace = singular_values.sum()

# standarised distance 
procrustes_distance = 1 - singular_value_trace**2

# optimum scaling of data
scaling_factor = singular_value_trace * ref_rcs / RCS(focalform)

Finally, we can apply the transformation. Note that the reference face was not changed in the process (e.g. not centered).

transformed_data =  Rotate(focalform, rotation_matrix) \
                    * scaling_factor \
                  + ref_centroid

PlotForms(data['face1'], transformed_data)

png

Believe it or not: This is the "perfect" superimposition of the two faces - iff you find the maths a perfect criterium, and given that you can ignore the overrepresention of the mouth part (think more about this before sputtering your scans with semi-landmarks).

In a Nutshell

All this, spiced up with the creation of a log of the involved steps, can be put into a function:

def Procrustes(focalform_raw, referenceform_raw, trafo_history = None):
    # http://stackoverflow.com/questions/18925181/procrustes-analysis-with-numpy
    # input forms must be numpy arrays of the shape [points, dimensions]

    ## work on copies, to be safe
    referenceform = referenceform_raw.copy()
    focalform = focalform_raw.copy()

    ## match nr of points and nr of dimensions
    if (not referenceform.shape[0] == focalform.shape[0]) or (not referenceform.shape[1] == focalform.shape[1]):
        raise IOError('forms must be of the same "shape" (Procrustes/numpy joke :)')

    ## store transformations
    if trafo_history is None:
        trafo_history = GetEmptyHistory()

    ## standardize
    referenceform = Standardize(referenceform)
    focalform = Standardize(focalform, trafo_history)

    ## rotation
    # optimum rotation matrix of focal form
    # cf. http://nghiaho.com/?page_id=671

    # calculate cross-dispersion (correlation) matrix (cf. http://www.kwon3d.com/theory/jkinem/rotmat.html)
    cross_dispersion = NP.dot(referenceform.T, focalform)

    # singular value decomposition of cross-dispersion matrix
    U,singular_values,V = NP.linalg.svd(cross_dispersion, full_matrices = False)

    # ... to get the rotation matrix R as follows:
    rotation_matrix = NP.dot(V.T, U.T)

    # if R has negative determinant, it is a reflection matrix and must be modified.
    if NP.linalg.det(rotation_matrix) < 0:
        rotation_matrix = NP.dot(NP.dot( V.T, NP.eye(3) * NP.array([1., 1., -1.]) ), U.T)
        singular_values[-1] *= -1

    ## further Procrustes measures
    singular_value_trace = singular_values.sum()

    # standarised distance 
    procrustes_distance = 1 - singular_value_trace**2

    ## scaling
    # optimum scaling of data
    scaling_factor = singular_value_trace * RCS(referenceform_raw) / RCS(focalform)

    ## transformation
    # ref_rcs = RCS(referenceform_raw)
    ref_centroid = Centroid(referenceform_raw)
    transformed_data =  Rotate(focalform, rotation_matrix) \
                        * scaling_factor \
                      + ref_centroid

    ## log
    MakeHistory(trafo_history \
                , rotate = rotation_matrix \
                , scale = scaling_factor \
                , translate = ref_centroid \
                , log = 'procrustes' \
                )

    # return transformed data, residual distance, and log
    return transformed_data, procrustes_distance, trafo_history

With this at hand, enjoy turning all your forms into shapes!

PLT.scatter(focalform[:2,0], -focalform[:2,1], s = 320, edgecolor = 'k', marker = 'o', facecolor = 'none')
PLT.plot(focalform[2:4,0], -focalform[2:4,1], color = 'k', lw = 2, ls = '-')
PLT.plot(focalform[4:,0], -focalform[4:,1], color = 'k', lw = 1, ls = '-')

PLT.gca().set_aspect('equal')
PLT.gca().get_xaxis().set_visible(False)
PLT.gca().get_yaxis().set_visible(False)
PLT.show();

png