RasmahRasmah

Differentiable Surface Meshing

println("Julia ", VERSION)
Julia 1.12.6

Overview

The volume-meshing tutorials freeze marching-tetrahedra topology and deform it. The surface path has the exact same shape: freeze_marching_cubes captures a triangle mesh’s connectivity once, and deform realizes it at new parameters as a smooth function of the field. This tutorial covers that frozen surface pipeline, plus the two surface tools that make it cheap and exact — icosphere for spheres and the smooth (blended) CSG operators for differentiable joints.

icosphere: the exact sphere mesh

For a sphere, marching cubes is overkill. icosphere(n) builds an exact, watertight triangle sphere by subdividing a regular icosahedron n times — the vertices already sit on the surface, so there is no staircase to begin with:

using Rasmah

for n in (1, 2, 3)
    m = icosphere(n)
    println("level ", n, ": ", size(m.vertices, 2), " vertices, ", size(m.faces, 2), " faces")
end
level 1: 42 vertices, 80 faces
level 2: 162 vertices, 320 faces
level 3: 642 vertices, 1280 faces

Its area converges to 4π from below as the subdivision deepens, and it is differentiable by construction:

area_r(r) = surface_area(TriangleMesh(r .* icosphere(3).vertices, icosphere(3).faces))
println("A(1)    = ", area_r(1.0), "   (4π = ", 4π, ")")
println("dA/dr   = ", derivative(area_r, 1.0), "   (8π = ", 8π, ")")
A(1)    = 12.506492733969862 m^2   (4π = 12.566370614359172)
dA/dr   = 25.012985467939725 m^2   (8π = 25.132741228718345)

Frozen marching cubes

For a general SDF, freeze_marching_cubes(f, lo, hi, n) stores the grid edges that the surface crossed and the face connectivity. deform(fm, f) then re-interpolates each crossing point from the new field — moving vertices, but never touching the faces:

sdf(r) = sphere(r)
lo, hi = bounds(sdf(1.0))
margin = 0.05 .* (hi .- lo)

fm = freeze_marching_cubes(x -> field(sdf(1.0), x), lo .- margin, hi .+ margin, 24)

m0 = deform(fm, x -> field(sdf(1.0), x))
m1 = deform(fm, x -> field(sdf(1.2), x))
println("same connectivity: ", size(m0.faces, 2), " == ", size(m1.faces, 2), " triangles")
same connectivity: 4472 == 4472 triangles

Because the connectivity is fixed and the vertex positions are smooth in the field, the surface area is differentiable through the mesh — the surface counterpart of the frozen-tet pipeline:

area(r) = surface_area(deform(fm, x -> field(sdf(r), x)))

println("dA/dr (frozen MC) = ", derivative(area, 1.0))
println("analytic 8π       = ", 8π)
dA/dr (frozen MC) = 25.14289096939937 m^2
analytic 8π       = 25.132741228718345

Smooth (blended) boolean operations

Exact booleans use min/max, which have kinks where the two fields are equal — those kinks are non-smooth, so their gradient is discontinuous and can stall an optimizer. The smooth CSG operators replace the kink with a rounded blend controlled by a parameter k (k = 0 recovers the exact operation, larger k rounds the joint more):

a = sphere(1.0)
b = translate(sphere(1.0), 1.2, 0.0, 0.0)

println("smooth union        at origin: ", field(smooth_union(a, b, 0.2), [0.0, 0.0, 0.0]))
println("smooth intersection at midpt : ", field(smooth_intersection(a, b, 0.2), [0.6, 0.0, 0.0]))
smooth union        at origin: -1.0
smooth intersection at midpt : -0.35000000000000003

(smooth_difference was shown in the CSG tutorial.) The smooth blend keeps the field differentiable through the joint, which is what lets a gradient flow through a boolean without the topology changing.

Next steps