Unit Sphere

A sphere is a useful model because, given its radius, the surface area and volume are known quantities. The error of any surface that approximates the analytical surface of the sphere can be easily quantified.

This page builds a unit sphere from a voxel segmentation. Voxels inside the sphere have value 1. Voxels outside have value 0. Marching cubes turns the segmentation into a triangulated surface. That surface approximates the isosurface of the segmentation at the value 0.5, halfway between 0 and 1. We measure the triangulated surface against the exact sphere.

Related Studies

Two other sections build unit spheres by other methods.

• Octa-Loop, on the Subdivision page, refines a unit octahedron into a sphere by Loop subdivision.
• Remesh: Unit Sphere remeshes a unit sphere of 1,088 facets.

Segmentation

The script unit_sphere_segmentation.py builds a sphere of radius $n$ voxels. The sphere sits at the center of a cube of $2n+1$ voxels per side. A voxel is inside when its center satisfies $x^2+y^2+z^2\le n^2$. An inside voxel has value 1, and every other voxel has value 0. The array axes are $(x,y,z)$, the order automesh reads from a .npy file.

uv run --with numpy unit_sphere_segmentation.py

The script writes five segmentations of increasing resolution: $n$ = 10, 20, 40, 80, and 160. Each voxel has side length $1/n$ once the sphere is scaled to radius 1. The volume column counts the inside voxels and multiplies by $1/n^3$. The error, $(\text{volume}-4\pi /3)/(4\pi /3)$, compares that volume with the exact volume $4\pi /3\approx 4.1888$.

The voxel volume falls short at every $n$ in the table, and the shortfall shrinks as $n$ grows.

Marching Cubes

Marching cubes turns a grid of values into a triangulated surface at one level of those values. The command automesh mesh tri --cubes marching runs the method of Lewiner et al.¹ For each segmentation, automesh

1. pads the segmentation by one voxel of zeros, so the surface closes,
2. runs marching cubes at level 0.5, halfway between outside (0) and inside (1),
3. winds every triangle outward, and
4. scales and translates the vertices to a unit sphere (see Transform).

The command below meshes the $n=10$ segmentation. The scale is $1/n$ and the translation is $-(n+0.5)/n$ in each direction.

automesh mesh tri -i unit_sphere_n010.npy -o unit_sphere_mc_n010.stl --cubes marching \
    --xscale 0.1 --yscale 0.1 --zscale 0.1 \
    --xtranslate -1.05 --ytranslate -1.05 --ztranslate -1.05

The loop below makes the surface for each of the five segmentations.

for n in 10 20 40 80 160; do
  s=$(python3 -c "print(1/$n)")
  t=$(python3 -c "print(-($n + 0.5)/$n)")
  automesh mesh tri -i unit_sphere_n$(printf %03d $n).npy \
    -o unit_sphere_mc_n$(printf %03d $n).stl --cubes marching \
    --xscale $s --yscale $s --zscale $s \
    --xtranslate $t --ytranslate $t --ztranslate $t
done

The surfaces below have the following counts. $v$, $e$, and $f$ count the vertices, edges, and faces. The volume is the signed volume that the faces enclose. The error compares that volume with the exact volume $4\pi /3\approx 4.1888$.

The volume created from marching cubes underestimates the true $4\pi /3$ value at every $n$ in the table. The surface also encloses less volume than the voxels it came from. This pattern holds at every $n$ from 4 to 160. Marching cubes cuts each corner of the voxel staircase with a flat triangle, and every cut removes a small volume.

Five integrity checks run on the five surfaces. All five surfaces give the same result on every check.

• Closed. No edge is open. An open edge belongs to one face only.
• Manifold. No edge is pinched. A pinched edge belongs to more than two faces.
• Oriented. No two faces traverse the same edge in the same direction.
• Connected. The surface is one piece.
• Euler characteristic. $\chi =v-e+f=2$.

For a closed, connected, orientable surface, $\chi =2-2g$, where $g$ is the genus, the number of handles. Here $\chi =2$, so $g=0$, and each surface is topologically a sphere.

Convergence

The script unit_sphere_convergence.py computes both volumes at every $n$ from 4 to 160, 157 values in all. It runs automesh mesh tri --cubes marching once per $n$.

Both volumes settle onto $4\pi /3$, but neither settles smoothly. The voxel error changes sign 46 times, and 25 of the 157 voxel volumes exceed $4\pi /3$. The marching cubes error changes sign 24 times, and 12 of its volumes exceed $4\pi /3$. The five values in the tables just happened to fall short.

