Surface Reconstruction

Surface reconstruction recovers a continuous surface from discrete measurements of it. The input is most often a point cloud from a laser scanner, a depth camera or multi-view stereo, sometimes with a normal per point; it can also be a set of range images or depth maps with known sensor positions, or a function, sampled on a grid or learned by a network, whose level set is the surface. The output is usually a triangle mesh. The task sits at the end of most 3D reconstruction pipelines, where it is also called meshing, and it is ill-posed: infinitely many surfaces pass through or near any finite set of points, so every method adds assumptions, which the survey of Berger et al. (2016) calls priors. This article covers the problem, the main families of methods and how results are measured. Poisson surface reconstruction, marching cubes, the signed distance function and TSDF fusion have their own articles.

Problem Definition

Hoppe, DeRose, Duchamp, McDonald and Stuetzle (1992) stated the general form of the problem. Given points $X = \{\mathbf{x}_1, \dots, \mathbf{x}_n\} \subset \mathbb{R}^3$ on or near an unknown surface $M$, a compact, connected, orientable 2-manifold that may have a boundary, find a surface $M'$ that approximates $M$, without being told its topology, whether it has boundaries, or its geometry (§1.1–§1.2). They model the measurements as $\mathbf{x}_i = \mathbf{y}_i + \mathbf{e}_i$ with $\mathbf{y}_i \in M$, and make two conditions precise:

• the sample is $\delta$-noisy if $\lVert \mathbf{e}_i \rVert \le \delta$ for every $i$;
• a noiseless sample is $\rho$-dense if every sphere of radius $\rho$ centred on $M$ contains a sample point, and a noisy sample is $\rho$-dense if it perturbs such a sample by at most $\delta$.

Two limits follow from these definitions, as the authors point out: features of $M$ smaller than $\delta$ cannot be recovered, nor features in regions sampled too sparsely; and for a surface with a boundary, such as a sphere with a disc removed, a hole in the sample cannot be told apart from a hole in the surface. Any method that closes holes therefore assumes the surface is closed, and any method that keeps them assumes the gaps are real.

Two variants of the goal are common. Interpolating methods return a mesh whose vertices are input points; approximating methods return a surface that passes near the points without having to touch them, which allows them to average out noise. Kazhdan, Bolitho and Hoppe (2006, §2) group the approximating methods into global fits, which consider all the data at once, and local fits, which consider subsets of nearby points at a time.

Input and Output

Input. The minimum is a set of 3D points. Berger et al. (2016, §2.2) list three kinds of extra input that change what is possible:

• Normals, either unoriented (an axis per point, for example from a plane fit to its neighbours, with no inside or outside) or oriented (pointing consistently outward or inward). Oriented normals let a method define an inside and an outside, which signed-distance and indicator-function methods require.
• Scanner information: the range images’ pixel lattice, which gives neighbours and a way to orient normals, and the lines of sight from the sensor to each point, which mark space in front of the surface as empty.
• Images registered to the depth data, such as the colour images of an RGB-D camera.

A depth map with known intrinsics and pose is an organized point cloud plus a line of sight per pixel, so methods built for range images, such as TSDF fusion, take it directly. An implicit function $f$ sampled on a grid, whether a fused signed distance, a solved indicator function or the output of a network, is turned into a mesh by isosurface extraction, which is what marching cubes does.

Output. A triangle mesh, which may be required to be watertight (closed, so that it bounds a volume) or allowed to have boundaries; sometimes an implicit field, or, for strongly prior-driven methods, primitives, skeletons or deformed template models instead of a mesh (Berger et al., 2016, Table 1). 3D representations compares these forms.

Variants

• From unorganized points, the general problem above; the methods differ in whether they need oriented normals.
• From range images or depth maps, where the lattice and lines of sight are known. Each range image is a height field and easy to mesh on its own (Bernardini et al., 1999, §2.1); the work is in merging many of them.
• From an implicit function: the surface is already defined as a level set and only has to be extracted.
• From images, through multi-view stereo or neural fields, where reconstruction and meshing are often combined (3D reconstruction).
• Objects or scenes. Closed objects scanned from all sides suit watertight methods; open scenes, such as a room or terrain scanned from one side, do not.
• With learned priors, where a network trained on a shape collection completes partial or sparse input of the same kind of shape.

Challenges

Berger et al. (2016, §2.1) characterize point clouds by five artifacts that matter most to reconstruction:

• Nonuniform sampling. Density varies with distance and angle to the scanner and with the shape. A neighbourhood must be large enough to describe the local surface and small enough to keep its features; with varying density no single radius works, so $k$-nearest neighbours are common.
• Noise, points scattered near the surface; for some scanners it lies along the line of sight. The usual goal is a surface near the points that does not fit the noise.
• Outliers, points far from the surface, sometimes in dense clusters, for example from false matches on specular surfaces in multi-view stereo; they should not be used at all.
• Misalignment of scans that were registered imperfectly (point cloud registration), including drift in sequential scanning.
• Missing data, where the density is zero: occluded, absorbing or out-of-range regions.

Two further difficulties are geometric. Orientation: normals from local plane fits have arbitrary signs, and choosing consistent signs over the whole cloud is the hard part of signed methods (below). Thin structures and sharp features: Amenta, Bern and Kamvysselis (1998, §3.2) observe that the middle of a thin plate has low curvature yet must be sampled densely so that its two sides are resolved as separate surfaces, and that at sharp edges and corners their sampling condition asks for infinitely dense samples.

Approaches

Combinatorial methods

Combinatorial methods build the mesh from the points themselves, as a subset of the faces of the 3D Delaunay triangulation, the decomposition into tetrahedra whose circumspheres contain no other sample. Bernardini et al. (1999, §2.2) divide these interpolating methods into sculpting methods, which compute a tetrahedralization and remove tetrahedra to carve out the shape, and region growing methods, which start from a seed triangle and attach points to its boundary. Amenta et al. (1998, §2) describe the earlier constructions:

