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.
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
After
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.
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_ecmcplt.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 = msh.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)
In a simulation, we need to carry out edge flips at every time step. The function flip_edge(hemesh: HeMesh, e: int) -> HeMesh does a single edge flip by modifying the connectivity arrays. Luckily, it is already JAX-compatible (we can JIT-compile it).
To carry out multiple flips, we must do the flips in sequence (otherwise, you risk leaving the mesh in an invalid state). To make things JAX-compatible, we do a jax.lax.scan scan over all half-edges.
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# let's detect all edges with negative dual length, and flip them.
dual_lengths = msh.get_signed_dual_he_length(geommesh.vertices, geommesh.face_positions, hemesh)
edges = jnp.where((dual_lengths < -0.05) & ~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)
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.
def flip_by_id(
hemesh:HeMesh, ids:Int[Array, 'flips'], to_flip:Bool[Array, 'flips']
)->HeMesh:
Flip half-edges from ids array if the to_flip is True. Wraps flip_edge.
to_flip = (dual_lengths < 0) & ~jnp.isnan(dual_lengths)
flipped_hemesh = flip_all(hemesh, to_flip=to_flip)flipped_hemesh = flip_all(hemesh, to_flip=(dual_lengths<0.02)) # no extra recompileflipped_geommesh = msh.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()
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)
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. Let’s implement two elementary operation, which are inverses of another: edge collapse and vertex split. Let’s start with edge collapse.
To split a half-edge e in a hemesh:
hemesh.heface[e], hemesh.heface[hemesh.twin[e]]hemesh.orig[e], hemesh.dest[e] We must be careful to preserve the manifold structure of the mesh, deal with edge cases. We can test the resulting half-edge mesh via plots, and use libigl to verify that the mesh is in a valid state.Finally, we also need a data structure to keep track of the map from the vertices/edges/faces of the initial mesh to that of the modified one.
WARNING: code not fully tested
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.
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().
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))
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 three arrays for remapping vertex, half-edge, and face indices from the original mesh to the new mesh.
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 129hemesh, 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_)40.6 μs ± 2.64 μ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))
Next, let’s create the opposite operation - splitting a vertex into two. To do so, we need to specify two half-edges, the “splitting axis”, which meet at a common vertex. Like before, we also need a map of how hold and new mesh elements are related. We append the newly created mesh elements at the end of the different arrays.
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 whose 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))
Looks good - the JAX-compatible triangular-mesh data structures seem to work. In particular, the tricky T1/edge-flip function. Next steps: linear operators on meshes and geometry processing.