Meshing Algorithms

Overview

Rasmah offers several meshing algorithms, each a different point on the surface-vs-volume and uniform-vs-adaptive spectrum. This tutorial surveys them, shows what each is for, and gives the “which one do I use?” rule of thumb.

The surface_mesh / tetrahedralize convenience functions dispatch over these: by default a sphere uses an exact icosphere, a box a 12-triangle box mesh, and a cylinder/cone a ruled analytic mesh — each both lighter and more accurate than marching cubes. A general SDF uses marching_cubes (or adaptive_surface_mesh, or dual_contouring via method = :dual), and a volume uses marching_tetrahedra (or conforming_tetmesh).

Throughout we mesh the unit sphere, whose exact surface area is 4π ≈ 12.566:

sdf = sphere(1.0)
f = x -> field(sdf, x)
lo, hi = [-1.5, -1.5, -1.5], [1.5, 1.5, 1.5]

([-1.5, -1.5, -1.5], [1.5, 1.5, 1.5])

Marching cubes (surface)

marching_cubes(f, lo, hi, n) extracts a watertight triangle mesh of the f = 0 level set on a uniform n × n × n grid. It resolves the ambiguous-face cases with the asymptotic decider, so the result is a closed 2-manifold:

mc = marching_cubes(f, lo, hi, 32)
println(size(mc.vertices, 2), " vertices, ", size(mc.faces, 2), " triangles")
println("surface area = ", round(surface_area(mc); digits = 3), "   (4π ≈ 12.566)")

8568 vertices, 4280 triangles
surface area = 12.532 m^2   (4π ≈ 12.566)

Cost is n³ — every cell is visited, even the ones far from the surface.

Adaptive octree (surface)

adaptive_surface_mesh(f, lo, hi; min_depth, max_depth) refines an octree only where the field is near zero, then extracts the surface with a dual scheme (surface nets): one vertex per surface cell, connected into fans around crossing edges. The far interior/exterior stay coarse, so the cost concentrates on the surface band:

ad = adaptive_surface_mesh(f, lo, hi; min_depth = 3, max_depth = 7)
println(size(ad.vertices, 2), " vertices, ", size(ad.faces, 2), " triangles")
println("surface area = ", round(surface_area(ad); digits = 3))

34328 vertices, 68652 triangles
surface area = 12.562 m^2

min_depth is the uniform base level; max_depth caps the near-surface refinement. The dual extraction is crack-free and watertight.

bracket = difference(box(2.0, 2.0, 2.0), sphere(0.5))
fine = surface_mesh(bracket; adaptive = true, max_depth = 6, method = :dual)
graded = surface_mesh(bracket; adaptive = true, max_depth = 6, method = :dual, curvature = true)
println(size(fine.faces, 2), " triangles -> ", size(graded.faces, 2), "   (graded)")
println("volume ", round(volume(fine); digits = 3), " -> ", round(volume(graded); digits = 3))

49752 triangles -> 19656   (graded)
volume 7.476 m^3 -> 7.476 m^3

The surface_mesh convenience

surface_mesh picks the algorithm for you. With n left at its default nothing, the standard primitives resolve to a lightweight analytic mesh — an exact icosphere for a sphere, a 12-triangle box_mesh for a box, a ruled cylinder/cone — while a general implicit uses marching_cubes at n = 32, or the adaptive octree when adaptive = true:

println(size(surface_mesh(sphere(1.0)).vertices, 2), " vertices   (icosphere)")
println(size(surface_mesh(box(1.0, 1.0, 1.0)).vertices, 2), " vertices   (box_mesh)")
println(size(surface_mesh(sdf; n = 32).vertices, 2), " vertices   (marching cubes)")

bracket = difference(box(2.0, 2.0, 2.0), sphere(0.5))
println(
    size(surface_mesh(bracket; n = 32).vertices, 2),
    " vertices   (marching cubes, non-primitive SDF)",
)

642 vertices   (icosphere)
8 vertices   (box_mesh)
642 vertices   (marching cubes)
6060 vertices   (marching cubes, non-primitive SDF)