• Alpha shapes (Edelsbrunner and colleagues; Edelsbrunner and Mücke tried them for surface reconstruction) include a simplex if it has a circumsphere of radius at most $\alpha$ with no sample inside. The parameter must be chosen by experiment, and with varying density there may be no good single value.
• Delaunay sculpting (Boissonnat) progressively eliminates tetrahedra from the Delaunay triangulation according to their circumspheres.

The crust (Amenta et al., 1998) computes the Voronoi diagram of the samples. For each sample it takes two poles, the farthest vertex of its Voronoi cell and the farthest vertex on the other side of the surface; the poles lie near the medial axis, the set of points with more than one nearest surface point. It then triangulates the samples together with the poles and keeps the Delaunay triangles whose three vertices are all samples (§4.2). The authors present it as the first algorithm for this problem with provable guarantees: for a good enough sample of a smooth surface, the crust contains a mesh topologically equivalent to the surface and converges to it as density increases. A good sample, an $r$-sample, places a sample within $r$ times the distance to the medial axis of every surface point, so density must be high where curvature is high or other parts of the surface are close, and can be low elsewhere (§3.2). Their proofs hold for small values such as $r \le 0.06$, but they observe that $r = 0.5$ generally suffices in practice; their Figure 4 reconstructs a surface from about 12,000 samples and from about 3,000 thinned to a 0.5-sample. The output interpolates the samples, needs no parameter, and is not necessarily a manifold (§4.2).

Ball pivoting (Bernardini, Mittleman, Rushmeier, Silva and Taubin, 1999) avoids computing the Delaunay triangulation, which they describe as prohibitively expensive and numerically unstable for millions of points (§2.2). Three points form a triangle if a ball of user-chosen radius $\rho$ touches them without containing another point. From a seed triangle, the ball pivots around each boundary edge, keeping contact with its two endpoints, until it touches a new point, which forms the next triangle; when no edge remains, a new seed is chosen (abstract, §3). The result is a manifold subset of an alpha shape, and so of the Delaunay triangulation (§1.1, §3). The method uses normals, assumed outward, to reject a triangle whose normal disagrees with the points’ (§3). Its failure modes follow from the ball: gaps wider than the ball leave holes, noisy samples below the surface are never touched and are dropped, and misaligned scans separated by about the ball’s size produce small spurious components (§3, Fig. 4). Repeating the pivoting with larger balls handles uneven density. They report linear time on datasets of millions of points and built a 14-million-triangle model of Michelangelo’s Florentine Pietà from more than 700 scans (§1, Fig. 1).

Kazhdan et al. (2006, §2) note the shared weakness: because these schemes interpolate all or most of the points, noisy data give jagged surfaces that are smoothed or refitted afterwards. Berger et al. (2016, §1.1) leave interpolating methods out of their survey for a related reason: their guarantees require sampling too dense for real scans with significant imperfections.

Signed distance from tangent planes

Hoppe et al. (1992) approximate instead. They estimate a signed distance $f$ to the unknown surface in a region near the data and contour its zero set (§3.1). They chose a signed rather than an unsigned distance because zero is a regular value of $f$ but not of $|f|$, so the implicit function theorem makes the zero set a manifold. The algorithm has four steps:

1. Tangent planes. For each point, the $k$ nearest points form its neighbourhood; the centre $\mathbf{o}_i$ is their centroid and the normal $\hat{\mathbf{n}}_i$ the eigenvector of their covariance with the smallest eigenvalue, a least-squares plane fit (point cloud describes the fit) (§3.2).
2. Orientation. The sign of each normal is chosen so that nearby planes agree. Maximizing $\sum \hat{\mathbf{n}}_i \cdot \hat{\mathbf{n}}_j$ over the edges of a proximity graph is NP-hard, by reduction to MAX-CUT, so the authors propagate signs instead (§3.3). Their graph is the Riemannian graph: the Euclidean minimum spanning tree of the centres, enriched with an edge between $i$ and $j$ whenever either centre is among the other’s $k$ nearest. Propagating along plain proximity went wrong at sharp features, the tips of a cat’s ears in their Figure 3b, so each edge gets the cost $1 - |\hat{\mathbf{n}}_i \cdot \hat{\mathbf{n}}_j|$, small for nearly parallel planes, and signs travel along the minimum spanning tree of these costs. This order tends to follow directions of low curvature. The plane with the highest centre is forced to point toward $+z$, and a depth-first traversal flips each plane that disagrees with its parent.
3. Signed distance. For a query point $\mathbf{p}$, find the plane whose centre $\mathbf{o}_i$ is nearest and take $f(\mathbf{p}) = (\mathbf{p} - \mathbf{o}_i) \cdot \hat{\mathbf{n}}_i$, the signed distance to that plane (§3.4). If the projection $\mathbf{z}$ of $\mathbf{p}$ onto the plane lies farther than $\rho + \delta$ from every sample, $f(\mathbf{p})$ is undefined: on a $\rho$-dense, $\delta$-noisy sample, such a point cannot be on the surface. The partition of space by nearest centre is a Voronoi diagram of the centres, so $f$ is discontinuous across cell boundaries, which the authors note does not affect the contouring.
4. Contouring. A variant of marching cubes that triangulates within tetrahedral decompositions of the cubes visits only cubes crossed by the zero set, and emits nothing in a cube with an undefined corner, which leaves boundaries where the data end. Cube edges shorter than $\rho + \delta$ estimate boundaries accurately. Edge collapses then remove badly shaped triangles (§3.5).

