3
$\begingroup$

Consider the following VoronoiMesh obtained after a few iterations of Lloyd's relaxation algorithm

L1 = 6; L2 = 10;
relaxed = 
  Nest[PropertyValue[{VoronoiMesh[#, {{-1, L2 + 2}, {-1, 
         L1 + 2}}], {2, All}}, 
     MeshCellCentroid] &, {RandomReal[L2 + 2, (L1 + 2) (L2 + 2)], 
     RandomReal[L1 + 2, (L1 + 2) (L2 + 2)]} // Transpose, 200];
mesh0 = VoronoiMesh[relaxed]

enter image description here

Now, inspired by the discussion in this question, I can first pick the "nicely shaped" interior cells, based on the cell area

mesh1 = MeshRegion[MeshCoordinates[mesh0], 
  With[{a = PropertyValue[{mesh0, 2}, MeshCellMeasure]}, 
   With[{m = 1.5}, Pick[MeshCells[mesh0, 2], UnitStep[a - m], 0]]]]

enter image description here

My goal is now to be able to easily manage the boundary (edges and corners independently) and interior, including/removing them in a "toggleable" way. This would define 9 regions in total (4 corners, 4 edges and 1 interior region). I'm already able to select either the boundary or the interior, as can be seen here

MeshRegion[MeshCoordinates[mesh1], MeshCells[mesh1, {2, "Frontier"}]]

enter image description here

and

MeshRegion[MeshCoordinates[mesh1], MeshCells[mesh1, {2, "Interior"}]]

enter image description here

Is there a good way to detect edges and/or corners of meshes like this? I'm aware of commands like IntersectingQ and Pick, that could prove to be useful. Any ideas?

Edit: For example, I would like to obtain meshes like the following (with the interior also "toggable")

enter image description here

Also, by edges I mean any of the four cell arrays that constitute the frontier/boundary of the mesh. Corners might be trickier to uniquely define, given that the mesh is slightly perturbed, but a solution for the edges would be great nonetheless.

$\endgroup$
3
$\begingroup$

We can identify the four corners by finding the coordinates of bounding rectangle of mesh1

cbounds = Tuples @ RegionBounds[mesh1];

and using the function Region`Mesh`MeshNearestCellIndex

fourcorners = Region`Mesh`MeshNearestCellIndex[mesh1] /@ cbounds]

To partition the frontier faces into 8 parts (four corners and four edges) we need to do some sorting:

