reconstruction

triwarp.reconstruction

Building a triangle mesh whose triangulation is the method's own, not any input's.

• From an oriented point cloud. screened_poisson fits an implicit indicator function and contours it, so it closes gaps and returns a watertight surface that need not pass through any input point; ball_pivoting instead interpolates the points, rolling a ball over them, and leaves a hole wherever the ball falls through. triangulate_point_cloud picks between them.
• From a mesh. resample_uniform samples an existing mesh's signed distance field on a uniform grid and re-extracts it, which is the way to get a clean, uniformly-sampled surface out of a self-intersecting or badly-triangulated one. It is shelved with the reconstructors rather than with the mesh edits because it shares their defining property: the output's triangulation comes from the grid, not from the input, so no vertex, edge or face survives the call.
• In the plane. delaunay_triangulation triangulates a 2D point set.

Extracting a level set of a field you already hold is levelset.marching_cubes, not a reconstruction: it lives with its consumers in triwarp.levelset, and screened_poisson and resample_uniform both end in it.

delaunay_triangulation

delaunay_triangulation(
    points: array[vec2], max_iter: int = 1000
) -> wp.array[wp.int32]

Delaunay triangulation of a 2D point set.

Ports igl::delaunay_triangulation: a sequential lexicographic incremental triangulation seeds the mesh, then interior edges are flipped until every one satisfies the empty-circumcircle (Delaunay) criterion. The flip loop reuses the parallel edge-flip machinery of flip_to_delaunay with a float64 in-circle predicate; ties (cocircular points) are left unflipped, which guarantees termination.

Parameters:

Name	Type	Description	Default
points	array[vec2]	(n,) 2D point positions on the target device.	required
max_iter	int	Maximum number of parallel flip passes.	1000

Returns:

Type	Description
array[int32]	Flat length-3 * n_faces buffer of CCW triangles over the input point indices, on points.device. Empty when the points are collinear.

Raises:

Type	Description
ValueError	If fewer than 3 points are given, or if a degenerate input (duplicate or collinear points) drives the seed past the 2 * n triangle bound a triangulation obeys.

See Also

flip_to_delaunay triangulate_point_cloud

Notes

The seed triangulation is an inherently sequential hull sweep, so it runs single-threaded on the CPU device and only the flips run on points.device -- a sequential sweep is genuinely faster on the CPU than as a single CUDA thread. A float64 in-circle determinant is used rather than exact predicates, so near-cocircular inputs may resolve either ambiguous diagonal.

Source code in triwarp/reconstruction.py

def delaunay_triangulation(points: wp.array[wp.vec2], max_iter: int = 1000) -> wp.array[wp.int32]:
    """
    Delaunay triangulation of a 2D point set.

    Ports ``igl::delaunay_triangulation``: a sequential lexicographic incremental triangulation
    seeds the mesh, then interior edges are flipped until every one satisfies the empty-circumcircle
    (Delaunay) criterion. The flip loop reuses the parallel edge-flip machinery of
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay] with a float64 in-circle predicate; ties
    (cocircular points) are left unflipped, which guarantees termination.

    Parameters
    ----------
    points
        ``(n,)`` 2D point positions on the target device.
    max_iter
        Maximum number of parallel flip passes.

    Returns
    -------
    wp.array[wp.int32]
        Flat length-``3 * n_faces`` buffer of CCW triangles over the input point indices, on
        ``points.device``. Empty when the points are collinear.

    Raises
    ------
    ValueError
        If fewer than 3 points are given, or if a degenerate input (duplicate or collinear points)
        drives the seed past the ``2 * n`` triangle bound a triangulation obeys.

    See Also
    --------
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]
    [`triangulate_point_cloud`][triwarp.reconstruction.triangulate_point_cloud]

    Notes
    -----
    The seed triangulation is an inherently sequential hull sweep, so it runs single-threaded on the
    **CPU** device and only the flips run on ``points.device`` -- a sequential sweep is genuinely
    faster on the CPU than as a single CUDA thread. A float64 in-circle determinant is used rather
    than exact predicates, so near-cocircular inputs may resolve either ambiguous diagonal.
    """
    device = points.device
    n = int(points.shape[0])
    if n < 3:
        raise ValueError(f"delaunay_triangulation requires at least 3 points, got {n}")

    faces_np = _lexicographic_triangulation(points)
    if faces_np.shape[0] == 0:
        return wp.empty(0, dtype=wp.int32, device=device)

    faces = wp.array(np.ascontiguousarray(faces_np.reshape(-1)), dtype=wp.int32, device=device)

    def launch(adjacency, adjacency_edges, unshared, sorted_keys, key_base, out_flip, out_quad):
        wp.launch(
            kernel_remesh.incircle_flip_candidates,
            dim=int(adjacency.shape[0]),
            inputs=[
                points,
                faces,
                adjacency,
                adjacency_edges,
                unshared,
                sorted_keys,
                key_base,
                out_flip,
                out_quad,
            ],
            device=device,
        )

    tw.remesh._flip_interior_edges(faces, n, launch, max_iter)
    return faces

triangulate_point_cloud