The examples in the paper used $k$ from 8 to 40, with $\rho + \delta$ set to infinity for surfaces known to be closed (§4, table). The authors found $k$ uncritical on data with little noise, and proposed adapting it to noise and density by monitoring the covariance eigenvalues (§5.1); they also noted that range data could orient the planes from the viewing directions. The method is local: Berger et al. (2016, §4.1) describe it as hugely influential on local smoothness methods, and note that noisy or inverted normals give inaccurate distances and that the nearest-plane rule gives a noisy output under nonuniform sampling. Other local methods blend the influence of several nearby points for continuity, for instance by moving least squares (Kazhdan et al., 2006, §2). The signed distance function article covers the field itself; the example below measures this pipeline.

Global implicit fits

Global methods fit one implicit function to all the points. Radial basis function methods, such as that of Carr and colleagues (2001), write the function as a low-degree polynomial plus a weighted sum of radially symmetric functions centred at the points, constrained to zero at the samples; off-surface points displaced along the oriented normals, with positive and negative values, rule out the trivial zero function (Berger et al., 2016, §4.2). The result is globally smooth and closes holes, but the ideal basis functions are globally supported, so the linear system is dense and ill-conditioned, and practical solutions on large datasets use adaptive RBF reduction and the fast multipole method (Kazhdan et al., 2006, §2).

Poisson reconstruction (Kazhdan et al., 2006) treats oriented points as samples of the gradient of the solid’s indicator function and solves for that function, which reduces to a Poisson equation; the surface is a level set of the solution. The authors argue that it combines the benefits of the two families: it is global, needing no heuristic neighbourhoods or blending, yet its basis functions live on an octree in space, are locally supported and give a sparse, well-conditioned system. Because the gradient is constrained everywhere, to zero away from the samples, spurious sheets seldom arise (§1). Screened Poisson reconstruction (Kazhdan and Hoppe, 2013) adds the points as interpolation constraints, which counters the over-smoothing several researchers had noted, keeps the system’s sparsity, and makes the solver linear in the number of points (abstract, §1). Both produce watertight surfaces and so invent surface over missing regions; the Poisson surface reconstruction article covers the formulation, parameters, trimming and an example with holes and flipped normals.

Volumetric fusion of range images

When the input is a set of aligned range images, the sensor positions are known. Curless and Levoy (1996) scan-convert each range image into a distance function, add it to a cumulative weighted signed distance on a voxel grid one image at a time, and extract the isosurface; under their assumptions it is optimal in the least-squares sense (abstract). Lines of sight also mark space as empty. Every voxel is classified as unseen, empty or near the surface, and a surface placed where empty space meets unseen space plugs holes that are hard to repair on the mesh, since holes are often non-planar and may need to join separate components (§4). TSDF fusion covers the method and KinectFusion its real-time form with a depth camera; Berger et al. (2016, §3.2) cite it as an example of visibility priors.

Bernardini et al. (1999, §2.1) describe the other route of the time, mesh stitching, which meshes each range image, finds overlaps, averages points there and stitches or zippers the meshes together, and they point out that the estimation of surface points need not be tied to the method that builds the mesh.

Learned implicit methods

Learned methods put the prior in a network trained on a collection of shapes. DeepSDF (Park et al., 2019) represents a class of shapes by one network that maps a latent code and a 3D point to a signed distance. Inference can use any number of SDF samples: completing a shape means solving for the code that best explains a partial observation, and the decoder supplies the rest (§6.3). From a single-view depth image, they create two samples per depth pixel, a small distance $\eta$ in front of and behind the surface along the estimated normal, and add free-space constraints between the surface and the camera. On partial range scans of unseen ShapeNet chairs, planes and sofas, their Table 4 shows DeepSDF better than 3D-EPN, a voxel-based completion network, on every listed measure; the depth must be given in the shape’s canonical frame.

Occupancy Networks (Mescheder et al., 2019) represent the surface as the decision boundary of a network that classifies points as inside or outside, conditioned on an encoding of the input. For reconstruction from point clouds, they sampled 300 points from each watertight ShapeNet model and added Gaussian noise of standard deviation 0.05. Their Table 2 reports a volumetric IoU of 0.778, a Chamfer-$L_1$ distance of 0.079 and a normal consistency of 0.895, against 0.674, 0.117 and 0.848 for Deep Marching Cubes and 0.565, 0.169 and 0.719 for 3D-R2N2, a voxel network (§4.3).

Both methods are meshed with marching cubes. Later work replaced the single global code with local features, fitted networks directly to raw points, and learned surfaces from images through differentiable rendering; neural implicit representations covers these, and neural radiance fields the density fields from which surfaces must be extracted afterwards.

Priors beyond smoothness

Berger et al. (2016, §3) organize the field by prior. Surface smoothness comes in local forms that smooth near the data and keep boundaries where sampling is missing (Hoppe et al.), global forms that seek smoothness at large scale and produce watertight surfaces of individual objects (radial basis functions, Poisson), and piecewise forms that keep sharp features of CAD models and man-made structures. Other priors are visibility (lines of sight), volumetric smoothness, geometric primitives, global regularity (symmetry, repetition, canonical relationships such as parallel and orthogonal planes), data-driven priors from shape collections, and user interaction. The survey argues that priors evolved with acquisition: smoothness for small objects from desktop scanners, then structural and data-driven priors once mobile real-time scanners captured whole scenes (§1).

How the field moved

In outline: the 1992 formulation removed assumptions about the data’s structure, such as contour stacks or single range images (Hoppe et al., §1). Range scanners of the 1990s made merging many aligned scans the practical problem; volumetric and mesh-stitching methods averaged redundant scans to correct error along the sensor’s line of sight, assuming negligible registration error (Bernardini et al., 1999, §2.1). Computational geometry supplied interpolating methods with guarantees, whose sampling conditions real scans rarely met (Berger et al., 2016, §1.1). Global implicit fits made watertight approximation of noisy, misregistered data robust and fast, and the survey singles out Poisson reconstruction as widely used because its code is stable and reliable (§11). The lines of sight of laser triangulation scanners brought visibility priors, and the Kinect’s real-time depth and colour brought data-driven methods for indoor scenes (§12); mobile scanning of whole scenes prompted structural priors (§1); and learned shape models extended data-driven priors to completion from sparse or partial input.

