Lloyd's algorithm is a method that repeatedly generates the Voronoi diagram of a given set of points, and then replaces each point with the centroid of its corresponding Voronoi cell. It is useful for evenly redistributing points around a rectangular region.
Code for the usual Euclidean version can be found in this Mathematica Stack Exchange post.
In this post, I will be demonstrating the spherical version of Lloyd's algorithm.
The code has its roots in a Wolfram Demonstration written by Maxim Rytin for rendering a spherical Voronoi diagram. I had modernized it in this Mathematica Stack Exchange post, and then had the notion that this code could be extended to implement the spherical Lloyd algorithm.
First, we need a few auxiliary routines. Older versions of Mathematica had an unstable implementation of VectorAngle[], so I wrote a stable version:
vecang[v1_?VectorQ, v2_?VectorQ] := Module[{n1 = Norm[v1], n2 = Norm[v2]},
2 ArcTan[Norm[v1 n2 + n1 v2], Norm[v1 n2 - n1 v2]]]
(You should be able to replace this with VectorAngle[] from version 11.2 onwards.)
Next, I needed a faster version of RotationMatrix[] for rotating one vector into another, so I wrote a routine for doing that, too:
vectorRotate[vv1_?VectorQ, vv2_?VectorQ] :=
Module[{v1 = Normalize[vv1], v2 = Normalize[vv2], c, d, d1, d2, t1, t2},
d = v1.v2;
If[TrueQ[Chop[1 + d] == 0],
c = UnitVector[3, First[Ordering[Abs[v1], 1]]];
t1 = c - v1; t2 = c - v2; d1 = t1.t1; d2 = t2.t2;
IdentityMatrix[3] - 2 (Outer[Times, t2, t2]/d2 -
2 t2.t1 Outer[Times, t2, t1]/(d2 d1) + Outer[Times, t1, t1]/d1),
c = Cross[v1, v2];
d IdentityMatrix[3] + Outer[Times, c, c]/(1 + d) - LeviCivitaTensor[3].c]]
These two are needed for computing the "spherical centroid" of a spherical polygon, as defined by Buss and Fillmore:
(* exponential map for sphere *)
sphereExp[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] + 1 :=
With[{n = Norm[p]}, vectorRotate[{0, 0, 1}, q].Append[p Sinc[n], Cos[n]]]
(* inverse of exponential map for sphere *)
sphereLog[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] :=
Most[vectorRotate[q, {0, 0, 1}].p]/Sinc[vecang[p, q]]
SphericalPolygonCentroid[pts_?MatrixQ] := Module[{k = 0, n = Length[pts], cp, h, pr},
cp = Normalize[Total[pts]];
pr = Internal`EffectivePrecision[pts];
While[cp = sphereExp[cp, h = Sum[sphereLog[cp, p], {p, pts}]/n];
k++; Norm[h] > 10^(-pr/2) && k <= 30];
cp]
Before I show the algorithm, let us first generate some starting points that will be "relaxed" by Lloyd's algorithm:
BlockRandom[SeedRandom[1337, Method -> "MersenneTwister"];
points = {2 ? #1, ArcCos[2 #2 - 1]} & @@@ RandomReal[1, {50, 2}]]
(* convert to Cartesian *)
sp = Function[{u, v}, {Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]}] @@@ points;
The code that implements Lloyd's algorithm is now remarkably compact:
With[{maxit = 40, (* maximum iterations *)
tol = 0.002 (* distance tolerance *)},
fr = FixedPointList[Function[pts,
Block[{ch, polys, verts, vor},
ch = ConvexHullMesh[pts];
verts = MeshCoordinates[ch];
polys = First /@ MeshCells[ch, 2];
vor = Normalize[Cross[verts[[#2]] - verts[[#1]],
verts[[#3]] - verts[[#1]]]] & @@@ polys;
SphericalPolygonCentroid[vor[[#]]] & /@ ch["VertexFaceConnectivity"]]],
sp, maxit,
SameTest -> (Max[MapThread[vecang, {#1, #2}]] < tol &)]];
I used FixedPointList[] so that I can visualize the progress of the algorithm in what follows. One can just use FixedPoint[], of course, if only the final result is wanted.
Now, we can visualize how the points move:
frames = Graphics3D[{{Opacity[.75], Sphere[]}, {Green, Sphere[#, 0.02]}}, Boxed -> False] & /@ fr;
ListAnimate[frames]
One might instead want to see how the Voronoi cells themselves shift around as Lloyd's algorithm proceeds. Visualizing that will require a bit more work.
First, we need a way to render individual spherical polygons, corresponding to the cells of a spherical Voronoi diagram. In this Mathematica Stack Exchange post, I gave routines for rendering a spherical polygon, which I reproduce here along with its auxiliary routines:
(* slerp for two points *)
slerp = Compile[{{q1, _Real, 1}, {q2, _Real, 1}, {f, _Real}},
Module[{n1 = Norm[q1], n2 = Norm[q2], omega, so},
omega = 2 ArcTan[Norm[q1 n2 + n1 q2], Norm[q1 n2 - n1 q2]];
If[Chop[so = Sin[omega]] == 0, q1, Sin[{1 - f, f} omega].{q1, q2}/so]]];
(* stripped down version of functions in
https://mathematica.stackexchange.com/a/10385 *)
sphericalLinearInterpolation[pts_?MatrixQ] := Module[{times},
times = Accumulate[Prepend[vecang @@@ Partition[pts, 2, 1, 1], 0]];
{Last[times], sphericalInterpolatingFunction[times, pts]}]
sphericalInterpolatingFunction[times_, vecs_][t_?NumericQ] :=
With[{k = GeometricFunctions`BinarySearch[times, t]},
slerp[vecs[[k]], vecs[[k + 1]], Rescale[t, times[[{k, k + 1}]], {0, 1}]]]
SphericalPolygon[pts_?MatrixQ, p_: 8] := SphericalPolygon[pts, {p, p}]
SphericalPolygon[pts_?MatrixQ, {p_, q_}] := Module[{ch, cp, en, pt, ql, sp, spl},
cp = SphericalPolygonCentroid[pts];
{en, sp} = sphericalLinearInterpolation[ArrayPad[pts, {{0, 1}}, "Periodic"]];
ch = Sin[? Range[p]/(2 p)]^2; (* rescaled Chebyshev points *)
spl = Most[First[Cases[ParametricPlot3D[sp[t], {t, 0, en}, PlotPoints -> q],
Line[l_] :> l, ?]]];
ql = Length[spl];
pt = Prepend[Apply[Join, Outer[Normalize[#1.{cp, #2}] &,
Transpose[{Append[Reverse[Most[ch]], 0], ch}],
spl, 1]], cp];
GraphicsComplex[pt, {EdgeForm[],
Polygon[PadLeft[Partition[Range[ql] + 1, 2, 1], {Automatic, 3}, 1]
~Join~
Flatten[Apply[Join[Reverse[#1], #2] &,
Partition[Partition[Range[p ql] + 1, ql],
{2, 2}, {1, 1},
{{1, 1}, {-1, 1}}],
{2}], 1]]}, VertexNormals -> pt]]
Here is a routine for generating the spherical Voronoi diagram itself (cf. the Lloyd implementation given above):
sphericalVoronoi[pts_?MatrixQ] := Module[{ch, polys, verts, vor},
ch = ConvexHullMesh[pts];
verts = MeshCoordinates[ch];
polys = First /@ MeshCells[ch, 2];
vor = Normalize /@ (Cross[verts[[#2]] - verts[[#1]],
verts[[#3]] - verts[[#1]]] & @@@ polys);
{vor, ch["VertexFaceConnectivity"]}]
The following code will then generate the cartoon at the beginning of this post:
frames = Function[pts,
Graphics3D[{MapIndexed[{ColorData[99] @@ #2, SphericalPolygon[#1, 12]} &,
(Function[idx, #1[[idx]]] /@ #2) & @@ sphericalVoronoi[pts]],
{Green, Sphere[pts, 0.02]}},
Boxed -> False, Lighting -> "Neutral"]] /@ fr;
ListAnimate[frames]