Under the oscillation, both errors fall toward zero. A least-squares fit of $\log |\text{error}|$ against $\log n$ gives a slope of −1.71 for the voxels and −1.76 for marching cubes. Both lie between the two guides. The marching cubes volume stays below the voxel volume at all 157 values.

Why Lewiner?

Lorensen and Cline introduced marching cubes in 1987.² Their method visits each cube of eight neighboring samples. It classifies each corner as inside or outside, which gives 256 cases. Symmetry reduces those to 15 base cases, and each base case has one fixed triangulation.

Some cases are ambiguous. On a face with two inside corners on one diagonal and two outside corners on the other, the corners alone cannot say whether the inside corners connect across the face. An interior ambiguity arises the same way inside a cube. A fixed table picks one answer without looking at the data. Neighboring cubes can then disagree, and the surface can take the wrong topology.

Chernyaev extended the table to 33 cases.³ His table resolves each ambiguity with the trilinear interpolant of the eight corner values, the smooth function the samples define inside the cube. Lewiner et al. gave an efficient and complete implementation of that table.¹ automesh uses the Lewiner method.

The Lewiner method is better than the original table because neighboring cubes always agree. Each ambiguity is resolved from the data, so two cubes that share a face choose the same connection across it. The surface then has the topology of the interpolant's level set.

The sphere reaches no ambiguous case, so the Lewiner method has nothing to resolve here, and no edge of the surface is pinched at any $n$. On data that do reach such cases, the original table can return a surface of different topology. Newman and Yi illustrate each ambiguous case and survey the methods that resolve them.⁴

A binary mask needs one further caution. On a mask of 0 and 1, the level 0.5 is exactly the saddle value of an ambiguous face. At that tie, the Lewiner method can pinch an edge between four faces. automesh uses the level 0.5. The sphere is a special case that never reaches the tie.

Transform

Marching cubes returns vertices in voxel units. Voxel $i$ spans $[i,i+1]$, so the voxel at index $n$ spans $[n,n+1]$. That voxel is the center of the sphere, so the center sits at its midpoint, $(n+0.5,n+0.5,n+0.5)$. Two operations map the surface to a unit sphere at the origin.

1. Scale by $1/n$, which makes the radius 1.
2. Translate by $-(n+0.5)/n$, which moves the center $(n+0.5)/n$ to the origin.

The --xscale and --xtranslate options of automesh scale first and translate second, in that order. The translation must use $n+0.5$. A translation of $-(n+1)/n$ treats the center as the upper face of voxel $n$, at $n+1$. It leaves the surface off center by $0.5/n$ in each direction.

The vertex radii show how close the surface comes to the unit sphere. The coefficient of variation, CoV $=\text{std}/\text{mean}$, measures how much the radii spread. A perfect sphere has CoV $=0$.

The mean radius stays within 0.1% of 1 at every $n$, so a scale of $1/n$ leaves no bias. The largest radius is always $1+0.5/n$, at the six poles. There a vertex sits halfway between the last inside voxel, at radius $n$, and the first outside voxel, at radius $n+1$.

The CoV halves each time $n$ doubles, so the spread of the radii falls as $1/n$. That is first order, slower than the volume error in Convergence. The volume averages over the whole surface, and the staircase steps partly cancel in that average. The radius spread sees every step.

Surface Comparison

Hausdorff Distance

The Hausdorff distance between two surfaces is the farthest that any point of either one lies from the other. For a surface $M$ and the unit sphere $S$,

$$d_H(M,S)=\max \left\{\underset{p\in M}{\sup }d(p,S),\underset{q\in S}{\sup }d(q,M)\right\}.$$

The script unit_sphere_comparison.py computes the first term exactly. A point $p$ lies $d(p,S)=||p|-1|$ from the unit sphere. Over a triangle, $|p|$ is convex. Its largest value sits at a vertex, and its smallest sits at the point of the triangle closest to the origin.⁵ The largest $|p|-1$ over all triangles gives the outward distance. The largest $1-|p|$ gives the inward distance.

The second term needs no computation. Each surface is closed and encloses the origin. The ray from the origin through a sphere point $q$ therefore crosses $M$ at some point $p$. Then $d(q,M)\le |q-p|=||p|-1|\le {\sup }_{p\in M}d(p,S)$. The second term never exceeds the first, so the Hausdorff distance is the larger of the outward and inward distances.

