algorithms: Geometry-processing algorithms

This notebook defines geometry processing algorithms building on the geometric primitives defined in the preceding modules.

Mesh and point-to-point alignment

Kabsch alignment finds the optimal (proper) rotation and translation between two matched sets of points. Here, we provide a JAX-compatible version.

---

source

kabsch_align

def kabsch_align(
    v:Float[Array, 'n_vertices dim'], # vertices to be aligned
    v_ref:Float[Array, 'n_vertices dim'], # reference vertices
)->tuple: # v optimally aligned to v_ref

Optimally rotate/translate v onto v_ref (differentiable Kabsch algorithm).

The resulting rotation R is proper (det(R) = 1), so no reflections are allowed. Works in any dimension, and is compatible with jax.jit and jax.grad. The convention is aligned == v @ R + d.

Smoothing and remeshing

Two operations to improve mesh quality without altering surface shape: 1. Tangential smoothing moves mesh vertices to the area-weighted average of the neighbors, along the tangential direction only 2. Delaunay flips improve triangle aspect ratio by edge flips. This modifies the mesh topology, but leaves the vertex set unchanged.

# Load test mesh
mesh = TriMesh.read_obj("../test_meshes/disk.obj", dim=2)
hemesh = msh.HeMesh.from_triangles(mesh.vertices.shape[0], mesh.faces)
vertices = mesh.vertices

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

Mesh quality assessment

Functions to evaluate triangle mesh quality: per-face maximum corner angle, summary statistics, and a human-readable quality report.

---

source

get_face_angles

def get_face_angles(
    vertices:Float[Array, 'n_vertices dim'], # Vertex positions.
    hemesh:HeMesh, # Half-edge mesh connectivity.
)->Float[Array, 'n_faces 3']: # The three interior angles of each triangle.

Get the three corner angles of every face (radians).

Uses the half-edge corner angles. Angles within a face sum to pi.

# Test get_face_max_angles
angles = get_face_angles(vertices, hemesh)
max_angles = jnp.max(angles, axis=-1)
print(f"max angles: min={jnp.rad2deg(max_angles.min()):.1f}°,max={jnp.rad2deg(max_angles.max()):.1f}°, mean={jnp.rad2deg(max_angles.mean()):.1f}°")
assert max_angles.shape == (hemesh.n_faces,)
assert jnp.all(max_angles >= jnp.pi / 3 - 1e-6)  # max angle >= 60° always

max angles: min=60.4°,max=106.6°, mean=75.0°

---

source

get_mesh_quality_stats

def get_mesh_quality_stats(
    vertices:Float[Array, 'n_vertices dim'], # Vertex positions.
    hemesh:HeMesh, # Half-edge mesh connectivity.
    degenerate_angle:float=5.0, # Threshold (degrees) for a triangle to be considered degenerate.
(I.e. if max angle > pi-threshold or min angle < threshold.)
    digits:int=5, # Number of decimal digits to round the statistics to.
)->dict:

Compute mesh quality statistics.

# Test mesh quality report
stats = get_mesh_quality_stats(vertices, hemesh)
stats

{'areas_min': 0.00615,
 'areas_max': 0.0218,
 'areas_cv': 0.18788,
 'max_angle': 106.60695,
 'min_angle': 31.15478,
 'angles_std': 13.73223,
 'n_degenerate': 0,
 'n_total_faces': 224}

Delaunay flipping

Flip edges to improve triangle quality based on the Delaunay criterion: for each interior edge, the sum of the two opposite angles should not exceed π. Based on GeometryCentral: Extrinsic Delaunay Flipping.

Note: extrinsic flipping is not guaranteed to produce a fully Delaunay mesh, but generally improves quality in practice.

---

source

is_locally_delaunay

def is_locally_delaunay(
    vertices:Float[Array, 'n_vertices dim'], # Vertex positions.
    hemesh:HeMesh, # Half-edge mesh connectivity.
)->Bool[Array, 'n_hes']: # True where the edge satisfies the Delaunay condition.

Check the local Delaunay condition for each edge.

An interior edge is locally Delaunay when the sum of the two opposite angles does not exceed π. Boundary edges are always considered Delaunay.

