Topological modifications in half-edge meshes

We often need to not only move the vertices of a mesh, but modify the connectivity. In half-edge meshes, there are several “elementary” mesh modifications. For triangulax, by far the most important one is the edge flip (see below). It is the only modification that preserves the number of all mesh elements, and is thus relatively easy to make compatible with JAX and differentiable programming.

Design note: for JIT-compatibility, none of the topology modification functions (flip_edge, collapse_edge, split_vertex) check in advance whether they will produce a valid mesh. Separate predicate functions (can_flip_edge, can_collapse_edge, can_split_vertex) are provided for this purpose: call them before the modification if you need to guard against invalid operations.

Edge flips / T1s

In our simulations, cells will exchange neighbors (T1-event). In the triangulation, this corresponds to an edge flip. We now implement the edge flip algorithm for HeMeshes. We basically edit the various connectivity arrays (in a JAX-compatible way).

The algorithm (and the naming conventions in flip_edge) are from here.

Before

image.png

After

image.png

---

source

flip_edge

def flip_edge(
    hemesh:HeMesh, e:Int[Array, ''], check_boundary:bool=False
)->HeMesh:

Flip half-edge e in a half-edge mesh.

See https://jerryyin.info/geometry-processing-algorithms/half-edge/. The algorithm is slightly modified since we keep track of the origin and destination of a half-edge, and use arrays instead of pointers. Returns a new HeMesh, does not modify in-place.

Does not check whether the flip produces a valid mesh. Use can_flip_edge to check first.

---

source

can_flip_edge

def can_flip_edge(
    hemesh:HeMesh, e:Int[Array, '']
)->Bool[Array, '']:

Check whether flipping half-edge e would produce a valid mesh.

An edge can be flipped if it is interior (not boundary) and the two opposite vertices are not already connected (which would create a duplicate edge).

mesh = TriMesh.read_obj("../test_meshes/disk.obj")
hemesh = HeMesh.from_triangles(mesh.vertices.shape[0], mesh.faces)
geommesh = GeomMesh(*hemesh.n_items, mesh.vertices, mesh.face_positions)

Warning: readOBJ() ignored non-comment line 3:
  o flat_tri_ecmc

plt.triplot(*geommesh.vertices.T, hemesh.faces)
ax = plt.gca()
p = msh.cellplot(hemesh, geommesh.face_positions,
                 cell_colors=np.array([0,0,0,0.1]), mpl_polygon_kwargs={"lw": 1, "ec": "k"})
plt.gca().add_collection(p)

plt.axis("equal")

(np.float64(-1.10003475),
 np.float64(1.09628575),
 np.float64(-1.09934025),
 np.float64(1.09050125))

# flip edge and recompute face positions

flipped_hemesh = flip_edge(hemesh, e=335)
flipped_geommesh = geom.set_voronoi_face_positions(geommesh, flipped_hemesh)

# connectivity is still valid

igl.is_edge_manifold(hemesh.faces)[0], igl.is_edge_manifold(flipped_hemesh.faces)[0], flipped_hemesh.iterate_around_vertex(100)

(True, True, Array([298, 299, 630, 632], dtype=int64))

# you can see the flipped edge between vertices 126-117 in the plot below (middle right)

fig = plt.figure(figsize=(8,8))

plt.triplot(*geommesh.vertices.T, hemesh.faces)
plt.triplot(*flipped_geommesh.vertices.T, flipped_hemesh.faces)

ax = plt.gca()
p1 = msh.cellplot(hemesh, geommesh.face_positions,
         cell_colors=np.array([0.,0.,0.,0.]), mpl_polygon_kwargs={"lw": 1, "ec": "k"})
p2 = msh.cellplot(flipped_hemesh, flipped_geommesh.face_positions,
              cell_colors=np.array([0.,0.,0.,0.]), mpl_polygon_kwargs={"lw": 1, "ec": "tab:orange"})
ax.add_collection(p1)
ax.add_collection(p2)
plt.axis("equal")

msh.label_plot(geommesh.vertices, hemesh.faces, fontsize=10, face_labels=False)