Passing n scales the sphere’s icosphere subdivision level (so n = 64 refines it further), which keeps the mesh exact while honouring the resolution hint; shapes with no analytic resolution fall back to marching cubes at n.

Dual contouring (sharp features)

marching_cubes rounds off sharp edges and corners: its vertices live on grid edges, so a box’s corner falls somewhere between grid crossings and comes out blunted. dual_contouring fixes this with the QEF (quadratic error function): one vertex per surface cell, placed at the minimizer of the Hermite error built from the cell’s crossing points and normals. Three mutually-orthogonal crossings (a box corner) solve to exactly the corner, so CAD edges survive:

bracket = difference(box(2.0, 2.0, 2.0), sphere(0.5))

mc = surface_mesh(bracket; n = 32)
dc = surface_mesh(bracket; method = :dual, n = 32)
println("marching cubes: ", size(mc.vertices, 2), " vertices")
println("dual contouring: ", size(dc.vertices, 2), " vertices")
println("dual volume = ", round(volume(dc); digits = 3), "   (box 8 − sphere ≈ 7.48)")

marching cubes: 6060 vertices
dual contouring: 6064 vertices
dual volume = 7.467 m^3   (box 8 − sphere ≈ 7.48)

Use method = :dual with adaptive = true to combine QEF vertices with the graded octree, and pass qef = true to adaptive_surface_mesh directly for the same effect. dual_contouring(f, lo, hi, n) is the uniform-grid form.

Decimation (simplify for display and export)

A fine marching_cubes mesh has more triangles than you need for a quick VTK.js preview or an STL export. decimate collapses edges in order of increasing quadric error (Garland–Heckbert QEM), so the surface shrinks to a fraction of its faces while its area and volume are preserved:

fine = dual_contouring(f, lo, hi, 48)
small = decimate(fine, 0.1)
println(size(fine.faces, 2), " triangles -> ", size(small.faces, 2))
println("area ", round(surface_area(fine); digits = 3), " -> ", round(surface_area(small); digits = 3))
println("volume ", round(volume(fine); digits = 3), " -> ", round(volume(small); digits = 3))

9564 triangles -> 956
area 12.566 m^2 -> 12.544 m^2
volume 4.187 m^3 -> 4.165 m^3

preserve_boundary = false lets boundary edges collapse too (for an open surface); the default keeps a mesh’s outer rim intact.

marching_cubes and icosphere emit one vertex per crossing edge, so interior edges carry duplicate vertices (a “triangle soup”). Topology-based operations — decimate included — want a welded mesh with shared vertices, which weld_vertices produces in one pass:

raw = marching_cubes(f, lo, hi, 24)
welded = weld_vertices(raw)
println(size(raw.vertices, 2), " soup vertices -> ", size(welded.vertices, 2), " welded")
println("area preserved: ", round(surface_area(welded); digits = 3), " (was ", round(surface_area(raw); digits = 3), ")")

4728 soup vertices -> 1158 welded
area preserved: 12.505 m^2 (was 12.505 m^2)

dual_contouring already welds by construction, so it needs no such step.

Local remeshing (repair a bad tetrahedral mesh)

remesh_optimize (§ mesh-quality) can only re-freeze a topology; it cannot fix one that is already badly distorted. local_remesh performs the actual repair: edge splits and collapses, 2-3 / 3-2 / 4-4 flips, and interior relaxation, each guarded by a scaled-Jacobian quality criterion:

good = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2)
nodes = copy(good.nodes); nodes[:, 14] .+= [0.3, 0.3, 0.3]
bad = TetMesh(nodes, good.elements)
println("min quality: ", round(min_quality(good); digits = 3), " -> ", round(min_quality(bad); digits = 3))

fixed = local_remesh(bad; target_quality = 0.3, iterations = 2)
println("after local_remesh: ", round(min_quality(fixed); digits = 3))
println("volume preserved: ", round(volume(fixed); digits = 6), " (was ", round(volume(bad); digits = 6), ")")

min quality: 0.408 -> 0.225
after local_remesh: 0.408
volume preserved: 1.0 m^3 (was 1.0 m^3)