# Test is_locally_delaunay
delaunay = is_locally_delaunay(vertices, hemesh)
print(f"Non-Delaunay edges: {(~delaunay & hemesh.is_unique & ~hemesh.is_bdry_edge).sum()}")
assert delaunay.shape == (hemesh.n_hes,)
assert jnp.all(delaunay[hemesh.is_bdry_edge])  # boundary edges are always Delaunay

Non-Delaunay edges: 4

---

source

fix_delaunay

def fix_delaunay(
    vertices:Float[Array, 'n_vertices dim'], # Vertex positions (unchanged by flips).
    hemesh:HeMesh, # Half-edge mesh connectivity.
    max_iters:int=2, # Maximum number of sweep iterations.
    max_flips:int=10, # Maximum number of flips at each iteration.
)->tuple:

Flip non-Delaunay edges iteratively until convergence.

Each iteration identifies non-Delaunay interior edges and flips them using topology.flip_by_score. Stops when no more flips are needed or max_iters is reached.

Important: when JIT-compiling, use jax.jit with static_argnames=[‘max_iters’, ‘max_flips’].

Tangential vertex smoothing

Vertex smoothing moves each vertex towards the average of its neighborhood. triangulax implements laplacian smoothing, which moves vertices towards the average of neighboring vertex positions.

In either case, for 3D meshes the displacement is projected tangentially (normal component removed). Boundary conditions can be 'fixed' (boundary vertices immobile) or 'free'.

Based on: GeometryCentral: Tangential Vertex Smoothing.

---

source

smooth_vertices_laplacian

def smooth_vertices_laplacian(
    vertices:Float[Array, 'n_vertices dim'], # Vertex positions.
    hemesh:HeMesh, # Half-edge mesh connectivity.
    step_size:float=1.0, # Fraction of displacement to apply (1 = full step).
    bc:str='fixed', # Boundary condition: 'fixed' freezes boundary vertices,
'free' allows them to move, and 'slide' allows them to
move only tangentially along the boundary.
)->Float[Array, 'n_vertices dim']: # Updated vertex positions.

One step of tangential Laplacian vertex smoothing.

Moves each vertex towards the mean position of its neighbours. For 3D meshes, the displacement is projected onto the tangent plane.

plt.triplot(noisy_v[:, 0], noisy_v[:, 1], hemesh.faces)
plt.triplot(smoothed_v[:, 0], smoothed_v[:, 1], hemesh.faces)
plt.axis('equal')

(np.float64(-1.10003475),
 np.float64(1.09628575),
 np.float64(-1.0993421965883812),
 np.float64(1.0905421283560033))

# Test 3D tangential smoothing on a sphere
mesh_3d = TriMesh.read_obj("../test_meshes/sphere.obj", dim=3)
hemesh_3d = msh.HeMesh.from_triangles(mesh_3d.vertices.shape[0], mesh_3d.faces)
verts_3d = mesh_3d.vertices

key = jax.random.PRNGKey(1)
noise_3d = 0.1 * jax.random.normal(key, shape=verts_3d.shape)
noisy_3d = verts_3d + noise_3d

print("Before 3D smoothing:")
print(get_mesh_quality_stats(noisy_3d, hemesh_3d))

smoothed_3d = noisy_3d
for _ in range(5):
    smoothed_3d = smooth_vertices_laplacian(smoothed_3d, hemesh_3d, step_size=0.3, bc='free')

print("\nAfter 5 Laplacian smoothing steps (3D):")
print(get_mesh_quality_stats(smoothed_3d, hemesh_3d))

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

Before 3D smoothing:
{'areas_min': 0.05567, 'areas_max': 0.31956, 'areas_cv': 0.35048, 'max_angle': 120.09799, 'min_angle': 19.84483, 'angles_std': 17.23441, 'n_degenerate': 0, 'n_total_faces': 80}

After 5 Laplacian smoothing steps (3D):
{'areas_min': 0.10323, 'areas_max': 0.20764, 'areas_cv': 0.13673, 'max_angle': 79.5562, 'min_angle': 45.60472, 'angles_std': 7.12533, 'n_degenerate': 0, 'n_total_faces': 80}

p = meshplot.plot(noisy_3d, hemesh_3d.faces, shading={"wireframe":True}, return_plot=True)

p.add_mesh(smoothed_3d + np.array([3, 0, 0]), hemesh_3d.faces, shading={"wireframe":True})

1