Marching Squares in 3D: Isosurface Scanner with Marching Tetrahedra, Euler Characteristic and Slice Monitor

Setting up the scanner...

Presets

Live readouts

Graphs

Marching-squares cases in the 2D slice: how many cells fall in each of the 16 corner patterns (0 and 15 are empty)

Contour length of the slice against the iso level (current field, current grid)

Convergence: relative error of the circle perimeter and of the sphere area against grid resolution N, with the 1/N^2 slope (dashed)

Try this

Metaballs, N = 12 then N = 64: thin bridges appear and disappear as the grid refines.
Turn off interpolation: the perimeter error jumps (check the Sphere preset).
Torus preset: V - E + F reads 0; move the iso level up until the hole closes and it becomes 2.

Numerical checks

running...

How to use

The scanner shows the isosurface of the 3D field inside the cube frame. The slice height z moves the copper scan ring and the translucent plane; the monitor shows the same 2D marching-squares slice, and the yellow line on the plane is that contour. Drag to orbit, wheel or pinch to zoom, Space pauses, R resets. Use the camera buttons for preset views; switch on the wireframe to see the individual triangles, or pick a preset such as the torus to watch the Euler characteristic go to 0. The turntable and cube frame were modelled in Blender (baked PBR textures); the mesh is generated live in JavaScript.

The science

The problem. A scalar field f(x, y) (or f(x, y, z)) is known only at the corners of a regular grid. The iso-contour f = c is a curve (a surface in 3D). Marching squares walks over every cell once and draws the piece of the curve inside it.

Marching squares. Each of the 4 corners is either inside (f > c) or outside, giving a 4-bit index and 16 cases. Cases 0 and 15 are empty; 1, 2, 3, 4, 6, 7, 8, 9, 11, 12, 13, 14 produce one segment; cases 5 and 10 (two opposite inside corners) produce two segments and are ambiguous: the corners alone do not say whether the two inside regions are joined. The usual fix is the asymptotic decider: use the average of the four corners as the value at the cell centre. Where an edge joins an inside and an outside corner, the crossing is placed by linear interpolation t = (c - f0) / (f1 - f0), P = P0 + t (P1 - P0) Without interpolation the crossing sits at the midpoint and the contour is visibly jagged.

Accuracy. With linear interpolation the position error of a smooth curve is O(h^2) with h = 2 / N, so doubling N cuts the error about four times. The convergence graph shows this slope and the checks measure it.

3D: marching tetrahedra. Classic marching cubes needs a 256-entry case table with its own ambiguities. Here every cube is split into 6 tetrahedra around its body diagonal; a tetrahedron has only a few cases (none, one corner cut off: 1 triangle, two corners: a quadrilateral = 2 triangles). The split is the same in every cube, so neighbouring cells agree on their faces and the surface is watertight. Vertices are shared via the grid edge they lie on, so V - E + F can be counted from the mesh: chi = V - E + F = 2 (1 - genus) per closed component (sphere 2, torus 0, two separate blobs 4). Normals come from the field gradient, n = -grad f / |grad f|, interpolated along the same edge.

Measured quantities. The surface area is the sum of triangle areas, the enclosed volume uses the divergence theorem V = (1/6) sum v0 . (v1 x v2). For the exact sphere (radius 0.6) and torus (R = 0.6, r = 0.25) these are compared with 4 pi r^2, 4/3 pi r^3 and 4 pi^2 R r.

Honest simplifications. (1) The fields are analytic functions sampled on a grid, not scanner data. (2) Tetrahedra make more (and thinner) triangles than marching cubes with the same N, but need no ambiguity table. (3) The 2D decider resolves case 5/10 from the centre average only; it is not guaranteed topologically correct for every field. (4) Features thinner than a cell are lost: raise N to see them return. (5) Linear interpolation of the gradient gives smooth shading but not exact normals.

Checks run in the page. Circle perimeter at N = 128 against 2 pi r; second-order convergence (error ratio 4 when N doubles) for both the perimeter and the sphere area; sphere volume; Euler characteristic of a sphere (2) and a torus (0) from the mesh; torus area; and chi = 2 x number of blobs for five metaballs.