It is discrete (topology is a fresh freeze, never differentiated), which is exactly what remesh_optimize needs as a fallback — pass remesh = :local to repair the mesh every time the quality check trips:

# remesh_optimize(build, θ0; remesh = :local)  # runs local_remesh before each check

Marching tetrahedra (volume)

marching_tetrahedra(f, lo, hi, n) fills the interior with tetrahedra, one per-cell decomposition of each cube into 6 tets, clipped to the f < 0 region. Boundary vertices lie on grid edges, so the surface is a staircase that only approximates the true surface:

mt = marching_tetrahedra(f, lo, hi, 16)
println(size(mt.nodes, 2), " nodes, ", size(mt.elements, 2), " tets")

bnd = unique(vec(boundary_faces(mt).faces))
println("max |sdf| on the boundary: ", round(maximum(abs(field(sdf, mt.nodes[:, i])) for i in bnd); digits = 4))

2205 nodes, 7860 tets
max |sdf| on the boundary: 0.0106

Conforming tetrahedra (volume)

conforming_tetmesh removes the staircase: it projects every boundary vertex onto the f = 0 level set (Newton via snap_to_zero) and relaxes the interior with the harmonic laplacian_deform, so boundary nodes sit on the surface and interior tets stay well-shaped:

ct = conforming_tetmesh(f, lo, hi, 16)
println(size(ct.nodes, 2), " nodes, ", size(ct.elements, 2), " tets")
println("max |sdf| on the boundary: ", round(maximum(abs(field(sdf, ct.nodes[:, i])) for i in bnd); digits = 4))

2205 nodes, 7860 tets
max |sdf| on the boundary: 0.0

Connectivity is frozen — only vertex positions move — so the conforming step stays differentiable. tetrahedralize(g, n; conform = true) is the convenience entry point:

println(size(tetrahedralize(sdf, 16).nodes, 2), " nodes   (plain)")
println(size(tetrahedralize(sdf, 16; conform = true).nodes, 2), " nodes   (conforming)")

4565 nodes   (plain)
4565 nodes   (conforming)

The frozen form freeze_conforming_tetmesh keeps this snap-and-relax pipeline differentiable: it freezes the marching-tetrahedra connectivity once, and each deform re-projects the boundary and relaxes the interior at the new field:

fc = freeze_conforming_tetmesh(f, lo, hi, 16)

mt2 = deform(fc, f)
println(size(mt2.nodes, 2), " nodes, ", size(mt2.elements, 2), " tets")

6499 nodes, 7860 tets

Structured box tetrahedralization

tetrahedralize_box(w, h, t, nx, ny, nz) fills an axis-aligned box with a structured grid split into 6 tets per cell — no SDF at all, just an exact background mesh:

tb = tetrahedralize_box(1.0, 1.0, 1.0, 3, 3, 3)
println(size(tb.nodes, 2), " nodes, ", size(tb.elements, 2), " tets")

64 nodes, 162 tets

Post-processing: smoothing and snapping

marching_cubes can leave slivers (long, thin triangles). smooth_surface applies Taubin’s two-pass λ|μ filter, which smooths without shrinking the way plain Laplacian smoothing does. The ImplicitGeometry overload re-projects every vertex back onto the surface each pass:

raw = marching_cubes(f, lo, hi, 24)
smoothed = smooth_surface(raw; iterations = 3)
println("smooth_surface: ", size(smoothed.vertices, 2), " vertices (unchanged topology)")

smooth_surface: 4728 vertices (unchanged topology)

The projection primitive behind all of this is snap_to_zero, a Newton step x ← x − f ∇f/‖∇f‖²:

p = snap_to_zero(f, [0.7, 0.0, 0.0]; iters = 5)
println("snap_to_zero: ", round.(p; digits = 4), "   (|sdf| = ", round(abs(field(sdf, p)); digits = 4), ")")

snap_to_zero: [1.0, 0.0, 0.0]   (|sdf| = 0.0)

snap_surface runs it over an entire surface mesh at once.

Native Delaunay meshing and exact predicates