uv run --with numpy unit_sphere_comparison.py

The Hausdorff distance is exactly $0.5/n$ at every $n$, half a voxel. The outward distance sets it, at the six poles (see Transform). The inward distance approaches it from below. The Hausdorff distance falls as $1/n$, first order, like the radius CoV. The triangle count $f$ grows as $n^2$, so in terms of triangles the distance falls only as $f^{-1/2}$.

Surface Area

The volume converges, as the Convergence section shows. The surface area does not converge. The script unit_sphere_comparison.py adds up the areas of the marching cubes triangles. The error compares that area with the exact area $4\pi \approx 12.5664$.

The voxel staircase itself would overstate the area by 50% at any $n$. Marching cubes cuts the voxel staircase corners, which removes most of that excess. Howver, the surface area still settles near 8.8% too large.

Smaller voxels make smaller steps, but the triangles keep the same few orientations relative to the true surface. So the Hausdorff distance falls to zero while the area error does not. A surface can come arbitrarily close to the sphere and still have the wrong area.

Triangle Quality

The script unit_sphere_comparison.py also runs automesh metrics on each surface. The Triangular Metrics page defines its measures.

A binary mask puts every vertex at the midpoint of a cube edge. So marching cubes can make only a few triangle shapes. Grouped by maximum edge ratio and minimum scaled Jacobian, the triangular surfaces created by marching cubes produce exactly five shapes. The table gives the fraction of triangles of each shape.

Refinement does not change the quality. The same five shapes appear at every $n$, in nearly the same proportions. More than a quarter of the triangles have the worst scaled Jacobian, $1/\sqrt{3}\approx 0.577$, with a smallest angle of 30°. No triangle is worse.

Smaller voxels shrink the triangles but never reshape them. A histogram of any measure would show only these five spikes, so the table is shown instead of a histogram.

Mesh

Sculpt

The Subdivision page meshes its Octa-Loop spheres with Sculpt. At levels 3 to 7, it reports a grid of 26×26×26 cells, 7,731 nodes, and 6,672 elements. Here Sculpt meshes the marching cubes surfaces on the same grid.

The Subdivision baseline runs Sculpt through Cubit's sculpt parallel command, which chooses the grid itself. For Octa-Loop level 3, Cubit places 26 cells in each direction across a box of $\pm 1.240409$, a cell size of about 0.0954. Standalone Sculpt with those options reproduces the baseline exactly, at 7,731 nodes and 6,672 elements. The same options mesh each marching cubes surface.

B=1.240409
GRID="-x 26 -y 26 -z 26 -t -$B -u -$B -v -$B -q $B -r $B -s $B -SS 2"
# octa_loop03.stl downloads from the Refinement table of the Subdivision page
sculpt -j 1 -stl octa_loop03.stl $GRID -e unit_sphere_sculpt_octa03
for n in 010 020 040 080 160; do
  sculpt -j 1 -stl unit_sphere_mc_n$n.stl $GRID -e unit_sphere_sculpt_n$n
done

Left to choose its own grid, Cubit sizes the cells to the small marching cubes facets. For $n=40$ it picks 73×73×73 cells and makes 180,201 elements. Fixing the grid keeps the meshes comparable with the baseline.

The quality is the minimum scaled Jacobian, as Sculpt reports it.

From $n=20$ on, every surface gives the same 7,827 nodes and 6,768 elements, 96 more of each than the baseline, or about 1.4%. The worst element is better than the baseline's at every $n$.

The difference between the middle and right panels comes from the size of a voxel step against the size of a Sculpt cell, about 0.095.

• At $n=10$, a voxel is 0.1 wide, so each step is about as large as a cell. The mesh boundary follows the steps, and the sphere looks lumpy.
• At $n=160$, a voxel is 0.00625 wide, so a cell spans about 15 voxels. The steps are far smaller than the elements, and the mesh looks like the baseline.
• From $n=20$ on, a cell spans about two voxels or more, and the mesh stops changing.

The run used Sculpt 16.08, from Cubit 16.08. The files unit_sphere_sculpt_*.e.1.0 are Sculpt outputs, and the repository does not contain them. The scripts of this page read them from this directory, so regenerating the Sculpt figures and table requires running the Sculpt commands above. automesh reads Sculpt's Exodus files directly, and automesh metrics gives the same minimum (0.343) and mean (0.888) as Sculpt for the baseline mesh.