Example

The script reconstructs a known shape with the tangent-plane method of Hoppe et al. and measures what orientation, density and noise do to the result. The shape is a thin rounded plate: every point at distance 0.08 from a box with half-sizes $(0.6, 0.4, 0.04)$, so 1.36 long, 0.96 wide and 0.24 thick, with an exact signed distance (negative inside) and outward normals. Points are drawn with a density that falls tenfold from $x = -0.68$ to $x = +0.68$, then perturbed by Gaussian noise of standard deviation $\sigma$ per coordinate. All orientation rules share the same tangent planes ($k = 15$) and Riemannian graph:

• Hoppe: minimum spanning tree with cost $1 - |\hat{\mathbf{n}}_i \cdot \hat{\mathbf{n}}_j|$;
• proximity: the same propagation with Euclidean edge length as the cost, our reading of the “geometric proximity” order whose failure the authors show in their Figure 3b;
• oracle: each fitted normal signed by the true normal, which separates errors of orientation from errors of the distance estimate.

The field $f$ is evaluated on the whole grid box around the plate (spacing 0.0125, about 0.12 of margin), with $\rho + \delta$ infinite as the authors set it for closed surfaces, and contoured with scikit-image’s marching_cubes at level 0 rather than the authors’ tetrahedral variant, without edge collapses. Each mesh is sampled with 100,000 points and scored against 100,000 samples of the true surface by the Chamfer-$L_1$ distance (the mean of the two mean nearest-neighbour distances), the F-score at $\tau = 0.01$ and $0.02$, and the normal consistency (NC), defined below. Two columns count tree edges, per seed: across, edges whose two true normals are more than 90° apart, such as edges from the top sheet to the bottom or around the rounded rim; ambiguous, edges whose fitted normals have $|\hat{\mathbf{n}}_i \cdot \hat{\mathbf{n}}_j| < 0.5$, along which a sign is copied between planes far from parallel. Each configuration runs three seeds (0, 1 and 2). For a given seed and $n$, the sample positions are the same at every $\sigma$, because positions are drawn before the noise.

import numpy as np
from scipy.spatial import cKDTree, Delaunay
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import minimum_spanning_tree, depth_first_order
from skimage.measure import marching_cubes

# Known shape: a thin rounded plate, all points at distance RAD from a box of half-sizes B.
B = np.array([0.6, 0.4, 0.04]); RAD = 0.08         # extents x +-0.68, y +-0.48, z +-0.12

def sdf(p):                                         # exact signed distance, negative inside
    q = np.abs(p) - B
    return np.linalg.norm(np.maximum(q, 0), axis=1) + np.minimum(q.max(axis=1), 0) - RAD

def grad(p, e=1e-5):                                # unit outward normal
    g = np.stack([(sdf(p + e * d) - sdf(p - e * d)) / (2 * e) for d in np.eye(3)], 1)
    return g / np.linalg.norm(g, axis=1, keepdims=True)

def sample_mesh(v, f, n, rng, weight=None):
    """n points on triangle mesh (v, f), with density proportional to area (times weight)."""
    a, b, c = v[f[:, 0]], v[f[:, 1]], v[f[:, 2]]
    area = 0.5 * np.linalg.norm(np.cross(b - a, c - a), axis=1)
    w = area if weight is None else area * weight((a + b + c) / 3)
    t = rng.choice(len(f), n, p=w / w.sum())
    r1, r2 = rng.uniform(size=(2, n)); s = np.sqrt(r1)
    return (1 - s)[:, None] * a[t] + (s * (1 - r2))[:, None] * b[t] + (s * r2)[:, None] * c[t], t

def project(p):                                     # onto the exact surface
    return p - sdf(p)[:, None] * grad(p)

# Fine mesh of the true surface, used only to draw samples (then projected exactly).
h = 0.005
ax = [np.arange(-0.72, 0.72 + h / 2, h), np.arange(-0.52, 0.52 + h / 2, h), np.arange(-0.16, 0.16 + h / 2, h)]
G = np.stack(np.meshgrid(*ax, indexing="ij"), -1)
fv, ff, _, _ = marching_cubes(sdf(G.reshape(-1, 3)).reshape(G.shape[:3]), 0.0, spacing=(h, h, h))
fv += np.array([a[0] for a in ax])
ramp = lambda x: 1 - 0.9 * (x[:, 0] + 0.68) / 1.36  # density falls 10:1 from x = -0.68 to +0.68

# --- Hoppe et al. (1992) -------------------------------------------------------------
def tangent_planes(x, k=15):
    """Centre o_i (centroid of the k-neighbourhood) and unoriented normal (PCA), section 3.2."""
    _, idx = cKDTree(x).query(x, k=k)
    nb = x[idx]; o = nb.mean(1)
    C = np.einsum("nki,nkj->nij", nb - o[:, None], nb - o[:, None])
    return o, np.linalg.eigh(C)[1][:, :, 0]