The Delaunay pipeline rests on exact geometric predicates — Shewchuk adaptive-precision arithmetic (error-free two-sum/two-product + expansions), so the cavity tests are exact even for nearly co-circular/co-spherical input. They accept points as vectors or NTuples and return ±1/0:

# orient2d: sign of the 2D cross product (positive = counter-clockwise)
println(orient2d([0.0, 0.0], [1.0, 0.0], [0.0, 1.0]))      # +1 (CCW)
println(orient2d([0.0, 0.0], [0.0, 1.0], [1.0, 0.0]))      # -1 (CW)
println(orient2d([0.0, 0.0], [1.0, 0.0], [2.0, 0.0]))      #  0 (collinear)

# incircle: sign of the in-circle test (positive = inside the circumcircle)
println(incircle([0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.2, 0.2]))   # +1 (inside)

# orient3d / insphere: the 3D analogues
println(orient3d([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]))
println(insphere([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.1, 0.1, 0.1]))

1
-1
0
1
-1
-1

delaunay_refine_2d / delaunay_refine insert circumcenters to enforce a radius–edge bound ratio, confined to inside, with optional Ruppert boundary encroachment splitting via boundary (edge pairs in 2D, face triples in 3D):

# 2D unit square with the boundary recovered exactly
pts = [0.0 1.0 1.0 0.0; 0.0 0.0 1.0 1.0]
bnd = [(1, 2), (2, 3), (3, 4), (4, 1)]
m2 = delaunay_refine_2d(
    pts; inside=p -> 0.0 <= p[1] <= 1.0 && 0.0 <= p[2] <= 1.0,
    boundary=bnd, ratio=0.7, max_steiner=100)
println(size(m2.faces, 2), " quality triangles")

# 3D refinement confined to the unit ball
p3 = hcat([0.6 * randn(3) for _ in 1:60]...)
m3 = delaunay_refine(p3; inside=p -> p[1]^2 + p[2]^2 + p[3]^2 < 1.0, ratio=2.0)
println(!isempty(m3.elements), "  volume = ", round(volume(m3); digits=3))

172 quality triangles
true  volume = 6.687 m^3

triangulate_region_2d grid-samples a planar region and Delaunay-refines it into a quality (z = 0) triangle mesh:

reg = triangulate_region_2d(p -> p[1]^2 + p[2]^2 < 1.0, [-1.0, -1.0], [1.0, 1.0]; n=16)
println(all(reg.vertices[3, :] .== 0))

true

Advancing-front triangulation

advancing_front_2d(boundary, target) triangulates a planar polygon by repeatedly taking the shortest front edge and forming a triangle with a reused or newly placed vertex:

sq = [0.0 1.0 1.0 0.0; 0.0 0.0 1.0 1.0]
af = advancing_front_2d(sq; target=0.3)
println(size(af.faces, 2), " triangles, all CCW: ",
    all(orient2d(af.vertices[:, af.faces[1, j]],
                 af.vertices[:, af.faces[2, j]],
                 af.vertices[:, af.faces[3, j]]) > 0 for j in axes(af.faces, 2)))

22 triangles, all CCW: true

Level-of-detail (LOD) pyramids

lod_pyramid builds a decimated LOD chain (each level a fraction of the last), and select_lod picks the level for a viewing distance:

fine = icosphere(4)
py = lod_pyramid(fine, 4)          # fine → coarse (py[1] === fine)
println(length(py), " levels; face counts ", [size(m.faces, 2) for m in py])

println(select_lod(py, 0.5))       # nearest → finest level
println(select_lod(py, 1e9))       # very far → coarsest level

4 levels; face counts [5120, 2560, 1280, 640]
1
4

The quad- and hexa- meshers return QuadMesh / HexMesh (quad_dominant, hexahedralize_box), each with its own vertex/element layout.

Point-cloud reconstruction

Given only a point cloud (a scan, or samples from a shape), Rasmah reconstructs a surface mesh. estimate_normals recovers per-point normals from the local neighbourhood, and poisson_reconstruction solves a screened Poisson equation to produce a watertight surface:

pts = surface_mesh(sdf; n = 40).vertices          # a noisy point cloud stand-in

