# 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[6];
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... Feb 16, 2018 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. Feb 16, 2018 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. Feb 16, 2018 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. Feb 16, 2018 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. Feb 17, 2018 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[[2]]/u[[1]]]]


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. Feb 16, 2018 at 18:56
• interesting solution too, +1 Feb 16, 2018 at 19:05
• @Szabolcs Updated accordingly. It is much faster with BoundingRegion. Feb 16, 2018 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[6];
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 Feb 16, 2018 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. Feb 16, 2018 at 19:15
• @anderstood, I do not agree. I have tested it also, and works. ;)) Feb 16, 2018 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 :)) Feb 16, 2018 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... Feb 16, 2018 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[2]&, 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[[2]]]}]


• I am puzzled, you are working on pts, not on rpts. This later set has no "exact" rotational symmetry. What am I missing? Feb 17, 2018 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. Feb 17, 2018 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. Feb 17, 2018 at 15:08
• @anderstood: Thank you for your explanation.With my wording in the comments&answer I intended to focus this issue! Feb 17, 2018 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. Feb 17, 2018 at 15:39