def riemannian_edges(o, k=15):
    """EMST of the centres (from Delaunay edges) plus symmetric k-NN edges, section 3.3."""
    tri = Delaunay(o).simplices
    e = np.unique(np.sort(np.vstack([tri[:, [i, j]] for i in range(4) for j in range(i + 1, 4)]), 1), axis=0)
    n = len(o)
    L = np.linalg.norm(o[e[:, 0]] - o[e[:, 1]], axis=1)
    emst = minimum_spanning_tree(coo_matrix((L, (e[:, 0], e[:, 1])), shape=(n, n))).tocoo()
    _, idx = cKDTree(o).query(o, k=k)
    knn = np.stack([np.repeat(np.arange(n), k - 1), idx[:, 1:].ravel()], 1)
    return np.unique(np.sort(np.vstack([np.stack([emst.row, emst.col], 1), knn]), 1), axis=0)

def orient(o, nrm, E, rule):
    """Propagate a sign along the MST of the Riemannian graph, depth first from the top."""
    n = len(o)
    if rule == "Hoppe":                             # cost 1 - |n_i . n_j|: nearly parallel first
        w = 1 - np.abs(np.sum(nrm[E[:, 0]] * nrm[E[:, 1]], 1))
    else:                                           # "proximity": Euclidean edge length
        w = np.linalg.norm(o[E[:, 0]] - o[E[:, 1]], axis=1)
    T = minimum_spanning_tree(coo_matrix((w + 1e-12, (E[:, 0], E[:, 1])), shape=(n, n)))
    root = int(np.argmax(o[:, 2]))                  # largest z, normal forced to +z
    nrm = nrm.copy()
    if nrm[root, 2] < 0: nrm[root] *= -1
    order, pred = depth_first_order(T, root, directed=False, return_predecessors=True)
    for j in order[1:]:
        if np.dot(nrm[pred[j]], nrm[j]) < 0: nrm[j] *= -1
    return nrm, np.stack([order[1:], pred[order[1:]]], 1)   # tree edges (child, parent)

H = 0.0125                                          # contouring grid, whole box around the plate
AX = [np.arange(-0.8, 0.8 + H / 2, H), np.arange(-0.6, 0.6 + H / 2, H), np.arange(-0.25, 0.25 + H / 2, H)]
GRID = np.stack(np.meshgrid(*AX, indexing="ij"), -1).reshape(-1, 3)
SHAPE = tuple(len(a) for a in AX)

def hoppe_mesh(o, nrm):
    """f(p) = (p - o_i) . n_i for the nearest centre o_i (section 3.4), zero set by marching cubes."""
    _, i = cKDTree(o).query(GRID)
    f = np.einsum("ij,ij->i", GRID - o[i], nrm[i]).reshape(SHAPE)
    v, fa, _, _ = marching_cubes(f, 0.0, spacing=(H, H, H))
    return v + np.array([a[0] for a in AX]), fa

# --- Evaluation ----------------------------------------------------------------------
R, _ = sample_mesh(fv, ff, 100000, np.random.default_rng(1000)); R = project(R); Rn = grad(R)

def scores(P, Pn):
    da, ia = cKDTree(R).query(P); db, _ = cKDTree(P).query(R)
    f = lambda t: 2 * (da < t).mean() * (db < t).mean() / max((da < t).mean() + (db < t).mean(), 1e-12)
    return 0.5 * (da.mean() + db.mean()), f(0.01), f(0.02), np.abs(np.sum(Pn * Rn[ia], 1)).mean()

def evaluate_mesh(v, fa):
    a, b, c = v[fa[:, 0]], v[fa[:, 1]], v[fa[:, 2]]
    tn = np.cross(b - a, c - a); tn /= np.linalg.norm(tn, axis=1, keepdims=True) + 1e-30
    P, t = sample_mesh(v, fa, 100000, np.random.default_rng(7))
    return scores(P, tn[t])

P2, _ = sample_mesh(fv, ff, 100000, np.random.default_rng(1001)); P2 = project(P2)
cd, f1, f2, nc = scores(P2, grad(P2))
print(f"floor (second exact sampling): CD-L1 {cd:.4f}  F@0.01 {f1:.3f}  F@0.02 {f2:.3f}  NC {nc:.3f}")

rules = ["Hoppe", "proximity", "oracle"]
for n in [1000, 3000, 10000]:
    y0, _ = sample_mesh(fv, ff, n, np.random.default_rng(0), ramp); y0 = project(y0)
    d15 = cKDTree(y0).query(y0, k=15)[0][:, -1]
    print(f"\nn = {n}: median distance to 15th neighbour, x < -0.4: {np.median(d15[y0[:, 0] < -0.4]):.3f}, "
          f"x > 0.4: {np.median(d15[y0[:, 0] > 0.4]):.3f}  (plate thickness 0.24)")
    print(f"{'sigma':>6s} {'rule':9s} {'orient. err % (range)':>22s} {'CD-L1 (range)':>22s} {'F@0.01':>7s} {'F@0.02':>7s} {'NC':>6s}"
          f"  {'across':>8s}  {'ambiguous':>11s}")
    for sigma in [0.0, 0.005, 0.01, 0.02]:
        res = {r: [] for r in rules}; tree = {r: [] for r in rules}
        for seed in range(3):
            rng = np.random.default_rng(seed)
            y, _ = sample_mesh(fv, ff, n, rng, ramp); y = project(y); ny = grad(y)
            x = y + rng.normal(scale=sigma, size=y.shape)
            o, nrm = tangent_planes(x)
            E = riemannian_edges(o)
            for r in rules:
                if r == "oracle":
                    on = nrm * np.sign(np.sum(nrm * ny, 1))[:, None]
                else:
                    on, T = orient(o, nrm, E, r)    # tree edges across the plate / with ambiguous sign
                    tree[r].append((np.sum(np.sum(ny[T[:, 0]] * ny[T[:, 1]], 1) < 0),
                                    np.sum(np.abs(np.sum(nrm[T[:, 0]] * nrm[T[:, 1]], 1)) < 0.5)))
                err = 100 * (np.sum(on * ny, 1) < 0).mean()
                res[r].append((err, *evaluate_mesh(*hoppe_mesh(o, on))))
        for r in rules:
            a = np.array(res[r])
            print(f"{sigma:6.3f} {r:9s} {a[:,0].mean():6.1f} ({a[:,0].min():4.1f}-{a[:,0].max():4.1f})"
                  f"      {a[:,1].mean():.4f} ({a[:,1].min():.4f}-{a[:,1].max():.4f})"
                  f" {a[:,2].mean():7.3f} {a[:,3].mean():7.3f} {a[:,4].mean():6.3f}"
                  + ("  " + "/".join(str(t[0]) for t in tree[r]).rjust(8) + "  " + "/".join(str(t[1]) for t in tree[r]).rjust(11)
                     if tree[r] else f"  {'-':>8s}  {'-':>11s}"))

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; about a minute on a laptop CPU):

