subdivide_to_size() does not split shared edges

[brummieb]

Summary

subdivide_to_size() does not appear to check if an edge being split is shared by another triangle; meaning you can end up with a vertex that, if subsequently moved, will result in tears/holes in a mesh.

This may be related to the "Fix trimesh.remesh.subdivide subdivide_to_size so that if you specify a subset of faces it splits their neighbors in a way to maintain watertightness." comment in #1557

Details

Using Python 3.9.2 and Trimesh 3.15.8 I'm loading a simple STL into a mesh, then subdividing it to make no edge longer than 5 units:

mesh = trimesh.load('part.stl')
new_v, new_f = trimesh.remesh.subdivide_to_size(mesh.vertices, mesh.faces, 5.0)
mesh = trimesh.Trimesh(vertices=new_v, faces=new_f)

Using the viewer I can check the original mesh (without subdivision) in both shaded and wireframe (the first and second images in the attached jpg). The third image shows the same mesh after the subdivide_to_size() call above, and the fourth is a close up of that mesh.

The circled areas show where triangles have been split, but the formerly shared edge with another triangle hasn't. If these vertices are then moved they will leave a hole in the mesh.

Is there any way to get trimesh to subdivide a mesh and ensure that edges shared with other triangles are also split?

[mikedh]

Yeah this is true, currently the only way to get a subdivision to maintain nice properties is vanilla subdivide for every face (i.e. without passing a subset of faces). I couldn't think of an easy way to do it at the time, and I was using subdivide_to_size for voxelization which doesn't need good connectivity. I think you'd have to do the stitching per-iteration? Or maybe you could do one pass using a tree of some sort? Happy to accept PR's!

[brummieb]

I have experimented with something like this before; whereby ever time a triangle is subdivided you search all other triangles to see if any have shared edges (i.e. an edge that is now split) and force the shared triangle to be split. It wasn't a particularly quick algorithm though, and certainly wasn't particularly "pythonic" code.

[johannes-lee]

Here's an implementation using numpy that subdivides the faces with long edges and splits the residual quads into new triangles (along the shortest diagonal) in order to create a consistent mesh (or at least, the consistency is inherited from the inputs). This should also preserve vertex ordering according to the inputted faces (as long as process=False doesn't change this). However, you should expect that some faces have low quality (small angles).

