Surfaces — meshes and implicit fields

Lessons 03 and 04 handed us aligned point clouds — the raw output of every depth sensor, back-projected and registered into one frame. But a point cloud is gappy: it samples the surface at scattered locations and says nothing about the continuous sheet between the samples, nothing about which side is "inside." A surface is exactly that missing structure. This lesson builds it two ways — explicit meshes you store directly and implicit fields whose level set is the surface — and the marching-cubes bridge that converts between them, making concrete the explicit↔implicit split we sketched abstractly in lesson 01.

The plan

Four moves. (1) Frame the two philosophies for turning points into surfaces — store the surface (explicit) vs store a function whose level set is the surface (implicit) — now with concrete objects instead of the lesson-01 abstraction. (2) Build the explicit side: triangle meshes, why triangles, normals as cross products, manifoldness, and Poisson surface reconstruction from the oriented normals lesson 03 gave us. (3) Build the implicit side: occupancy and signed-distance functions, the eikonal property ‖∇s‖ = 1, and why implicit fields are the natural home for learned shape (DeepSDF, Occupancy Networks). (4) The bridge: marching cubes, the algorithm that extracts an explicit mesh from an implicit field by linear interpolation on a grid — the shape-world counterpart of the \exp/\log dictionary from lesson 02.

1 · From points to surfaces — two philosophies

Lesson 01 split every 3D representation by one question: do you store the shape, or a function that tests the shape? With aligned point clouds now in hand (lessons 03–04), that split stops being philosophy and becomes a concrete engineering fork. A point cloud is neither answer — it is a bag of samples {(x,y,z)} with holes between them and no notion of interior. To use it — render it, measure it, collide against it, learn from it — you must commit to a surface. Two ways to commit:

Explicit — store the surface directly. Write down the surface as primitives with coordinates: a mesh of vertices joined into triangles. The surface is the data; there is no function to evaluate, you just read off the geometry. This is what graphics hardware draws, so it renders in real time.
Implicit — store a function whose level set is the surface. Keep a scalar field f:ℝ³→ℝ and define the surface as the set where f crosses a threshold (occupancy at ½, signed distance at 0). You never list the surface; you query the function at any point. This is smooth, differentiable, and — the property that made it eat modern shape learning — free of any fixed resolution or topology.

The trade-off, previewed here and paid off across sections 2–4:

