# Finding a visually pleasing rotation of point sets with approximate symmetries

Suppose we have an arrangement of points in 2D that are visually symmetric, and therefore they have a rotation that seems "natural".

For example, consider these points:

pts = {{0, 0}, {2, 0}, {0, 1}, {2, 1}};
rpts = pts.RotationMatrix[RandomReal[2 Pi]];
rpts = # + RandomReal[0.01 {-1, 1}, 2] & /@ rpts;

Graphics[{PointSize[0.05], Point[rpts]}] It would be nicer to draw them like this: In this particular case, this can be achieved with principal component analysis:

Graphics[{PointSize[0.05], Point[PrincipalComponents@rpts]}]


However, this approach fails when the points have an (approximate) rotational symmetry:

pts = {{0, 0}, {1, 0}, {0, 1}, {1, 1}};
rpts = pts.RotationMatrix[RandomReal[2 Pi]];
rpts = # + RandomReal[0.01 {-1, 1}, 2] & /@ rpts;

Graphics[{PointSize[0.05], Point[PrincipalComponents@rpts]}] pts = CirclePoints;
rpts = pts.RotationMatrix[RandomReal[2 Pi]];
rpts = # + RandomReal[0.01 {-1, 1}, 2] & /@ rpts;

Graphics[{PointSize[0.05], Point[PrincipalComponents@rpts]}] Question: What method would work better in these nearly rotationally symmetric cases, while also being able to handle any other case?

Application: Rotating graph layouts obtained with force-directed methods (which often produce symmetric results if the graph has symmetries).

Here's a more complicated point set for testing:

pts = Uncompress[
"1:eJwBUQGu/iFib1JlAgAAABQAAAACAAAA9NeXQqZd6T8vLQSgIiOyP3FgrKVIVum/\
jpcK2AsTs78cwn5lTgbUP1LqwmkTaOe/gfuthKTXxj+NFQSPAcToP9aS8KD+K+O/\
V586UvPa4L/llJnBFw7mv7b6zL0pJ9o/kyJOEEMbxT/RsQ4owyDZv55vXyqkybY/\
qEKnmNK92j9MEeOEUf7XP/XZ8IOQDcy/HL1oVXsV1T8FPrI3VVPSP0CYVbBz/tS/u3P7y/\
f+0b84R93SiCfYv54206OWhMw/fo+y4Bvs2r+2kJ5NEImkv6K1ic6w9eU/9UOI0wkR2r/\
50WjzpEzjP6aCSvzg2+A/7wIJEZ8T2z+M/dGqofKiP4PFpBk7are/\
1B39elLA2r9LFR1DriPFv8BotjoaNNk/LvMBAYhsxr9tb8rXvt7ov4sIM9IORNS/\
zLGHTttd5z9Yyq9Z"]

• in a hurry, my thought is that your points comprise a convex hull object with high rotational symmetry. So, you can obtain the inertia axes of the region, and make them vertical and horizontal... – José Antonio Díaz Navas Feb 16 '18 at 12:52
• @José Antonio Díaz Navas: Your approach yields eigenvectors ~{1,0},{0,1}. Because there are many symmetries rotation doesn't change the inertia matrix. – Ulrich Neumann Feb 16 '18 at 16:08
• I should have noted that the reason why I think this is possible to do is that Mathematica's own force directed graph layouts are almost always rotated in a pleasing way. – Szabolcs Feb 16 '18 at 18:08
• @UlrichNeumann what you say it is not true. The moment of inertia is sensitive to the distribution of the points delimiting the region, so their inertia´s axes (check it !). Anyway, I have checked this is not a good approach, as the inertia axes are not aligned with points of interest in the boundary. – José Antonio Díaz Navas Feb 16 '18 at 18:19
• @José Antonio Díaz Navas: What I want to point out is that the inertia tensor of the given points, with obviously 10 symmetry axes, doesn't have two unique principal axes. Just look at my answer. – Ulrich Neumann Feb 17 '18 at 14:17

You could find the rotation angle that minimizes the height of the bounding box:

thetaOpt = theta /. Last@Minimize[Differences@
MinMax@(pts.RotationMatrix[theta])[[All, 2]], theta]
Graphics[{PointSize[0.05], Point[pts.RotationMatrix[thetaOpt]]}] Another approach using BoundingRegion as suggested by Szabolcs (u is the first direction of the base of the bounding parallelogram):