floor (second exact sampling): CD-L1 0.0029  F@0.01 1.000  F@0.02 1.000  NC 1.000

n = 1000: median distance to 15th neighbour, x < -0.4: 0.098, x > 0.4: 0.198  (plate thickness 0.24)
 sigma rule       orient. err % (range)          CD-L1 (range)  F@0.01  F@0.02     NC    across    ambiguous
 0.000 Hoppe       14.5 ( 0.0-43.5)      0.0190 (0.0132-0.0300)   0.579   0.741  0.839     0/1/0        0/0/0
 0.000 proximity    0.0 ( 0.0- 0.0)      0.0133 (0.0129-0.0138)   0.633   0.808  0.900     0/0/0        0/3/1
 0.000 oracle       0.0 ( 0.0- 0.0)      0.0133 (0.0129-0.0138)   0.633   0.808  0.900         -            -
 0.005 Hoppe       14.3 ( 0.0-42.9)      0.0192 (0.0136-0.0299)   0.570   0.741  0.836     0/1/0        0/0/0
 0.005 proximity    1.8 ( 0.0- 5.4)      0.0157 (0.0134-0.0201)   0.603   0.784  0.862     0/0/0        1/2/2
 0.005 oracle       0.0 ( 0.0- 0.0)      0.0137 (0.0134-0.0142)   0.621   0.806  0.893         -            -
 0.010 Hoppe       16.3 ( 0.0-49.0)      0.0204 (0.0149-0.0312)   0.554   0.728  0.820     0/1/0        0/0/0
 0.010 proximity    0.0 ( 0.0- 0.0)      0.0147 (0.0139-0.0151)   0.608   0.796  0.880     0/1/0        2/2/2
 0.010 oracle       0.0 ( 0.0- 0.0)      0.0147 (0.0139-0.0151)   0.608   0.796  0.880         -            -
 0.020 Hoppe       18.7 ( 0.0-56.1)      0.0218 (0.0165-0.0318)   0.482   0.708  0.800     0/1/0        0/0/0
 0.020 proximity    0.0 ( 0.0- 0.0)      0.0165 (0.0159-0.0171)   0.528   0.769  0.850     0/1/0        4/4/6
 0.020 oracle       0.0 ( 0.0- 0.0)      0.0165 (0.0159-0.0171)   0.528   0.769  0.850         -            -

n = 3000: median distance to 15th neighbour, x < -0.4: 0.054, x > 0.4: 0.120  (plate thickness 0.24)
 sigma rule       orient. err % (range)          CD-L1 (range)  F@0.01  F@0.02     NC    across    ambiguous
 0.000 Hoppe        0.0 ( 0.0- 0.0)      0.0068 (0.0068-0.0069)   0.844   0.922  0.967     0/0/0        0/0/0
 0.000 proximity    0.0 ( 0.0- 0.0)      0.0068 (0.0068-0.0069)   0.844   0.922  0.967     0/0/0        0/0/0
 0.000 oracle       0.0 ( 0.0- 0.0)      0.0068 (0.0068-0.0069)   0.844   0.922  0.967         -            -
 0.005 Hoppe        0.0 ( 0.0- 0.0)      0.0070 (0.0069-0.0071)   0.842   0.921  0.962     0/0/0        0/0/0
 0.005 proximity    0.0 ( 0.0- 0.0)      0.0070 (0.0069-0.0071)   0.842   0.921  0.962     0/0/0        0/0/0
 0.005 oracle       0.0 ( 0.0- 0.0)      0.0070 (0.0069-0.0071)   0.842   0.921  0.962         -            -
 0.010 Hoppe        0.0 ( 0.0- 0.0)      0.0077 (0.0075-0.0078)   0.825   0.921  0.945     0/0/0        0/0/0
 0.010 proximity    0.0 ( 0.0- 0.0)      0.0077 (0.0075-0.0078)   0.825   0.921  0.945     0/0/0        0/2/2
 0.010 oracle       0.0 ( 0.0- 0.0)      0.0077 (0.0075-0.0078)   0.825   0.921  0.945         -            -
 0.020 Hoppe        0.5 ( 0.5- 0.6)      0.0120 (0.0118-0.0125)   0.646   0.874  0.810     0/0/0        4/0/0
 0.020 proximity   29.7 (19.1-47.4)      0.0244 (0.0211-0.0280)   0.485   0.682  0.586     0/0/0   87/103/108
 0.020 oracle       0.0 ( 0.0- 0.0)      0.0116 (0.0113-0.0119)   0.650   0.879  0.818         -            -