Axis	Explicit (mesh)	Implicit (occupancy / SDF)
What is stored	the surface itself (verts + faces)	a function; surface = a level set
Rendering	fast — GPU rasterizes triangles natively	slow — ray-march / extract a mesh first
Optimization	hard — connectivity is discrete, non-differentiable	easy — smooth in inputs & weights
Topology change	rigid — merging/splitting needs re-meshing	free — any topology, zero bookkeeping
Memory	∝ surface area (#verts)	fixed (network weights), resolution-free
Query "am I inside?"	awkward (ray casting, winding number)	trivial (sign of f)

Neither wins outright — that is the whole point. You render explicit, you optimize implicit, and you convert between them constantly. Section 4's marching cubes is that converter; the widget makes it tangible in 2D. But first, each side on its own terms.

2 · Explicit: triangle meshes

A mesh is two lists: vertices V = {v_i ∈ ℝ³} and faces F, each face an ordered tuple of vertex indices. Almost universally the faces are triangles, and the reasons are not aesthetic:

A triangle is always planar. Three points define a plane exactly; there is no way to make a triangle non-flat. A quad with four points can be non-coplanar (a "bowtie" or a warped saddle), which makes its normal and its rasterization ambiguous. Triangles have no such degeneracy.
The GPU is built for them. The entire fixed-function and shader pipeline (lesson 06) rasterizes triangles; it is the hardware primitive. Anything else is triangulated first.
Barycentric interpolation is trivial. Any point inside a triangle is a convex combination p = αv₀ + βv₁ + γv₂ with α+β+γ = 1. Those same weights (α,β,γ) interpolate any per-vertex attribute — color, normal, texture coordinate, depth — across the face for free. Interpolation over a quad or n-gon has no such unique, linear answer.

Connectivity and adjacency

The raw vertex/face lists answer "where are the triangles?" but not "which triangles touch this edge?" — a query you need constantly (smoothing, subdivision, collapse, normal propagation). A half-edge (doubly-connected edge list) structure fixes this: each undirected edge is split into two oppositely-oriented half-edges, each storing its origin vertex, its face, its next half-edge around that face, and its twin on the neighboring face. From those pointers, "all faces around a vertex," "the two faces sharing an edge," and "walk the boundary of a hole" become O(1)-per-step traversals instead of scans.

Normals

A surface normal is the direction perpendicular to the surface — needed for shading, for orientation (which way is "out"), and as the sampled gradient Poisson reconstruction will consume below. On a mesh there are two:

Face normal. For a triangle with vertices v₀, v₁, v₂ (counter-clockwise as seen from outside), take two edge vectors and cross them:

n_f = \dfrac{(v₁ − v₀) × (v₂ − v₀)}{‖(v₁ − v₀) × (v₂ − v₀)‖} .

The cross product is perpendicular to both edges — hence to the triangle's plane — and the winding order (CCW) fixes its sign so it points outward. Its un-normalized length is twice the triangle's area, a fact we exploit next.

Vertex normal. A single vertex is shared by several faces; for smooth (Gouraud/Phong) shading we want one normal per vertex, a weighted average of the incident face normals:

n_v = \dfrac{\sum_{f ∋ v} w_f\, n_f}{‖\sum_{f ∋ v} w_f\, n_f‖},\qquad w_f = \text{face area or incident angle at } v .

Area-weighting falls out for free (use the un-normalized cross product, whose length is 2·\text{area}); angle-weighting (the interior angle of the face at v) is more robust to irregular tessellation, since it doesn't let one large triangle dominate a vertex it barely touches. Either way, averaging normals is what turns a faceted polyhedron into a visually smooth surface without adding geometry.

Manifoldness

A mesh is manifold if it locally looks like a flat sheet everywhere: every edge is shared by at most two faces, and the faces around any vertex form a single fan (or a single open fan on a boundary). Violations — an edge shared by three faces, or two cones meeting at a single "pinch" vertex (a non-manifold junction) — break the assumptions of nearly every downstream algorithm: normals become ill-defined, "inside" stops making sense, and half-edge traversal (which assumes exactly one twin) fails. Reconstruction and processing pipelines work hard to output watertight, manifold meshes for exactly this reason.

The trap: winding order and non-manifold output

Two silent mesh bugs bite constantly. (1) Inconsistent winding: if some triangles are CCW and others CW, their face normals point in random directions, so vertex-normal averaging cancels and shading goes black-and-blotchy — the mesh looks "inside-out" in patches. (2) Non-manifold junctions from naive stitching: an edge with three incident faces has no well-defined normal and no valid half-edge twin, so it crashes or corrupts every adjacency query. Always verify: consistent orientation, and every edge with ≤2 faces.

Surface reconstruction from an oriented point cloud

Now the concrete bridge from lessons 03–04: given a point cloud with oriented normals (lesson 03 estimated per-point normals; lesson 04 aligned the scans), how do we build a mesh? The dominant global method is Poisson surface reconstruction, and its idea is beautiful.

Define the indicator function χ(x) that is 1 inside the solid and 0 outside. Its gradient ∇χ is zero everywhere except at the surface, where it spikes in the direction of the surface normal. So the oriented point normals {n_i} are precisely samples of ∇χ — a noisy measurement of the gradient of the very indicator we want. Build a smooth vector field \vec V from those normal samples and solve for the scalar χ whose gradient best matches it, in the least-squares sense. Taking the divergence of both sides of ∇χ = \vec V turns it into a Poisson equation:

∇·∇χ = ∇·\vec V,\qquad \text{i.e.}\quad Δχ = ∇·\vec V ,

a single sparse linear system for χ over the domain. Solve it, then extract the surface as the level set χ = ½ (the halfway crossing between inside and outside) — which you do with marching cubes from section 4. Because it is one global solve driven by all the normals at once, Poisson reconstruction is smooth, watertight by construction (a level set of a continuous function has no holes), and robust to noise and non-uniform sampling — each noisy normal contributes only a little to a big system, so errors average out.

The local alternative is ball-pivoting: roll a ball of fixed radius ρ over the points; whenever it rests on three points without containing any other, emit that triangle, then pivot the ball around each edge to the next. It is fast, memory-light, and interpolates the actual samples exactly (no smoothing) — but being purely local it leaves holes where sampling is sparser than ρ and is sensitive to the radius choice, whereas Poisson fills gaps globally. Rule of thumb: Poisson for clean watertight models from noisy scans, ball-pivoting for faithfully interpolating dense, low-noise data.

3 · Implicit: occupancy and signed distance functions

The implicit side never writes down the surface. It stores a scalar field f:ℝ³→ℝ and defines the surface as a level set — the set of points where f hits a chosen iso-value. Two fields dominate.

Occupancy. o(x) = P(x \text{ is inside}) ∈ [0,1]. The surface is the decision boundary o(x) = ½: inside where o > ½, outside where o < ½. It answers "inside?" directly, but it is flat away from the boundary — deep inside and just-inside both read ≈1, so o tells you nothing about how far the surface is.

Signed distance function (SDF). s(x) = the signed distance to the nearest surface point: |s(x)| is the Euclidean distance, and the sign is negative inside, positive outside. The surface is the zero level set s(x) = 0. This is strictly richer than occupancy: it encodes distance, and that distance is what makes two things work. First, sphere-traced rendering (lesson 09) — standing at a point you know you can safely march a full |s(x)| along the ray without overshooting the surface, so you reach it in a handful of steps instead of tiny fixed increments. Second, normals for free: the gradient of a distance field points away from the surface along the shortest path, so

n(x) = \dfrac{∇s(x)}{‖∇s(x)‖} \Big|_{s=0} .

The eikonal property

A true distance function obeys a defining constraint almost everywhere:

‖∇s(x)‖ = 1 .

The reason is direct: if s measures distance to the surface, then moving one unit of length through space (in the direction away from the surface) changes the distance by exactly one unit — so the gradient's magnitude is 1. This is the eikonal equation. It holds everywhere except on the medial axis / surface itself, where the nearest-point is non-unique and s is non-differentiable. When you learn an SDF with a neural network there is nothing forcing the output to be a genuine distance — a raw MLP will fit the zero set but have arbitrary slope elsewhere, breaking sphere tracing and normals. So you add an eikonal loss, a soft penalty \mathbb E_x\big[(‖∇_x s_θ(x)‖ − 1)^2\big] that regularizes the field toward unit-gradient, i.e. toward behaving like a real SDF.

Why implicit is the natural home for learned shape

Four properties, each dodging a specific pain from the explicit side and from lesson 01's curse of dimension:

Continuous, resolution-free. The field is a function you can query at any real coordinate. There is no grid, so there is no resolution to pick and no aliasing — you evaluate exactly where you need the surface, at arbitrary precision.
Any topology for free. A shape that merges two parts into one, or grows a handle, is just a different level set of a slightly different field — no re-meshing, no connectivity edits. This is the implicit advantage lesson 01's widget dramatized, and the one this lesson's widget makes you feel: slide two blobs together and the surface merges with zero bookkeeping.
Differentiable by construction. An MLP is smooth in its inputs and its weights, so gradients flow — the property lesson 01 §4 argued the whole modern era needs. You can backprop an image or geometry loss straight into the shape.
Fixed memory. The shape is the network's weights, whose size is independent of resolution. This is the direct escape from the O(N³) voxel curse of lesson 01: a 1024³ grid is gigabytes, but an MLP that represents the same surface at unbounded resolution is a few MB.

Two landmark instantiations, both circa 2019:

DeepSDF. An MLP s_θ(x, z) mapping a query point x and a per-shape latent code z to a signed distance. It is trained as an auto-decoder: there is no encoder — each training shape owns a latent z that is optimized jointly with the shared weights θ, so the latent space becomes a learned prior over shapes you can interpolate and complete within.
Occupancy Networks. The same recipe with occupancy: o_θ(x, z) ∈ [0,1], trained as a binary classification of inside/outside at sampled points, surface read off at o = ½.

In both, the surface normal comes from autodiff: differentiate the network output with respect to its input coordinate, ∇_x s_θ — no finite differences, exact to machine precision, and (with the eikonal loss) already close to unit length.

4 · The bridge — marching cubes

An implicit field is ideal to optimize but you cannot hand a raster GPU a function — to display or export the surface you must extract an explicit mesh from the level set. Marching cubes is the standard algorithm, and it is a model of how a global geometric problem becomes a tiny per-cell lookup.

The recipe:

Sample the field f on a regular 3D grid of corner values.
Classify each cube. For each cube of 8 corners, label each corner inside (f < \text{iso}) or outside (f ≥ \text{iso}). Eight binary labels give a case index in 2^8 = 256 configurations — reduced by rotational and complementary symmetry to just 15 base cases, precomputed in a lookup table.
Locate crossings. The surface crosses an edge exactly when its two endpoints straddle the iso-value (opposite signs). Find where along the edge by linear interpolation: for endpoint values f₀, f₁ the iso-crossing sits at parameter

t = \dfrac{\text{iso} − f₀}{f₁ − f₀},\qquad p = p₀ + t\,(p₁ − p₀) .

Emit triangles. The case index selects, from the table, which crossing points to connect into triangles inside that cube. March to the next cube; shared edges reuse the same crossing, so the emitted triangles knit into a watertight mesh.

Two facts to hold onto. Resolution ↔ cost/smoothness: a finer grid resolves detail and produces a smoother polyline/mesh but multiplies cells (and, in 3D, at O(N³) — the lesson-01 curse returns for extraction). Ambiguous cases: certain corner patterns — the saddle configurations, where a face has two diagonally-opposite inside corners — admit two topologically different triangulations; if adjacent cubes resolve the same shared face inconsistently, the mesh gets holes or cracks. Robust variants (asymptotic decider, marching tetrahedra) enforce a consistent choice.

The reframe

Marching cubes is the explicit↔implicit dictionary — the shape-world analog of the \exp/\log maps from lesson 02. There, \exp took the flat, optimizer-friendly Lie-algebra coordinate onto the curved rotation manifold and \log read it back. Here, marching cubes is the \exp-like map from the optimizer-friendly implicit field to the render-friendly explicit mesh; Poisson reconstruction (§2) and sampling a mesh into an SDF run the other direction. You optimize in the representation that is smooth, and render in the one that is fast, converting across this bridge as needed.

The widget below is marching cubes in 2D — marching squares — where a cube is a 4-corner cell and a triangle is a line segment. Everything transfers: same corner classification, same linear-interpolation crossings, same case table, same saddle ambiguity.

Marching squares — the 2D marching cubes

An implicit field over the canvas: the signed-distance-like field of a union of circles (a min-of-SDFs / metaball), drawn as a diverging heatmap — blue inside (negative), red outside (positive), pale at the boundary. A uniform N×N grid overlays it, and the black polyline is the iso-contour extracted by true marching squares: classify each cell's 4 corners, then draw segments through the linearly-interpolated crossings t = (\text{iso} − f_a)/(f_b − f_a). Watch three things: the polyline sharpens as N grows (it only approximates the smooth level set); sliding the gap until the blobs touch merges two contours into one with zero bookkeeping — topology for free, the implicit advantage from lesson 01; and moving iso shows the field is a whole family of surfaces at once.

iso-level: 0 grid N: 20 gap between blobs: 150

Grid N

Segments emitted

—

Topology

—

iso value

Show the core JS

// implicit field: signed distance to a UNION of circles = min of per-circle SDFs
function field(px, py){
  var d = Infinity;
  for (var c = 0; c < circles.length; c++){
    var dc = Math.hypot(px - circles[c].x, py - circles[c].y) - circles[c].r;
    d = Math.min(d, dc);            // union of solids  ->  min of SDFs
  }
  return d;                          // <0 inside, magnitude ~ distance to surface
}

// linear interpolation of the iso-crossing on an edge between corner values fa, fb
function cross(fa, fb, pa, pb){
  var t = (iso - fa) / (fb - fa);    // where along the edge f == iso
  return [ pa[0] + t*(pb[0]-pa[0]), pa[1] + t*(pb[1]-pa[1]) ];
}

// MARCHING SQUARES: per cell, build a 4-bit case from the corners, draw the segment(s)
for (var i = 0; i < N; i++) for (var j = 0; j < N; j++){
  var f0 = F[i][j],   f1 = F[i+1][j],       // corner values (TL, TR,
      f2 = F[i+1][j+1], f3 = F[i][j+1];      //                 BR, BL)
  var code = (f0 < iso?1:0) | (f1 < iso?2:0) | (f2 < iso?4:0) | (f3 < iso?8:0);
  // edge crossings (only computed for edges whose endpoints straddle iso):
  var top = cross(f0,f1, TL,TR), rgt = cross(f1,f2, TR,BR),
      bot = cross(f2,f3, BR,BL), lft = cross(f3,f0, BL,TL);
  // the 16 cases -> which crossings to join. Saddle cases 5 & 10 are AMBIGUOUS;
  // resolve them the SAME way everywhere (here: split into the two "corner" arcs)
  // so neighbouring cells agree and the contour never cracks.
  segments.push.apply(segments, CASES[code]([top,rgt,bot,lft]));
}

Where this points next

We can now manufacture a continuous surface from gappy points — explicitly as a watertight mesh (Poisson from oriented normals), implicitly as an occupancy/SDF field you can learn and query anywhere — and cross between the two with marching cubes. But a representation is inert until you can turn it into an image, and every lesson so far has quietly assumed a renderer. Lesson 06 builds the one the whole field runs on: the classical rasterization pipeline that projects 3D triangles through the camera, decides which pixels each triangle covers, and resolves occlusion with a z-buffer — how a mesh becomes pixels, in real time. Lesson 07 then makes those pixels look right (perspective-correct interpolation, texture mapping and its aliasing, shading) and closes on the catch that reshapes the rest of the track: the hard rasterizer's coverage test is a discrete decision with no useful gradient (lesson 01 §4), so you cannot invert it by gradient descent. That single limitation is exactly what forces the soft, volumetric representations that follow — the radiance field of lesson 09, where a scene becomes colored semi-transparent fog rendered by an integral you can backpropagate. And it inherits directly from this lesson: NeRF's density field is a cousin of the occupancy field, and the eikonal-regularized SDF returns as the geometry backbone of the best NeRF-SDF hybrids.

Takeaway

A point cloud is gappy and has no notion of surface or inside; a surface supplies both, two ways. Explicit = store it as a triangle mesh (triangles because they are always planar, GPU-native, and barycentric-interpolable); normals are the normalized cross product n_f = (v₁−v₀)×(v₂−v₀) per face and an area/angle-weighted average per vertex; meshes must stay manifold (≤2 faces per edge); and you reconstruct one from oriented normals via Poisson — treat the normals as samples of ∇χ and solve Δχ = ∇·\vec V for a global, watertight χ = ½ level set. Implicit = store a field whose level set is the surface: occupancy at ½, or an SDF at 0 that also gives distance (sphere tracing) and normals (∇s) and obeys the eikonal constraint ‖∇s‖ = 1 (enforced as a training loss). Implicit fields are continuous, any-topology, differentiable, and fixed-memory — why DeepSDF and Occupancy Networks put learned shape here. Marching cubes bridges the two: sample on a grid, classify corners (256→15 cases), find edge crossings by linear interpolation t = (\text{iso}−f₀)/(f₁−f₀), emit triangles — the shape-world \exp/\log. Optimize implicit, render explicit, convert as needed.

Interview prompts

Contrast explicit meshes and implicit fields as surface representations, with the render/optimize/topology trade-offs. (§1 — explicit stores the surface (fast render, hard to optimize, rigid topology); implicit stores a function/level set (easy to optimize, any topology, resolution-free, slow to render).)
Why are meshes made of triangles specifically, and how do you compute a face vs a vertex normal? (§2 — triangles are always planar, GPU-native, and barycentric-interpolable; face normal = normalized (v₁−v₀)×(v₂−v₀), vertex normal = area/angle-weighted average of incident face normals.)
Explain Poisson surface reconstruction from an oriented point cloud. (§2 — treat oriented normals as samples of ∇χ of the indicator, solve the Poisson equation Δχ = ∇·\vec V globally, extract the χ = ½ level set; watertight and noise-robust. Ball-pivoting is the local alternative.)
What is a signed distance function, why is it richer than occupancy, and what is the eikonal property? (§3 — SDF gives signed distance (neg inside), enabling sphere tracing and normals via ∇s; occupancy is flat away from the boundary; eikonal: ‖∇s‖=1 a.e., enforced as an eikonal loss on learned SDFs.)
Why are implicit fields the natural representation for learned shape (DeepSDF / Occupancy Networks)? (§3 — continuous/resolution-free, any topology for free, differentiable by construction, fixed memory in the weights — dodging the O(N³) voxel curse; DeepSDF is an auto-decoder s_θ(x,z), normals via autodiff.)
Walk through marching cubes and name where it can fail. (§4 — sample on a grid, classify 8 corners (256→15 cases), interpolate edge crossings at t=(\text{iso}−f₀)/(f₁−f₀), emit triangles from the table; ambiguous saddle cases produce holes if neighboring cells resolve them inconsistently, and cost/detail scale with resolution.)