How to find symmetry axis of set of 2D-points?

Hi,

My program generates sets of points with highly visible symmetry like this image bellow.

I would like to get symmetry axis of this set in some form (for example I would like to get coefficients of its equation y = a*x + b) and then rotate the data so the axis would be horizontal. Rotation part is easy but I really can't wrap my head around getting those coefficients.

Big thanks for anyone that can help.

[Torsten]

Maybe minimzing the sum of orthogonal distances of all points to the line ?

[John D'Errico]

@Torsten - effectively correct. Just use an SVD, which allows you to compute that orthogonal regression. But you don't know which of the lines identified is "the" axis, since the axis of symmetry need not be the one associated with the smallest or largest singular value.

Answer 1

[John D'Errico]

Interesting question. If there is a true exact line, one that forms a perfect line of symmetry, then the SVD should identify it. (I'm pretty sure this is true but I've not proven it. That should not be difficult though if my gut is correct.) On your data, for example:

load Data_Points.mat
mu = mean(pts,1);
[U,S,V] = svd(pts - mu,0);

Now I'll plot two lines, in this case, I'll just use a tool I wrote recently to plot a half plane on an existing figure.

plot(pts(:,1),pts(:,2),'o')
H1 = plot2dHalfPlane(V(:,1),dot(V(:,1),mu),xlim,ylim);
H2 = plot2dHalfPlane(V(:,2),dot(V(:,2),mu),xlim,ylim);
axis equal
axis tight
H1.FaceColor = 'r';
H2.FaceColor = 'none';

Effectively, there were two lines defined. One of them, in this case, the line that defines the red half plane is the required axis of symmetry. It need not always be the first singular vector though. The line of symmetry here is:

syms x y
vpa(dot(V(:,1),[x,y]) == dot(V(:,1),mu),16)

ans = 

If the data is not truly symmetrical around any line at all, then the above scheme using an SVD will not be exactly correct, but in that case there is no exactly correct solution.

Finally, in some cases there may be TWO lines of symmetry. Of course the SVD will identify them both in that case. (And, yes, I know, your followup question will be how to identify which of those perpendicular lines is the line of symmetry? The simple answer should be to try reflecting the data around each of the lines identified. The one that reproduces the original data perfectly will tell you which is the line of symmetry, or if both axes were perfect lines of symmetry.)

Lets try it out now, to see if we were successful. First, if the first vector defines the line of symmetry:

ref1 = mu + (pts-mu)*V*diag([-1 1])*V';
plot(pts(:,1),pts(:,2),'bo')
hold on
plot2dHalfPlane(V(:,1),dot(V(:,1),mu),xlim,ylim)
plot(ref1(:,1),ref1(:,2),'r+')
hold off
axis equal
axis tight

Instead, try reflecting around the alternative line.

ref2 = mu + (pts-mu)*V*diag([1 -1])*V';
plot(pts(:,1),pts(:,2),'bo')
hold on
plot(ref2(:,1),ref2(:,2),'r+')
plot2dHalfPlane(V(:,2),dot(V(:,2),mu),xlim,ylim)
hold off
axis equal
axis tight

Clearly it is the first one that is the line of symmetry.

[Fryderyk Kukowski]

Thank you very much!

I managed to set up conversion from line parameters to roation and all is well:

"And, yes, I know, your followup question will be how to identify which of those perpendicular lines is the line of symmetry? The simple answer it to try reflecting the data around each of the lines identified. The one that reproduces the original data perfectly will tell you which is the line of symmetry [...]"

I solved this problem by comparing bounding boxes (so smallest rectangular area circumscribed on the data) of the data after both symmetries to original bounding box. Then I chose that symmetry that had less bounding box error and that's it. It is not super robust but so far after testing it on several other sets of data it worked perfectly every time.

Again, thank you very much!

[John D'Errico]

I'm happy to have helped. I have to admit, all I can think of is Rorschach test when I look at that data. :)

[Fryderyk Kukowski]

About those comparison plot that you've just added...

Very nice, thank you!

But unfortunately I need a programmatical way to determine which line is the one because I'm going to set up automatic processing (so automatic reorientation of the data) of many, many (probably a few thousands) similar sets of points that first part of my program generates. The bounding box method worked for me just fine. But if you have a better idea that is just as fast (this method takes 0.000730 seconds) I'd be glad to read about it.

Cheers!

[Fryderyk Kukowski]

"Rorschach test"

Right, thanks, I'll look into it.

[John D'Errico]

My thought was to use knnsearch.

load Data_Points.mat
mu = mean(pts,1);
[U,S,V] = svd(pts - mu,0);

ref1 = mu + (pts-mu)*V*diag([-1 1])*V';
ref2 = mu + (pts-mu)*V*diag([1 -1])*V';

The problem is, each of the points in the original set should have a companion after reflection that is close. But you don't know which point. knnsearch is very fast though.

[~,D1] = knnsearch(pts,ref1);
[~,D2] = knnsearch(pts,ref2);
norm(D1)

ans = 1.3862e-13

norm(D2)

ans = 4.4693

The choice with the smallest norm, should, if there is a viable axis of symmetry identified, be the one you want. It seems pretty clear here since norm(D1) is virtually zero.

timeit(@() knnsearch(pts,ref1))

ans = 4.5982e-04

That seems decent to me in terms of time. It should also be a little more robust than a bounding box. I could be wrong though.

[Fryderyk Kukowski]

That's pretty clever. Thanks, I'll check it out.