pythonmeshtrimesh

How do I use Python Trimesh to get boundary vertex indices?

How do I use the Trimesh Python library to retrieve the indices of the vertices that make up the boundary of a mesh?

For example, for a planar mesh, I expect only the vertices that are on the outer edges. If the planar mesh has a hole, I also expect the vertices that mark the edges of the hole. For an open cylindrical mesh, I expect only the vertices that line the two openings.

All of my meshes are open, like pipes or boxes without tops and bottoms. They are not watertight.

For a given mesh instance, Neither the edges (which returns more entries than I have vertices!), nor the edges_unique properties return what I expect. The facets_boundary works for planar meshes, but fails spectacularly for the cylindrical mesh. In reading the project API documentation, I find it difficult to understand what I should expect from these properties.

Do I have to find the edges myself, e.g. using vertex_counts = np.bincount(faces.flatten()) ?

For my meshes, this produces results as follows:

Planar mesh (4x5)
vertex_counts:  [2 3 3 3 3 1 3 6 6 6 6 3 3 6 6 6 6 3 3 6 6 6 6 3 1 3 3 3 3 2]

Planar mesh with hole (8x8 with a 3x3 hole in the middle)
vertex_counts:  [2 3 3 3 3 3 3 3 1 3 6 6 6 6 6 6 6 3 3 6 4 3 3 3 5 6 3 3 6 3 0 0 0 3 6 3 3
 6 3 0 0 0 3 6 3 3 6 3 0 0 0 3 6 3 3 6 5 3 3 3 4 6 3 3 6 6 6 6 6 6 6 3 1 3
 3 3 3 3 3 3 2]

Cylindrical mesh (8 circular divisions, 6 longitudinal divisions)
vertex_counts:  [3 3 3 3 3 3 3 3 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
 6 6 6 6 6 6 6 6 6 6 6 3 3 3 3 3 3 3 3]
Solution

  • For a given mesh instance, neither the edges (which returns more entries than I have vertices!), nor the edges_unique properties return what I expect.

    Let's start from scratch. A (triangular) mesh is a list of 3D coordinates (the vertices) plus a list of triplets of vertices (the triangles). Let's consider a super-simple mesh, representing a square:

    0 --- 1
| \   |
|   \ |
2 --- 3

    The vertices 0, 1, 2 and 3 are connected thanks to two triangles: [0, 2, 3] and [0, 3, 1]. In this context, the edges of a mesh are the individual segments in all triangles, that is, the list of couples (0, 2), (2, 3), (3,0), (0, 3), (3, 1) and (1, 0). As you can see, the diagonal edge is repeated since it appears in both triangles. The edges_unique property makes sure that this repetition is not included. In general, for a triangular mesh with $V$ vertices you can expect to have roughly $2V$ faces and $3V$ edges (see here).

    The facets_boundary works for planar meshes, but fails spectacularly for the cylindrical mesh. In reading the project API documentation, I find it difficult to understand what I should expect from these properties.

    From the documentation: "Return a list of face indices for coplanar adjacent faces". I agree that it is not super clear, but what it does is to group faces (triangles) that are both adjacent and co-planar. Take as an example a cylinder. The top and bottom circles are both made of a set of triangles that touch each other and lie on the same planes. The two sets of triangles are grouped as facets. Similarly, a box has six facets (though each facet could be divided into many triangles!). Another way of thinking about it: a facet is a group of triangles that could be replaced by a planar polygon.

    Now, back to your original question: how to find the boundary? The answer is to keep all edges that belong to only one triangle. If an unordered list of edges is enough for you, you can do it simply with:

    unique_edges = mesh.edges[trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)]

    (see this GitHub issue). If you want 3D coordinates, consider using Trimesh.outline() instead.

    If you need an ordered list of vertices, then it's a bit more involved. First of all, this makes sense only if your mesh is a manifold (no edge can belong to more than two triangles). Then, as explained in this answer, after getting the list of boundary edges you need to traverse them in order. Here is a function that can do so:

    def boundary(mesh, close_paths=True):
    # Set of all edges and of boundary edges (those that appear only once).
    edge_set = set()
    boundary_edges = set()

    # Iterate over all edges, as tuples in the form (i, j) (sorted with i < j to remove ambiguities).
    # For each edge, three cases are possible:
    # 1. The edge has never been visited before. In this case, we can add it to the edge set and as a boundary
    #    candidate as well.
    # 2. The edge had already been visited once. We want to keep it into the set of all edges but remove it from the
    #    boundary set.
    # 3. The edge had already been visited at least twice. This is generally an indication that there is an issue with
    #    the mesh. More precisely, it is not a manifold, and boundaries are not closed-loops.
    for e in map(tuple, mesh.edges_sorted):
        if e not in edge_set:
            edge_set.add(e)
            boundary_edges.add(e)
        elif e in boundary_edges:
            boundary_edges.remove(e)
        else:
            raise RuntimeError(f"The mesh is not a manifold: edge {e} appears more than twice.")

    # Given all boundary vertices, we create a simple dictionary that tells who are their neighbours.
    neighbours = defaultdict(lambda: [])
    for v1, v2 in boundary_edges:
        neighbours[v1].append(v2)
        neighbours[v2].append(v1)

    # We now look for all boundary paths by "extracting" one loop at a time. After obtaining a path, we remove its edges
    # from the "boundary_edges" set. The algorithm terminates when all edges have been used.
    boundary_paths = []

    while len(boundary_edges) > 0:
        # Given the set of remaining boundary edges, get one of them and use it to start the current boundary path.
        # In the sequel, v_previous and v_current represent the edge that we are currently processing.
        v_previous, v_current = next(iter(boundary_edges))
        boundary_vertices = [v_previous]

        # Keep iterating until we close the current boundary curve (the "next" vertex is the same as the first one).
        while v_current != boundary_vertices[0]:
            # We grow the path by adding the vertex "v_current".
            boundary_vertices.append(v_current)

            # We now check which is the next vertex to visit.
            v1, v2 = neighbours[v_current]
            if v1 != v_previous:
                v_current, v_previous = v1, v_current
            elif v2 != v_previous:
                v_current, v_previous = v2, v_current
            else:
                # This line should be un-reachable. I am keeping it only to detect bugs in case I made a mistake when
                # designing the algorithm.
                raise RuntimeError(f"Next vertices to visit ({v1=}, {v2=}) are both equal to {v_previous=}.")

        # Close the path (by repeating the first vertex) if needed.
        if close_paths:
            boundary_vertices.append(boundary_vertices[0])

        # "Convert" the vertices from indices to actual Cartesian coordinates.
        boundary_paths.append(mesh.vertices[boundary_vertices])

        # Remove all boundary edges that were added to the last path.
        boundary_edges = set(e for e in boundary_edges if e[0] not in boundary_vertices)

    # Return the list of boundary paths.
    return boundary_paths

    Note that it returns a list of 3D paths. If you want list of vertices, just change the line boundary_paths.append(mesh.vertices[boundary_vertices]) to boundary_paths.append(boundary_vertices). Also note that if the parameter close_paths is True, the first element of a path is repeated at the end as well to "make it" a closed path.

    Finally, there may be ways to do the same that are much more efficient than what I just wrote!