Yet another take at drawing trimesh meshes :)

[villares]

I have tweaked a bit my draw_mesh() function using numpy to draw meshes removing the "coplanar face edges", drawing only "facet boundaries".

import py5, py5_tools
import shapely
import trimesh
import numpy as np

def setup():
    global subtracted_box, crescent_mesh
    py5.size(500, 500, py5.P3D)
    c = shapely.Point(0, 0).buffer(50)  # circle from point, buffer is radius
    d = shapely.affinity.translate(c, 50, 0)  # new translated circle
    crescent = c - d   # same as c.difference(d)
    crescent_mesh = trimesh.creation.extrude_polygon(crescent, 250)  # extrusion
    crescent_mesh.apply_translation((0, 0, -125))  # single argument, tuple
    apply_rotation(crescent_mesh, py5.PI/2, direction=(0, 1, 0))  # my rotation helper
    hole = trimesh_box(0, 0, 0, 180, 300, 180)
    square_tube = trimesh_box(0, 0, 0, 200).difference(hole)  
    plus = trimesh_box(0, 0, 0, 100, 50, 300).union(trimesh_box(0, 0, 0, 50, 100, 300))
    subtracted_box = square_tube.difference(crescent_mesh).difference(plus)
    # to export gif
    py5_tools.animated_gif('out.gif', duration=0.05, frame_numbers=range(1, 361, 3))
   
def draw():
    py5.background(0, 100, 100)
    py5.lights()
    py5.translate(py5.width / 2, py5.height / 2)
    py5.rotate_x(py5.PI / 8)
    py5.rotate_y(py5.radians(py5.frame_count))
    py5.fill(200, 200, 0)
    draw_mesh(subtracted_box)
    #py5.shape(subtracted_box)  # for comparison
    #py5.fill(0, 0, 100)       # to debug
    #draw_mesh(crescent_mesh)

def draw_mesh(m):
    vs = m.vertices
    bs = m.facets_boundary
    # draw the triangulated faces without stroke
    py5.push_style()  # this is so I can revert to previous stroke settings
    py5.no_stroke()
    with py5.begin_closed_shape(py5.TRIANGLES):
        py5.vertices(vs[np.concatenate(m.faces)])
    py5.pop_style()
    # now draw just the facet boundary lines
    a, b = np.vstack(bs).T
    py5.lines(np.column_stack((vs[a], vs[b])))
  
def apply_rotation(obj, angle, direction=[1, 0, 0], center=[0, 0, 0]):
    rot_matrix = trimesh.transformations.rotation_matrix(angle, direction, center)
    obj.apply_transform(rot_matrix)

def trimesh_box(x, y, z, w, h=None, d=None):
    h = h or w
    d = d or h
    c = trimesh.creation.box((w, h, d))
    c.apply_translation((x, y, z))  # in place!
    return c

py5.run_sketch(block=False)

For comparison, with all edges visible:

[hx2A]

This is nice, and I like your implementation better than the one I came up with. This is simpler.

[villares]

I found a problem with this strategy: It looks like single triangle faces are not considered "facets" and don't have their edges listed in mesh.facets_boundary. 😢

I'd love to have a (simple) solution for this...

import py5
import py5_tools 
import trimesh
import numpy as np

def setup():
    global a, b, c, d, e
    py5.size(500, 500, py5.P3D)
    a = box_mesh(0, 0, 0, 200)
    b = a.copy()
    apply_rotation(b, py5.radians(45), direction=(0, 0, 1))
    c = a.copy()
    apply_rotation(c, py5.radians(45), direction=(0, 1, 0))
    apply_rotation(c, py5.radians(45), direction=(1, 0, 0))
    d = c.intersection(b)
    # Export animation, skip empty frame 0 created on setup
    py5_tools.animated_gif('demo_mesh.gif', duration=0.08,
                           frame_numbers=range(1, 361, 5)) 

def draw():
    py5.background(200, 0, 200)
    py5.lights()
    py5.translate(py5.width / 2, py5.height / 2)
    py5.rotate_x(py5.PI / 8)
    py5.rotate_y(py5.radians(py5.frame_count))
    py5.stroke_weight(2)
    py5.fill(0, 200, 200)
    py5.stroke(0)
    draw_mesh(d)
    if py5.is_key_pressed:
        py5.scale(1.01)
        py5.no_fill()
        py5.stroke(255)
        draw_mesh(b)
        draw_mesh(c)
    
def box_mesh(x, y, z, w, h=None, d=None):
    h = h or w
    d = d or h
    mesh = trimesh.creation.box((w, h, d))
    mesh.apply_translation((x, y, z))
    return mesh
    
def apply_rotation(mesh, angle, direction=[1, 0, 0], center=[0, 0, 0]):
    rot_matrix = trimesh.transformations.rotation_matrix(angle, direction, center)
    mesh.apply_transform(rot_matrix)
    
def draw_mesh(m):
    vs = m.vertices
    bs = m.facets_boundary
    py5.push_style()  # To be able to get the stroke settings back
    py5.no_stroke()   # Turn off strokes, to draw without edges
    # Draw triangulated faces without edges
    with py5.begin_closed_shape(py5.TRIANGLES):
        py5.vertices(vs[np.concatenate(m.faces)])
    py5.pop_style() # Get previous stroke settings
    # Draw only edges that are facet boundaries
    a, b = np.vstack(bs).T
    py5.lines(np.column_stack((vs[a], vs[b])))
    
py5.run_sketch()

[hx2A]

I found a problem with this strategy: It looks like single triangle faces are not considered "facets" and don't have their edges listed in mesh.facets_boundary. 😢

I've faced this problem also. The way to address it is to get the list of faces in each facet and subtract that from all of the faces. The remaining faces are leftovers and are single face facets.

[villares]

Perhaps something likenp.setdiff1d(mesh.faces, np.concatenate(mesh.facets)) ?

[hx2A]

Perhaps something likenp.setdiff1d(mesh.faces, np.concatenate(mesh.facets)) ?

Yes, I think that could work and would be fast and efficient

[villares]

Another rabbit hole: are you de-reduplicating edges for better pen-plotter output?

[hx2A]

Yes, if you are making SVG output. A tool like vpype can also do that for you.

[villares]

I'm trying coplanar_edges = mesh.face_adjacency_edges[mesh.face_adjacency_angles < tolerance] but it has a few problems...

• It seems to miss some edges
• The inverse doesn't give you all non-coplanar edges...
• It is hard to combine this with the other provided "edge collections", like using it to filter out co-planar edges.

[villares]

Experimenting with this, assuming edges have the indices ordered...:

coplanar_edges = mesh.face_adjacency_edges[mesh.face_adjacency_angles < tol]
mask = ~np.any(np.all(coplanar_edges[:, None] == mesh.edges_unique, axis=2), axis=0)
a, b = np.vstack(mesh.edges_unique[mask]).T
py5.lines(np.column_stack((vs[a], vs[b])))

Otherwise:

coplanar_edges = np.sort(mesh.face_adjacency_edges[mesh.face_adjacency_angles < tol], axis=1)
mask = ~np.any(np.all(coplanar_edges[:, None] == mesh.edges_sorted, axis=2), axis=0)
a, b = np.vstack(mesh.edges_sorted[mask]).T
py5.lines(np.column_stack((vs[a], vs[b])))

But I see no difference between the visual results of either.

My problem is that in this example I still get one co-planar edge shown:
https://github.com/villares/sketch-a-day/tree/main/2025/sketch_2025_08_03

[GitHub has been acting up for me, this page never ends loading so the JS that enables the edit & reply functions are unavailable]