automesh

automesh mesh hex meshes the $n=10$ marching cubes surface in two ways. The uniform mesh uses the Sculpt cell size, $2\times 1.240409/26\approx 0.095416$, so its lattice matches the Sculpt grid. The adaptive mesh uses the defaults.

automesh mesh hex -i unit_sphere_mc_n010.stl -o uniform.inp \
  -u 0.095416 --metrics uniform.csv
automesh mesh hex -i unit_sphere_mc_n010.stl -o adaptive.inp \
  --metrics adaptive.csv

Taubin Smoothing

Taubin smoothing rounds the voxel steps of the $n=10$ surface. This section measures how the number of smoothing iterations changes the Minimum Scaled Jacobian of the uniform and adaptive meshes.

automesh smooth applies Taubin smoothing with its defaults. The pass band is 0.1 and the scale is $\lambda =0.6307$, which gives $\mu =-0.6732$. Each iteration is one Laplacian step. Odd-numbered iterations shrink the surface with $\lambda$. Even-numbered iterations inflate it with $\mu$. The sweep runs every integer from 0 to 100 iterations, every tenth from 110 to 300, and every fiftieth from 350 to 1,000. Each smoothed surface gets a uniform mesh with the Sculpt cell size and an adaptive mesh with the defaults.

automesh smooth -i unit_sphere_mc_n010.stl -o smooth.stl --iterations 260
automesh mesh hex -i smooth.stl -o uniform.inp -u 0.095416 --metrics uniform.csv
automesh mesh hex -i smooth.stl -o adaptive.inp --metrics adaptive.csv

Uniform. The minimum rises with scatter from 0.161 at 0 iterations to 0.613 at 260 iterations. The mean rises from 0.827 to 0.901. Up to 300 iterations the hex count stays between 5,747 and 6,150, against 6,024 without smoothing. The peak is a plateau: 0.613, 0.612, and 0.612 at 260, 270, and 280 iterations, against 0.574 and 0.575 at 240 and 250. The minimum then falls to 0.483 at 300 iterations and scatters between 0.459 and 0.510 up to 1,000.

Adaptive. The hex count falls from 8,243 at 0 iterations to 243 at 1 iteration. From 34 to 700 iterations it stays at 111 hexes. The minimum peaks at 0.564 at 11 iterations, on 207 hexes. The peak is a spike, not a plateau. The minimum is 0.367 at 10 iterations and 0.395 at 12. The next highest minimum is 0.544 at 51 iterations, on 111 hexes. The 111-hex meshes start at 0.323 at 34 iterations and scatter between 0.299 and 0.544 up to 700.

Surface drift. The defaults inflate this surface, because $|\mu |$ exceeds $\lambda$. The mean vertex radius grows from 0.9992 at 0 iterations to 1.0207 at 260 and 1.0862 at 1,000. The enclosed volume is 4.126 at 11 iterations, 1.5% below $4\pi /3\approx 4.1888$. It is 4.428 at 260 iterations, 5.7% above, and 5.322 at 1,000, 27.1% above. The peak search therefore stops at 300 iterations.

Odd and even counts. An odd count ends on a shrink step and an even count on an inflate step, so neighboring counts give different surfaces. The enclosed volume is 4.115 at 1 iteration and 4.149 at 2. The adaptive minimum is 0.328 at 1 iteration and 0.085 at 2. The uniform minimum is 0.223 at 3 iterations and 0.266 at 4.

None of the 270 meshes has an inverted hex.

Quality After Smoothing

The figure below compares the two best meshes of the sweep with the same meshes before smoothing and with the Sculpt baseline. The uniform mesh comes from the surface smoothed 260 times, with a minimum of 0.613. The adaptive mesh comes from the surface smoothed 11 times, with a minimum of 0.564. Without smoothing, the minima are 0.161 and 0.086. The Sculpt minimum is 0.388.

The quality is the minimum scaled Jacobian, as in the Sculpt table. The first two rows, the Sculpt meshes of the Sculpt section, are the reference. Smoothing raises the uniform minimum from 0.161 to 0.613 and the adaptive minimum from 0.086 to 0.564. Of the two best meshes, the uniform mesh has the higher minimum, by 0.049, and the higher mean, by 0.108.

The adaptive mesh has 207 hexes, 3.4% of the uniform count, so its histograms sit about a factor of 30 lower. Compare the shapes of the curves, not their heights.

