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

56 次查看（过去 30 天）
Fryderyk Kukowski
Fryderyk Kukowski 2023-1-13
评论： Fryderyk Kukowski 2023-1-14
Hi,
My program generates sets of points with highly visible symmetry like this image bellow.
Data Points
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.
  2 个评论
Torsten
Torsten 2023-1-13
Maybe minimzing the sum of orthogonal distances of all points to the line ?
John D'Errico
John D'Errico 2023-1-13
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2023-1-13
编辑：John D'Errico 2023-1-13
Ran in:
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.
  6 个评论
显示 4更早的评论
John D'Errico
John D'Errico 2023-1-13
Ran in:
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by