## Procrustes Analysis

- Tags:
python
# Procrustes Analysis

When trying to find the location of an object on an image, one method matches points of an object that you have found with a template of points that sufficiently define the object. This post is a short description of a method for matching called the Procrustes distance using Python.

The Procrustes distance is named in honor of a story from Greek mythology. Procrustes was a son of Poseidon, a smith and bandit who lived in a cave. As a smith, he had made an iron bed. He would invite travelers to stay with him, and travelers could sleep in the bed.

However, Procrustes wanted the visitor to be the exact length of the bed. Fortunately, being a resourceful smith, he made tools to stretch the hapless visitors limbs if they were too short, or he would lop off the extra length of legs if they were too long.

Apparently, there was also a possibility that he actually had two beds, so that he could select the bed that was most ill-fitting. This does fly in the face of my OCD theory. In any event, Theseus eventually came along and put Procrustes out of his misery.

However, we are left with a motif of a process which attempts to match an object by adjusting the location, size, and orientation of two objects such that they are superimposed. Once superimposed, a degree of difference or distance between the two objects is calculated to help decide if there is a match.

This description will cover a 2D form of Procrustes. We will implement this step by step.

First, some housekeeping needs to be set up. We will import the tools we need from the math and numpy modules.

`from math import cos, sin, atan from numpy import array, dot`

The functions that we set up will act on a numpy 2-column array. There will be one array called the template, which we want to match, and another called points, which is the array to be tested.

`if __name__ == '__main__': # Define the template, just a triangle template = array([[-60., 0.], [60., 20.], [60., -20.], [-60., 0.]]) # Points that are tested points = array([[50., 170.], [70., 50.], [30., 50.], [50., 170.]]) # This is will be the function that we will call print 'procrustes distance:', procrustes(template, points)`

## Superimposition

The first step to superimpostion entails transferring the set of points that comprise the object for comparison as well as the template we are attempting to match to a zero origin.

`def translate(points): """ This function translates the points to center around the origin. """ return points - sum(points) / points.shape[0]`

The next step scales the set of points to a common size. When both the points and template are a common size, then any deviations is shape become apparent.

`def scale(points): """ This function scales the points. Assumes that points are already centered around the origin. """ scale = sum(pow(sum(pow(points, 2.0) / float(points.shape[0])), .5)) return points / scale`

Once both the template and points are a common size, the only difference aside from deviations in shape is the rotation. This next function calculates the angle that best minimizes the differences.

Note that this function assumes that there is a one-to-one corrspondence between points in the template and the collection of points to be tested.

`def get_rotation(template, points): numerator = sum(points[:, 0] * template[:, 1] - \ points[:, 1] * template[:, 0]) divisor = sum(points[:, 0] * template[:, 0] + \ points[:, 1] * template[:, 1]) # Avoiding dividing by zero if divisor == 0.0: divisor = 0.00000000001 return atan(numerator / divisor)`

Now, the only piece left is a function that rotates the points to the angle that represents the best fit.

`def rotate(points, theta, center_point=(0, 0)): """ Rotates the points around the center point. """ new_array = array(points) new_array[0, :] -= center_point[0] new_array[1, :] -= center_point[1] new_array = dot(rotation_matrix_2D(theta), new_array.transpose()).transpose() new_array[0, :] += center_point[0] new_array[1, :] += center_point[1] return new_array def rotation_matrix_2D(theta): return array([[cos(theta), -sin(theta)], [sin(theta), cos(theta)]])`

At this point, we have all the elements necessary for superimpostion. We can translate to a zero origin, scale, discover the best angle, and rotate. The last step entails calculating the distance between the two sets of adjusted point.

`def procrustes(template, points): """ This function computes the minizimed distance between the template and a test set of point once the test set has been scaled, translated, and rotated to the best fit. """ tmp_points = scale(translate(points)) tmp_template = scale(translate(template)) theta = get_rotation(tmp_template, tmp_points) r_points = rotate(tmp_points, theta) return pow(pow(tmp_template - r_points, 2.0).sum(), .5)`

At this point, we have all the elements to successfully calculate the Procrustes distance. So, to repeat from above:

`if __name__ == '__main__': # Define the template, just a triangle template = array([[-60., 0.], [60., 20.], [60., -20.], [-60., 0.]]) # Points that are tested points = array([[50., 170.], [70., 50.], [30., 50.], [50., 170.]]) # This is will be the function that we will call print 'procrustes distance:', procrustes(template, points)`

If we run this, we get:

`$ python procrustes.py procrustes distance: 6.01373026076e-12`

In this case, it is, aside from rounding, zero.