n = 10000: median distance to 15th neighbour, x < -0.4: 0.030, x > 0.4: 0.065  (plate thickness 0.24)
 sigma rule       orient. err % (range)          CD-L1 (range)  F@0.01  F@0.02     NC    across    ambiguous
 0.000 Hoppe        0.0 ( 0.0- 0.0)      0.0039 (0.0039-0.0040)   0.954   0.996  0.991     0/0/0        0/0/0
 0.000 proximity    0.0 ( 0.0- 0.0)      0.0039 (0.0039-0.0040)   0.954   0.996  0.991     0/0/0        0/0/0
 0.000 oracle       0.0 ( 0.0- 0.0)      0.0039 (0.0039-0.0040)   0.954   0.996  0.991         -            -
 0.005 Hoppe        0.0 ( 0.0- 0.0)      0.0042 (0.0042-0.0042)   0.950   0.995  0.982     0/0/0        0/0/0
 0.005 proximity    0.0 ( 0.0- 0.0)      0.0042 (0.0042-0.0042)   0.950   0.995  0.982     0/0/0        1/0/0
 0.005 oracle       0.0 ( 0.0- 0.0)      0.0042 (0.0042-0.0042)   0.950   0.995  0.982         -            -
 0.010 Hoppe        0.2 ( 0.1- 0.2)      0.0061 (0.0054-0.0066)   0.919   0.983  0.915     0/0/0        5/1/2
 0.010 proximity   43.4 (28.3-51.9)      0.0235 (0.0210-0.0248)   0.615   0.696  0.563     0/0/0  238/196/197
 0.010 oracle       0.0 ( 0.0- 0.0)      0.0058 (0.0053-0.0061)   0.922   0.986  0.920         -            -
 0.020 Hoppe        3.7 ( 3.3- 4.1)      0.0157 (0.0153-0.0162)   0.570   0.831  0.605     0/0/0        5/7/4
 0.020 proximity   49.9 (47.3-54.8)      0.0297 (0.0292-0.0300)   0.334   0.546  0.409     0/0/0  1574/1590/1590
 0.020 oracle       0.0 ( 0.0- 0.0)      0.0131 (0.0126-0.0137)   0.605   0.871  0.638         -            -

Clean, dense data. At $n = 10{,}000$ without noise the three rules orient every plane alike and give the same mesh. Its Chamfer distance, 0.0039, sits above the floor of 0.0029 that two exact samplings of the surface reach against each other, and its F-score at $\tau = 0.01$ is 0.954. The gap is the method’s own error; in our reading it comes from the piecewise-planar field, which jumps between Voronoi cells of the centres, and from centroids that lie inside the curved rim.

Density alone. With correct signs (the oracle), the Chamfer distance grows by a factor of about 1.7 from $n = 10{,}000$ to 3,000 and about 2 from 3,000 to 1,000 (0.0039, 0.0068, 0.0133), and the normal consistency drops from 0.991 to 0.900. At $n = 1{,}000$ the 15 nearest neighbours at the sparse end reach 0.198, close to the plate’s thickness of 0.24, so in our reading the fitted planes there average over the rim and the opposite sheet.

Each propagation rule fails in its own way. On the sparse cloud, Hoppe’s cost failed on one seed of three at every noise level (43.5% of the normals flipped without noise), and the across column shows one tree edge between planes whose true normals point apart, on seed 1, at every noise level. The sample positions are the same at every $\sigma$, so in our reading the same layout fails each time. In our reading, the graph’s neighbourhoods at the sparse end reach the opposite sheet; a plane on the top and a plane on the bottom are nearly parallel, so the cost $1 - |\hat{\mathbf{n}}_i \cdot \hat{\mathbf{n}}_j|$ is near zero and the tree takes the edge. Copying the sign across it gives the bottom sheet the top’s upward normal, and the error spreads to the whole subtree behind it. This is the thin-plate case that Amenta et al. single out. The proximity rule prefers the shortest edges, which stay on one sheet, and was correct there. On dense, noisy data the roles reverse. At $n = 10{,}000$ and $\sigma = 0.02$ the proximity rule flips about half the normals (49.9%) with no edge across the plate, but with about 1,600 ambiguous edges per seed: in our reading the nearest neighbours are now about as close as the noise, their fitted normals scatter, and a sign copied between planes far from parallel is close to a guess. Hoppe’s cost avoids those edges (4 to 7 per seed) and keeps the error at 3.7%, which is the reason the authors give for their ordering: propagate between nearly parallel planes. A single across edge need not cause an error (proximity, $n = 1{,}000$, $\sigma \ge 0.01$, seed 1); in our reading such an edge can lie on the rounded rim, where the two true normals are just over 90° apart while the fitted normals are less than 90° apart, so the copied sign is still right.

Orientation errors show up in every measure. Where the proximity rule flips half the normals, the Chamfer distance more than doubles against the oracle’s (0.0297 against 0.0131), the F-score at 0.01 falls from 0.605 to 0.334 and the normal consistency from 0.638 to 0.409. Where two oppositely oriented planes meet, $f$ changes sign across their Voronoi boundary, so the zero set can gain sheets that run away from the surface; since $f$ is evaluated on the whole grid box, such sheets can extend to its walls. Berger et al. (2016, §2.2.1) make a related distinction, illustrated with Poisson reconstruction: scattered flips act like noise, while flips clustered over large regions are hard to handle.

Noise and the neighbourhood size. Even with correct signs, noise raises the error, and more steeply at the strict threshold: at $n = 10{,}000$ the oracle’s F-score at 0.01 falls from 0.954 to 0.605 between $\sigma = 0$ and 0.02, while at 0.02 it falls from 0.996 to 0.871. The nearest-plane rule fits each plane to only its $k$ neighbours and averages no further. At $\sigma = 0.02$ the oracle does better at $n = 3{,}000$ than at 10,000 (Chamfer 0.0116 against 0.0131, NC 0.818 against 0.638): in our reading, with $k$ fixed, denser sampling shrinks the neighbourhood, and a plane fitted to 15 noisy points over a smaller patch is less reliable. Berger et al. (2016, §2.2.1) cite a bound on the error of such fitted normals with this trade-off: small neighbourhoods suit clean data, larger ones noisy data. Hoppe et al. (§5.1) proposed adapting $k$ for the same reason.