Repeated flips

In a simulation, we need to carry out edge flips at every time step. The function flip_edge does a single edge flip by modifying the connectivity arrays, and is already JIT-compatible.

To carry out multiple flips, we must do them in sequence (otherwise, you risk leaving the mesh in an invalid state). The simplest approach is flip_all, which does a jax.lax.scan over all half-edges. This is JIT-compatible because the scan length is fixed (= number of half-edges), but can be slow for large meshes since it visits every edge even if only a few need flipping.

A more efficient alternative is flip_n_shortest, which sorts edges by length, selects the max_flips shortest candidates, and scans only over those. This is significantly faster (e.g., 100–110 μs for 10 flips vs. 600 μs for a full scan on a typical mesh). The max_flips parameter is a static argument: changing it triggers recompilation, but within a simulation it is typically constant. See tutorials/03_vertex_models for a full usage example with per-edge cooldowns.

mesh = TriMesh.read_obj("../test_meshes/disk.obj")
hemesh = HeMesh.from_triangles(mesh.vertices.shape[0], mesh.faces)
geommesh = GeomMesh(*hemesh.n_items, mesh.vertices, mesh.face_positions)

Warning: readOBJ() ignored non-comment line 3:
  o flat_tri_ecmc

from jaxtyping import Float

def get_oriented_dual_he_length(vertices: Float[jax.Array, "n_vertices 2"],
                                face_positions: Float[jax.Array, "n_faces 2"],
                                hemesh: msh.HeMesh) -> Float[jax.Array, " n_hes"]:
    """Compute lengths of dual edges. Boundary dual edges get length 1. Negative sign = flipped edge."""
    dual_edges = face_positions[hemesh.heface]-face_positions[hemesh.heface[hemesh.twin]]

    edges = vertices[hemesh.orig]-vertices[hemesh.dest]
    edges_normalized = (edges.T / jnp.linalg.norm(edges, axis=-1)).T
    signed_dual_length = jnp.einsum('vi,vi->v', edges_normalized,
                                    dual_edges @ trig.get_rot_mat(-jnp.pi/2))
    signed_dual_length = jnp.where(hemesh.is_bdry_edge, 1, signed_dual_length)
    return signed_dual_length

from importlib import reload
reload(geom)

<module 'triangulax.geometry' from '/Users/nc1333/Documents/Princeton/Coding/triangulax/triangulax/geometry.py'>

# let's detect all edges with negative dual length, and flip them.

dual_lengths = geom.get_oriented_dual_he_length(geommesh.vertices, geommesh.face_positions, hemesh)
edges = jnp.where((dual_lengths < 0.0) & ~hemesh.is_bdry_edge & hemesh.is_unique)[0]
# we only want to flip unique hes!
edges, edges.size

(Array([  9, 185, 191, 335], dtype=int64), 4)

---

source

flip_all

def flip_all(
    hemesh:HeMesh, to_flip:Bool[Array, 'n_hes']
)->HeMesh:

Flip all (unique) half-edges where to_flip is True in a half-edge mesh. Wraps flip_edge.

Note: scans over all half-edges, which can be slow for large meshes. See flip_n_shortest for a more efficient alternative that only scans over a fixed number of candidate edges.

---

source

flip_by_id

def flip_by_id(
    hemesh:HeMesh, ids:Int[Array, 'flips'], to_flip:Bool[Array, 'flips']
)->HeMesh:

Flip half-edges from ids array where to_flip is True. Wraps flip_edge.

---

source

flip_n_shortest

def flip_n_shortest(
    hemesh:HeMesh, # The half-edge mesh.
    edge_lengths:Int[Array, 'n_hes'], # Per-half-edge edge lengths (e.g., dual/Voronoi edge lengths).
    threshold:float, # Edges shorter than this are flipped.
    max_flips:int=10, # Maximum number of edges to consider. Static argument (changing it triggers recompilation).
)->tuple: # The mesh after flipping.

Flip up to max_flips shortest edges below threshold.