nrm = estimate_normals(pts; k = 12)
rec = poisson_reconstruction(pts; normals = nrm, n = 40)
println("reconstructed surface: ", size(rec.faces, 2), " faces")

reconstructed surface: 22588 faces

The companion reconstruction methods are alpha_shape (ball-based), ball_pivot, and power_crust/voronoi_surface (medial-axis), with mls_projection / mls_smooth for moving-least-squares denoising. For scan pipelines, statistical_outlier_removal, radius_outlier_removal, and voxel_downsample clean the cloud, and icp / align_scans / fuse_scans register and merge multiple scans.

Mesh booleans and offsets

Surface meshes support boolean operations through a signed-distance round trip — mesh_union, mesh_intersection, mesh_difference, or the generic mesh_boolean(a, b, op):

b1 = box_mesh([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
b2 = box_mesh([0.5, 0.5, 0.5], [1.5, 1.5, 1.5])

println("union:        ", size(mesh_union(b1, b2).faces, 2), " faces")
println("intersection: ", size(mesh_intersection(b1, b2).faces, 2), " faces")
println("difference:   ", size(mesh_difference(b1, b2).faces, 2), " faces")

union:        32168 faces
intersection: 41768 faces
difference:   41516 faces

mesh_offset(m, d) and mesh_hollow(m, t) grow/shell a mesh by a distance.

Frame fields, curvature, and feature extraction

For structured re-meshing, cross_field / field_aligned_quads align a quad mesh to a cross field, and frame_field / frame_field_hex extend the idea to hexahedral meshing. Differential properties come from mesh_curvatures / vertex_gaussian_curvature / vertex_mean_curvature, and sharp features from feature_edges / feature_vertices.

Which one do I use?

• Triangle surface of a sphere → surface_mesh(sphere(r)) (exact icosphere).
• Triangle surface of a box → surface_mesh(box(a, b, c)) (12-triangle box_mesh).
• Triangle surface of any SDF → surface_mesh(sdf; n = …) (marching_cubes).
• Sharp edges and corners → surface_mesh(sdf; method = :dual) (dual contouring).
• Fine surface, coarse everywhere else → adaptive_surface_mesh.
• Graded surface (flat coarse, curved fine) → adaptive_surface_mesh(…; curvature = true).
• Tetrahedral volume for FEM → tetrahedralize(sdf, n).
• Volume whose boundary must lie on the surface → tetrahedralize(sdf, n; conform = true).
• Exact background box mesh → tetrahedralize_box.
• Shrink a mesh for preview or export → decimate(m, ratio) (weld first with weld_vertices).
• Repair a badly distorted tet mesh → local_remesh(m).
• Clean up marching-cubes slivers → smooth_surface(m, sdf).
• Delaunay triangulate a point set (2D/3D) → delaunay_triangulation / delaunay_tetrahedralization (exact predicates, near-linear).
• Quality tet mesh from an SDF → delaunay_tetmesh(sdf, n) / tetrahedralize(sdf, n; method = :delaunay).
• Quality triangle mesh of a planar region → triangulate_region_2d.
• Recover a boundary exactly → conforming_delaunay(points, boundary) (and flip_recover for flip-based edge recovery).
• Remove Delaunay slivers → exude_slivers(points).
• Optimize vertex positions → odt_smooth / cvt_smooth (Lloyd).
• Uniformize a triangle mesh → isotropic_remesh(m).
• Anisotropic / CFD refinement → anisotropic_refine_2d / anisotropic_refine.
• Prismatic / boundary layers → extrude_prisms / boundary_layers.
• Quad-dominant surface → quad_dominant(m); hexahedral box → hexahedralize_box(lo, hi; nx, ny, nz).
• Curve a surface to the exact geometry → curve_surface(m, g).
• Voronoi cell volumes → voronoi_volumes(m).

Next steps

• Watertight Marching Cubes — why marching_cubes is watertight and how to verify it.
• Differentiable Surface Meshing — frozen marching cubes and how deform keeps gradients flowing.
• Mesh Quality — per-element quality and the threshold that drives re-meshing.
• Linear Elasticity with FEM — tetrahedral meshing straight into a solve.