triangulate_point_cloud(
    points: array[vec3],
    normals: array[vec3] | None = None,
    *,
    num_neighbours: int = 0,
    radius: float = 0.0,
    crit_angle: float = math.pi / 2.0,
    boundary_angle: float = 0.9 * math.pi,
    max_neighbours: int = 32,
    crit_hole_length: float = -1.0,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Reconstruct a triangle mesh from an (optionally oriented) point cloud.

Local-triangulation reconstruction, on the GPU. Each point independently builds a local triangle fan over its nearest neighbours in the tangent plane, greedily optimised toward a Delaunay-like fan (build_local_triangulations); triangles that recur in two or three of these local triangulations are kept and assembled into a triangle-soup mesh (vertices are the input points). Non-manifold and degenerate faces are dropped and small boundary holes are filled.

Orientation relies on trusted normals rather than on a sequential propagation pass, which is what keeps the assembly parallel: supply normals for arbitrary geometry, or leave them None to estimate them by PCA oriented outward from the cloud centroid (valid for star-shaped clouds only — see estimate_normals).

Parameters:

Name	Type	Description	Default
points	array[vec3]	(n,) point positions on the target device.	required
normals	array[vec3] | None	Optional (n,) oriented unit normals. When None, normals are estimated.	None
num_neighbours	int	Number of nearest neighbours used per point when radius <= 0. Defaults to max_neighbours when neither num_neighbours nor radius is set.	0
radius	float	When positive, only neighbours within this distance are used (still capped at max_neighbours nearest). Mutually exclusive with num_neighbours.	0.0
crit_angle	float	Maximum dihedral angle (radians) tolerated in a fan triangle before an edge is flipped.	pi / 2.0
boundary_angle	float	A point is a fan boundary when its neighbour ring has an angular gap wider than this.	0.9 * pi
max_neighbours	int	Per-point neighbour cap. Must not exceed kernels.reconstruction.MAX_NEIGHBOURS (the compile-time scratch size).	32
crit_hole_length	float	Boundary loops with perimeter at most this value are filled. When negative, defaults to 0.1 * the point-cloud bounding-box diagonal.	-1.0

Returns:

Name	Type	Description
vertices	array[vec3]	The referenced input points, compacted from index zero (unreferenced points dropped), on points.device.
faces	array[int32]	Flat 3 * n_faces triangle index buffer into vertices.

Raises:

Type	Description
ValueError	If both num_neighbours and radius are positive, or max_neighbours exceeds the compile-time cap.
RuntimeError	If points and normals are not all on one device.

Notes

The per-point fan is bounded by max_neighbours and the search radius is never grown automatically, so very sparse or highly non-uniform clouds may leave extra boundary holes.

See Also

fill_min_weight triwarp.repair.remove_degenerate_faces

Source code in triwarp/reconstruction.py

def triangulate_point_cloud(
    points: wp.array[wp.vec3],
    normals: wp.array[wp.vec3] | None = None,
    *,
    num_neighbours: int = 0,
    radius: float = 0.0,
    crit_angle: float = math.pi / 2.0,
    boundary_angle: float = 0.9 * math.pi,
    max_neighbours: int = 32,
    crit_hole_length: float = -1.0,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Reconstruct a triangle mesh from an (optionally oriented) point cloud.

    Local-triangulation reconstruction, on the GPU. Each point independently builds a local
    triangle fan over its nearest neighbours in the tangent plane, greedily
    optimised toward a Delaunay-like fan (``build_local_triangulations``); triangles that recur in
    two or three of these local triangulations are kept and assembled into a triangle-soup mesh
    (vertices are the input points). Non-manifold and degenerate faces are dropped and small
    boundary holes are filled.

    Orientation relies on **trusted normals** rather than on a sequential propagation pass, which
    is what keeps the assembly parallel: supply ``normals`` for arbitrary geometry, or leave them
    ``None`` to estimate them by PCA oriented outward from the cloud centroid (valid for
    star-shaped clouds only — see [`estimate_normals`][triwarp.points.estimate_normals]).

    Parameters
    ----------
    points
        ``(n,)`` point positions on the target device.
    normals
        Optional ``(n,)`` oriented unit normals. When ``None``, normals are estimated.
    num_neighbours
        Number of nearest neighbours used per point when ``radius <= 0``. Defaults to
        ``max_neighbours`` when neither ``num_neighbours`` nor ``radius`` is set.
    radius
        When positive, only neighbours within this distance are used (still capped at
        ``max_neighbours`` nearest). Mutually exclusive with ``num_neighbours``.
    crit_angle
        Maximum dihedral angle (radians) tolerated in a fan triangle before an edge is flipped.
    boundary_angle
        A point is a fan boundary when its neighbour ring has an angular gap wider than this.
    max_neighbours
        Per-point neighbour cap. Must not exceed
        ``kernels.reconstruction.MAX_NEIGHBOURS`` (the compile-time scratch size).
    crit_hole_length
        Boundary loops with perimeter at most this value are filled. When negative, defaults to
        ``0.1 *`` the point-cloud bounding-box diagonal.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        The referenced input points, compacted from index zero (unreferenced points dropped),
        on ``points.device``.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer into ``vertices``.

    Raises
    ------
    ValueError
        If both ``num_neighbours`` and ``radius`` are positive, or ``max_neighbours`` exceeds the
        compile-time cap.
    RuntimeError
        If ``points`` and ``normals`` are not all on one device.

    Notes
    -----
    The per-point fan is bounded by ``max_neighbours`` and the search radius is never grown
    automatically, so very sparse or highly non-uniform clouds may leave extra boundary holes.

    See Also
    --------
    [`fill_min_weight`][triwarp.holes.fill_min_weight]
    [`triwarp.repair.remove_degenerate_faces`][]
    """
    require_same_device(points=points, normals=normals)
    if num_neighbours > 0 and radius > 0.0:
        raise ValueError("pass at most one of num_neighbours and radius")
    if max_neighbours > kernel_reconstruction.MAX_NEIGHBOURS:
        raise ValueError(
            "max_neighbours must be <= "
            f"{kernel_reconstruction.MAX_NEIGHBOURS}, got {max_neighbours}"
        )

    device = points.device
    n = int(points.shape[0])
    if n == 0:
        # Cloned, like every other exit from this module: the returned vertices are the function's
        # own buffer on the populated path (`_assemble_faces`, `_clean_reconstruction`), so handing
        # back the caller's array here would make the degenerate input the one case where mutating
        # the result mutates the input.
        return wp.clone(points), wp.empty(0, dtype=wp.int32, device=device)

    k = num_neighbours if num_neighbours > 0 else max_neighbours
    k = min(k, max_neighbours)

    # Dense (n, k+1) nearest-neighbour table; slot 0 is the point itself and is skipped in-kernel.
    neighbor_idx, neighbor_dist = tw.neighbors.query_nearest(points, points, k=k + 1, backend="bvh")

    if normals is None:
        normals = tw.points.estimate_normals(points, neighbor_idx)

    # Per-point local fan triangulation.
    out_tris = wp.empty((n, k, 3), dtype=wp.int32, device=device)
    out_valid = wp.zeros((n, k), dtype=wp.bool, device=device)
    wp.launch(
        kernel_reconstruction.build_local_triangulations,
        dim=n,
        inputs=[
            points,
            normals,
            neighbor_idx,
            neighbor_dist,
            wp.float32(radius),
            wp.float32(crit_angle),
            wp.float32(boundary_angle),
            out_tris,
            out_valid,
        ],
        device=device,
    )

    # Compact the emitted candidate triangles.
    kept = tw.array.flatnonzero(out_valid.reshape(-1))
    if int(kept.shape[0]) == 0:
        return wp.clone(points), wp.empty(0, dtype=wp.int32, device=device)
    candidates = twt.as_array2d(tw.array.gather(out_tris.reshape((n * k, 3)), kept), wp.int32)

    # Repeated oriented triangles: t3 (3 reps) preferred, then t2 (2 reps).
    t3 = _repeated_oriented_triangles(candidates, n, 3)
    t2 = _repeated_oriented_triangles(candidates, n, 2)

    return _assemble_faces(points, t3, t2, crit_hole_length)

screened_poisson

screened_poisson(
    points: array[vec3],
    normals: array[vec3],
    *,
    depth: int = 8,
    full_depth: int = 5,
    scale: float = 1.1,
    point_weight: float = 4.0,
    solver_iterations: int = 100,
    solver_tolerance: float = 1e-06,
    confidence: bool = False,
    method: Literal["dense", "adaptive"] = "dense",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Screened-Poisson surface reconstruction from an oriented point cloud (Kazhdan PoissonRecon).

GPU-native reimplementation of the screened-Poisson filter (PyMeshLab's generate_surface_reconstruction_screened_poisson, open3d's create_from_point_cloud_poisson). A dense node-centered grid over a padded bounding cube is used instead of PoissonRecon's octree: the oriented normals are trilinearly splatted into a smoothed vector field V (with a density weight W), the indicator function x is recovered by solving the screened-Poisson system (L_N + point_weight * W) x = -div V matrix-free with a Jacobi-preconditioned conjugate gradient, and the iso-surface is extracted with warp.MarchingCubes at the iso-value given by the average of x over the input points. A cascadic coarse-to-fine schedule solves the system from full_depth up to depth, prolonging each level's solution as the next level's initial guess (mirroring PoissonRecon's multigrid hierarchy).

normals are required and must be globally consistently oriented (all pointing outward or all inward): the reconstruction encodes the surface orientation in the sign of the vector field, and estimate_normals performs no global orientation.

Parameters:

Name	Type	Description	Default
points	array[vec3]	(n,) point positions on the target device (n >= 3).	required
normals	array[vec3]	(n,) oriented unit normals (their magnitude is used as a confidence weight when confidence is set).	required
depth	int	Finest grid depth: the finest grid has 2**depth + 1 nodes per axis. Memory grows as the cube of this; depth=8 (a 257**3 grid) is a safe default and depth is capped at 10. This is the single biggest lever on both cost and quality — see Notes before lowering it.	8
full_depth	int	Coarsest depth of the cascade (3 <= full_depth <= depth). The system is solved at every depth from full_depth to depth.	5
scale	float	Ratio between the reconstruction cube's side and the point cloud's largest bounding-box extent (>= 1, so the cube always contains the cloud); 1.1 pads the cloud by 10 %.	1.1
point_weight	float	Screening weight alpha tying the iso-surface to the input samples. 0 recovers the unscreened Poisson reconstruction (a tiny epsilon is still added to keep the operator SPD).	4.0
solver_iterations	int	Maximum conjugate-gradient iterations per cascade level.	100
solver_tolerance	float	Relative residual tolerance for the conjugate-gradient solve.	1e-06
confidence	bool	When True, weight each sample's splat by |normals[i]| (treating the normal magnitude as a per-sample confidence), matching PoissonRecon's confidence flag.	False
method	Literal['dense', 'adaptive']	Solver backend. "dense" (default) uses the dense node-centered grid described above. "adaptive" uses a warp.fem adaptive Nanogrid refined only near the samples with a variational (finite-element) assembly: fewer degrees of freedom for the same finest resolution, so it reaches higher depth on the same memory budget. The two backends agree up to discretization. A point-source weak form rings when the cells are much finer than the sampling, so the adaptive backend caps the finest near-surface cell at roughly the mean sample spacing (PoissonRecon-style); a depth above what the sampling supports then only refines the extraction lattice. It lazily imports warp.fem (a one-time codegen cost of tens of seconds on first use).	'dense'

Returns:

Name	Type	Description
vertices	array[vec3]	Iso-surface vertices on points.device.
faces	array[int32]	Flat 3 * n_faces triangle index buffer, free of zero-area triangles (remove_degenerate_faces) and oriented outward (make_normals_outward).

Raises:

Type	Description
ValueError	If points has fewer than 3 points, or the depth/scale parameters are out of range.
RuntimeError	If points and normals are not all on one device.

See Also

triangulate_point_cloud ball_pivoting estimate_normals

Notes

The whole system is solved in index space (grid spacing 1); the world scale is reapplied by marching cubes via the cube bounds, so the extracted level set is invariant to that scale. The per-node density normalization V /= max(W, eps) is cruder than PoissonRecon's kernel-density estimate and can bulge on strongly anisotropic clouds. The screening weight point_weight is not area-normalized as in PoissonRecon, so its numeric scale differs; 4 is a reasonable default for unit-scale clouds.

Choosing depth. The useful depth follows the cloud's sampling density, not its size, and past that point extra depth costs time without buying fidelity -- once the grid outruns the samples, accuracy can get worse rather than better. The default of 8 matches PoissonRecon's own; for a known cloud density it is worth lowering. The adaptive backend is much flatter in depth (it already caps near-surface refinement at the sample spacing) and is the better choice when the depth wanted for the extraction lattice exceeds what the sampling supports.

point_weight=0 is ill-conditioned, and its output is not reproducible. Screening is what conditions the operator; at 0 only the 1e-4 floor above keeps it SPD, so the conjugate gradient stops on a solution whose level set is genuinely uncertain, and the face count (and how many faces come out zero-area before the cleanup described under Returns) varies between runs on the same input; the surface is also not watertight. The default point_weight=4 is stable across runs and watertight. So 0 is for comparing against a screened reconstruction, not for producing one -- and do not pin a count taken from it.

Source code in triwarp/reconstruction.py

def screened_poisson(
    points: wp.array[wp.vec3],
    normals: wp.array[wp.vec3],
    *,
    depth: int = 8,
    full_depth: int = 5,
    scale: float = 1.1,
    point_weight: float = 4.0,
    solver_iterations: int = 100,
    solver_tolerance: float = 1e-6,
    confidence: bool = False,
    method: Literal["dense", "adaptive"] = "dense",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Screened-Poisson surface reconstruction from an oriented point cloud (Kazhdan PoissonRecon).

    GPU-native reimplementation of the screened-Poisson filter (PyMeshLab's
    ``generate_surface_reconstruction_screened_poisson``, ``open3d``'s
    ``create_from_point_cloud_poisson``).
    A dense node-centered grid over a padded bounding cube is used instead of PoissonRecon's octree:
    the oriented normals are trilinearly splatted into a smoothed vector field ``V`` (with a density
    weight ``W``), the indicator function ``x`` is recovered by solving the screened-Poisson system
    ``(L_N + point_weight * W) x = -div V`` matrix-free with a Jacobi-preconditioned conjugate
    gradient, and the iso-surface is extracted with ``warp.MarchingCubes`` at the iso-value given by
    the average of ``x`` over the input points. A **cascadic** coarse-to-fine schedule solves the
    system from ``full_depth`` up to ``depth``, prolonging each level's solution as the next level's
    initial guess (mirroring PoissonRecon's multigrid hierarchy).

    ``normals`` are **required** and must be globally consistently oriented (all pointing outward or
    all inward): the reconstruction encodes the surface orientation in the sign of the vector field,
    and [`estimate_normals`][triwarp.points.estimate_normals] performs no global orientation.

    Parameters
    ----------
    points
        ``(n,)`` point positions on the target device (``n >= 3``).
    normals
        ``(n,)`` oriented unit normals (their magnitude is used as a confidence weight when
        ``confidence`` is set).
    depth
        Finest grid depth: the finest grid has ``2**depth + 1`` nodes per axis. Memory grows as the
        cube of this; ``depth=8`` (a ``257**3`` grid) is a safe default and ``depth`` is capped at
        ``10``. This is the single biggest lever on both cost and quality — see Notes before
        lowering it.
    full_depth
        Coarsest depth of the cascade (``3 <= full_depth <= depth``). The system is solved at every
        depth from ``full_depth`` to ``depth``.
    scale
        Ratio between the reconstruction cube's side and the point cloud's largest bounding-box
        extent (``>= 1``, so the cube always contains the cloud); ``1.1`` pads the cloud by 10 %.
    point_weight
        Screening weight ``alpha`` tying the iso-surface to the input samples. ``0`` recovers the
        unscreened Poisson reconstruction (a tiny epsilon is still added to keep the operator SPD).
    solver_iterations
        Maximum conjugate-gradient iterations per cascade level.
    solver_tolerance
        Relative residual tolerance for the conjugate-gradient solve.
    confidence
        When ``True``, weight each sample's splat by ``|normals[i]|`` (treating the normal magnitude
        as a per-sample confidence), matching PoissonRecon's ``confidence`` flag.
    method
        Solver backend. ``"dense"`` (default) uses the dense node-centered grid described above.
        ``"adaptive"`` uses a ``warp.fem`` adaptive Nanogrid refined only near the samples with a
        variational (finite-element) assembly: fewer degrees of freedom for the same finest
        resolution, so it reaches higher ``depth`` on the same memory budget. The two backends agree
        up to discretization. A point-source weak form rings when the cells are much finer than the
        sampling, so the adaptive backend caps the finest near-surface cell at roughly the mean
        sample spacing (PoissonRecon-style); a ``depth`` above what the sampling supports then only
        refines the extraction lattice. It lazily imports ``warp.fem`` (a one-time codegen cost of
        tens of seconds on first use).

    Returns
    -------
    vertices : wp.array[wp.vec3]
        Iso-surface vertices on ``points.device``.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer, free of zero-area triangles
        ([`remove_degenerate_faces`][triwarp.repair.remove_degenerate_faces]) and oriented outward
        ([`make_normals_outward`][triwarp.repair.make_normals_outward]).

    Raises
    ------
    ValueError
        If ``points`` has fewer than 3 points, or the depth/scale parameters are out of range.
    RuntimeError
        If ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`triangulate_point_cloud`][triwarp.reconstruction.triangulate_point_cloud]
    [`ball_pivoting`][triwarp.reconstruction.ball_pivoting]
    [`estimate_normals`][triwarp.points.estimate_normals]

    Notes
    -----
    The whole system is solved in index space (grid spacing ``1``); the world scale is reapplied by
    marching cubes via the cube bounds, so the extracted level set is invariant to that scale. The
    per-node density normalization ``V /= max(W, eps)`` is cruder than PoissonRecon's kernel-density
    estimate and can bulge on strongly anisotropic clouds. The screening weight ``point_weight`` is
    not area-normalized as in PoissonRecon, so its numeric scale differs; ``4`` is a reasonable
    default for unit-scale clouds.

    **Choosing ``depth``.** The useful depth follows the cloud's *sampling density*, not its size,
    and past that point extra depth costs time without buying fidelity -- once the grid outruns the
    samples, accuracy can get worse rather than better. The default of ``8`` matches PoissonRecon's
    own; for a known cloud density it is worth lowering. The ``adaptive`` backend is much flatter in
    depth (it already caps near-surface refinement at the sample spacing) and is the better choice
    when the depth wanted for the extraction lattice exceeds what the sampling supports.

    **``point_weight=0`` is ill-conditioned, and its output is not reproducible.** Screening is what
    conditions the operator; at ``0`` only the ``1e-4`` floor above keeps it SPD, so the conjugate
    gradient stops on a solution whose level set is genuinely uncertain, and the face count (and
    how many faces come out zero-area before the cleanup described under ``Returns``) varies
    between runs on the same input; the surface is also not watertight. The default
    ``point_weight=4`` is stable across runs and watertight. So ``0`` is for comparing *against* a
    screened reconstruction, not for producing one -- and do not pin a count taken from it.
    """
    require_same_device(points=points, normals=normals)
    if not (3 <= full_depth <= depth <= 10):
        raise ValueError(
            "screened_poisson requires 3 <= full_depth <= depth <= 10, got "
            f"full_depth={full_depth}, depth={depth}."
        )
    if scale < 1.0:
        # A cube smaller than the cloud's own bounding box leaves points outside
        # ``[cube_lower, cube_upper]``; the grid samplers clamp their cell index into range rather
        # than raising, so those points would silently splat onto the boundary node instead of
        # being rejected or the cube being grown.
        raise ValueError(f"screened_poisson requires scale >= 1, got {scale}.")
    if method not in ("dense", "adaptive"):
        raise ValueError(f"screened_poisson method must be 'dense' or 'adaptive', got {method!r}.")

    device = points.device
    n = int(points.shape[0])
    if n < 3:
        raise ValueError(f"screened_poisson requires at least 3 points, got {n}.")

    cube_lower, cube_upper, cube_size = _poisson_cube(points, scale)
    # Effective screening weight: a floor keeps the operator SPD even at point_weight == 0.
    screen = max(float(point_weight), 1e-4)

    if method == "adaptive":
        vertices, faces = _screened_poisson_adaptive(
            points,
            normals,
            cube_lower,
            cube_upper,
            cube_size,
            depth=depth,
            full_depth=full_depth,
            screen=screen,
            solver_iterations=solver_iterations,
            solver_tolerance=solver_tolerance,
            confidence=confidence,
        )
    else:
        solution, res = _poisson_dense_solve(
            points,
            normals,
            cube_lower,
            cube_size,
            depth,
            full_depth,
            screen,
            confidence,
            solver_iterations,
            solver_tolerance,
            device,
        )
        inv_cell = float(res - 1) / cube_size
        sampled = wp.empty(n, dtype=wp.float32, device=device)
        wp.launch(
            kernel_reconstruction.sample_field_trilinear,
            dim=n,
            inputs=[solution, res, cube_lower, inv_cell, points, sampled],
            device=device,
        )
        iso = _poisson_iso_value(sampled, normals, confidence)
        vertices, faces = _extract_poisson_surface(solution, res, iso, cube_lower, cube_upper)

    # Drop zero-area triangles before orienting. Marching cubes emits one wherever the level set
    # grazes a lattice node, and such a face has no normal for ``make_normals_outward`` to orient
    # and hands the caller a NaN out of any closest-point query -- trimesh's ``closest_point``
    # divides by the squared length of the zero-length edge. This is deliberately *not* the full
    # ``_clean_reconstruction`` tail the other three reconstructions use: welding the coincident
    # vertices as well preserves the boundary-edge count but manufactures non-manifold edges, and
    # dedup would change the default path's output. As written it is **byte-identical** on a
    # well-screened reconstruction -- the default ``point_weight`` emits no degenerate face at all.
    vertices, faces = tw.repair.remove_degenerate_faces(vertices, faces)
    if int(faces.shape[0]) > 0:
        faces = tw.repair.make_normals_outward(vertices, faces)
    return vertices, faces

resample_uniform

resample_uniform(
    vertices: array[vec3],
    faces: array[int32],
    *,
    voxel_size: float | None = None,
    offset: float = 0.0,
    sign_mode: Literal["parity", "winding"] = "winding",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Rebuild a mesh by sampling its signed distance field on a uniform grid and re-extracting it.

MeshLab's generate_resampled_uniform_mesh. Three uses, all the same call:

• Regularize. The output triangulation comes from the grid, not from the input, so every pathology of the input topology — self-intersections, non-manifold edges, duplicated or inverted faces, a soup — is simply not carried over. This is the bluntest repair there is and the one that always works.
• Offset. A positive offset extracts the level set outside the surface (a dilation) and a negative one inside (an erosion): how shells, clearances and tool paths are built.
• Shrink-wrap. A coarse voxel_size with offset=0 gives a watertight envelope of a complicated input.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
voxel_size	float | None	Grid spacing. Defaults to 1 % of the bounding-box diagonal, which is a ~100-cell grid across the mesh. Cost is cubic in the reciprocal, so halving it is eight times the field evaluation; and no feature thinner than a voxel survives.	None
offset	float	Level set to extract, in world units. 0 (the default) reproduces the surface.	0.0
sign_mode	Literal['parity', 'winding']	How the sign of the distance field is decided, forwarded to signed_distance_on_mesh. The default here is "winding" rather than that function's "parity": resampling is usually applied to a mesh that is broken, and ray parity has no principled answer through a hole, while the generalized winding number degrades gracefully.	'winding'

Returns:

Name	Type	Description
vertices	array[vec3]	Resampled vertex positions on vertices.device.
faces	array[int32]	Flat 3 * n_faces triangle index buffer.

Raises:

Type	Description
ValueError	If voxel_size is not positive.
RuntimeError	If vertices and faces are not all on one device.

See Also

marching_cubes triwarp.proximity.signed_distance_on_mesh triwarp.remesh.cluster_decimate triwarp.voxels.voxelize_mesh

Notes

The grid is padded by three voxels beyond the bounding box plus offset, so a dilated level set is never clipped by the lattice boundary and the extracted surface is always closed. That padding is why the memory cost is a little above (extent / voxel_size) ** 3.

Source code in triwarp/reconstruction.py

def resample_uniform(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    voxel_size: float | None = None,
    offset: float = 0.0,
    sign_mode: Literal["parity", "winding"] = "winding",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Rebuild a mesh by sampling its signed distance field on a uniform grid and re-extracting it.

    MeshLab's ``generate_resampled_uniform_mesh``. Three uses, all the same call:

    - **Regularize.** The output triangulation comes from the grid, not from the input, so every
      pathology of the input topology — self-intersections, non-manifold edges, duplicated or
      inverted faces, a soup — is simply not carried over. This is the bluntest repair there is and
      the one that always works.
    - **Offset.** A positive ``offset`` extracts the level set *outside* the surface (a dilation)
      and a negative one inside (an erosion): how shells, clearances and tool paths are built.
    - **Shrink-wrap.** A coarse ``voxel_size`` with ``offset=0`` gives a watertight envelope of a
      complicated input.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    voxel_size
        Grid spacing. Defaults to ``1 %`` of the bounding-box diagonal, which is a ~100-cell grid
        across the mesh. **Cost is cubic in the reciprocal**, so halving it is eight times the field
        evaluation; and no feature thinner than a voxel survives.
    offset
        Level set to extract, in world units. ``0`` (the default) reproduces the surface.
    sign_mode
        How the sign of the distance field is decided, forwarded to
        [`signed_distance_on_mesh`][triwarp.proximity.signed_distance_on_mesh]. The default here is
        ``"winding"`` rather than that function's ``"parity"``: resampling is usually applied to a
        mesh that is *broken*, and ray parity has no principled answer through a hole, while the
        generalized winding number degrades gracefully.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        Resampled vertex positions on ``vertices.device``.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer.

    Raises
    ------
    ValueError
        If ``voxel_size`` is not positive.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`marching_cubes`][triwarp.levelset.marching_cubes]
    [`triwarp.proximity.signed_distance_on_mesh`][triwarp.proximity.signed_distance_on_mesh]
    [`triwarp.remesh.cluster_decimate`][triwarp.remesh.cluster_decimate]
    [`triwarp.voxels.voxelize_mesh`][triwarp.voxels.voxelize_mesh]

    Notes
    -----
    The grid is padded by three voxels beyond the bounding box *plus* ``offset``, so a dilated level
    set is never clipped by the lattice boundary and the extracted surface is always closed. That
    padding is why the memory cost is a little above ``(extent / voxel_size) ** 3``.
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_faces = int(faces.shape[0]) // 3
    if n_faces == 0:
        return wp.clone(vertices), wp.clone(faces)

    lower, upper = tw.bounds.aabb(vertices)
    diagonal = float(wp.length(upper - lower))
    if voxel_size is None:
        voxel_size = 0.01 * diagonal
    if voxel_size <= 0.0:
        raise ValueError(f"resample_uniform requires voxel_size > 0, got {voxel_size}")

    # Three voxels of slack beyond the box and beyond the offset, so a dilated level set closes
    # inside the lattice instead of being cut off by it.
    pad = 3.0 * voxel_size + max(offset, 0.0)
    grid_lower = wp.vec3(lower[0] - pad, lower[1] - pad, lower[2] - pad)
    grid_upper = wp.vec3(upper[0] + pad, upper[1] + pad, upper[2] + pad)
    resolution = wp.vec3i(
        *(
            max(2, math.ceil((grid_upper[axis] - grid_lower[axis]) / voxel_size) + 1)
            for axis in range(3)
        )
    )
    spacing = wp.vec3(
        *((grid_upper[axis] - grid_lower[axis]) / float(resolution[axis] - 1) for axis in range(3))
    )

    n_points = int(resolution[0]) * int(resolution[1]) * int(resolution[2])
    points = wp.empty(n_points, dtype=wp.vec3, device=device)
    wp.launch(
        kernel_reconstruction.lattice_points,
        dim=(int(resolution[0]), int(resolution[1]), int(resolution[2])),
        inputs=[resolution, grid_lower, spacing, points],
        device=device,
    )
    field = tw.proximity.signed_distance_on_mesh(vertices, faces, points, sign_mode=sign_mode)
    out_vertices, out_faces = tw.levelset.marching_cubes(
        twt.as_array3d(
            field.reshape((int(resolution[0]), int(resolution[1]), int(resolution[2]))), wp.float32
        ),
        iso=offset,
        bounds=(grid_lower, grid_upper),
    )
    # Warp's extractor emits a vertex per crossing per cell, so coincident duplicates are normal
    # rather than exceptional; welding them is what makes the result a closed surface. Hole filling
    # is off (``crit_hole_length=0``): a level set of a signed field is already closed, so any hole
    # here would be a symptom worth surfacing rather than patching over.
    return _clean_reconstruction(out_vertices, out_faces, 0.0)

ball_pivoting

ball_pivoting(
    points: array[vec3],
    normals: array[vec3] | None = None,
    *,
    radius: float = 0.0,
    clustering: float = 0.2,
    crease_angle: float = math.pi / 2.0,
    max_waves: int = 0,
    crit_hole_length: float = 0.0,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Ball-pivoting surface reconstruction from an oriented point cloud (interpolating).

GPU-native, wave-parallel port of the ball-pivoting algorithm (PyMeshLab's generate_surface_reconstruction_ball_pivoting, open3d's create_from_point_cloud_ball_pivoting, vcglib BPA). A ball of the given radius is rolled over the oriented points; wherever it rests on three points without containing another point it forms a triangle, and it then pivots around each boundary edge of the growing front to add the next triangle. The output mesh interpolates the input points (its vertices are a subset of points), unlike the implicit screened_poisson.

The advancing front is a persistent device-side structure — an edge hash table plus a compacted list of its boundary edges, mutated in place as triangles commit. Each wave pivots every live front edge independently through a hash-grid ball query and commits a conflict-free independent set via a two-phase vertex claim. Manifoldness is protected by an interior-edge guard (an edge already shared by two triangles never gains a third). The result is passed through the standard cleanup tail (deduplicate, drop degenerate / non-manifold faces, orient outward, optionally fill small holes).

Parameters:

Name	Type	Description	Default
points	array[vec3]	(n,) point positions on the target device (n >= 3).	required
normals	array[vec3] | None	Optional (n,) oriented unit normals. When None, they are estimated by PCA oriented outward from the centroid (estimate_normals; valid for star-shaped clouds only).	None
radius	float	Ball radius. When <= 0, it is auto-guessed as 1.5 x the mean nearest-neighbour spacing — a device reduction, so the guess is not bit-reproducible between calls and neither is the triangulation built from it (see Notes).	0.0
clustering	float	Candidate points closer than clustering * radius to an edge endpoint are rejected (vcglib clustering, as a fraction of radius).	0.2
crease_angle	float	Maximum dihedral angle (radians) between a new triangle and the one it pivots from; larger folds are rejected. Use >= pi to disable the crease guard.	pi / 2.0
max_waves	int	Maximum number of parallel waves, where each wave either pivots the current front or seeds orphan triangles. 0 picks a generous safety cap from the point count.	0
crit_hole_length	float	Boundary loops with perimeter at most this value are filled at the end. 0 disables hole filling; a negative value defaults to 0.1 x the bounding-box diagonal.	0.0

Returns:

Name	Type	Description
vertices	array[vec3]	The referenced input points, compacted from index zero, on points.device.
faces	array[int32]	Flat 3 * n_faces triangle index buffer, oriented outward.

Raises:

Type	Description
ValueError	If fewer than 3 points are given.
RuntimeError	If points and normals are not all on one device.

See Also

screened_poisson triangulate_point_cloud estimate_normals

Notes

A single radius is used (the classic multi-radius schedule is a documented follow-up), and the triangle budget grows on demand rather than raising.

The triangulation is reproducible at an explicit radius: repeated runs on one device and build return the same set of triangles, because a wave resolves competing proposals by a key packed from their own vertices and the wave counter rather than by the order they reached an atomic counter — device state the loop advances deterministically, so the ordering varies between waves (which keeps a front edge from being starved by a globally fixed key) without varying between runs. Three caveats on how far that reaches:

• the default radius <= 0 is not covered. The auto-guess averages the nearest-neighbour spacings with a device reduction whose summation order is not fixed, so two calls on one cloud can guess radii a few ULP apart — enough to tip a borderline pivot and change the triangle set. Pass radius explicitly wherever the result has to be reproducible;
• the face buffer's row order is not pinned — commit_triangles appends with a wp.atomic_add — so compare reconstructions as a set of triangles, not buffer-to-buffer;
• neither is the winding of a component that the cleanup tail cannot orient by volume, since make_winding_consistent seeds each connected component from an arbitrary face. Compare with the row entries sorted, or orient both sides first.

Results are not comparable across devices, where floating-point contraction differs.

The result is interpolating and edge-manifold, and on a densely, uniformly sampled closed surface it is watertight in practice: a subdivided icosphere reconstructs to exactly 2 v - 4 faces with no boundary edge at all, and bunny (a real scan, so genuinely open in places) comes out at 1.97 faces per referenced vertex with 1.7% boundary edges. Watertightness is not guaranteed — the wave-parallel front has no per-vertex fan structure to align colliding sheets with, so pathological sampling can still leave a seam. Regions sampled too sparsely for the ball to rest are left as holes by construction (widen radius or set crit_hole_length); for an implicit surface that is watertight by construction use screened_poisson.

Source code in triwarp/reconstruction.py

def ball_pivoting(
    points: wp.array[wp.vec3],
    normals: wp.array[wp.vec3] | None = None,
    *,
    radius: float = 0.0,
    clustering: float = 0.2,
    crease_angle: float = math.pi / 2.0,
    max_waves: int = 0,
    crit_hole_length: float = 0.0,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Ball-pivoting surface reconstruction from an oriented point cloud (interpolating).

    GPU-native, wave-parallel port of the ball-pivoting algorithm (PyMeshLab's
    ``generate_surface_reconstruction_ball_pivoting``, ``open3d``'s
    ``create_from_point_cloud_ball_pivoting``, vcglib BPA). A ball of the given ``radius`` is rolled
    over the oriented points; wherever it rests on three points without containing another point it
    forms a triangle, and it then pivots around each boundary edge of the growing front to add the
    next triangle. The output mesh **interpolates** the input points (its vertices are a subset of
    ``points``), unlike the implicit
    [`screened_poisson`][triwarp.reconstruction.screened_poisson].

    The advancing front is a **persistent** device-side structure — an edge hash table plus a
    compacted list of its boundary edges, mutated in place as triangles commit. Each wave pivots
    every live front edge independently through a hash-grid ball query and commits a conflict-free
    independent set via a two-phase vertex claim. Manifoldness is protected by an interior-edge
    guard (an edge already shared by two triangles never gains a third). The result is passed
    through the standard cleanup tail (deduplicate, drop degenerate / non-manifold faces, orient
    outward, optionally fill small holes).

    Parameters
    ----------
    points
        ``(n,)`` point positions on the target device (``n >= 3``).
    normals
        Optional ``(n,)`` oriented unit normals. When ``None``, they are estimated by PCA oriented
        outward from the centroid ([`estimate_normals`][triwarp.points.estimate_normals]; valid for
        star-shaped clouds only).
    radius
        Ball radius. When ``<= 0``, it is auto-guessed as ``1.5 x`` the mean nearest-neighbour
        spacing — a device reduction, so the guess is not bit-reproducible between calls and neither
        is the triangulation built from it (see ``Notes``).
    clustering
        Candidate points closer than ``clustering * radius`` to an edge endpoint are rejected
        (vcglib clustering, as a fraction of ``radius``).
    crease_angle
        Maximum dihedral angle (radians) between a new triangle and the one it pivots from; larger
        folds are rejected. Use ``>= pi`` to disable the crease guard.
    max_waves
        Maximum number of parallel waves, where each wave either pivots the current front or seeds
        orphan triangles. ``0`` picks a generous safety cap from the point count.
    crit_hole_length
        Boundary loops with perimeter at most this value are filled at the end. ``0`` disables hole
        filling; a negative value defaults to ``0.1 x`` the bounding-box diagonal.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        The referenced input points, compacted from index zero, on ``points.device``.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer, oriented outward.

    Raises
    ------
    ValueError
        If fewer than 3 points are given.
    RuntimeError
        If ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`screened_poisson`][triwarp.reconstruction.screened_poisson]
    [`triangulate_point_cloud`][triwarp.reconstruction.triangulate_point_cloud]
    [`estimate_normals`][triwarp.points.estimate_normals]

    Notes
    -----
    A single radius is used (the classic multi-radius schedule is a documented follow-up), and the
    triangle budget grows on demand rather than raising.

    The triangulation is **reproducible at an explicit ``radius``**: repeated runs on one device and
    build return the same set of triangles, because a wave resolves competing proposals by a key
    packed from their own vertices and the wave counter rather than by the order they reached an
    atomic counter — device state the loop advances deterministically, so the ordering varies
    between waves (which keeps a front edge from being starved by a globally fixed key) without
    varying between runs. Three caveats on how far that reaches:

    * the default ``radius <= 0`` is **not** covered. The auto-guess averages the nearest-neighbour
      spacings with a device reduction whose summation order is not fixed, so two calls on one cloud
      can guess radii a few ULP apart — enough to tip a borderline pivot and change the triangle
      set. Pass ``radius`` explicitly wherever the result has to be reproducible;
    * the face buffer's *row order* is not pinned — ``commit_triangles`` appends with a
      ``wp.atomic_add`` — so compare reconstructions as a set of triangles, not buffer-to-buffer;
    * neither is the *winding* of a component that the cleanup tail cannot orient by volume, since
      [`make_winding_consistent`][triwarp.repair.make_winding_consistent] seeds each connected
      component from an arbitrary face. Compare with the row entries sorted, or orient both sides
      first.

    Results are not comparable across devices, where floating-point contraction differs.

    The result is **interpolating and edge-manifold**, and on a densely, uniformly sampled closed
    surface it is watertight in practice: a subdivided icosphere reconstructs to exactly ``2 v - 4``
    faces with no boundary edge at all, and ``bunny`` (a real scan, so genuinely open in places)
    comes out at 1.97 faces per referenced vertex with 1.7% boundary edges. Watertightness is not
    *guaranteed* — the wave-parallel front has no per-vertex fan structure to align colliding
    sheets with, so pathological sampling can still leave a seam. Regions sampled too sparsely for
    the ball to rest are left as holes by construction (widen ``radius`` or set
    ``crit_hole_length``); for an implicit surface that is watertight by construction use
    [`screened_poisson`][triwarp.reconstruction.screened_poisson].
    """
    require_same_device(points=points, normals=normals)
    device = points.device
    n = int(points.shape[0])
    if n < 3:
        return wp.clone(points), wp.empty(0, dtype=wp.int32, device=device)

    # Nearest-neighbour table drives both the radius auto-guess and (if needed) normal estimation.
    neighbor_idx, neighbor_dist = tw.neighbors.query_nearest(points, points, k=7, backend="bvh")
    if radius <= 0.0:
        # The whole table, not columns ``1:``: the self-distance in slot 0 is exactly zero and the
        # positive-finite filter drops it, so flattening costs nothing and keeps the reduction on
        # the device -- a column slice would be strided and could not be flattened at all.
        mean_spacing = _mean_positive_finite(neighbor_dist.flatten())
        radius = 1.5 * (mean_spacing if mean_spacing is not None else 1.0)
    if normals is None:
        normals = tw.points.estimate_normals(points, neighbor_idx)

    # Cell width equal to the ball radius, not to the ``2 * radius`` neighbourhood the pivot
    # searches: the inner empty-ball test is by far the most frequent query, and a cell twice its
    # radius makes it enumerate roughly 8x the points it needs. Going finer than this loses more to
    # cell-probe overhead than it saves in point tests.
    grid = tw.neighbors.hashgrid_from_points(points, radius)
    # The pivot search walks this cooperatively (one warp an edge); the empty-ball test inside
    # it stays on the hash grid, which is the better structure for a *serial* per-lane query.
    bvh = tw.neighbors.bvh_from_points(points)
    crease_cos = math.cos(crease_angle) if crease_angle < math.pi else -1.0

    state = _BpaState(points, normals, grid, bvh, radius, clustering, crease_cos, 4 * n + 16)
    _bpa_run(state, max_waves if max_waves > 0 else 16 * n)

    count = int(read_scalar(state.counters, kernel_bpa.CNT_FACE))
    faces = wp.clone(state.all_faces[: count * 3])
    return _clean_reconstruction(points, faces, crit_hole_length)