Given a set of points, how do I approximate the major axis of its shape?

Question

Given a "shape" drawn by the user, I would like to "normalize" it so they all have similar size and orientation. What we have is a set of points. I can approximate the size using bounding box or circle, but the orientation is a bit more tricky.

The right way to do it, I think, is to calculate the majoraxis of its bounding ellipse. To do that you need to calculate the eigenvector of the covariance matrix. Doing so likely will be way too complicated for my need, since I am looking for some good-enough estimate. Picking min, max, and 20 random points could be some starter. Is there an easy way to approximate this?

Edit: I found Power method to iteratively approximate eigenvector. Wikipedia article. So far I am liking David's answer.

Thrum

Answer 1

You'd be calculating the eigenvectors of a 2x2 matrix, which can be done with a few simple formulas, so it's not that complicated. In pseudocode:

// sums are over all points
b = -(sum(x * x) - sum(y * y)) / (2 * sum(x * y))
evec1_x = b + sqrt(b ** 2 + 1)
evec1_y = 1
evec2_x = b - sqrt(b ** 2 + 1)
evec2_y = 1

You could even do this by summing over only some of the points to get an estimate, if you expect that your chosen subset of points would be representative of the full set.

Edit: I think x and y must be translated to zero-mean, i.e. subtract mean from all x, y first (eed3si9n).

Redman

If this calculates the eigenvector of the covariance matrix, it's great. Is there a link that points to this method? – Thrum

Check out number-none.com/product/My%20Friend,%20the%20Covariance%20Body/… and the sample code at gdmag.com/code.htm (sep02.zip) – Stair

I found a quick cheat for 2x2 matrix's eigenvector: tinyurl.com/cd998p. For C = |E(xx) E(xy)||E(xy) E(yy)|, L1-d/c comes out to be (sum(xx) - sum(yy))/(2 * sum(xy))+sqrt(((sum(xx) - sum(yy))/(2 * sum(xy)))^2 + 1), which is exactly b + sqrt(b^2 + 1). – Thrum

As I added to the answer, I think x and y need to be adjusted to zero-mean. – Thrum

You're right, they do... I guess I was implicitly assuming that when I wrote out the formulas. – Redman

Answer 2

Here's a thought... What if you performed a linear regression on the points and used the slope of the resulting line? If not all of the points, at least a sample of them.

The r^2 value would also give you information about the general shape. The closer to 0, the more circular/uniform the shape is (circle/square). The closer to 1, the more stretched out the shape is (oval/rectangle).

Frasco

Answer 3

The ultimate solution to this problem is running PCA
I wish I could find a nice little implementation for you to refer to...

Quest

I was going to say PCA. This is the ultimate correct answer. – Nnw

Answer 4

Here you go! (assuming x is a nx2 vector)

def majAxis(x):
    e,v = np.linalg.eig(np.cov(x.T)); return v[:,np.argmax(e)]

Frequentative

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Elyseelysee