Sorts edges by length, selects the max_flips shortest unique, non-boundary candidates, and flips those below threshold. Much faster than flip_all for large meshes.

to_flip = (dual_lengths < 0) & ~jnp.isnan(dual_lengths)

flipped_hemesh = flip_all(hemesh, to_flip=to_flip)

igl.is_edge_manifold(flipped_hemesh.faces)[0], igl.is_vertex_manifold(flipped_hemesh.faces)[0]

(True, np.True_)

flipped_hemesh = flip_all(hemesh, to_flip=(dual_lengths<0.02)) # no extra recompile

flipped_geommesh = geom.set_voronoi_face_positions(geommesh, flipped_hemesh)

fig = plt.figure(figsize=(8,8))

plt.triplot(*geommesh.vertices.T, hemesh.faces)
plt.triplot(*flipped_geommesh.vertices.T, flipped_hemesh.faces)

ax = plt.gca()
p1 = msh.cellplot(hemesh, geommesh.face_positions,
         cell_colors=np.array([0.,0.,0.,0.]), mpl_polygon_kwargs={"lw": 1, "ec": "k"})
p2 = msh.cellplot(flipped_hemesh, flipped_geommesh.face_positions,
              cell_colors=np.array([0.,0.,0.,0.]), mpl_polygon_kwargs={"lw": 1, "ec": "tab:orange"})
ax.add_collection(p1)
ax.add_collection(p2)
plt.axis("equal")

msh.label_plot(geommesh.vertices, hemesh.faces, fontsize=10, face_labels=False)

Splitting and collapsing vertices

The edge flip is the only topological modification of a half-edge mesh that leaves the number of vertices, edges, and faces constant. This makes it especially easy, and compatible with JAX’s “static array size” paradigm.

However, we may also want to simulate processes (like cell division or death) where the number of cells does change. We implement two elementary operations, which are inverses of one another: edge collapse and vertex split.

To collapse a half-edge e in a hemesh:

1. Delete faces hemesh.heface[e], hemesh.heface[hemesh.twin[e]]
2. Delete all the half-edges in those faces.
3. Glue the “gap” back together.
4. Merge the vertices hemesh.orig[e], hemesh.dest[e]

We must be careful to preserve the manifold structure of the mesh and deal with edge cases. We test the resulting half-edge mesh via plots and use libigl to verify that the mesh is in a valid state.

We also need a data structure (MeshReindexMap) to keep track of how vertices/edges/faces of the initial mesh map to those of the modified one.

---

source

remap_inds_removal_reverse

def remap_inds_removal_reverse(
    N:int, removed:Int[Array, 'n_removed']
)->Int[Array, 'N-n_removed']:

Remap indices after removal. Reverse of remap_inds_removal_forward.

---

source

remap_inds_removal_forward

def remap_inds_removal_forward(
    N:int, removed:Int[Array, 'n_removed']
)->Int[Array, 'N']:

Remap indices after removal. Returns array arr[i] = i - (removed < i).sum().

---

source

MeshReindexMap

def MeshReindexMap(
    v_forward:Int[Array, 'n_vertices_old'], v_reverse:Int[Array, 'n_vertices_new'],
    f_forward:Int[Array, 'n_faces_old'], f_reverse:Int[Array, 'n_faces_new'], he_forward:Int[Array, 'n_hes_old'],
    he_reverse:Int[Array, 'n_hes_new'], info:dict=<factory>
)->None:

Old↔︎new index maps produced by topology-changing operations.

N = 10
removed = jnp.array([5, 2])
forward = remap_inds_removal_forward(N, removed)
reverse = remap_inds_removal_reverse(N, removed)

forward[6], reverse[2], jnp.allclose(forward[reverse], jnp.arange(N - removed.shape[0]))

(Array(4, dtype=int64), Array(3, dtype=int64), Array(True, dtype=bool))

---

source

can_collapse_edge

def can_collapse_edge(
    hemesh:HeMesh, e:Int[Array, '']
)->Bool[Array, '']:

Check whether collapsing half-edge e would produce a valid mesh (link condition).

