## 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.
### 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
[`HeMesh`](https://nikolas-claussen.github.io/triangulax/src/halfedge_datastructure.html#hemesh)es.
We basically edit the various connectivity arrays (in a JAX-compatible
way).
The algorithm (and the naming conventions in
[`flip_edge`](https://nikolas-claussen.github.io/triangulax/src/topological_modifications.html#flip_edge))
are from
[here](https://jerryyin.info/geometry-processing-algorithms/half-edge/).
**Before**
image.png
**After**
image.png
------------------------------------------------------------------------
source
### flip_edge
``` python
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.
``` python
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
``` python
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))
``` python
# flip edge and recompute face positions
flipped_hemesh = flip_edge(hemesh, e=335)
flipped_geommesh = msh.set_voronoi_face_positions(geommesh, flipped_hemesh)
```
``` python
# 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))
``` python
# 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(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.
``` python
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
``` python
# 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)
------------------------------------------------------------------------
source
### flip_all
``` python
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.*
------------------------------------------------------------------------
source
### flip_by_id
``` python
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.*
``` python
to_flip = (dual_lengths < 0) & ~jnp.isnan(dual_lengths)
flipped_hemesh = flip_all(hemesh, to_flip=to_flip)
```
``` python
flipped_hemesh = flip_all(hemesh, to_flip=(dual_lengths<0.02)) # no extra recompile
```
``` python
flipped_geommesh = msh.set_voronoi_face_positions(geommesh, flipped_hemesh)
```
``` python
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)
```
### 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. 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`:
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, 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
------------------------------------------------------------------------
source
### remap_inds_removal_reverse
``` python
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
``` python
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
``` python
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=
)->None:
```
*Old↔new index maps produced by topology-changing operations.*
``` python
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
### collapse_edge
``` python
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.
``` python
# 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
``` python
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))
``` python
# 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_)
``` python
```
40.6 μs ± 2.64 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
``` python
# 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”)
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.
------------------------------------------------------------------------
source
### split_vertex
``` python
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.
``` python
# 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)
``` python
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
``` python
# 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
``` python
# 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]
```
``` python
# 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))
### Next steps
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.
