How can we apply Lloyd's relaxation algorithm to a VoronoiMesh?



Thanks a lot KennyColnago for your answer. Based on Simon Woods suggestion of using PropertyValue[{vmesh, 2}, MeshCellCentroid], I think this code can be simplified a bit, e.g.

VoronoiAdaptive2[p_, iter_] := Block[{vmesh = VoronoiMesh[p]},
vmesh = VoronoiMesh[PropertyValue[{vmesh, 2}, MeshCellCentroid]],
{i, iter}];

Now I have a question about the algorithm itself: According to Wikipedia,

Each time a relaxation step is performed, the points are left in a slightly more even distribution: closely spaced points move farther apart, and widely spaced points move closer together.

But using a large number of iterations with the above code results in a mesh with a concentration of very small cells in the center and large cells on the boundary (omitting the line of KennyColnago's code that excludes cells beyond a certain radius).

For example, with 200 points and the same random seed as KennyColnago and 200 iterations, the following mesh is produced:

VoronoiAdaptive2[RandomReal[{-1, 1}, {200, 2}], 200]

200 iters

Going up to 1200 iterations produces a visually identical image. Is this the expected output? I was hoping to see something more like on this website

  • 4
    $\begingroup$ It would help potential answerers if you included some code to generate a VoronoiMesh to work with, and any progress you have made in solving the problem yourself. Did you try implementing the description of the algorithm in the linked Wikipedia article? $\endgroup$ – Simon Woods Oct 3 '15 at 17:19
  • 1
    $\begingroup$ For a neat application (v10) of adaptive Voronoi meshes to images, see the Help page for VoronoiMesh > Applications > Image Processing. With iteration, the cells become nearly circular, and their sizes adapt locally to the image structure. $\endgroup$ – KennyColnago Dec 25 '15 at 1:35
  • $\begingroup$ You might be interested in section 3 of my post on the Wolfram Community. $\endgroup$ – Silvia Mar 9 '16 at 5:12

To address the update to the question, you can use the second argument of VoronoiMesh to set a rectangular boundary which will let the algorithm converge to a uniform spacing. It looks like the linked animation is also inserting additional points at the centre, or perhaps it starts with a high density of points near the centre of the mesh. Here is something similar - I keep inserting a point at {0,0} until there are 200 points:

vmesh = VoronoiMesh[RandomReal[{-1, 1}, {20, 2}]];
 pts = PropertyValue[{vmesh, 2}, MeshCellCentroid];
 If[Length@pts < 200, AppendTo[pts, {0, 0}]];
 vmesh = VoronoiMesh[pts, {{-1, 1}, {-1, 1}}]]

Here's my take on Lloyd's algorithm. The code I present here should not be too hard to encapsulate into a routine; I have only elected to present it in this way so that I can animate the progress of the relaxation method:

BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
            pts = RandomReal[{-1, 1}, {50, 2}]];

pl = With[{maxit = 50, (* maximum iterations *)
           tol = 0.005 (* distance tolerance *)}, 
          FixedPointList[Function[pts, Block[{cells},
                         cells = MeshPrimitives[VoronoiMesh[pts,
                                                {{-1, 1}, {-1, 1}}], "Faces"];
                         RegionCentroid /@ cells[[SparseArray[
                         Outer[#2 @ #1 &, pts, RegionMember /@ cells, 1],
                         Automatic, False]["NonzeroPositions"][[All, 2]]]]]],
                         pts, maxit, 
                         SameTest -> (Max[MapThread[EuclideanDistance,
                                                    {#1, #2}]] < tol &)]];

MapThread[Show, {VoronoiMesh /@ Rest[pl], 
                 MapThread[Graphics[{AbsolutePointSize[3], Line[{##}],
                                    {Black, Point[#1]}, {Red, Point[#2]}}] &,
                           #] & /@ Partition[pl, 2, 1]}] // ListAnimate

Lloyd relaxation

Note: I elected not to use PropertyValue[{(* mesh *), 2}, MeshCellCentroid], since the centroids are not returned in the same order as the points that generated their corresponding cells, thus necessitating a complicated termination criterion.

In this answer, ilian shows an undocumented function that can be used to simplify the implementation of Lloyd's algorithm. Here's how it goes:

pl = With[{maxit = 50, (* maximum iterations *)
           tol = 0.005 (* distance tolerance *)}, 
          FixedPointList[Function[pts, Block[{cells, ci, vm},
                         vm = VoronoiMesh[pts, {{-1, 1}, {-1, 1}}];
                         cells = MeshPrimitives[vm, "Faces"]; 
                         ci = Region`Mesh`MeshMemberCellIndex[vm];
                         RegionCentroid /@ cells[[ci[pts][[All, -1]]]]]], pts, maxit, 
                         SameTest -> (Max[MapThread[EuclideanDistance,
                                                    {#1, #2}]] < tol &)]];
  • $\begingroup$ This is nice. Thanks for pointing me to this. I still think it's crazy there's no easy way to get corresponding cells for the sites. $\endgroup$ – RunnyKine Mar 10 '16 at 19:35
  • $\begingroup$ It annoys me a lot that I had to complicate my implementation here because of this scrambling. Maybe if enough of us complain… $\endgroup$ – J. M.'s ennui Mar 10 '16 at 19:43

The following imperfect code is what I play with. Given an input set of points p, run VoronoiMesh iter times, each time replacing the points with the centroids of the cells.

VoronoiAdaptive[p_, iter_] :=
   Block[{points = p, vmesh, coord, poly, centroids, 
          error = ConstantArray[0, iter]},
         vmesh = VoronoiMesh[points];
         coord = MeshCoordinates[vmesh];
         poly = MeshCells[vmesh, 2];
         centroids = Map[Mean[coord[[#[[1]]]]] &, poly];
            points = centroids;
            vmesh = VoronoiMesh[points];
            coord = MeshCoordinates[vmesh];
            poly = MeshCells[vmesh, 2];
            centroids = Map[Mean[coord[[#[[1]]]]] &, poly];
            error[[i]] = Mean[Abs[Flatten[points - centroids]]],
            {i, 1, iter}];
         {coord, poly}

Here is a plot of the original mesh for 200 random points. Cells with vertices outside an arbitrary radius of 1.5 are removed.

Block[{p, coord, poly, iter = 0},
   p = RandomReal[{-1, 1}, {20, 2}];
   {coord, poly} = VoronoiAdaptive[p, iter];
   poly = DeleteCases[poly, _?(Max[Map[Norm, coord[[#[[1]]]]]] > 1.5 &)];
      EdgeForm[{Thickness[0.001], White}],
      GraphicsComplex[coord, poly],
      Red, Point[p],
      Frame -> True, Background -> Black]

no iterations

Setting iter=10 gives the following.

10 iterations

  • 1
    $\begingroup$ You could use PropertyValue[{vmesh, 2}, MeshCellCentroid] to get the centroids $\endgroup$ – Simon Woods Oct 3 '15 at 18:33
  • $\begingroup$ Thanks a lot, @KennyColnago. I'm wondering though why using a large number of iterations doesn't produce a very uniform mesh like in the animation here. $\endgroup$ – nasosev Oct 4 '15 at 0:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.