The two tree-edge columns were added after the first run of the script, prompted by its results, to show where each rule goes wrong; nothing else in the set-up or the printed values changed.

Datasets and Benchmarks

Ground truth for surface reconstruction is hard to obtain from physical objects, so evaluations commonly use digital shapes as ground truth and sample or virtually scan them in place of a real scanner (Berger et al., 2016, §11), as the example does. Hoppe et al. (1992, §4) sampled points from existing meshes and ray-traced constructive solid geometry models from eight viewpoints to simulate multi-view range data. Kazhdan and Hoppe (2013, §6.1) used the benchmark of Berger and colleagues (cited there as 2011; the survey lists it as a 2013 paper in ACM Transactions on Graphics), virtual scans of five implicit surfaces (the Anchor, Dancing Children, Daratech, Gargoyle and Quasimodo models) that simulate nonuniform sampling, noise and misalignment.

For learned methods, ShapeNet supplies object meshes, from which point clouds, depth views and inside or outside labels are generated; Mescheder et al. (2019) evaluated against watertight versions of the models. Scene-scale benchmarks score whole pipelines from images, with laser or structured-light scans as reference: DTU for objects, Tanks and Temples for large scenes from video, and ETH3D; their protocols, including the alignment, cropping and observability masks they apply before scoring, are covered in their articles.

Evaluation Metrics

Reconstructed surfaces are compared with reference geometry as point sets sampled from both. Berger et al. (2016, §11) name the Hausdorff distance, the mean distance and the error in normals as the common geometric measures, and add topological accuracy, such as the correct genus, and the recovery of structure.

• Chamfer distance. The mean nearest-neighbour distance from reconstruction to reference (accuracy) and back (completeness), combined. Occupancy Networks average the two means (Chamfer-$L_1$, Mescheder et al., 2019, §4), measured in a unit of a tenth of the longest edge of each object’s bounding box. Definitions differ between papers; the Chamfer distance article covers the conventions and their floor at a finite number of samples, which the example’s first line shows.
• F-score. Tanks and Temples defines precision $P(d)$ as the percentage of reconstructed points within $d$ of the reference, recall $R(d)$ as the percentage of reference points within $d$ of the reconstruction, and $F(d) = 2P(d)R(d)/(P(d) + R(d))$ (Knapitsch et al., 2017, Eq. 4–7 of the authors’ preprint). Because $F$ goes to zero when either term does, a high score at a strict threshold needs a surface that is both accurate and complete (F1 score, precision and recall).
• Normal consistency. The mean absolute dot product between the normal at each point of one mesh and the normal at its nearest neighbour in the other, which Mescheder et al. (2019, §4) defined to measure how well a method captures higher-order information. By its definition the absolute value ignores orientation, so a flipped but well-placed patch still scores well; in the example NC fell with orientation errors, in our reading because the flips create misplaced sheets.
• Volumetric IoU. For closed surfaces, intersection over union of the enclosed volumes, which Mescheder et al. estimate from 100,000 random points in the bounding volume. It needs watertight meshes.
• Distances to held-out points. Without a reference surface for real scans, Kazhdan and Hoppe (2013, §6.1) split each scan’s points into two equal halves, reconstructed from one and reported the RMS distance from the other half to the surface; they found the maximum error too sensitive to individual outliers to be informative.

The measures disagree in useful ways. The Chamfer distance grows with how far wrong parts lie; the F-score counts only whether points lie within $\tau$; normal consistency responds to the direction of the surface, ignoring its orientation, and so to roughness on a surface that lies close to the reference. In the example, noise lowers the F-score at $\tau = 0.01$ far more than at 0.02, so a threshold must be stated with every score.

Applications

• 3D scanning and digitization. Merging range scans of statues and artefacts was the original driver: ball pivoting built the Pietà model from hundreds of scans (Bernardini et al., 1999), and volumetric fusion integrated as many as 70 range images into models of up to 2.6 million triangles (Curless and Levoy, abstract).
• Image-based modelling. Meshing is the last geometric step after structure from motion and multi-view stereo in photogrammetry pipelines (3D reconstruction).
• Medical and scientific data. Hoppe et al. (1992, §4) reconstructed a femur from points on 39 planar contours of a CT scan, a case that methods built for contour stacks handle separately.
• Manufacturing and design. Rapid prototyping needs watertight models, one motivation for hole filling (Curless and Levoy, 1996, §4); Hoppe et al. (1992, §1) also list surfaces sketched by moving a tracked stylus.
• Mapping with depth cameras. Real-time fusion of depth frames, as in KinectFusion and dense SLAM, produces meshes of rooms and objects.

Open Problems

• Orientation. Signed and indicator methods depend on consistent normals, and propagation can fail in clusters, as the example shows; Berger et al. (2016, §2.2.1) describe methods that remain sensitive to imperfect data. Their survey mentions a later reweighting that favours graph edges lying in the tangent plane, so that propagation does not jump between opposite sides.
• Missing data. Watertight methods invent surface over holes and open methods leave them; which is right depends on whether the gap is real, the ambiguity Hoppe et al. stated in 1992. Visibility, structural and data-driven priors narrow it for particular kinds of shapes and scenes (Berger et al., 2016, §2.1, §3).
• Sharp features and thin parts. Both demand sampling far denser than smooth regions (Amenta et al., 1998, §3.2), and smoothness priors round edges off.
• Scale and generality of learned priors. Learned methods complete shapes from little data but, as DeepSDF’s completion experiments show, may assume a canonical frame and a known class; neural implicit representations discusses the limits of a single global code.
• Comparable evaluation. Scores depend on the metric’s definition, the number of sampled points, the threshold and the alignment, so results are comparable only under the same protocol (Chamfer distance).