The two minima differ in how robust they are. The uniform value belongs to a plateau: 0.613, 0.612, and 0.612 at 260, 270, and 280 iterations. The adaptive value belongs to a spike: 0.367 at 10 iterations, 0.564 at 11, and 0.395 at 12. A small change in the iteration count would lose the adaptive result.

Control Study

automesh mesh hex's adaptive octree combines two refinement signals:

1. curvature-driven sizing, gated by a chord-error tolerance (--tolerance, disabled by default), and
2. local-thickness sizing, from the shape diameter function of Shapira, Shamir, and Cohen-Or.⁶

No script on this page passes --tolerance, so curvature-driven sizing never activates here. Local thickness alone (not curvature) sets the octree's refinement target.

Taubin smoothing raised the Minimum Scaled Jacobian of the marching cubes meshes, from 0.161 to 0.613 for the uniform mesh and from 0.086 to 0.564 for the adaptive mesh. The control study asks whether a surface that is already smooth does better than that. The control is Octa-Loop level 4, a Loop subdivision of an octahedron from the Subdivision page. Its file, octa_loop04.stl, downloads from the Refinement table. The surface is closed, has 2,048 triangles, and its 1,026 vertices lie within 0.3% of the unit sphere.

Level 4 is the level whose triangle count is closest to that of the smoothed $n=10$ surface. Smoothing moves the vertices and leaves the count at 3,800. Level 3 has 512 triangles, level 4 has 2,048, and level 5 has 8,192. The mean triangle edge points the same way: 0.118 for level 4 against 0.092 for the smoothed surface, and about 0.059 for level 5, since each level halves the edge.

automesh mesh hex -i octa_loop04.stl -o octa_uniform.inp \
  -u 0.095416 --metrics octa_uniform.csv
automesh mesh hex -i octa_loop04.stl -o octa_adaptive.inp \
  --metrics octa_adaptive.csv

The smooth surface does not beat the best smoothed surface. The uniform mesh of Octa-Loop level 4 has a minimum of 0.518, which is 0.095 below the best smoothed marching cubes mesh, and a mean of 0.903, which is 0.002 above it. The adaptive mesh has a minimum of 0.558, 0.006 below, and a mean of 0.775, 0.018 below. The adaptive mesh has 111 hexes, against 207.

The comparison favors the smoothed surface, because its result is the maximum of 135 iteration counts and the Octa-Loop result is one run. Only 17 of the 135 counts give a uniform minimum above 0.518. Only 1 gives an adaptive minimum above 0.558, the best one itself. From 100 iterations on, the median minimum is 0.503 for the uniform mesh and 0.492 for the adaptive mesh, and the Octa-Loop values are above both. Against the unsmoothed marching cubes meshes, the Octa-Loop minima are 0.518 against 0.161 for the uniform mesh and 0.558 against 0.086 for the adaptive mesh.

A smooth input surface reaches, in one run and with no search over iteration counts, a minimum within 0.006 of the best smoothed adaptive result, and 0.095 below the best smoothed uniform result. None of the four meshes has an inverted hex.

1. Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira, and Geovan Tavares. "Efficient implementation of Marching Cubes' cases with topological guarantees." Journal of Graphics Tools 8(2) (2003) 1–15. https://doi.org/10.1080/10867651.2003.10487582 ↩ ↩2

2. William E. Lorensen and Harvey E. Cline. "Marching cubes: A high resolution 3D surface construction algorithm." ACM SIGGRAPH Computer Graphics 21(4) (1987) 163–169. https://doi.org/10.1145/37402.37422 ↩

3. Evgeni V. Chernyaev. Marching Cubes 33: Construction of topologically correct isosurfaces. Technical Report CN/95-17, CERN, 1995. ↩

4. Timothy S. Newman and Hong Yi. "A survey of the marching cubes algorithm." Computers & Graphics 30(5) (2006) 854–879. https://doi.org/10.1016/j.cag.2006.07.021 ↩

5. Christer Ericson. Real-Time Collision Detection. Morgan Kaufmann, 2005. Section 5.1.5, the closest point on a triangle to a point. ↩

6. Lior Shapira, Ariel Shamir, and Daniel Cohen-Or. "Consistent mesh partitioning and skeletonisation using the shape diameter function." The Visual Computer 24 (2008) 249–259. https://doi.org/10.1007/s00371-007-0197-5 ↩