serra

Analytical multi-material meshes from voxelized segmentations. It is named after the artist Richard Serra, who is known for making beautiful smooth geometric forms out of rusted metal.

serra turns a 3-D array of integer labels into one triangle mesh per label. It is built for connectomics-scale data, where a single chunk may contain hundreds of thousands of distinct objects, so it makes one pass over the volume regardless of how many labels are present.

Documentation at [https://alleninstitute.github.io/serra/]

pip install serra-mesh

Building from source needs a Rust toolchain: see the developer guide.

import serra_mesh

mesher = serra_mesh.Mesher(voxel_resolution=[4, 4, 40])
mesher.mesh(cutout)          # a 3-D array of integer labels
mesh = mesher.get(504)       # mesh.vertices, mesh.faces

What makes it different

serra uses multi-label surface nets (dual contouring) rather than marching cubes, following Frisken (2022) — see docs/references.md. Vertices sit inside cells at a position determined by where the label boundary crosses the cell, instead of being pinned to voxel-edge midpoints. This removes the staircase artifact, which shows up most clearly in surface area.

Measured on analytically-known spheres (isotropic voxels), against zmesh:

Marching-cubes area error does not shrink as resolution rises — it is a systematic bias, not a sampling error. With the optional relaxation pass enabled (relaxation=3) serra's area error drops to +0.38%, while volume stays within 0.2% of analytic in every case.

It also includes simplification and smoothing routines which maintain watertight and produce no non-manifold vertices or edges, while maintaining reproducible vertices at chunk boundaries to facilitate large scale mesh generation via chunking. Dual contouring reads two cell layers per face, so one voxel of halo is not enough — see docs/chunked.md.)

Smoothing is optional, bounded by max_deviation, and safe to run per chunk because seam vertices are pinned:

serra_mesh.Mesher(fairing=20, fairing_taubin=True)   # recommended
serra_mesh.Mesher(taubin=20)                         # same filter, per label

Prefer fairing. It is Frisken's surface fairing: one position per cell, shared by every label present there, rather than a private copy per label. The accuracy is identical — 0.068 voxels of mean surface error against an analytic tube either way — but smoothing each label separately pulls the two copies of a wall between touching objects apart, by up to 2.2 voxels on real neuropil, and the segmentation stops being a partition of space. Sharing the cell makes that impossible: the copies are one number, not two that started equal.

There is also relaxation=k, a plain Laplacian. It is not recommended: it shrinks, losing 2.5% of an object's volume at k=3, 7% at k=10, and 41% of a two-voxel-radius tube. It is kept because it is cheap and harmless on large objects. See docs/accuracy.md.

How it looks

Three objects from the 512³ MICrONS test volume, meshed by zmesh and by serra and rendered from an identical camera with flat shading. Flat shading is deliberate — it makes individual triangles visible, which is exactly what distinguishes a staircased surface from a smooth one.

The objects are sampled around the 80th size percentile rather than taken from the top: the largest object in a cutout is a cell body or a trunk crossing the whole box, and says little about the surfaces most objects get. The top row of each figure is the whole object, the bottom row a close-up spanning 1400 nm of the same surface.

	38K voxels, 47K faces
	42K voxels, 54K faces
	47K voxels, 57K faces

Marching cubes produces axis-aligned terraces because its vertices are pinned to voxel-edge midpoints. serra's are placed inside each cell from where the label boundary actually crosses it, so the terracing is gone even before relaxation; relaxation=3 removes the remaining faceting. Face counts are within 1% across all three, so this is not a resolution difference.

The same objects decimated 10× — the regime PyChunkedGraph actually stores, and the one where the input surface matters most, since a quadric simplifier keeps whatever the extractor gave it:

	47K → 4.7K faces
	54K → 5.4K faces
	57K → 5.7K faces

Against the mesh MICrONS publishes

A 5 µm cutout around segment 864691136144674612 at 32×32×40 nm, with both meshers decimated to the face count of the LOD-0 mesh the dataset actually serves for that segment — so all four panels are drawn on the same budget:

Regenerate with:

python bench/render_comparison.py --zmesh ../zmesh --out docs/images
python bench/render_comparison.py --zmesh ../zmesh --out docs/images --simplify 10
python bench/render_segment.py --zmesh ../zmesh --out docs/images

Performance

On the 512³ connectomics volume (2524 objects), Apple M4 Pro (14 cores), against zmesh 1.15.0 from PyPI. Each figure is the median of three runs, one process per implementation so peak RSS is not contaminated:

zmesh has no threading, so the honest single-threaded comparison is the first column, and zmesh traverses the volume about 2.2× faster there. serra needs more than one core to match it, reaching 0.47 s on 14. What serra buys for that is a vertex placement that marching cubes cannot express — see accuracy — and about 5× on extraction.

The extraction gap is architectural rather than a matter of tuning. Marching cubes emits a triangle soup: every triangle carries its own three vertices with no sharing. Turning that into an indexed mesh means deduplicating them, and zmesh does it per object with a hash map keyed on packed coordinates. On this volume that is 268.6M soup vertices collapsing to 45.0M unique — a 6× redundancy. The soup also drives the memory: 268.6M packed vertices is 2.1 GB held live while the map is being built.

serra never creates the duplicates. The extractor assigns one vertex per cell per connected component up front and quads reference those indices directly, so get() is a coordinate conversion and a triangulation.

Both therefore finish with about 45M unique vertices, which is the output row above. They differ in how many they touch on the way: zmesh generates 268.6M and collapses them, serra emits 44.8M and is done. Per vertex processed the two are within 15% of each other — roughly 20 ns for zmesh against 23 ns for serra — so the ~5× in the table is the 6× of intermediate vertices that only one of them ever creates, not a slower inner loop.

Reproduce with:

python bench/compare_zmesh.py serra
python bench/compare_zmesh.py zmesh

Controlling parallelism

serra_mesh.Mesher(threads=0)   # default: every core
serra_mesh.Mesher(threads=1)   # fully sequential
serra_mesh.Mesher(threads=4)   # exactly four

Set threads=1 if you are already parallelising at a higher level — one chunk per process in a pipeline, say — otherwise every process tries to claim every core and they fight each other.

Any value above 1 gets a private thread pool, so the setting is honoured exactly, is not overridden by RAYON_NUM_THREADS, and does not disturb other users of rayon in the same process. Only threads=0 defers to RAYON_NUM_THREADS. mesher.effective_threads reports what will actually be used.

Scaling on the volume above, and peak memory measured in a fresh process each time:

Output is byte-identical at every thread count, which the test suite checks directly rather than assuming. The volume is split into bands along one axis; a band cannot emit its own first cell layer's quads, since those read the layer below, so a short serial pass produces them afterwards and splices them into each label's face list in the position a single traversal would have put them.

Chunked meshing

Each chunk owns a disjoint range of voxels and is passed to mesh() with a 1-voxel halo on every side — so neighbouring input arrays overlap by 2 voxels. This is what makes the dual cells along a seam shared between both chunks, and therefore what makes their vertices identical.

One voxel of halo is enough at any relaxation setting. Iterative smoothing normally propagates one cell per iteration, which would mean k iterations need k + 1 voxels of halo. serra instead holds the outermost layer of cells fixed — precisely the vertices whose one-ring the chunk does not fully contain — so relaxation never reads past the halo. A chunk's mesh is therefore reproducible from that chunk's own array alone, whatever k is.

The trade-off is deliberate: a chunk's interior smooths slightly more than the band around its seams, so a stitched surface is self-consistent and watertight, but not identical to the same volume meshed in one piece.

Naming

The distribution is serra-mesh and the module is serra_mesh, because the name serra on PyPI was taken.

Status

Under active development. See docs/ for the user guide and the developer guide to the module layout.

Citing

The method is Frisken, S. F. (2022), SurfaceNets for Multi-Label Segmentations with Preservation of Sharp Boundaries, Journal of Computer Graphics Techniques 11(1), 34-54. For this implementation see CITATION.cff; the full reference list is in docs/references.md.

License

MIT