14
|
12742 Views
|
3 Replies
|
15 Total Likes
View groups...
Share
GROUPS:

# The spherical Lloyd algorithm

Posted 6 years ago
 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 = InternalEffectivePrecision[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 = GeometricFunctionsBinarySearch[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] 
3 Replies
Sort By:
Posted 6 years ago
 - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!
Posted 6 years ago
 This is quite useful for many practical problems.
Posted 6 years ago
 Nice! Thanks for sharing!