def subdivide_to_size_iter(vertices, faces, max_edge, return_index=False):
    # vertices ((n, 3) float) – Vertices in space
    # faces ((m, 3) int) – Indices of vertices which make up triangles
    # max_edge (float) – Maximum length of any edge in the result
    # return_index (bool) – If True, return index of original face for new faces
    
    mesh = trimesh.Trimesh(vertices, faces, process=False)
    
    # Each face falls into one of 4 cases:
    #   1. Preserved: No edges need to be halved. Preserve this face.
    #   2. Split: One edge needs to be halved. Split this face into 2 new faces along the midpoint and the opposite vertex
    #   3. Triquad: Two edges need to be halved. Create a face from the shared vertex and the two midpoints, resulting in a tri and a quad
    #        Create two triangles from quad along the shortest diagonal
    #   4. Subdivide: Three edges need to be halved. Subdivide normally.
    
    # By construction, all faces must have 3 edges. -> shape (num_faces, 3)
    face_to_edge = np.arange(mesh.edges.shape[0]).reshape((-1, 3))
    face_to_edge_unique = mesh.edges_unique_inverse[face_to_edge] # which unique edge each of the ordered edges in each face corresponds to.
    
    edge_criterion_individual = (mesh.edges_unique_length > max_edge) # whether each edge is too long
    face_criterion = np.any(edge_criterion_individual[face_to_edge_unique], axis=-1) # whether the face has any edges that are too long
    edge_criterion = np.zeros(mesh.edges_unique.shape[0], dtype='bool') # whether each edge should be halved (based on self and adjacent faces)
    edge_criterion[np.unique(face_to_edge_unique[face_criterion])] = True
    face_edge_criterion = edge_criterion[face_to_edge_unique] # shape (num_faces, 3 faces_per_edge). Whether each edge of each face should be halved
    
    
    halved_midpoints = mesh.vertices[mesh.edges_unique[edge_criterion]].mean(1)
    
    # new_midpoint_inverse maps from unique edge index to new_vert idx of its midpoint
    new_midpoint_inverse = np.cumsum(edge_criterion) - 1
    new_midpoint_inverse[~edge_criterion] = len(edge_criterion) # this should an invalid index when referencing and is used as a sanity check
    new_midpoint_inverse += mesh.vertices.shape[0]
    
    
    new_verts = np.concatenate([mesh.vertices, halved_midpoints], axis=0)
    face_split_type = face_edge_criterion.sum(-1)
    
    preserve_inds = np.nonzero(face_split_type == 0)[0]
    split_inds = np.nonzero(face_split_type == 1)[0]
    triquad_inds = np.nonzero(face_split_type == 2)[0]
    subdivide_inds = np.nonzero(face_split_type == 3)[0]
    ########################
    preserve_faces = mesh.faces[preserve_inds]
    ########################
    # either 0, 1, or 2. Referrering to which edge is split in face_edge_criterion
    # split_edge[i] = j means that edge j (vertices split_faces[i, [j, (j+1)%3]]) is split and vertex (j + 2)%3 = 2 is opposite
    split_edge = np.argmax(face_edge_criterion[split_inds], axis=1)
    opposite_vertex = (split_edge + 2) % 3

    # new faces are [opposite vertex, (opp+1)%3, split midpoint]
    #           and [(opp+2)%3, opp, midpoint]
    face_to_edge_unique_split = face_to_edge_unique[split_inds]
    split_faces = np.array([
        [mesh.faces[split_inds, opposite_vertex], mesh.faces[split_inds, (opposite_vertex+1)%3], new_midpoint_inverse[face_to_edge_unique_split[np.arange(len(split_inds)), split_edge]]],
        [mesh.faces[split_inds, (opposite_vertex+2)%3], mesh.faces[split_inds, opposite_vertex], new_midpoint_inverse[face_to_edge_unique_split[np.arange(len(split_inds)), split_edge]]],
    ]).transpose(2, 0, 1).reshape((-1, 3))
    
    ########################
    # either 0, 1, or 2. Referrering to which edge is not split in face_edge_criterion for triquad faces
    # triquad_edge[i] = j means that edge j (vertices split_faces[i, [j, (j+1)%3]]) is split and vertex (j + 2)%3 = 2 is opposite
    triquad_edge = np.argmin(face_edge_criterion[triquad_inds], axis=1)
    triquad_opposite_vertex = (triquad_edge + 2) % 3
    face_to_edge_unique_triquad = face_to_edge_unique[triquad_inds]

    # triquad_tri is the triangle between the shared vertex and the two midpoints
    triquad_tri = np.array([
        [mesh.faces[triquad_inds, triquad_opposite_vertex], new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge-1)%3]], new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge+1)%3]]]
    ]).transpose(2, 0, 1).reshape((-1, 3))
    triquad_edge_diag1 = np.array([new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge-1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex-1)%3]]).T
    triquad_edge_diag2 = np.array([new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge+1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex+1)%3]]).T
    triquad_diag1_length = np.linalg.norm(np.diff(new_verts[triquad_edge_diag1], axis=1), axis=-1).flatten()
    triquad_diag2_length = np.linalg.norm(np.diff(new_verts[triquad_edge_diag2], axis=1), axis=-1).flatten()

    # faces created from drawing edge from midpoint following shared vertex with vertex preceding shared vertex 
    triquad_faces1 = np.array([
        [new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge-1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex+1)%3], mesh.faces[triquad_inds, (triquad_opposite_vertex-1)%3]],
        [new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge-1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex-1)%3], new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge+1)%3]]],
    ]).transpose(2, 0, 1)[triquad_diag1_length <= triquad_diag2_length].reshape((-1, 3))
    # faces created from drawing edge from midpoint preceding shared vertex with vertex following shared vertex
    triquad_faces2 = np.array([
        [new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge+1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex+1)%3], mesh.faces[triquad_inds, (triquad_opposite_vertex-1)%3]],
        [new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge+1)%3]], new_midpoint_inverse[face_to_edge_unique_triquad[np.arange(len(triquad_inds)), (triquad_edge-1)%3]], mesh.faces[triquad_inds, (triquad_opposite_vertex+1)%3], ],
    ]).transpose(2, 0, 1)[triquad_diag1_length >  triquad_diag2_length].reshape((-1, 3))

    triquad_faces = np.concatenate([triquad_tri, triquad_faces1, triquad_faces2], axis=0)
    
    ########################
    # faces from subdivision
    subdivide_faces = np.array([
        [mesh.faces[subdivide_inds, 0], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 0]], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 2]]],
        [mesh.faces[subdivide_inds, 1], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 1]], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 0]]],
        [mesh.faces[subdivide_inds, 2], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 2]], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 1]]],
        [new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 0]], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 1]], new_midpoint_inverse[face_to_edge_unique[subdivide_inds, 2]]],
    ]).transpose(2, 0, 1).reshape((-1, 3))
    ########################
    new_faces = np.concatenate([preserve_faces, split_faces, triquad_faces, subdivide_faces], axis=0)
    
    if return_index:
        index = np.concatenate([
            preserve_inds,
            np.repeat(split_inds, 2),
            triquad_inds,
            np.repeat(triquad_inds[triquad_diag1_length >= triquad_diag2_length], 2),
            np.repeat(triquad_inds[triquad_diag1_length <  triquad_diag2_length], 2),
            np.repeat(subdivide_inds, 4),
        ])
        return new_verts, new_faces, index
    
    return new_verts, new_faces
def subdivide_to_size(vertices, faces, max_edge, max_iter=10, return_index=False):
    # vertices ((n, 3) float) – Vertices in space
    # faces ((m, 3) int) – Indices of vertices which make up triangles
    # max_edge (float) – Maximum length of any edge in the result
    # max_iter (int) – The maximum number of times to run subdivision. A non-positive value will use as many iterations as needed.
    # return_index (bool) – If True, return index of original face for new faces


    max_length = trimesh.Trimesh(vertices, faces, process=False).edges_unique_length.max()
    n_iter = int(np.ceil(np.log2(max_length/max_edge)))
    n_iter = min(max_iter, n_iter) if max_iter > 0 else n_iter
    index_maps = [np.arange(faces.shape[0])]
    for _ in range(n_iter):
        if not return_index:
            vertices, faces = subdivide_to_size_iter(vertices, faces, max_edge)
        else:
            vertices, faces, index_iter = subdivide_to_size_iter(vertices, faces, max_edge, return_index=True)
            index_maps.append(index_iter)
    if not return_index:
        return vertices, faces
    index = index_maps[-1]
    for index_prev in reversed(index_maps[0:-1]):
        index = index_prev[index]
    return vertices, faces, index
def subdivide_to_size_mesh(mesh, max_edge):
    # mesh: trimesh.Trimesh object
    # max_edge (float) – Maximum length of any edge in the result
    new_mesh = trimesh.Trimesh(*subdivide_to_size(mesh.vertices, mesh.faces, max_edge, max_iter=-1), process=False)
    return new_mesh

[brandonrwin]

@johannes-lee, could you add a license to this code? Otherwise, it's very hard to use.