Message Boards Message Boards

GROUPS:

The spherical Lloyd algorithm

Posted 6 months ago
994 Views
|
3 Replies
|
14 Total Likes
|

visualization of Lloyd relaxation with a Voronoi diagram

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]

progress of points being relaxed by the Lloyd algorithm


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]
POSTED BY: J. M.
Answer
3 Replies

Nice! Thanks for sharing!

This is quite useful for many practical problems.

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract