Poisson Surface Reconstruction

Poisson surface reconstruction turns points with normals into a closed triangle mesh. Its idea, due to Kazhdan, Bolitho and Hoppe (2006), is that the normals of a solid’s surface sample the gradient of the solid’s indicator function, 1 inside and 0 outside. Finding the function whose gradient best matches the normals is a Poisson equation, and the surface is a level set of the solution. It is a common last step of 3D reconstruction pipelines, after laser scanning or multi-view stereo, and PoissonRecon and Open3D implement its screened variant (Kazhdan and Hoppe, 2013).

Problem

A point cloud records where surfaces were measured, not the surface between the samples. Scans are noisy, unevenly sampled, incomplete and, when merged, misregistered (Kazhdan et al., 2006, §2); the surface should approximate the points without fitting their noise and fill gaps plausibly.

Inputs and Outputs

• Input: points $\mathbf{p}_i$ with consistently oriented normals $\mathbf{n}_i$. Normals from local plane fits have arbitrary signs and need an orientation step first (point cloud).
• Output: a watertight triangle mesh and, in implementations, a per-vertex value measuring how well the samples supported each part of it.

Intuition

Blur the indicator function slightly. Far from the surface it is constant, so its gradient is zero; across the surface it rises from 0 to 1, so its gradient points along the inward normal. Oriented points are therefore samples of this gradient field, and integrating the field recovers the indicator up to a constant.

The fit is global, and the gradient is constrained everywhere: to zero far from the data. Kazhdan et al. (2006, §1) give this as the reason spurious sheets seldom appear away from the samples, unlike in methods that constrain an implicit function only near them.

Algorithm

As described by Kazhdan et al. (2006, §4):

1. Build an octree: the smallest one in which every sample falls in a leaf at the maximum depth $D$.
2. Form the vector field $\vec V$: spread each normal over the eight nearest depth-$D$ nodes with trilinear weights (Eq. 3); for non-uniform sampling, weight samples by the area they represent and splat sparse ones at a coarser depth (§4.5).
3. Solve $\Delta\chi = \nabla\cdot\vec V$, a sparse symmetric system, depth by depth in a multigrid scheme with conjugate gradients at each depth (§4.3).
4. Choose the iso-value $\gamma$, the average of $\chi$ at the samples (§4.4).
5. Extract $\chi = \gamma$ with marching cubes adapted to the octree (§4.4).

Mathematical Formulation

Let $M$ be the solid, $\chi_M$ its indicator, $\tilde F$ a smoothing filter and $\mathbf{N}_{\partial M}$ the inward normal. By the divergence theorem, the gradient of the smoothed indicator is the smoothed normal field (Kazhdan et al., 2006, Eq. 1):

$$\nabla\left(\chi_M * \tilde F\right)(\mathbf{q}) = \int_{\partial M} \tilde F(\mathbf{q} - \mathbf{p})\, \mathbf{N}_{\partial M}(\mathbf{p})\, d\mathbf{p}.$$

The surface is unknown, so the integral becomes a sum over samples $s$, each standing for a surface patch $\mathcal{P}_s$ (Eq. 2):

$$\vec V(\mathbf{q}) = \sum_{s} |\mathcal{P}_s|\, \tilde F(\mathbf{q} - s.\mathbf{p})\, s.\mathbf{N}.$$

The filter approximates a Gaussian whose variance is on the order of the sampling resolution (§3). $\vec V$ is generally not curl-free, so no $\chi$ has exactly this gradient; the least-squares fit

$$\min_\chi \int \left\lVert \nabla\chi(\mathbf{p}) - \vec V(\mathbf{p}) \right\rVert^2 d\mathbf{p}$$

has the Euler–Lagrange equation $\Delta\chi = \nabla\cdot\vec V$ (Kazhdan and Hoppe, 2013, Eq. 1 and §3).

Discretization. Write $\chi = \sum_i x_i B_i$ over functions attached to octree nodes and scaled to their width. In the 2006 paper the base function is the convolution of three box filters, a piecewise quadratic B-spline supported on $[-1.5, 1.5]^3$ cell widths, so each function overlaps at most 124 others at its depth (§4.1). The Galerkin method gives a sparse system $A\mathbf{x} = \mathbf{b}$ (Kazhdan and Hoppe, 2013, Eq. 2):

$$A_{ij} = \langle \nabla B_i, \nabla B_j \rangle, \qquad b_i = \langle \vec V, \nabla B_i \rangle.$$

Iso-value. $\chi$ is fixed only up to a constant and is never exactly 0 and 1, so the 2006 paper takes

$$\gamma = \frac{1}{|S|} \sum_{s \in S} \chi(s.\mathbf{p}).$$

Scaling $\chi$ does not move this level set, so $\vec V$ is needed only up to a constant factor (§4.4).

Parameters

Names and defaults follow the PoissonRecon README (version 18.76, as of October 2026).

• Depth $D$ (--depth, default 8): the maximum octree depth, at most a $2^D$-per-side grid. In the 2006 paper, depths 6, 8 and 10 (64³ to 1024³ grids) capture progressively finer detail (§5.1).
• Samples per node (--samplesPerNode, default 1.5): the minimum number of samples per node, which stops refinement where data is sparse; 1 to 5 for noise-free and 15 to 20 for noisy data, per the README. Raising it reduces noise but also detail (Kazhdan and Hoppe, 2013, §6.3).
• Point weight $\alpha$ (--pointWeight, default twice the B-spline degree, which defaults to 1): the screening weight (see Variants); 0 gives the original method.
• Boundary type (--bType, default Neumann; see Variants).

Open3D’s TriangleMesh.create_from_point_cloud_poisson wraps Kazhdan’s code and returns the mesh with a per-vertex density. Version 0.20.0, the latest release as of October 2026, adds samples_per_node (default 1.5), point_weight (2.0) and full_depth to the depth, width, scale, linear_fit and n_threads arguments of 0.19.0.

Density trimming. Open3D’s tutorial (unchanged in 0.20.0) notes that the method extrapolates surface into areas of low point density, and its example removes vertices below the 0.01 density quantile. PoissonRecon’s --density writes each vertex’s estimated sampling depth, and its SurfaceTrimmer discards the parts below a --trim value. The DTU benchmark meshed the points of the three multi-view stereo methods it evaluated at depth 11 and trimmed areas with depth below 8 (Aanæs et al., 2016, §5 of the authors’ manuscript).

Complexity

A regular grid of resolution $r$ needs $O(r^3)$ memory and, with the FFT, $O(r^3 \log r)$ time, while the surface occupies only $O(r^2)$ cells; hence the octree (Kazhdan et al., 2006, §2, §4.1). Table 1 of the 2006 paper reports 6, 26, 126 and 633 s and 19, 75, 155 and 699 MB at depths 7 to 10 on the dragon. The text calls the growth about fourfold per depth; the table’s memory ratios are 3.9, 2.1 and 4.5. Kazhdan and Hoppe (2013, §5) make the solver linear in the number of points; in their Table I (David, depth 10), the screened method took 182 s against 412 s for the original implementation.

Implementation

The example reconstructs a torus from 20,000 noisy oriented samples on a dense $128^3$ grid rather than an octree. It splats inward normals with trilinear weights, smooths them and takes the divergence in the Fourier domain, and solves $-|\mathbf{k}|^2 \hat\chi = \widehat{\nabla\cdot\vec V}$, assuming periodic boundaries in a padded box. Errors are exact distances to the torus from the sign changes of $\chi - \gamma$ along grid edges, where marching cubes would place vertices. Three corrupted clouds follow: a patch removed from the top of the tube, 20% of normals flipped at random, and the patch’s normals flipped. A vertical line through the patch’s centre shows where the surface lies.

import numpy as np
from scipy.ndimage import map_coordinates

R, r = 0.6, 0.25          # torus: ring radius R, tube radius r, axis along z
N, L = 128, 1.25          # grid of N^3 cells covering the box [-L, L]^3
h = 2 * L / N             # cell size

def sample_torus(n, rng):
    """n area-uniform points on the torus, with outward unit normals."""
    u, v = rng.uniform(0, 2 * np.pi, (2, 4 * n))
    keep = rng.uniform(0, R + r, 4 * n) < R + r * np.cos(v)   # area element is proportional to R + r cos v
    u, v = u[keep][:n], v[keep][:n]
    nrm = np.stack([np.cos(v) * np.cos(u), np.cos(v) * np.sin(u), np.sin(v)], axis=1)
    pts = np.stack([R * np.cos(u), R * np.sin(u), np.zeros(n)], axis=1) + r * nrm
    return pts, nrm

def dist_to_torus(x):
    """Exact distance from points x to the torus surface."""
    return np.abs(np.hypot(np.hypot(x[:, 0], x[:, 1]) - R, x[:, 2]) - r)

def to_grid(x):
    return ((x + L) / h - 0.5).T        # world -> continuous cell index (cell centres at integers)

def reconstruct(pts, nrm, area, sigma=1.0):
    """Indicator chi (1 inside, 0 outside) and the iso-value gamma."""
    # 1. Splat each sample's inward normal (the gradient of chi) into the grid with
    #    trilinear weights, scaled by the surface area the sample stands for.
    g = to_grid(pts).T
    i0 = np.floor(g).astype(int)
    f = g - i0
    V = np.zeros((3, N, N, N))
    for corner in np.ndindex(2, 2, 2):
        w = np.prod(np.where(corner, f, 1 - f), axis=1)
        idx = tuple((i0 + corner).T)
        for c in range(3):
            np.add.at(V[c], idx, -area * w * nrm[:, c] / h**3)
    # 2. Smooth with a Gaussian (sigma cells) and take the divergence, in the Fourier domain.
    k = 2 * np.pi * np.fft.fftfreq(N, d=h)
    K = np.stack(np.meshgrid(k, k, k, indexing="ij"))
    K2 = (K**2).sum(axis=0)
    div_hat = 1j * sum(K[c] * np.fft.fftn(V[c]) for c in range(3))
    div_hat *= np.exp(-0.5 * (sigma * h) ** 2 * K2)
    # 3. Solve the Poisson equation lap(chi) = div(V) with periodic boundaries:
    #    -|k|^2 chi_hat = div_hat. The constant (k = 0) term is free; set it to 0.
    K2[0, 0, 0] = 1.0
    chi_hat = -div_hat / K2
    chi_hat[0, 0, 0] = 0.0
    chi = np.fft.ifftn(chi_hat).real
    # 4. Iso-value: the average of chi at the input samples.
    gamma = map_coordinates(chi, to_grid(pts), order=1).mean()
    return chi, gamma

def surface_points(chi, gamma):
    """Where chi - gamma changes sign along grid edges (the vertices marching cubes would place)."""
    s, out = chi - gamma, []
    for ax in range(3):
        a = np.moveaxis(s, ax, 0)
        m = (a[:-1] > 0) != (a[1:] > 0)
        idx = np.argwhere(m).astype(float)
        idx[:, 0] += a[:-1][m] / (a[:-1][m] - a[1:][m])   # linear interpolation along the edge
        xyz = np.empty_like(idx)
        xyz[:, [ax] + [j for j in range(3) if j != ax]] = idx
        out.append(xyz)
    return np.vstack(out) * h - L + 0.5 * h

def crossings_on_line(chi, gamma, x, y):
    """z values where the vertical line through (x, y) crosses the surface."""
    z = np.linspace(-L + h / 2, L - h / 2, 4 * N)
    vals = map_coordinates(chi, to_grid(np.stack([np.full_like(z, x), np.full_like(z, y), z], 1)), order=1) - gamma
    i = np.where(np.diff(np.sign(vals)) != 0)[0]
    return z[i] - vals[i] * (z[i + 1] - z[i]) / (vals[i + 1] - vals[i])

def make_cloud(seed, case, n=20000, noise=0.005):
    rng = np.random.default_rng(seed)
    pts, nrm = sample_torus(n, rng)
    pts = pts + rng.normal(scale=noise, size=pts.shape)
    area = 4 * np.pi**2 * R * r / n                       # torus area / number of samples
    patch = (np.abs(np.arctan2(pts[:, 1], pts[:, 0])) < 0.5) & (pts[:, 2] > 0.05)  # top of the tube near +x
    if case == "missing patch":
        pts, nrm = pts[~patch], nrm[~patch]
    elif case == "20% flipped at random":
        nrm = np.where(rng.uniform(size=(n, 1)) < 0.2, -nrm, nrm)
    elif case == "patch flipped":
        nrm = np.where(patch[:, None], -nrm, nrm)
    return pts, nrm, area

cases = ["clean", "missing patch", "20% flipped at random", "patch flipped"]
print(f"grid {N}^3, cell size h = {h:.4f}; torus R = {R}, r = {r}; true top of tube at z = {r}")
print(f"seed 4\n{'case':24s} {'gamma':>6s} {'mean dist':>10s} {'max dist':>9s}   line x=R, y=0 crosses at z =")
for case in cases:
    pts, nrm, area = make_cloud(4, case)          # seed 4; the loop below gives the range over seeds 0-9
    chi, gamma = reconstruct(pts, nrm, area)
    d = dist_to_torus(surface_points(chi, gamma))
    zc = crossings_on_line(chi, gamma, R, 0.0)
    print(f"{case:24s} {gamma:6.3f} {d.mean():10.4f} {d.max():9.4f}   {np.array2string(zc, precision=3)}")
    if case == "clean":
        inside, outside = map_coordinates(chi, to_grid(np.array([[R, 0, 0], [1.1, 1.1, 1.1]])), order=1)
        print(f"{'':24s} chi = {inside:.3f} at the tube's centre, {outside:.3f} near the box corner")

print("\nseeds 0-9: range of mean / max distance to the true torus")
for case in cases:
    stats = []
    for seed in range(10):
        pts, nrm, area = make_cloud(seed, case)
        chi, gamma = reconstruct(pts, nrm, area)
        d = dist_to_torus(surface_points(chi, gamma))
        stats.append((d.mean(), d.max()))
    s = np.array(stats)
    print(f"{case:24s} mean {s[:, 0].min():.4f}-{s[:, 0].max():.4f}   max {s[:, 1].min():.4f}-{s[:, 1].max():.4f}")

Output (NumPy 2.3.5, SciPy 1.16.3):

grid 128^3, cell size h = 0.0195; torus R = 0.6, r = 0.25; true top of tube at z = 0.25
seed 4
case                      gamma  mean dist  max dist   line x=R, y=0 crosses at z =
clean                     0.434     0.0006    0.0028   [-0.25   0.251]
                         chi = 0.947 at the tube's centre, -0.047 near the box corner
missing patch             0.412     0.0043    0.1095   [-0.242  0.142]
20% flipped at random     0.263     0.0009    0.0055   [-0.249  0.25 ]
patch flipped             0.373     0.0153    0.2447   [-0.235 -0.035  0.273  0.347]

seeds 0-9: range of mean / max distance to the true torus
clean                    mean 0.0005-0.0008   max 0.0024-0.0041
missing patch            mean 0.0041-0.0049   max 0.1041-0.1168
20% flipped at random    mean 0.0008-0.0011   max 0.0044-0.0083
patch flipped            mean 0.0151-0.0158   max 0.2447-0.2496

On clean data, $\chi$ jumps by almost exactly 1 across the surface, offset by a constant (the Fourier solve makes its mean zero) that the iso-value absorbs. The mean error, 0.0006, is far below the per-point noise of 0.005.

Without the patch, the surface stays closed, but the line meets its top at $z = 0.142$ instead of 0.25: an invented lid sagging five to six cells. Flipping 20% of the normals at random scales the expected field by $1 - 2(0.2) = 0.6$; $\gamma$ drops by about that factor and the surface barely moves, as the scale invariance of the iso-value predicts. Flipping a coherent patch does damage: the line crosses four times ($-0.235$, $-0.035$, 0.273, 0.347) instead of twice, with errors up to the tube radius. In our reading, coherent orientation errors, not scattered ones, are the dangerous kind.

Properties and Behavior

• Watertight. The level set separates the region where $\chi > \gamma$ from the rest, so it bounds a volume unless it reaches the domain boundary, which Dirichlet conditions prevent; no hole filling is needed (our reasoning; the 2006 paper sets a watertight surface as its goal, §3).
• Equivalent to the FFT method on uniform data. The basic formulation is mathematically equivalent to Kazhdan’s earlier Fourier-based method (2006, Appendix A); the octree adds adaptivity and a smoother fit where sampling is sparse (§5.2).

Limitations

• Normals must be consistently oriented. As the example shows, a wrongly oriented region creates spurious sheets.
• It always closes the surface. Unobserved regions receive invented surface. Without line-of-sight information, the 2006 reconstruction joined the two feet of the Happy Buddha statue, where no samples lay between them (§5.2). Open scenes, such as a room or a one-sided scan, need trimming.
• Over-smoothing of detail, noted by several researchers (Kazhdan and Hoppe, 2013, §1), which motivated screening.
• Non-uniform density. Sparse regions get a coarser, smoother fit (2006, §4.5), losing detail there.

Variants

Screened Poisson reconstruction. Kazhdan and Hoppe (2013) define $\chi$ as $\tfrac12$ inside and $-\tfrac12$ outside, so the surface is the zero level set, and add a term that pulls $\chi$ toward zero at the samples (Eq. 3):

$$E(\chi) = \int \left\lVert \vec V(\mathbf{p}) - \nabla\chi(\mathbf{p}) \right\rVert^2 d\mathbf{p} + \alpha\, \frac{\mathrm{Area}(P)}{\sum_{\mathbf{p} \in P} w(\mathbf{p})} \sum_{\mathbf{p} \in P} w(\mathbf{p})\, \chi^2(\mathbf{p}),$$

with unit weights $w$ in their experiments. The minimizer satisfies a screened Poisson equation $(\Delta - \alpha \tilde I)\chi = \nabla\cdot\vec V$ (Eq. 6), where $\tilde I$ applies the screening only at the samples (Appendix A). The matrix gains a term $\alpha \langle B_i, B_j \rangle$ summed over the points (Eq. 5) with the same sparsity, so the multigrid solver still applies, and the weight doubles with each depth to keep the terms balanced (§4.3). With $\alpha = 4$, on scans of five models with half the points held out, the authors report lower RMS distances to the held-out points than the original method in every case (Fig. 4b), and similar results on clean synthetic samples (Fig. 5b). On the Berger et al. benchmark, Calakli and Taubin’s smooth signed distance method was more accurate, which they attribute to better extrapolation into missing regions (§6.1). Screening also fits misaligned scans and noise too closely (§6.3).

Boundary conditions. The original solver fixes $\chi = -\tfrac12$ on the domain boundary (Dirichlet); the 2013 version adds Neumann conditions, a zero normal derivative there, which the authors find less restrictive when data is missing. On an angel scanned only from the front, the Dirichlet surface closes off before the boundary, while the Neumann surface extends out to it (§4.4, Fig. 3).

Alternatives

TSDF fusion (Curless and Levoy, 1996) fuses signed distances from range images taken from known sensor positions, so it can use lines of sight. Kazhdan et al. (2006, §5.2) found that it separated the Buddha’s feet but showed “lipping” at sharp creases, because its distances follow the viewing direction rather than the surface normal. Poisson reconstruction needs only oriented points.

Combinatorial methods build the mesh from the points themselves. Delaunay, alpha-shape and Voronoi methods typically interpolate all or most of the points, so noisy data gives jagged surfaces that need smoothing or refitting (Kazhdan et al., 2006, §2). Ball pivoting (Bernardini et al., 1999) rolls a ball of chosen radius over the points, forming a triangle wherever it touches three points without containing another, with larger balls for sparse regions; missing points leave holes the ball cannot fill (§3). Leaving holes rather than inventing surface may suit a measurement application (our reading).