An edge can be collapsed if it is interior and the two endpoint vertices share exactly two common neighbors (the opposite vertices of the two adjacent faces). This is the discrete “link condition” that ensures the collapse preserves manifoldness.

---

source

collapse_edge

def collapse_edge(
    hemesh:HeMesh, e:int, check_boundary:bool=False
)->tuple:

Collapse half-edge e in a half-edge mesh. Keeps the origin vertex of e.

Returns a new HeMesh (does not modify in-place), and a MeshReindexMap for remapping vertex, half-edge, and face indices from the original mesh to the new mesh.

Does not check whether the collapse produces a valid mesh. Use can_collapse_edge to check first.

JIT-compatible, but calling with different numbers of vertices/edges/faces will cause recompilation.

# test on the existing example mesh. pick some interior, unique half-edge

candidates = np.where(np.asarray(hemesh.is_unique & (~hemesh.is_bdry_edge)))[0]
e_collapse = candidates[40]
print("Collapsing half-edge", e_collapse, "with vertices", int(hemesh.orig[e_collapse]), int(hemesh.dest[e_collapse]))

hemesh_collapsed, remap = collapse_edge(hemesh, e_collapse,)
vertices_collapsed =  geommesh.vertices[remap.v_reverse]

Collapsing half-edge 42 with vertices 8 129

hemesh, hemesh_collapsed # removes 1 vertex, 6 half-edges, and 2 faces

(HeMesh(N_V=131, N_HE=708, N_F=224), HeMesh(N_V=130, N_HE=702, N_F=222))

# stil valid mesh
(msh.test_mesh_validity(hemesh_collapsed), igl.is_edge_manifold(hemesh_collapsed.faces)[0],
 igl.is_vertex_manifold(hemesh_collapsed.faces)[0])

(True, True, np.True_)

44.8 μs ± 5.34 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

# visualize before/after
fig, ax = plt.subplots(1, 2, figsize=(8, 4))

plt.sca(ax[0])
v1, v2 = (hemesh.orig[e_collapse], hemesh.dest[e_collapse])
plt.scatter(*geommesh.vertices[v1], c="tab:orange")
plt.scatter(*geommesh.vertices[v2], c="tab:orange")
plt.triplot(np.asarray(geommesh.vertices)[:, 0], np.asarray(geommesh.vertices)[:, 1], np.asarray(hemesh.faces), lw=0.5, color="k")
plt.title("before")
plt.axis("equal")
plt.axis("off")

plt.sca(ax[1])
plt.scatter(*geommesh.vertices[v1], c="tab:orange")
plt.triplot(np.asarray(vertices_collapsed)[:, 0], np.asarray(vertices_collapsed)[:, 1],
            np.asarray(hemesh_collapsed.faces), lw=0.5, color="k")
plt.title("after collapse")
plt.axis("equal")
plt.axis("off")

(np.float64(-1.10003475),
 np.float64(1.09628575),
 np.float64(-1.09934025),
 np.float64(1.09050125))

Split vertex (“cell division”)

The opposite of edge collapse: splitting a vertex into two. We specify two half-edges (the “splitting axis”) that originate at a common vertex. Like before, we need a MeshReindexMap tracking how old and new mesh elements are related. New elements are appended at the end of the arrays.

---

source

can_split_vertex

def can_split_vertex(
    hemesh:HeMesh, e1:int, e2:int
)->bool:

Check whether splitting a vertex along half-edges e1 and e2 is valid.

Both half-edges must originate from the same vertex, and both must be on interior faces.

---

source

split_vertex

def split_vertex(
    hemesh:HeMesh, e1:int, e2:int, check_args:bool=False
)->tuple:

Split a vertex into two along a “splitting axis” given by two half-edges originating at that vertex.

New vertex inserted at origin of e2. The new vertex will be the final one in the array.

This function is not JIT-compatible, since it depends on iterating around the vertex to update origins/destinations.

# test split on the existing example mesh