u = BoundingRegion[pts, "MinOrientedRectangle"][[2, 1]];
pts2 = pts.RotationMatrix[ArcTan[u[]/u[]]]


returns the same result as above:

Graphics[{Gray, BoundingRegion[pts, "MinOrientedRectangle"],
LightGray, BoundingRegion[pts2, "MinOrientedRectangle"],
Red, Point@pts, Blue, Point@pts2}] • BoundingRegion seems to be able to do something similar too. – Szabolcs Feb 16 '18 at 18:56
• interesting solution too, +1 – José Antonio Díaz Navas Feb 16 '18 at 19:05
• @Szabolcs Updated accordingly. It is much faster with BoundingRegion. – anderstood Feb 16 '18 at 19:13

Based on what José mentioned (plus an extra step to tweak the alignment):

alignPoints[pts_List, vec_List] :=
Module[{obj = DelaunayMesh[pts], pts2, line},
pts2 = MeshCoordinates@
TransformedRegion[obj,
Composition[AffineTransform[-Eigenvectors[MomentOfInertia[obj]]],
TranslationTransform[-RegionCentroid[obj]]]];
line = First@MaximalBy[Subsets[pts2, {2}], EuclideanDistance @@ # &];
RotationTransform[{-Subtract @@ line, vec}] /@ pts2
]


The second argument is which "direction" you find more visually appealing to align with:

For

pts = {{0, 0}, {1, 0}, {0, 1}, {1, 1}};
rpts = pts.RotationMatrix[RandomReal[2 Pi]];
rpts = # + RandomReal[0.01 {-1, 1}, 2] & /@ rpts;

GraphicsRow[
Graphics[{PointSize[0.1], Point[#]}] & /@ {rpts,
alignPoints[rpts, {0, 1}]}, 100] GraphicsRow[
Graphics[{PointSize[0.1], Point[#]}] & /@ {rpts,
alignPoints[rpts, {1, 1}]}, 100] pts = CirclePoints;
rpts = pts.RotationMatrix[RandomReal[2 Pi]];
rpts = # + RandomReal[0.01 {-1, 1}, 2] & /@ rpts;

GraphicsRow[
Graphics[{PointSize[0.1], Point[#]}] & /@ {rpts,
alignPoints[rpts, {0, 1}]}, 100] GraphicsRow[
Graphics[{PointSize[0.1], Point[#]}] & /@ {rpts,
alignPoints[rpts, {1, 0}]}, 100] pts = Uncompress[
"1:eJwBUQGu/iFib1JlAgAAABQAAAACAAAA9NeXQqZd6T8vLQSgIiOyP3FgrKVIVum/\
jpcK2AsTs78cwn5lTgbUP1LqwmkTaOe/gfuthKTXxj+NFQSPAcToP9aS8KD+K+O/\
V586UvPa4L/llJnBFw7mv7b6zL0pJ9o/kyJOEEMbxT/RsQ4owyDZv55vXyqkybY/\
qEKnmNK92j9MEeOEUf7XP/XZ8IOQDcy/HL1oVXsV1T8FPrI3VVPSP0CYVbBz/tS/u3P7y/\
f+0b84R93SiCfYv54206OWhMw/fo+y4Bvs2r+2kJ5NEImkv6K1ic6w9eU/9UOI0wkR2r/\
50WjzpEzjP6aCSvzg2+A/7wIJEZ8T2z+M/dGqofKiP4PFpBk7are/\
1B39elLA2r9LFR1DriPFv8BotjoaNNk/LvMBAYhsxr9tb8rXvt7ov4sIM9IORNS/\
zLGHTttd5z9Yyq9Z"]

GraphicsRow[
Graphics[{PointSize[0.1], Point[#]}] & /@ {pts,
alignPoints[pts, {0, 1}]}, 100] • It is nice that someone check some ideas work !! +1 – José Antonio Díaz Navas Feb 16 '18 at 19:04

I have tried another workaround based on ConvexHull. These are the points: Let us obtain its convex hull:

conv = ConvexHullMesh[rpts] Take one of their lines delimiting its boundary, and calculate its angle with the horizontal axis. This will serve to rotate the original points to obtain a more pleased distribution:

orientations = MeshPrimitives[conv, 1];
ang = VectorAngle[Subtract @@ orientations[[1, 1]], {0, 1}];
Graphics[{PointSize[0.02], Point@(rpts.RotationMatrix[ang]),
{PointSize[0.01], Gray,Point@ rpts}}] which seems to be a not too bad workaround either.

• That works for rotational symmetry, but might fail for a parallelogram for instance. – anderstood Feb 16 '18 at 19:15
• @anderstood, I do not agree. I have tested it also, and works. ;)) – José Antonio Díaz Navas Feb 16 '18 at 19:18
• Try rpts = {{0, 0}, {1, 0}, {1.2, .1}, {.2, .1}}. The result will depend on the boundary segment that you choose. (note that I upvoted your answer all the same :)) – anderstood Feb 16 '18 at 19:22
• In this extreme cases, the solution proposed by @Szabolcs by using PrincipalComponents is quite useful. Remember that the premise is a set of points exhibiting some rotational symmetry... Oh, I appreciate your comments, of course... – José Antonio Díaz Navas Feb 16 '18 at 19:27

Perhaps I misunderstood the hint from @José Antonio Díaz Navas but I tried to verify his suggestion. The inertia matrix of the points pts

m = Chop@Mean[pts];
M = Total@Map[Outer[Times, # - m, # - m] - (# - m).(# - m) IdentityMatrix&, pts]


ews = Eigensystem[M]
Graphics[{Point[pts], Red, Line[{m, m + ews[[2, 1]]}],Line[{m, m + ews[[2, 2]]}]}] which unfortunately doesn't show the espected symmetry properties. I believe the reason why this approach doesn't work is because of the 10 symmetriy axes!

appendix Here I want to give an approach, which directly tries to find the mirror symmetry axis depending on \[CurlyPhi]. The underlying idea is to minimize the distance between the points pts and the mirrored points, thereby considering only the Nearest neighbors.

JJ[\[CurlyPhi]_?NumericQ, punkte_] :=
Block[{p\[CurlyPhi], p\[CurlyPhi]S, nb, J},
p\[CurlyPhi] =Map[( RotationMatrix[\[CurlyPhi]].#) &, pts];(*Punkte gedreht...*)
p\[CurlyPhi]S =Map[{-1, 1} ( RotationMatrix[\[CurlyPhi]].#) &,
pts] ;(*...und gespiegelt*)

nb = Flatten[Map[Nearest[p\[CurlyPhi]S, #, 1] &, p\[CurlyPhi]],1]; (* Spiegelnachbar*)
J = Total[(nb - p\[CurlyPhi])^2, 2] (* Symmetrie: J\[Equal]0*)
]


Minimization gives one of the possible solutions

opt = NMinimize[{JJ[\[CurlyPhi], pts] ,0 < \[CurlyPhi] < .5}, \[CurlyPhi]]
Graphics[{Point[pts], Red,Line[{-{-Cos[\[CurlyPhi]], +Sin[\[CurlyPhi]]}, {-Cos[\[CurlyPhi]],Sin[\[CurlyPhi]]}} /. opt[]]}] • I am puzzled, you are working on pts, not on rpts. This later set has no "exact" rotational symmetry. What am I missing? – José Antonio Díaz Navas Feb 17 '18 at 14:33
• A cross section that has more than 2 axis of symmetry has a spherical inertia tensor (all the direction are principal directions) so I don't think it can be used to find the axes of symmetry. – anderstood Feb 17 '18 at 14:43
• @José Antonio Díaz Navas: I took the points given in the final question pts=Uncompress[...] . My argumentation is based on these symmetry properties. – Ulrich Neumann Feb 17 '18 at 15:08
• @anderstood: Thank you for your explanation.With my wording in the comments&answer I intended to focus this issue! – Ulrich Neumann Feb 17 '18 at 15:13
• @UlrichNeumann Yes I meant to confirm that with some mathematical justification. But I'm starting to believe it's maybe a terminology issue: José might be talking about the moment of inertia while you and me are talking of the second moment of area. My previous comment was about second moment of area. – anderstood Feb 17 '18 at 15:39