I have developed some code for VoronoiMesh objects, and I am able to find the polygon with the highest number of adjacent polygons. How do I then find the indices of the adjacent polygons?

  • 3
    $\begingroup$ Please post your code. $\endgroup$ Jan 29, 2016 at 21:47

4 Answers 4


Finding the points which correspond to the Voronoi cells with maximum number of neighboring cells is easy by sorting the VertexDegree on underlying DelanunayMesh.

Lets find the indexes of the pts with maximum vertex degree which assures that they will be surrounded by most number of Voronoi cells.

pts = RandomReal[{-1, 1}, {250, 2}];
Vmesh = VoronoiMesh[pts];
graph = With[{m = DelaunayMesh[pts]}, 
    MeshCells[m, 1] /. Line[{start_, end_}] -> {start, end}, 
    VertexCoordinates -> MeshCoordinates[m], 
    GraphLayout -> "PlanarEmbedding", 
    EdgeStyle -> Directive[[email protected], Dashed, White], 
    VertexStyle -> White, VertexSize -> 3]];
highestdeg = Flatten@Position[VertexDegree[graph], VertexDegree[graph] // Max]

{76, 91, 247}

So we have our VoronoiMesh and the point indices highestdeg whose neighboring cells we will search with the following function. This function will return the cell index which contains the point with given index and the cell indices of it's neighboring mesh cells.

FindAdjPoly[ptIndex_, {pts_, Vmesh_}] := 
 Block[{pt, mpt, regs, lines, poly, loc, all},
  mpt = MeshCoordinates[Vmesh];
  (* Get the mesh polygons *)
  regs = MeshCells[Vmesh, 2];
  pt = pts[[ptIndex]];
  (* Select the polygon containing pt *)
  poly = SelectFirst[(regs /. Polygon[a__] :> {a, Polygon@mpt[[a]]}), 
    RegionMember[Last@#, pt] &];
  (* Find the index of the polygon in mesh *)
  loc = MeshCellIndex[Vmesh, {Polygon[First@poly]}][[1, 2]];
  (*Get the edges of the polygon *)
  lines = Partition[First@poly, 2, 1, 1];
  (* Find polygons that share the above edges *)
  all = DeleteDuplicates@
    Flatten[With[{ed = #}, 
         regs /. Polygon[a__] :> (MemberQ[a, #] & /@ ed), {True,True}]] & /@ lines];
  {loc, Complement[all, {loc}]}

Calling the function:

polys = FindAdjPoly[#, {pts, Vmesh}] & /@ highestdeg

{{248, {11, 28, 35, 99, 149, 198, 202, 210, 234, 244}}, {249, {8, 14, 16, 59, 62, 106, 119, 169, 186, 188}}, {250, {75, 77, 80, 84, 85, 138, 147, 172, 176, 192}}}

Why not check if everything is working as expected... enter image description here

Code to visualize:

  Table[Prepend[#, 2] & /@ Transpose@{polys[[i, 2]]}, {i,Length@polys}]], graph, 
 Graphics[{White, Disk[#, 0.023] & /@ pts[[highestdeg]], Red,
 [email protected], Point /@ pts[[highestdeg]]}]]
  • $\begingroup$ Hi, Thanks a lot for the code. It's a very nice code. But, when I try to visualize them, I found an error. May be because of version problem. However, I have a question. Why you use Table and transpose those values. And, It would be great for me if you describe a bit of the "FindAdjPoly" in the code. I want to learn how you do it? Thanks $\endgroup$
    – Odrisso
    Jan 31, 2016 at 0:20
  • $\begingroup$ @Odrisso I have added comments to the code for FindAdjPoly. Hope this helps you. About the visualization I am working on version 10.3.1 and I see no errors. The Table and Transpose is used to process the output polys so that HighlightMesh can make sense of it. $\endgroup$ Jan 31, 2016 at 11:03

Here are the points at the centers of the Voronoi cells:

myCenters = RandomReal[1, {100, 2}];

This creates a DelanunayMesh of the center points, i.e., the graph linking centers that is mathematically dual to the Voronoi cell representation:

myDelaunayMesh =  DelaunayMesh[myCenters];

This extracts the links that include the most surrounded point:

myLinkstomostSurroundedPoint = (Last@
      Select[myLinks, (#[[1, 1]] == myCenters[[i]] || #[[1, 2]] == 
           myCenters[[i]])  &], {i, Length[myCenters]}]) /. 
  Line[{a__, b__}] -> {a, b};

This finds the center that appears the most often in that list, i.e., is the most-surrounded cell center:

mostSurroundedPoint = 
 Commonest@Flatten[myLinkstomostSurroundedPoint, 1][[1]];

Here are the other points in the list, i.e., the "surrounding" points:

theSurroundingPoints = 
    Flatten[myLinkstomostSurroundedPoint, 1]), {mostSurroundedPoint}];

Here is the position in the list of the most surrounded point:

Position[myCenters, mostSurroundedPoint[[1]]]




Here are the indexes of the other, surrounding points:

Flatten@(Position[myCenters, #] & /@ theSurroundingPoints)


{78, 80, 43, 8, 86, 10, 84, 30, 46}


  • $\begingroup$ Hi. Thanks. But,I still found problem. Here is my code: v = VoronoiMesh[mydata]; p1 = MeshPrimitives[v, 2] dv = Length[p1] n = Length /@ Level[pv, {2}] sv = Max[n] nv = Count[n, sv] i1 = Level[Position[n, sv], {2}] In this way, I determine the region with the highest number of adjacent regions. But, I can't determine the adjacent regions now. Please let me know about it. I run your code. But, can't understand it. $\endgroup$
    – Odrisso
    Jan 29, 2016 at 23:50
  • $\begingroup$ The adjacent regions have as centers given by the last line in my code. $\endgroup$ Jan 30, 2016 at 0:03
  • $\begingroup$ When I use the last line, it choose all the polygons. Not the surrounding ones. I think I have to Set up a For loop to find the Union of neighbouring adjacent polygons. Can you please kindly tell me about it. $\endgroup$
    – Odrisso
    Jan 30, 2016 at 0:37
  • $\begingroup$ The last line gives the indexes in myCenters of the neighbors. So the two dimensional points are myCenters[[78]], myCenters[[80]], etc. Clearly not "all" the polygons. $\endgroup$ Jan 30, 2016 at 0:47

You can use properties of a MeshRegion object to find adjacent faces/meshes. The only problem is that MeshRegion properties return face indices, and unfortunately, face labels don't match up with face indices. So, here are two functions to convert between face labels and face indices:

fromLabel[vm_, labels_] := Ordering[Length /@ vm["Faces"]][[labels]]
toLabel[vm_, indices_] := ReplaceAll[
    Thread[Ordering[Length /@ vm["Faces"]]->Range[Length[vm["Faces"]]]]

Then, the following will highlight adjacent faces:

highlightAdjacentFaces[vm_, label_] := HighlightMesh[
    Thread[{2, toLabel[vm, vm["FaceFaceConnectivity"][[fromLabel[vm, label]]]]}]


pts = RandomReal[10, {250, 2}];

vm = VoronoiMesh[pts];
highlightAdjacentFaces[vm, 200]
highlightAdjacentFaces[vm, {100, 150, 200}]

enter image description here

enter image description here


It has never been easier as with the new IGMeshCellAdjacencyGraph tool in IGraph/M, (version 0.3.93 or later), which is based on this answer. The following method works for an arbitrary MeshRegion.

Here is a mesh:

pts = RandomReal[{-1, 1}, {250, 2}];
Vmesh = VoronoiMesh[pts];

This is how you find a polygon with the maxmimal number of neighbors:

g = IGMeshCellAdjacencyGraph[Vmesh, 2];
  Position[VertexDegree[g], Max[VertexDegree[g]]]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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