# choose an interior vertex (avoid boundary) and two outgoing half-edges for split axis
v_split = jnp.where(~hemesh.is_bdry)[0][10]
v_new = hemesh.n_vertices  # new vertex index

ring = hemesh.iterate_around_vertex(v_split)
h1 = int(ring[0])
h2 = int(ring[len(ring)//2])
print("Splitting vertex", v_split, "axis hes", h1, h2)

hemesh_split, smap = split_vertex(hemesh, h1, h2)
print("Old:", hemesh, "new:", hemesh_split)

Splitting vertex 13 axis hes 56 407
Old: HeMesh(N_V=131, N_HE=708, N_F=224) new: HeMesh(N_V=132, N_HE=714, N_F=226)

assert hemesh_split.n_vertices == hemesh.n_vertices + 1
assert hemesh_split.n_hes == hemesh.n_hes + 6
assert hemesh_split.n_faces == hemesh.n_faces + 2

F_split = np.asarray(hemesh_split.faces, dtype=np.int32)
print("edge manifold:", igl.is_edge_manifold(F_split)[0])
print("vertex manifold:", igl.is_vertex_manifold(F_split)[0])
print("Valid HE mesh:", msh.test_mesh_validity(hemesh_split))

edge manifold: True
vertex manifold: True
Valid HE mesh: True

# inverse consistency check: split then collapse the inserted edge
e_join = hemesh.n_hes+1 # error for  2*hemesh.n_vertices+1 ?
hemesh_back, back_map = collapse_edge(hemesh_split, e_join)

F0 = msh._canonical_faces_np(hemesh.faces)
F_back = msh._canonical_faces_np(hemesh_back.faces)
print("back to original counts?", hemesh_back.n_items == hemesh.n_items)
print("back to original faces?", np.array_equal(F0, F_back))
print("Valid HE mesh after collapse:", msh.test_mesh_validity(hemesh_back))

back to original counts? True
back to original faces? True
Valid HE mesh after collapse: True

# offset the new vertex slightly for visibility

vertices_split = np.concatenate([geommesh.vertices, geommesh.vertices[v_split][None, :]], axis=0)
d = trig.get_perp_2d(geommesh.vertices[hemesh.dest[h1]] - geommesh.vertices[hemesh.dest[h2]])

eps = 0.2
vertices_split[v_new] = vertices_split[v_new] +  eps * d
vertices_split[v_split] = vertices_split[v_split] - eps * d

vertices_collapsed =  vertices_split[back_map.v_reverse]

# quick visualization (triangulation plot)

fig, ax = plt.subplots(1, 3, figsize=(12, 4))

plt.sca(ax[0])
plt.triplot(*geommesh.vertices.T, hemesh.faces, lw=0.5, color="k")
plt.scatter(*geommesh.vertices[v_split], c="tab:orange")
plt.title("before split")
plt.axis("equal"); plt.axis("off")

plt.sca(ax[1])
plt.triplot(*vertices_split.T, hemesh_split.faces, lw=0.5, color="k")
plt.scatter(*vertices_split[v_split], c="tab:orange")
plt.scatter(*vertices_split[v_new], c="tab:red")
plt.title("after split (new vertex in red)")
plt.axis("equal"); plt.axis("off")

plt.sca(ax[2])
plt.triplot(*vertices_collapsed.T, hemesh_back.faces, lw=0.5, color="r")
plt.scatter(*vertices_split[v_split], c="tab:orange")
plt.title("Remerged connectivity")
plt.axis("equal"); plt.axis("off")

(np.float64(-1.10003475),
 np.float64(1.09628575),
 np.float64(-1.09934025),
 np.float64(1.09050125))

Not yet implemented

The following topological operations are not yet available in triangulax:

• Edge split: insert a new vertex on an existing edge, splitting it and both adjacent faces. (Distinct from vertex split above.)
• Edge contraction with boundary support: the current collapse_edge only handles interior edges.
• Batch collapse / split: JIT-compatible routines for performing multiple collapses or splits per time step, analogous to flip_all / flip_n_shortest for edge flips.
• Vertex removal: remove a vertex and re-triangulate the resulting hole.