I have a MeshRegion that defines a 2D surface embedded in 3D. Is there a function that returns all the boundaries of the mesh? More specifically, I need for each boundary loop, a list of its vertices (or edges) ordered along the boundary. The ordering is important, i.e. subsequent vertices (or edges) in the list should be connected in the mesh.

I have a rather large obj file, but lets say I had the MeshRegion defined below that gives half of a sphere (I can imagine more elegant code to get this mesh). Then I would like to extract all the vertices on the boundary (in this case only one loop), but in order of traversal along the boundary.

s = BoundaryDiscretizeGraphics[Graphics3D[{Ball[]}]];
pos = Flatten[Position[MeshCoordinates[s], {x_, _, _} /; x < 0]];
f = DeleteCases[MeshCells[s, 2][[All, 1]], 
x_ /; Intersection[pos, x] != {}];
mesh = MeshRegion[MeshCoordinates[s], Polygon /@ f]
  • $\begingroup$ I added some more detail to the question $\endgroup$ – Beginner Jul 21 '17 at 16:27

Here is an attempt.

What distinguishes an interior vertex from an exterior vertex? For an interior vertex, if you take all the neighboring vertices and the edges connecting them, they form a cycle (with the original vertex at its center). So we can use NeighborhoodGraph and FindCycle to do the work here.

First convert your mesh region to a Graph

graph = Graph[Range@Length@MeshCoordinates@mesh, 
  MeshCells[mesh, 1] /. Line[{a_, b_}] :> UndirectedEdge[a, b]]

enter image description here

Now a utility for finding exterior vertices

exteriorVertexQ[graph_, vertex_] := SameQ[{},

and a check to see that it works

exteriorPoints = Select[VertexList[graph], exteriorVertexQ[graph, #] &];

enter image description here

One way to make sure these points are in order, is to again use FindCycle

orderedExteriorPoints = Subgraph[graph,
       ] // FindCycle // First // 
    ReplaceAll[UndirectedEdge -> Sequence] // DeleteDuplicates;
GraphicsComplex[MeshCoordinates@mesh, Line@orderedExteriorPoints] // Graphics3D

enter image description here

You can visualize the result in your original mesh as well

boundaryEdges = Position[
    MeshCells[mesh, 1],
    Line[{Alternatives @@ exteriorPoints, 
      Alternatives @@ exteriorPoints}]
    ] // Flatten;

HighlightMesh[mesh, {1, boundaryEdges}]

enter image description here


1. "ConnectivityMatrix"

We can process the array returned by mesh["ConnectivityMatrix"[1, 2]] (a SparseArray where entry $ij$ is 1 if the 1-dimensional element with index $i$ is connected to the 2-dimensional element with index $j$) to select rows with a single non-zero entry (i.e. edges connected to a single face):

boundaryedgeindices = Flatten @ Position[
    Length /@ mesh["ConnectivityMatrix"[1, 2]]["AdjacencyLists"], 1];

HighlightMesh[mesh, Style[{1, boundaryedgeindices}, Thick, Red]]

enter image description here

To get the boundary lines, use boundaryedgeindices with MeshPrimitives:

boundarylines = MeshPrimitives[mesh, {1, boundaryedgeindices}];


enter image description here

2. "EdgeFaceConnectivityRules"

If you create the mesh region in an alternative way, you can also use the property "EdgeFaceConnectivityRules" to get the indices of boundary edges:

mesh2 = DiscretizeGraphics @
   Select[And @@ NonNegative[#[[1, All, 1]]] &] @ MeshPrimitives[s, 2];

be2 = Keys @ Select[#[[1]] == 0 &] @  Association[mesh2["EdgeFaceConnectivityRules"]];

be2 == boundaryedgeindices
HighlightMesh[mesh2, Style[{1, be2}, Thick, Red]]

enter image description here


Here is a procedure to obtain the oriented boundary from your mesh region: the idea is to find out how many neighbours each point of the mesh has and retain only those points whose number of neighbours is the less. These will give you the boundary. Start by computing the number of neighbours for each point:

  MeshPrimitives[mesh, 2];
  Level[%, {2}];
  Flatten[%, 1];
  neighbours =Tally@%;

The output is a list of the points of the mesh and next to each point we have the number of its neighbours. Next we select those points with less than 5 neighbours. In this particular mesh these are the boundary points:

  boundarypoints = First /@Select[neighbours, Last@# < 5 &]

We confirm that they are actually the boundary points:

  Graphics3D[Point /@ boundarypoints, Axes -> True, AxesLabel -> {"x", "y", "z"}]

enter image description here

Finally we need to sort these points to get the oriented boundary. A look at the graph suggests the following procedure:

   boundarypoints /. {x_, y_, z_} :> {{x, y, z}, ArcTan[z, y]};
   SortBy[%, Last];
   sortedpoints = First /@ %;

We close the oriented curve to get the loop:

   AppendTo[sortedpoints, First@sortedpoints];

enter image description here

Perhaps this procedure can be generalized to other more complicated meshes.

  • $\begingroup$ Thanks, but this solution will only work for the simplistic mesh I provided as an example. In general, I need a procedure that works for arbitrary meshes, i.e. where boundary points can have any number of neighbors and where the sorting cannot be done based on 3D coordinate positions, but needs to be based on the connectivity. I was thinking of a solution that is based on the connectivity graph, e.g. converts the boundary points to a graph based on the connectivity and then applies some loop search. But I couldn't find an efficient or elegant solution yet. $\endgroup$ – Beginner Jul 23 '17 at 20:40
  • $\begingroup$ To think of a more general solution for your problem it would be useful if you could post an example with a more complex mesh whose oriented boundary you wish to compute. $\endgroup$ – Alfonso Jul 24 '17 at 17:02
  • $\begingroup$ It's an obj file, how can you post this here? $\endgroup$ – Beginner Jul 25 '17 at 18:01
  • $\begingroup$ I think you need to provide a link where the .obj file can be downloaded. $\endgroup$ – Alfonso Jul 27 '17 at 7:17

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.