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]