ClearAll[sortCounterClockwise]
sortCounterClockwise = SortBy[#2, 
    Function[x, ArcTan @@ (PropertyValue[{#, x}, MeshCellCentroid] - RegionCentroid[#])]] &;

sortedCorners = sortCounterClockwise[mesh1, fourcorners];

interiorCells = MeshCellIndex[mesh1, {2, "Interior"}];

sortedFrontier = sortCounterClockwise[mesh1, MeshCellIndex[mesh1, {2, "Frontier"}]];

sortedFrontier = RotateLeft[sortedFrontier, 
   -1 + Flatten[Position[sortedFrontier, First @ sortedCorners]]];

ClearAll[cornerCells, edgeCells]
Evaluate[Riffle[Array[cornerCells, 4], Array[edgeCells, 4]]] = 
  TakeList[sortedFrontier, 
   Riffle[-1 + Differences@ 
    Append[Flatten[Position[sortedFrontier, #] & /@ sortedCorners], 
       1 + Length@sortedFrontier], 1, {1, -2, 2}]];

Style each corner and edge with a different color and add legend:

cornerColors = {Red, Blue, Green, Cyan};
edgeColors = ColorData[97] /@ Range[4];

legend = SwatchLegend[Join[cornerColors, edgeColors, {LightBlue}], 
   Join[Subscript["corner", #] & /@ Range[4], 
     Subscript["edge", #] & /@ Range[4], {"interior"}], 
   LegendMarkerSize -> {20, 20}];


Legended[HighlightMesh[mesh1, 
  Join[MapThread[Style, {cornerCells /@ Range[4], cornerColors}], 
   MapThread[Style, {edgeCells /@ Range[4], edgeColors}], 
   {Style[interiorCells, LightBlue]}], 
 ImageSize -> Large], 
 Placed[legend, Right]]

enter image description here

Update: Use FlipView to toggle visibility of face groups:

ClearAll[fliP]
fliP = FlipView[{Style[#, #2], Style[#, Opacity[0]]}] &; 

cornerPrims = MeshPrimitives[mesh1, #] & /@ cornerCells /@ Range[4];
edgePrims = MeshPrimitives[mesh1, #] & /@ edgeCells /@ Range[4];
interiorPrims = MeshPrimitives[mesh1, #] & /@ interiorCells;

Graphics[Join[MapThread[fliP, {cornerPrims, cornerColors}], 
   MapThread[fliP, {edgePrims, edgeColors}],
   {fliP[interiorPrims, LightBlue]}],
   ImageSize -> Large] // Deploy

enter image description here

| improve this answer | |
$\endgroup$
4
$\begingroup$

You could use ToBoundaryMesh to extract the points and edges like so:

Needs["NDSolve`FEM`"]
L1 = 6; L2 = 10;
relaxed = 
  Nest[PropertyValue[{VoronoiMesh[#, {{-1, L2 + 2}, {-1, 
         L1 + 2}}], {2, All}}, 
     MeshCellCentroid] &, {RandomReal[L2 + 2, (L1 + 2) (L2 + 2)], 
     RandomReal[L1 + 2, (L1 + 2) (L2 + 2)]} // Transpose, 200];
mesh0 = VoronoiMesh[relaxed]
mesh1 = MeshRegion[MeshCoordinates[mesh0], 
  With[{a = PropertyValue[{mesh0, 2}, MeshCellMeasure]}, 
   With[{m = 1.5}, Pick[MeshCells[mesh0, 2], UnitStep[a - m], 0]]]]
mr1 = MeshRegion[MeshCoordinates[mesh1], 
  MeshCells[mesh1, {2, "Frontier"}]]
mr2 = MeshRegion[MeshCoordinates[mesh1], 
  MeshCells[mesh1, {2, "Interior"}]]
(b0 = ToBoundaryMesh[RegionBoundary[mesh1], 
    "MeshOrder" -> 1])["Wireframe"]
(b1 = ToBoundaryMesh[RegionBoundary[mr1], 
    "MeshOrder" -> 1])["Wireframe"]
(b2 = ToBoundaryMesh[RegionBoundary[mr2], 
    "MeshOrder" -> 1])["Wireframe"]
lp1 = ListPlot[{b0["Coordinates"], b2["Coordinates"]}, 
  PlotLegends -> {"Outer", "Inner"}]
fc0 = Flatten@
    FindCycle@Graph[ElementIncidents[b0["BoundaryElements"]][[1]]] /. 
   UndirectedEdge -> List;
fc2 = Flatten@
    FindCycle@Graph[ElementIncidents[b2["BoundaryElements"]][[1]]] /. 
   UndirectedEdge -> List;
Graphics[{Red, GraphicsComplex[b0["Coordinates"], Line[fc0]], Blue, 
  GraphicsComplex[b2["Coordinates"], Line[fc2]]}]

Extract Points and Edges

Update Due To Confusion About Edges

When the term "edges" was used, I presumed that the 1d edges were desired to be extracted. The OP was edited with a new clarified definition. The following approach might be helpful.

First, I extracted the connectivity from the polygons to construct a graph and then used VertexDegree to find the nodes with the highest degree of sharing.

mc1 = MeshCells[mr1, 2];
conn = Delete[0] /@ MeshCells[mr1, 2];
edges = MapThread[List, {#, RotateLeft[#]}] & /@ conn;
g = Graph[Flatten[edges, 1]];
HighlightVertexDegree[g_, vd_] := 
  HighlightGraph[g, 
   Table[Style[VertexList[g][[i]], 
     ColorData["TemperatureMap"][vd[[i]]/Max[vd]]], {i, 
     VertexCount[g]}]];
vd = VertexDegree[g];
HighlightVertexDegree[g, vd]

Vertex Degree

In my case, VertexDegree identified three corner point that had a high degree of sharing with other Polygons.

The following code extracts the high degree vertices and finds the polygons containing those vertices. Those polygons are the corners and the complement are the "edges".

threenodes = Flatten@Position[vd, n_ /; n > 4];
cas = ContainsAny[#] & /@ conn;
polyswiththree = Through[cas[{#}]] & /@ threenodes;
notpolyswiththree = Map[Not, polyswiththree, 2] /. Not -> Identity;
poly = Pick[mc1, #] & /@ polyswiththree;
notpoly = Pick[mc1, #] & /@ notpolyswiththree;
Graphics[{Blue, 
  GraphicsComplex[MeshCoordinates[mr1], 
   Complement[Flatten@Union@notpoly, Flatten@poly]]}]
Graphics[MapIndexed[{Hue[#2/4], 
    GraphicsComplex[MeshCoordinates[mr1], #1]} &, poly]]

Edges and Corners

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you! However, I'm not entirely following your construction. Maybe I wasn't clear enough. I've edited my question to exemplify the sort of meshes I'm looking for. Please take a look. $\endgroup$ – sam wolfe Jan 14 at 13:28
  • 1
    $\begingroup$ I was thinking of "edges" in the typical 1d sense found in FEM, so it is different from the way you are thinking edges. I will keep my answer up for a while in case I think of anything, but I may delete it. $\endgroup$ – Tim Laska Jan 14 at 14:57
  • $\begingroup$ Ah! Now I see what you are doing. I will clarify the definition I'm using. Nonetheless, this is indeed a very neat construction, could be very useful for me later, but please do as you like. $\endgroup$ – sam wolfe Jan 14 at 15:06
  • 1
    $\begingroup$ I updated the post with my understanding of edges and corners. Perhaps it is helpful. There are probably more streamlined approaches, but it seems to work. $\endgroup$ – Tim Laska Jan 15 at 4:40

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.