proximity

triwarp.proximity

Queries that ask where a point stands relative to a triangle mesh.

Three answers, in increasing order of what they need from the mesh: closest_point_on_mesh and normals_at_closest_faces need only a surface; signed_distance_on_mesh needs a consistent winding to give the distance a sign; and winding_number needs neither watertightness nor manifoldness, which is why it is the robust inside test on damaged input. containing_faces_2d is the planar case -- point location in a 2D triangulation -- and query_mesh_aabb_with_offsets is the low-level box query the others are built over.

Everything here takes raw (vertices, faces) buffers. The queries phrased the other way round -- how far away the surface is rather than where it is, all taking a prebuilt wp.Mesh and a set of points to measure at -- live in triwarp.visibility: ambient occlusion and obscurance outward, shape diameter and thickness inward, and the maximal tangent sphere in every direction at once.

Everything signed here, plus contains_points, follows Warp's SDF sign convention: outside positive, inside negative. Trimesh's signed_distance uses the opposite sign.

Point-set acceleration structures (wp.Bvh / wp.HashGrid) and raw neighbor queries live in triwarp.neighbors; axis-aligned bounding boxes in triwarp.bounds.

ITEMS_PER_QUERY_SLICE module-attribute

ITEMS_PER_QUERY_SLICE = 128

closest_point_on_mesh

closest_point_on_mesh(
    vertices: array[vec3],
    faces: array[int32],
    points: array[vec3],
    *,
    max_dist: float | None = None,
    mesh: Mesh | None = None,
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.float32],
    wp.array[wp.int32],
]

For each query point, find the closest point on any triangle of the mesh.

Uses wp.mesh_query_point_no_sign via wp.Mesh. Distances are unsigned Euclidean lengths in float32.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n,) mesh vertex positions as wp.vec3.	required
faces	array[int32]	(f * 3,) flat triangle index array as wp.int32. The internally built wp.Mesh aliases these buffers rather than copying them; do not mutate them for the duration of the call.	required
points	array[vec3]	(m,) query positions in space as wp.vec3.	required
max_dist	float | None	Maximum search radius per query. Faces farther than this are ignored. When None, derived from the axis-aligned box enclosing mesh vertices and query points.	None
mesh	Mesh | None	A wp.Mesh already built over vertices and faces, to spare the clone and BVH build this otherwise pays on every call. Purely an optimization: the answer is identical either way, and it is not checked against vertices / faces -- passing a mesh over different geometry silently answers for that geometry, since only mesh is queried. Trimesh.warp_mesh is a cached property and is what to pass.	None

Returns:

Type	Description
closest	(m, 3) closest point on the mesh surface for each query.
distance	(m,) unsigned distance from each query to its closest surface point.
triangle_id	(m,) index of the triangle containing each closest point, or -1 when no face lies within max_dist.

Raises:

Type	Description
RuntimeError	If vertices, faces, points and mesh are not all on one device.

Source code in triwarp/proximity.py

def closest_point_on_mesh(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    points: wp.array[wp.vec3],
    *,
    max_dist: float | None = None,
    mesh: wp.Mesh | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.float32], wp.array[wp.int32]]:
    """
    For each query point, find the closest point on any triangle of the mesh.

    Uses ``wp.mesh_query_point_no_sign`` via ``wp.Mesh``. Distances are unsigned
    Euclidean lengths in ``float32``.

    Parameters
    ----------
    vertices
        ``(n,)`` mesh vertex positions as ``wp.vec3``.
    faces
        ``(f * 3,)`` flat triangle index array as ``wp.int32``.
        The internally built ``wp.Mesh`` aliases these buffers rather than copying them;
        do not mutate them for the duration of the call.
    points
        ``(m,)`` query positions in space as ``wp.vec3``.
    max_dist
        Maximum search radius per query. Faces farther than this are ignored.
        When ``None``, derived from the axis-aligned box enclosing mesh
        vertices and query points.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the clone and BVH
        build this otherwise pays on every call. Purely an optimization: the answer is identical
        either way, and it is not checked against ``vertices`` / ``faces`` -- passing a mesh over
        *different* geometry silently answers for that geometry, since only ``mesh`` is queried.
        [`Trimesh.warp_mesh`][triwarp.mesh.Trimesh.warp_mesh] is a cached property and is what to
        pass.

    Returns
    -------
    closest
        ``(m, 3)`` closest point on the mesh surface for each query.
    distance
        ``(m,)`` unsigned distance from each query to its closest surface point.
    triangle_id
        ``(m,)`` index of the triangle containing each closest point, or ``-1``
        when no face lies within ``max_dist``.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces``, ``points`` and ``mesh`` are not all on one device.
    """
    require_same_device(vertices=vertices, faces=faces, points=points, mesh=mesh)
    device = vertices.device
    m = int(points.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if m == 0:
        return (
            wp.empty(0, dtype=wp.vec3, device=device),
            wp.empty(0, dtype=wp.float32, device=device),
            wp.empty(0, dtype=wp.int32, device=device),
        )
    if n_faces == 0:
        nan = float("nan")
        out_closest = wp.full(m, wp.vec3(nan, nan, nan), dtype=wp.vec3, device=device)
        out_distance = wp.full(m, float("inf"), dtype=wp.float32, device=device)
        out_face = wp.full(m, -1, dtype=wp.int32, device=device)
        return out_closest, out_distance, out_face

    if mesh is None:
        require_nonempty_mesh(faces, "closest_point_on_mesh")
        # The mesh aliases the caller's buffers and is discarded here, so it needs no copy.
        mesh = wp.Mesh(points=vertices, indices=faces)
    if max_dist is None:
        max_dist = tw.bounds.enclosing_diagonal(mesh.points, points)

    out_closest = wp.empty(m, dtype=wp.vec3, device=device)
    out_distance = wp.empty(m, dtype=wp.float32, device=device)
    out_face = wp.empty(m, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.closest_point_on_mesh,
        dim=m,
        inputs=[mesh.id, points, wp.float32(max_dist), out_closest, out_distance, out_face],
        device=device,
    )
    return out_closest, out_distance, out_face

closest_point_on_edges

closest_point_on_edges(
    vertices: array[vec3],
    edges: Array2dInt32,
    queries: array[vec3],
    *,
    max_dist: float | None = None,
    bvh: Bvh | None = None,
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.float32],
    wp.array[wp.int32],
]

For each query point, find the closest point on any edge of an edge set.

The wireframe counterpart of closest_point_on_mesh: same return triple, same max_dist semantics, but the geometry is a set of segments rather than a surface. That is the query a crease set, a seam, a boundary rim or a feature curve wants -- all four are already produced as an edge array by triwarp.seams, triwarp.boundary and triwarp.edges.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n,) positions the edges index, as wp.vec3.	required
edges	Array2dInt32	(n_edges, 2) wp.int32 vertex-index pairs. Order within a pair is irrelevant, and edges may share vertices or repeat.	required
queries	array[vec3]	(m,) query positions in space as wp.vec3.	required
max_dist	float | None	Maximum search distance per query; an edge farther than this is ignored and the query reports a miss. When None, derived from the box enclosing vertices and queries, which no real query can exceed.	None
bvh	Bvh | None	A wp.Bvh already built over this edge set's per-edge boxes, to spare the bounds pass and the build. Purely an optimization -- and, as with closest_point_on_mesh's mesh, it is not checked against vertices / edges: a BVH over a different edge set silently answers for the boxes it holds while the distances are computed from these vertices.	None

Returns:

Type	Description
closest	(m,) closest point on the edge set for each query, as wp.vec3.
distance	(m,) unsigned distance from each query to that point.
edge_id	(m,) row of edges the closest point lies on, or -1 when no edge lies within max_dist.

Raises:

Type	Description
TypeError	If edges is not a rank-2 wp.int32 array.
ValueError	If edges does not have two columns.
RuntimeError	If vertices, edges, queries and bvh are not all on one device.

Notes

A miss reports the query point itself and max_dist, alongside the -1 id -- the same convention closest_point_on_mesh uses, so the two are interchangeable in a caller that only reads the id.

The traversal is iterative deepening over a per-edge-box BVH and is exact, not a broad-phase approximation: a scan of the cube of half-extent r about a query enumerates every edge whose closest point lies within r, so a best distance under r certifies the answer. The tempting shortcut -- one degenerate (a, b, b) triangle per edge, queried with wp.mesh_query_point_no_sign -- does not work: Warp's mesh BVH rejects a zero-area triangle outright, on both devices.

See Also

closest_point_on_mesh The surface form. On a closed mesh its answer is never farther than this one's. triwarp.polyline.polyline_point_distance The same computation for an ordered chain, where the segments are consecutive vertices and no index structure is built.

Source code in triwarp/proximity.py

def closest_point_on_edges(
    vertices: wp.array[wp.vec3],
    edges: twt.Array2dInt32,
    queries: wp.array[wp.vec3],
    *,
    max_dist: float | None = None,
    bvh: wp.Bvh | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.float32], wp.array[wp.int32]]:
    """
    For each query point, find the closest point on any edge of an edge set.

    The **wireframe** counterpart of
    [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]: same return triple, same
    ``max_dist`` semantics, but the geometry is a set of segments rather than a surface. That is the
    query a crease set, a seam, a boundary rim or a feature curve wants -- all four are already
    produced as an edge array by [`triwarp.seams`][triwarp.seams],
    [`triwarp.boundary`][triwarp.boundary] and [`triwarp.edges`][triwarp.edges].

    Parameters
    ----------
    vertices
        ``(n,)`` positions the edges index, as ``wp.vec3``.
    edges
        ``(n_edges, 2)`` ``wp.int32`` vertex-index pairs. Order within a pair is irrelevant, and
        edges may share vertices or repeat.
    queries
        ``(m,)`` query positions in space as ``wp.vec3``.
    max_dist
        Maximum search distance per query; an edge farther than this is ignored and the query
        reports a miss. When ``None``, derived from the box enclosing ``vertices`` and ``queries``,
        which no real query can exceed.
    bvh
        A ``wp.Bvh`` already built over this edge set's per-edge boxes, to spare the bounds pass and
        the build. Purely an optimization -- and, as with
        [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]'s ``mesh``, it is not
        checked against ``vertices`` / ``edges``: a BVH over a *different* edge set silently answers
        for the boxes it holds while the distances are computed from these vertices.

    Returns
    -------
    closest
        ``(m,)`` closest point on the edge set for each query, as ``wp.vec3``.
    distance
        ``(m,)`` unsigned distance from each query to that point.
    edge_id
        ``(m,)`` row of ``edges`` the closest point lies on, or ``-1`` when no edge lies within
        ``max_dist``.

    Raises
    ------
    TypeError
        If ``edges`` is not a rank-2 ``wp.int32`` array.
    ValueError
        If ``edges`` does not have two columns.
    RuntimeError
        If ``vertices``, ``edges``, ``queries`` and ``bvh`` are not all on one device.

    Notes
    -----
    A miss reports the query point itself and ``max_dist``, alongside the ``-1`` id -- the same
    convention [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh] uses, so the two
    are interchangeable in a caller that only reads the id.

    The traversal is iterative deepening over a per-edge-box BVH and is **exact**, not a broad-phase
    approximation: a scan of the cube of half-extent ``r`` about a query enumerates every edge whose
    closest point lies within ``r``, so a best distance under ``r`` certifies the answer. The
    tempting shortcut -- one degenerate ``(a, b, b)`` triangle per edge, queried with
    ``wp.mesh_query_point_no_sign`` -- does **not** work: Warp's mesh BVH rejects a zero-area
    triangle outright, on both devices.

    See Also
    --------
    [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]
        The surface form. On a closed mesh its answer is never farther than this one's.
    [`triwarp.polyline.polyline_point_distance`][triwarp.polyline.polyline_point_distance]
        The same computation for an *ordered* chain, where the segments are consecutive vertices and
        no index structure is built.
    """
    require_same_device(vertices=vertices, edges=edges, queries=queries, bvh=bvh)
    device = vertices.device
    twt.ensure_edge_pairs(edges, "edges")
    m = int(queries.shape[0])
    n_edges = int(edges.shape[0])

    if m == 0:
        return (
            wp.empty(0, dtype=wp.vec3, device=device),
            wp.empty(0, dtype=wp.float32, device=device),
            wp.empty(0, dtype=wp.int32, device=device),
        )
    if n_edges == 0:
        return (
            wp.clone(queries),
            wp.full(m, float("inf"), dtype=wp.float32, device=device),
            wp.full(m, -1, dtype=wp.int32, device=device),
        )

    if bvh is None:
        lower = wp.empty(n_edges, dtype=wp.vec3, device=device)
        upper = wp.empty(n_edges, dtype=wp.vec3, device=device)
        wp.launch(
            kernel_edges.edge_aabb_bounds,
            dim=n_edges,
            inputs=[vertices, edges, lower, upper],
            device=device,
        )
        bvh = tw.neighbors.bvh_from_bounds(lower, upper)
    if max_dist is None:
        max_dist = tw.bounds.enclosing_diagonal(vertices, queries)
    # The scene box bounds each query's *complete* search radius, so a query outside the geometry
    # still terminates exactly rather than growing to ``max_dist``.
    min_bound, max_bound = tw.bounds.aabb(vertices)

    out_closest = wp.empty(m, dtype=wp.vec3, device=device)
    out_distance = wp.empty(m, dtype=wp.float32, device=device)
    out_edge = wp.empty(m, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.closest_point_on_edges,
        dim=m,
        inputs=[
            vertices,
            edges,
            queries,
            bvh.id,
            wp.float32(max_dist),
            wp.float32(_EDGE_INITIAL_RADIUS_SCALE * max_dist),
            min_bound,
            max_bound,
            out_closest,
            out_distance,
            out_edge,
        ],
        device=device,
    )
    return out_closest, out_distance, out_edge

mesh_to_mesh_distance

mesh_to_mesh_distance(
    vertices_a: array[vec3],
    faces_a: array[int32],
    vertices_b: array[vec3],
    faces_b: array[int32],
    *,
    upper_bound: float | None = None,
) -> tuple[float, int, int]

Smallest distance between two triangle meshes, with the pair of faces that achieves it.

The clearance between two parts, and zero when they touch or overlap. Unlike a vertex-to-mesh query this is the true minimum over the surfaces: two boxes edge to edge realise their clearance between edge interiors, and every vertex of each is further from the other than that.

Two phases. An upper bound comes first -- the smallest distance from a sample of A's vertices to B, which is a real distance between the surfaces and therefore an upper bound on their minimum. Then every face of A queries a BVH over B's faces with its own bounding box grown by that bound, and each candidate pair gets the exact triangle-triangle distance. The bound is what makes the broad phase sound rather than heuristic: the true minimum is at most the bound, so the pair achieving it has boxes within that distance and cannot be culled -- and that argument needs an upper bound rather than a tight one, which is why sampling is sound and why supplying your own coarse upper_bound is too. The returned distance is exact either way; a looser bound only leaves more candidates for the narrow phase to reject.

Parameters:

Name	Type	Description	Default
vertices_a	array[vec3]	First mesh: (n_vertices,) positions and a length-3 * n_faces index buffer.	required
faces_a	array[vec3]	First mesh: (n_vertices,) positions and a length-3 * n_faces index buffer.	required
vertices_b	array[vec3]	Second mesh, same layout.	required
faces_b	array[vec3]	Second mesh, same layout.	required
upper_bound	float | None	A distance known to be at least the answer, which prunes the broad phase. Supply one when you have it -- from a previous frame, or from a bounding-volume gap -- and the vertex query that would otherwise derive it is skipped. Too small a bound gives a wrong answer, not a slow one: it culls the pair that would have won. None derives a sound bound.	None

Returns:

Name	Type	Description
distance	float	The minimum distance. 0.0 exactly when some pair of faces crosses.
face_a	int	The face of the first mesh achieving it, or -1 if either mesh is empty.
face_b	int	The face of the second mesh achieving it, or -1.

Raises:

Type	Description
ValueError	If upper_bound is negative.
RuntimeError	If vertices_a, faces_a, vertices_b and faces_b are not all on one device.

See Also

closest_point_on_mesh Point-to-mesh, which is the query this derives its bound from. mesh_with_mesh The zero-distance case in detail: every intersecting pair and the segments they cross on. face_self_intersecting_mask The one-mesh analogue of that.

Notes

Witness faces are ambiguous under ties

Two parallel plates have a continuum of closest pairs and any of them is a correct answer, and ties are the common case rather than the exotic one: whenever the closest approach is realised at a vertex, every face around that vertex achieves the minimum exactly. The tie-break on face_a is the lowest index; face_b is unspecified among the faces attaining it, and which one comes back can differ between runs of the identical input even though the squared distance is bit-identical. Compare distances against another implementation, and faces only where the configuration is generic.

Source code in triwarp/proximity.py

def mesh_to_mesh_distance(
    vertices_a: wp.array[wp.vec3],
    faces_a: wp.array[wp.int32],
    vertices_b: wp.array[wp.vec3],
    faces_b: wp.array[wp.int32],
    *,
    upper_bound: float | None = None,
) -> tuple[float, int, int]:
    """
    Smallest distance between two triangle meshes, with the pair of faces that achieves it.

    The clearance between two parts, and zero when they touch or overlap. Unlike a vertex-to-mesh
    query this is the true minimum over the *surfaces*: two boxes edge to edge realise their
    clearance between edge interiors, and every vertex of each is further from the other than
    that.

    Two phases. An upper bound comes first -- the smallest distance from a *sample* of ``A``'s
    vertices to ``B``, which is a real distance between the surfaces and therefore an upper bound on
    their minimum. Then every face of ``A`` queries a BVH over ``B``'s faces with its own bounding
    box grown by that bound, and each candidate pair gets the exact triangle-triangle distance. The
    bound is what makes the broad phase sound rather than heuristic: the true minimum is at most the
    bound, so the pair achieving it has boxes within that distance and cannot be culled -- and that
    argument needs an upper bound rather than a *tight* one, which is why sampling is sound and why
    supplying your own coarse ``upper_bound`` is too. The returned distance is exact either way; a
    looser bound only leaves more candidates for the narrow phase to reject.

    Parameters
    ----------
    vertices_a, faces_a
        First mesh: ``(n_vertices,)`` positions and a length-``3 * n_faces`` index buffer.
    vertices_b, faces_b
        Second mesh, same layout.
    upper_bound
        A distance known to be at least the answer, which prunes the broad phase. Supply one when
        you have it -- from a previous frame, or from a bounding-volume gap -- and the vertex query
        that would otherwise derive it is skipped. **Too small a bound gives a wrong answer**, not a
        slow one: it culls the pair that would have won. ``None`` derives a sound bound.

    Returns
    -------
    distance : float
        The minimum distance. ``0.0`` exactly when some pair of faces crosses.
    face_a : int
        The face of the first mesh achieving it, or ``-1`` if either mesh is empty.
    face_b : int
        The face of the second mesh achieving it, or ``-1``.

    Raises
    ------
    ValueError
        If ``upper_bound`` is negative.
    RuntimeError
        If ``vertices_a``, ``faces_a``, ``vertices_b`` and ``faces_b`` are not all on one device.

    See Also
    --------
    [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]
        Point-to-mesh, which is the query this derives its bound from.
    [`mesh_with_mesh`][triwarp.intersection.mesh_with_mesh]
        The zero-distance case in detail: every intersecting pair and the segments they cross on.
    [`face_self_intersecting_mask`][triwarp.validation.face_self_intersecting_mask]
        The one-mesh analogue of that.

    Notes
    -----
    !!! note "Witness faces are ambiguous under ties"
        Two parallel plates have a continuum of closest pairs and any of them is a correct answer,
        and ties are the common case rather than the exotic one: whenever the closest approach is
        realised at a *vertex*, every face around that vertex achieves the minimum exactly. The
        tie-break on ``face_a`` is the lowest index; ``face_b`` is **unspecified** among the faces
        attaining it, and which one comes back can differ between runs of the identical input even
        though the squared distance is bit-identical. Compare *distances* against another
        implementation, and faces only where the configuration is generic.
    """
    require_same_device(
        vertices_a=vertices_a, faces_a=faces_a, vertices_b=vertices_b, faces_b=faces_b
    )
    if upper_bound is not None and upper_bound < 0.0:
        raise ValueError(f"upper_bound must be non-negative, got {upper_bound}")
    device = faces_a.device
    n_faces_a = int(faces_a.shape[0]) // 3
    n_faces_b = int(faces_b.shape[0]) // 3
    if n_faces_a == 0 or n_faces_b == 0:
        return math.inf, -1, -1

    # One structure over B, shared by both phases. It is bound to a name for the whole call because
    # the kernels read its BVH by id through ``wp.mesh_get_bvh``, and a collected ``wp.Mesh`` would
    # leave them a dangling id.
    require_nonempty_mesh(faces_b, "mesh_to_mesh_distance")
    mesh_b = wp.Mesh(points=vertices_b, indices=faces_b)

    if upper_bound is None:
        # A vertex-to-surface distance is a distance between the surfaces, so its minimum bounds the
        # answer from above. One readback, and it is what lets the broad phase cull at all.
        #
        # Over a *subsample* of A's vertices, because a minimum over a subset is still an upper
        # bound and this query can otherwise dominate the whole call on a large mesh -- see
        # ``_BOUND_SAMPLE_TARGET``. A looser limit prunes less, so the traversal examines a
        # superset of the candidates it did before and its minimum is the same value.
        stride = max(1, int(vertices_a.shape[0]) // _BOUND_SAMPLE_TARGET)
        probe = vertices_a if stride == 1 else vertices_a[::stride]
        _points, distances, _faces = closest_point_on_mesh(vertices_b, faces_b, probe, mesh=mesh_b)
        upper_bound = float(tw.reduce.min(distances))

    # The per-face AABBs stay: the kernel's box-gap prune reads them, so they are not merely the
    # input to a build. There is no second acceleration structure built over them -- the kernels
    # read the ``wp.Mesh`` built above's own BVH directly with ``wp.mesh_get_bvh``.
    lower = wp.empty(n_faces_b, dtype=wp.vec3, device=device)
    upper = wp.empty(n_faces_b, dtype=wp.vec3, device=device)
    wp.launch(
        kernel_triangles.face_aabb_bounds,
        dim=n_faces_b,
        inputs=[vertices_b, faces_b, lower, upper],
        device=device,
    )
    distance_sq = wp.empty(n_faces_a, dtype=wp.float32, device=device)
    witness = wp.empty(n_faces_a, dtype=wp.int32, device=device)
    # Seeded at the bound the vertex query already paid for, so every thread prunes against it from
    # its first candidate instead of waiting for some other thread to publish one.
    #
    # Seeded at *exactly* ``upper_bound ** 2`` this is wrong, and that is why it used to be ``inf``:
    # the prune skips a candidate whose box gap is ``>=`` the limit, so when the bound *is* the
    # answer -- two spheres whose closest points are vertices -- the very pair achieving it is
    # skipped and the result comes back ``inf``. The relative bump is what keeps that pair, and
    # ``1e-4`` rather than an ulp because ``upper_bound`` comes from ``mesh_query_point_no_sign``,
    # which is documented off by up to 2.1e-5; the margin is ~5x that and ~800 float32 ulps, and it
    # weakens the prune by nothing measurable. ``max`` covers touching meshes, where the bound is
    # ``0`` and any positive limit keeps the exactly-zero-gap pair.
    global_best_sq = wp.full(
        1,
        max(upper_bound * upper_bound * (1.0 + 1e-4), _MIN_POSITIVE_FLOAT32),
        dtype=wp.float32,
        device=device,
    )
    # The broad phase is wildly unbalanced -- most query faces return no candidate at all, while a
    # few carry most of the traversal -- so the walk runs in two passes on CUDA: a thread per face,
    # capped, then a *block* per face that exceeded the cap. See ``_QUERY_CANDIDATE_CAP``. On the
    # cpu device ``wp.launch_tiled`` runs one lane per block, so the second pass would be a serial
    # re-walk; there the cap is disabled and the first pass settles every face instead.
    tiled = prefers_tiled_reduction(device)
    candidate_cap = _QUERY_CANDIDATE_CAP if tiled else INT32_MAX
    overflow = wp.empty(n_faces_a if tiled else 1, dtype=wp.int32, device=device)
    counter = wp.zeros(1, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.face_to_mesh_distance,
        dim=n_faces_a,
        inputs=[
            vertices_a,
            faces_a,
            vertices_b,
            faces_b,
            lower,
            upper,
            mesh_b.id,
            wp.float32(upper_bound),
            wp.int32(candidate_cap),
            global_best_sq,
            distance_sq,
            witness,
            counter,
            overflow,
        ],
        device=device,
    )
    if tiled:
        # One readback, and it is what sizes the second launch. Skipping it by launching
        # ``n_faces_a`` blocks would put an empty block on nearly all of them.
        n_overflow = int(read_scalar(counter, 0))
        if n_overflow > 0:
            wp.launch_tiled(
                kernel_proximity.face_to_mesh_distance_tiled,
                dim=n_overflow,
                inputs=[
                    vertices_a,
                    faces_a,
                    vertices_b,
                    faces_b,
                    lower,
                    upper,
                    mesh_b.id,
                    wp.float32(upper_bound),
                    overflow,
                    global_best_sq,
                    distance_sq,
                    witness,
                ],
                device=device,
                block_dim=_QUERY_TILE_WIDTH,
            )
    keys = wp.empty(n_faces_a, dtype=wp.int64, device=device)
    wp.launch(
        kernel_proximity.face_distance_keys,
        dim=n_faces_a,
        inputs=[distance_sq, keys],
        device=device,
    )
    # The key's low 32 bits are the winning face, so one reduction and one 8-byte read give both the
    # distance and the argmin -- no second pass over the candidates.
    best_face_a = int(tw.reduce.min(keys)) & 0xFFFFFFFF
    # Two 4-byte reads, not two ``.numpy()`` calls: each of those copies the *whole* per-face array
    # to the host to index one element, which matters increasingly on a large mesh where the copy
    # would dominate.
    return (
        math.sqrt(float(read_scalar(distance_sq, best_face_a))),
        best_face_a,
        int(read_scalar(witness, best_face_a)),
    )

normals_at_closest_faces

normals_at_closest_faces(
    mesh: Mesh,
    points: array[vec3],
    *,
    max_dist: float | None = None,
    face_normals: array[vec3] | None = None,
) -> wp.array[wp.vec3]

Return unit face normals at the closest mesh triangle for each query point.

For each position in points, runs an unsigned closest-point query on mesh and returns the normal of the hit triangle.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Warp mesh (BVH built by caller).	required
points	array[vec3]	(m,) query positions as wp.vec3.	required
max_dist	float | None	Maximum search radius per query. When None, derived from the axis-aligned box enclosing mesh vertices and query points.	None
face_normals	array[vec3] | None	Optional length-n_faces unit face normals of mesh (face_normals_and_areas); recomputed when None. This function only gathers from them, so a caller issuing several query batches against one mesh should pass them -- Trimesh.face_normals has them cached.	None

Returns:

Type	Description
array[vec3]	(m,) face normals at the closest triangle for each query. A query with no face within max_dist reports the first face's normal; use closest_point_on_mesh directly, whose triangle_id is -1 there, when a miss has to be detected. A mesh with no faces reports an all-NaN normal for every query, matching closest_point_on_mesh's own zero-face convention.

Raises:

Type	Description
RuntimeError	If mesh, points and face_normals are not all on one device.

See Also

closest_point_on_mesh

Source code in triwarp/proximity.py

def normals_at_closest_faces(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    max_dist: float | None = None,
    face_normals: wp.array[wp.vec3] | None = None,
) -> wp.array[wp.vec3]:
    """
    Return unit face normals at the closest mesh triangle for each query point.

    For each position in ``points``, runs an unsigned closest-point query on
    ``mesh`` and returns the normal of the hit triangle.

    Parameters
    ----------
    mesh
        Warp mesh (BVH built by caller).
    points
        ``(m,)`` query positions as ``wp.vec3``.
    max_dist
        Maximum search radius per query. When ``None``, derived from the
        axis-aligned box enclosing mesh vertices and query points.
    face_normals
        Optional length-``n_faces`` unit face normals of ``mesh``
        ([`face_normals_and_areas`][triwarp.triangles.face_normals_and_areas]); recomputed when
        ``None``. This function only gathers from them, so a caller issuing several query batches
        against one mesh should pass them --
        [`Trimesh.face_normals`][triwarp.mesh.Trimesh.face_normals] has them cached.

    Returns
    -------
    wp.array[wp.vec3]
        ``(m,)`` face normals at the closest triangle for each query. A query with no face
        within ``max_dist`` reports the first face's normal; use
        [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh] directly, whose
        ``triangle_id`` is ``-1`` there, when a miss has to be detected. A mesh with no faces
        reports an all-NaN normal for every query, matching
        [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]'s own zero-face
        convention.

    Raises
    ------
    RuntimeError
        If ``mesh``, ``points`` and ``face_normals`` are not all on one device.

    See Also
    --------
    [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]
    """
    require_same_device(mesh=mesh, points=points, face_normals=face_normals)
    device = points.device
    m = int(points.shape[0])
    if m == 0:
        return wp.empty(0, dtype=wp.vec3, device=device)
    if int(mesh.indices.shape[0]) == 0:
        # closest_point_on_mesh's own n_faces == 0 branch leaves out_face at -1 for every query,
        # which the miss-clamp below would map to face 0 of a face_normals array that has no face
        # 0 -- an out-of-bounds Python-scope gather (array.gather bounds-checks nothing) that
        # segfaults rather than raising. Match closest_point_on_mesh's own convention instead.
        nan = float("nan")
        return wp.full(m, wp.vec3(nan, nan, nan), dtype=wp.vec3, device=device)

    _closest, _distance, out_face = closest_point_on_mesh(
        mesh.points, mesh.indices, points, max_dist=max_dist, mesh=mesh
    )
    # A miss leaves ``-1``, which would gather out of bounds; clamping to face 0 costs one map
    # over ``m`` and keeps the read in range.
    hit_face = wp.empty(m, dtype=wp.int32, device=device)
    wp.map(wp.max, out_face, wp.int32(0), out=hit_face)
    if face_normals is None:
        face_normals, _areas = face_normals_and_areas(mesh.points, mesh.indices)
    return tw.array.gather(face_normals, hit_face)

signed_distance_on_mesh

signed_distance_on_mesh(
    vertices: array[vec3],
    faces: array[int32],
    points: array[vec3],
    *,
    max_dist: float | None = None,
    sign_mode: Literal["parity", "winding"] = "parity",
    n_sample: int = 5,
    perturbation_scale: float = 0.1,
    accuracy: float = 2.0,
    winding_threshold: float = 0.5,
    mesh: Mesh | None = None,
) -> wp.array[wp.float32]

Signed distance from each query point to a triangle mesh (Warp SDF convention).

Distances follow Warp's signed-distance field convention:

• Points outside the mesh have positive distance.
• Points inside have negative distance.
• Points within TOLERANCE_MERGE of the surface return positive unsigned distance.

Trimesh signed_distance uses the opposite sign; negate its output to compare. See also contains_points (inside iff signed distance is negative, except on the on-surface tolerance band).

The unsigned distance is identical in both sign_mode values — only the sign differs.

Choosing a sign_mode

"parity" (default) uses wp.mesh_query_point_sign_parity: it casts n_sample perturbed rays and votes on the crossing parity. Exact on a watertight mesh, cheap, but it has no principled answer on an open or holed surface — a ray that escapes through a hole flips the verdict.

"winding" uses wp.mesh_query_point_sign_winding_number, which evaluates the generalized winding number on the mesh BVH (a Barnes-Hut style traversal governed by accuracy) and compares it against winding_threshold. This is the Jacobson et al. robust inside/outside criterion and it degrades gracefully on non-watertight input, which is why it is the mode to reach for on raw scan data.

The two modes agree on watertight meshes, but on a surface with an open patch "winding" reproduces the exact generalized winding number's sign far more often than "parity" does. The cost is a slower query and a larger wp.Mesh: support_winding_number=True stores a solid-angle expansion per BVH node, using noticeably more device memory.

"winding" is still much cheaper than thresholding winding_number yourself, because that sums the exact solid angle over every face for every query rather than using the BVH's Barnes-Hut traversal, and the gap widens with the face count. Reach for winding_number only when you need the winding value — Warp exposes no builtin for the approximated value, only its sign.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n,) mesh vertex positions as wp.vec3.	required
faces	array[int32]	(f * 3,) flat triangle index array as wp.int32. The internally built wp.Mesh aliases these buffers rather than copying them; do not mutate them for the duration of the call.	required
points	array[vec3]	(m,) query positions in space as wp.vec3.	required
max_dist	float | None	Maximum search radius per query. When None, derived from the axis-aligned box enclosing mesh vertices and query points.	None
sign_mode	Literal['parity', 'winding']	"parity" (default) for ray-parity sign, "winding" for the generalized winding-number sign. See the note above.	'parity'
n_sample	int	Perturbed rays for parity voting (off-triangle sign branch). "parity" only.	5
perturbation_scale	float	Uniform perturbation scale for parity rays. "parity" only.	0.1
accuracy	float	Barnes-Hut accuracy for the winding-number traversal: a node is expanded when the query point is within accuracy times the node's radius, so larger values are more accurate and slower. "winding" only; Warp's default is 2.0.	2.0
winding_threshold	float	Winding number above which a point counts as inside. "winding" only; 0.5 is the standard choice for a once-wound closed surface.	0.5
mesh	Mesh | None	A wp.Mesh already built over vertices and faces, to spare the clone and BVH build. Only valid with sign_mode="parity": the winding mode needs a mesh built with support_winding_number=True, and wp.Mesh exposes no way to read that flag back, so a supplied mesh cannot be checked and is refused rather than silently degraded to parity.	None

Returns:

Type	Description
array[float32]	(m,) signed distances in float32.

Raises:

Type	Description
ValueError	If sign_mode is not "parity" or "winding", or if mesh is supplied together with sign_mode="winding".
RuntimeError	If vertices, faces, points and mesh are not all on one device.

See Also

winding_number contains_points

Source code in triwarp/proximity.py

def signed_distance_on_mesh(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    points: wp.array[wp.vec3],
    *,
    max_dist: float | None = None,
    sign_mode: Literal["parity", "winding"] = "parity",
    n_sample: int = 5,
    perturbation_scale: float = 0.1,
    accuracy: float = 2.0,
    winding_threshold: float = 0.5,
    mesh: wp.Mesh | None = None,
) -> wp.array[wp.float32]:
    """
    Signed distance from each query point to a triangle mesh (Warp SDF convention).

    Distances follow Warp's signed-distance field convention:

    * Points **outside** the mesh have **positive** distance.
    * Points **inside** have **negative** distance.
    * Points within [`TOLERANCE_MERGE`][triwarp.constants.TOLERANCE_MERGE] of the surface
      return positive unsigned distance.

    Trimesh ``signed_distance`` uses the opposite sign; negate its output to compare.
    See also [`contains_points`][triwarp.ray.contains_points] (inside iff signed distance is
    negative, except on the on-surface tolerance band).

    The **unsigned** distance is identical in both ``sign_mode`` values — only the sign differs.

    !!! note "Choosing a `sign_mode`"

        ``"parity"`` (default) uses ``wp.mesh_query_point_sign_parity``: it casts ``n_sample``
        perturbed rays and votes on the crossing parity. Exact on a watertight mesh, cheap, but
        it has no principled answer on an open or holed surface — a ray that escapes through a
        hole flips the verdict.

        ``"winding"`` uses ``wp.mesh_query_point_sign_winding_number``, which evaluates the
        *generalized winding number* on the mesh BVH (a Barnes-Hut style traversal governed by
        ``accuracy``) and compares it against ``winding_threshold``. This is the
        Jacobson et al. robust inside/outside criterion and it degrades gracefully on
        non-watertight input, which is why it is the mode to reach for on raw scan data.

        The two modes agree on watertight meshes, but on a surface with an open patch ``"winding"``
        reproduces the exact generalized winding number's sign far more often than ``"parity"``
        does. The cost is a slower query and a larger ``wp.Mesh``: ``support_winding_number=True``
        stores a solid-angle expansion per BVH node, using noticeably more device memory.

        ``"winding"`` is still much cheaper than thresholding
        [`winding_number`][triwarp.proximity.winding_number] yourself, because that sums the exact
        solid angle over *every* face for *every* query rather than using the BVH's Barnes-Hut
        traversal, and the gap widens with the face count. Reach for
        [`winding_number`][triwarp.proximity.winding_number] only when you need the winding
        *value* — Warp exposes no builtin for the approximated value, only its sign.

    Parameters
    ----------
    vertices
        ``(n,)`` mesh vertex positions as ``wp.vec3``.
    faces
        ``(f * 3,)`` flat triangle index array as ``wp.int32``.
        The internally built ``wp.Mesh`` aliases these buffers rather than copying them;
        do not mutate them for the duration of the call.
    points
        ``(m,)`` query positions in space as ``wp.vec3``.
    max_dist
        Maximum search radius per query. When ``None``, derived from the
        axis-aligned box enclosing mesh vertices and query points.
    sign_mode
        ``"parity"`` (default) for ray-parity sign, ``"winding"`` for the generalized
        winding-number sign. See the note above.
    n_sample
        Perturbed rays for parity voting (off-triangle sign branch). ``"parity"`` only.
    perturbation_scale
        Uniform perturbation scale for parity rays. ``"parity"`` only.
    accuracy
        Barnes-Hut accuracy for the winding-number traversal: a node is expanded when the query
        point is within ``accuracy`` times the node's radius, so larger values are more accurate
        and slower. ``"winding"`` only; Warp's default is ``2.0``.
    winding_threshold
        Winding number above which a point counts as inside. ``"winding"`` only; ``0.5`` is the
        standard choice for a once-wound closed surface.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the clone and BVH
        build. **Only valid with ``sign_mode="parity"``**: the winding mode needs a mesh built with
        ``support_winding_number=True``, and ``wp.Mesh`` exposes no way to read that flag back, so a
        supplied mesh cannot be checked and is refused rather than silently degraded to parity.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` signed distances in ``float32``.

    Raises
    ------
    ValueError
        If ``sign_mode`` is not ``"parity"`` or ``"winding"``, or if ``mesh`` is supplied together
        with ``sign_mode="winding"``.
    RuntimeError
        If ``vertices``, ``faces``, ``points`` and ``mesh`` are not all on one device.

    See Also
    --------
    [`winding_number`][triwarp.proximity.winding_number]
    [`contains_points`][triwarp.ray.contains_points]
    """
    require_same_device(vertices=vertices, faces=faces, points=points, mesh=mesh)
    if sign_mode not in ("parity", "winding"):
        raise ValueError(f"sign_mode must be 'parity' or 'winding', got {sign_mode!r}")

    device = vertices.device
    m = int(points.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if m == 0:
        return wp.empty(0, dtype=wp.float32, device=device)
    if n_faces == 0:
        return wp.full(m, float("inf"), dtype=wp.float32, device=device)

    if mesh is not None and sign_mode == "winding":
        # wp.Mesh exposes no way to read back support_winding_number, so a supplied mesh cannot be
        # checked for the per-node solid-angle expansion the winding builtin needs -- and without it
        # the builtin silently degrades to ray parity. Refusing is the only safe answer; the parity
        # mode has no such requirement and accepts any mesh.
        raise ValueError(
            "sign_mode='winding' cannot use a supplied mesh: it needs "
            "wp.Mesh(support_winding_number=True), which cannot be verified after construction. "
            "Omit mesh=, or use sign_mode='parity'."
        )
    if mesh is None:
        require_nonempty_mesh(faces, "signed_distance_on_mesh")
        # The winding-number builtin silently degrades to ray parity unless the mesh carries the
        # per-node solid-angle expansion, so the flag is bound to sign_mode here rather than
        # exposed.
        # The mesh aliases the caller's buffers and is discarded here, so it needs no copy.
        mesh = wp.Mesh(
            points=vertices, indices=faces, support_winding_number=sign_mode == "winding"
        )
    if max_dist is None:
        max_dist = tw.bounds.enclosing_diagonal(mesh.points, points)
    out_distance = wp.empty(m, dtype=wp.float32, device=device)
    if sign_mode == "winding":
        wp.launch(
            kernel_proximity.signed_distance_on_mesh_winding,
            dim=m,
            inputs=[
                mesh.id,
                points,
                wp.float32(max_dist),
                wp.float32(accuracy),
                wp.float32(winding_threshold),
                out_distance,
            ],
            device=device,
        )
        return out_distance
    wp.launch(
        kernel_proximity.signed_distance_on_mesh,
        dim=m,
        inputs=[
            mesh.id,
            points,
            wp.float32(max_dist),
            wp.int32(n_sample),
            wp.float32(perturbation_scale),
            out_distance,
        ],
        device=device,
    )
    return out_distance

signed_distance_grid

signed_distance_grid(
    vertices: array[vec3],
    faces: array[int32],
    voxel_size: float | None = None,
    *,
    bounds: tuple[vec3, vec3] | None = None,
    pad: int = 2,
    sign_mode: Literal["parity", "winding"] = "parity",
    mesh: Mesh | None = None,
) -> tuple[twt.Array3dFloat32, tuple[wp.vec3, wp.vec3]]

Sample the signed distance to a mesh on a regular lattice, as a field and the box it spans.

The bridge from a surface to a level set, and the missing half of the implicit round trip: triwarp.voxels.to_field already gives triwarp.levelset.marching_cubes an occupancy lattice, but occupancy is 0 or 1 and thresholding it at anything other than 0.5 does not move the surface anywhere. A distance field does, which is what makes marching_cubes(*signed_distance_grid(...), iso=d) an offset surface at distance d -- see triwarp.levelset.offset_mesh, the named entry point for it.

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	Lattice spacing, isotropic. None takes triwarp.voxels.resolve_voxel_grid's default of 1 % of the bounding-box diagonal, which is this package's one definition of an unspecified grid.	None
bounds	tuple[vec3, vec3] | None	(lower, upper) box to sample, before padding. None uses the mesh's own axis-aligned box, which is what an offset wants -- an inward offset needs no more, and an outward one needs pad to cover it.	None
pad	int	Cells of margin added on every side, so the lattice extends pad * voxel_size beyond the box. Two is enough for the surface itself to be enclosed; an outward offset of d needs pad >= d / voxel_size + 1 or its level set is clipped by the lattice boundary.	2
sign_mode	Literal['parity', 'winding']	How the sign is decided, forwarded to signed_distance_on_mesh: "parity" counts ray crossings, "winding" sums solid angles and is the one that survives a mesh with open rims.	'parity'
mesh	Mesh | None	A wp.Mesh already built over vertices and faces, to spare the build. Forwarded as-is, so signed_distance_on_mesh's rule applies unchanged: it is usable with sign_mode="parity" only, since the winding sign needs a mesh built with support_winding_number=True and that cannot be verified after construction.	None

Returns:

Name	Type	Description
field	Array3dFloat32	(nx, ny, nz) signed distances, negative inside. Exactly the first argument triwarp.levelset.marching_cubes takes.
bounds	tuple[vec3, vec3]	The (lower, upper) corners the lattice actually spans, padded and snapped so the spacing is exactly voxel_size on every axis. Pass it straight through as that function's bounds.

Raises:

Type	Description
ValueError	If voxel_size is not positive, pad is negative, or faces is empty.
RuntimeError	If vertices, faces and mesh are not all on one device.

Examples:

field, box = tw.proximity.signed_distance_grid(v, f, voxel_size=0.05, pad=4)
shell_v, shell_f = tw.levelset.marching_cubes(field, 0.1, bounds=box)

Notes

The whole lattice is sampled, so size it deliberately: the field costs 4 * nx * ny * nz bytes and one closest-point query per sample, which is 16.7 M queries and 67 MB at 256 ** 3. There is no narrow band, and that is deliberate rather than missing -- signed_distance_on_mesh reports +max_dist for a query that finds no face within the limit, so a banded field would carry a positive value deep inside the solid and silently invert the level set. Reduce the resolution instead.

The lattice is a corner lattice: field[0, 0, 0] sits exactly on the returned lower. That is triwarp.voxels.grid_points's convention and marching_cubes's, and it is not the voxel-centre convention the rest of triwarp.voxels uses.

See Also

signed_distance_on_mesh The per-query form this samples, and where the sign conventions are documented. triwarp.levelset.offset_mesh What to call instead when the answer wanted is the offset surface rather than the field. triwarp.voxels.to_field The occupancy lattice, when a binary inside test is all that is needed.

Source code in triwarp/proximity.py

def signed_distance_grid(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    voxel_size: float | None = None,
    *,
    bounds: tuple[wp.vec3, wp.vec3] | None = None,
    pad: int = 2,
    sign_mode: Literal["parity", "winding"] = "parity",
    mesh: wp.Mesh | None = None,
) -> tuple[twt.Array3dFloat32, tuple[wp.vec3, wp.vec3]]:
    """
    Sample the signed distance to a mesh on a regular lattice, as a field and the box it spans.

    The bridge from a surface to a **level set**, and the missing half of the implicit round trip:
    [`triwarp.voxels.to_field`][triwarp.voxels.to_field] already gives
    [`triwarp.levelset.marching_cubes`][triwarp.levelset.marching_cubes] an
    *occupancy* lattice, but occupancy is ``0`` or ``1`` and thresholding it at anything other than
    ``0.5`` does not move the surface anywhere. A distance field does, which is what makes
    ``marching_cubes(*signed_distance_grid(...), iso=d)`` an offset surface at distance ``d`` --
    see [`triwarp.levelset.offset_mesh`][triwarp.levelset.offset_mesh], the named entry point for
    it.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    voxel_size
        Lattice spacing, isotropic. ``None`` takes
        [`triwarp.voxels.resolve_voxel_grid`][triwarp.voxels.resolve_voxel_grid]'s default of 1 % of
        the bounding-box diagonal, which is this package's one definition of an unspecified grid.
    bounds
        ``(lower, upper)`` box to sample, **before** padding. ``None`` uses the mesh's own
        axis-aligned box, which is what an offset wants -- an inward offset needs no more, and an
        outward one needs ``pad`` to cover it.
    pad
        Cells of margin added on every side, so the lattice extends ``pad * voxel_size`` beyond the
        box. Two is enough for the surface itself to be enclosed; an **outward offset of ``d``
        needs ``pad >= d / voxel_size + 1``** or its level set is clipped by the lattice boundary.
    sign_mode
        How the sign is decided, forwarded to
        [`signed_distance_on_mesh`][triwarp.proximity.signed_distance_on_mesh]: ``"parity"`` counts
        ray crossings, ``"winding"`` sums solid angles and is the one that survives a mesh with open
        rims.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the build. Forwarded
        as-is, so ``signed_distance_on_mesh``'s rule applies unchanged: it is usable with
        ``sign_mode="parity"`` only, since the winding sign needs a mesh built with
        ``support_winding_number=True`` and that cannot be verified after construction.

    Returns
    -------
    field : twt.Array3dFloat32
        ``(nx, ny, nz)`` signed distances, negative inside. Exactly the first argument
        [`triwarp.levelset.marching_cubes`][triwarp.levelset.marching_cubes] takes.
    bounds : tuple[wp.vec3, wp.vec3]
        The ``(lower, upper)`` corners the lattice actually spans, padded and snapped so the spacing
        is exactly ``voxel_size`` on every axis. Pass it straight through as that function's
        ``bounds``.

    Raises
    ------
    ValueError
        If ``voxel_size`` is not positive, ``pad`` is negative, or ``faces`` is empty.
    RuntimeError
        If ``vertices``, ``faces`` and ``mesh`` are not all on one device.

    Examples
    --------
    ```python
    field, box = tw.proximity.signed_distance_grid(v, f, voxel_size=0.05, pad=4)
    shell_v, shell_f = tw.levelset.marching_cubes(field, 0.1, bounds=box)
    ```

    Notes
    -----
    **The whole lattice is sampled, so size it deliberately**: the field costs
    ``4 * nx * ny * nz`` bytes and one closest-point query per sample, which is 16.7 M queries and
    67 MB at ``256 ** 3``. There is no narrow band, and that is deliberate rather than missing --
    ``signed_distance_on_mesh`` reports ``+max_dist`` for a query that finds no face within the
    limit, so a banded field would carry a **positive** value deep inside the solid and silently
    invert the level set. Reduce the resolution instead.

    The lattice is a *corner* lattice: ``field[0, 0, 0]`` sits exactly on the returned ``lower``.
    That is [`triwarp.voxels.grid_points`][triwarp.voxels.grid_points]'s convention and
    ``marching_cubes``'s, and it is **not** the voxel-centre convention the rest of
    [`triwarp.voxels`][triwarp.voxels] uses.

    See Also
    --------
    [`signed_distance_on_mesh`][triwarp.proximity.signed_distance_on_mesh]
        The per-query form this samples, and where the sign conventions are documented.
    [`triwarp.levelset.offset_mesh`][triwarp.levelset.offset_mesh]
        What to call instead when the answer wanted is the offset surface rather than the field.
    [`triwarp.voxels.to_field`][triwarp.voxels.to_field]
        The occupancy lattice, when a binary inside test is all that is needed.
    """
    require_same_device(vertices=vertices, faces=faces, mesh=mesh)
    if pad < 0:
        raise ValueError("pad must be non-negative")
    if int(faces.shape[0]) == 0:
        raise ValueError("signed_distance_grid needs at least one face")
    device = vertices.device
    spacing, _origin = tw.voxels.resolve_voxel_grid(
        vertices, voxel_size, None, caller="signed_distance_grid"
    )

    lower, upper = bounds if bounds is not None else tw.bounds.aabb(vertices)
    margin = float(pad) * spacing
    lower = wp.vec3(lower[0] - margin, lower[1] - margin, lower[2] - margin)
    # One sample per spacing, and at least the two marching cubes needs to have a cell at all.
    shape = tuple(
        max(2, math.floor((float(upper[axis]) + margin - float(lower[axis])) / spacing) + 1)
        for axis in range(3)
    )
    snapped_upper = wp.vec3(
        *(float(lower[axis]) + (shape[axis] - 1) * spacing for axis in range(3))
    )

    samples = tw.voxels.grid_points(shape, bounds=(lower, snapped_upper), device=device)
    distances = signed_distance_on_mesh(vertices, faces, samples, sign_mode=sign_mode, mesh=mesh)
    return twt.as_array3d(distances.reshape(shape), wp.float32), (lower, snapped_upper)

winding_number

winding_number(
    vertices: array[vec3],
    faces: array[int32],
    points: array[vec3],
    *,
    tiled: bool = True,
) -> wp.array[wp.float32]

Generalized winding number at each query point.

Sums the signed solid angle subtended by each oriented triangle. For a closed, consistently oriented watertight mesh, interior points have winding number near 1 and exterior points near 0.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n,) mesh vertex positions as wp.vec3.	required
faces	array[int32]	(f * 3,) flat triangle index array as wp.int32.	required
points	array[vec3]	(m,) query positions in space as wp.vec3.	required
tiled	bool	When True (default), sum solid angles with the face list partitioned across threads: one thread per (query, face slice) walks a strided slice of ITEMS_PER_QUERY_SLICE faces and accumulates one wp.atomic_add per slice, so the summation order is nondeterministic and the result can differ in the last float32 digits between runs. When False, each query thread loops over all faces serially — orders of magnitude slower on large meshes, but the fixed left-to-right summation makes it the exact-sum reference.	True

Returns:

Type	Description
array[float32]	(m,) winding numbers in float32.

Raises:

Type	Description
RuntimeError	If vertices, faces and points are not all on one device.

Source code in triwarp/proximity.py

def winding_number(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    points: wp.array[wp.vec3],
    *,
    tiled: bool = True,
) -> wp.array[wp.float32]:
    """
    Generalized winding number at each query point.

    Sums the signed solid angle subtended by each oriented triangle. For a
    closed, consistently oriented watertight mesh, interior points have
    winding number near ``1`` and exterior points near ``0``.

    Parameters
    ----------
    vertices
        ``(n,)`` mesh vertex positions as ``wp.vec3``.
    faces
        ``(f * 3,)`` flat triangle index array as ``wp.int32``.
    points
        ``(m,)`` query positions in space as ``wp.vec3``.
    tiled
        When ``True`` (default), sum solid angles with the face list partitioned across
        threads: one thread per ``(query, face slice)`` walks a strided slice of
        [`ITEMS_PER_QUERY_SLICE`][triwarp.proximity.ITEMS_PER_QUERY_SLICE] faces and accumulates
        one ``wp.atomic_add`` per slice, so the summation order is nondeterministic and the result
        can differ in the last float32 digits between runs. When ``False``, each query thread
        loops over all faces serially — orders of magnitude slower on large meshes,
        but the fixed left-to-right summation makes it the exact-sum reference.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` winding numbers in ``float32``.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces`` and ``points`` are not all on one device.
    """
    require_same_device(vertices=vertices, faces=faces, points=points)
    device = points.device
    n_queries = int(points.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if n_queries == 0:
        return wp.empty(0, dtype=wp.float32, device=device)
    if n_faces == 0:
        return wp.zeros(n_queries, dtype=wp.float32, device=device)

    out_winding = (
        wp.zeros(n_queries, dtype=wp.float32, device=device)
        if tiled
        else wp.empty(n_queries, dtype=wp.float32, device=device)
    )
    if tiled:
        n_face_slices = max(1, (n_faces + ITEMS_PER_QUERY_SLICE - 1) // ITEMS_PER_QUERY_SLICE)
        wp.launch(
            kernel_proximity.winding_number_tiled,
            dim=(n_queries, n_face_slices),
            inputs=[
                vertices,
                faces,
                wp.int32(n_faces),
                wp.int32(n_face_slices),
                points,
                out_winding,
            ],
            device=device,
        )
    else:
        wp.launch(
            kernel_proximity.winding_number,
            dim=n_queries,
            inputs=[vertices, faces, wp.int32(n_faces), points, out_winding],
            device=device,
        )
    return out_winding

query_mesh_aabb_with_offsets

query_mesh_aabb_with_offsets(
    mesh: Mesh,
    query_lower: array[vec3],
    query_upper: array[vec3],
    *,
    max_hits: int = 16,
) -> tuple[
    wp.array[wp.int32],
    wp.array[wp.int32],
    wp.array[wp.int32],
]

Low-level mesh AABB query with per-query axis-aligned bounds.

For each query primitive k, tests intersection of [query_lower[k], query_upper[k]] against every triangle in mesh via wp.mesh_query_aabb. At most max_hits candidate face indices are recorded per query.

Requires the default Warp mesh BVH backend; bvh_constructor="cubql" meshes do not support AABB queries.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Target warp.Mesh built with the default BVH backend.	required
query_lower	array[vec3]	Length-m lower corners of the query boxes, on the target device.	required
query_upper	array[vec3]	Length-m upper corners of the query boxes.	required
max_hits	int	Maximum candidate faces recorded per query. Hits past this cap are dropped, so the result is a bounded sample rather than the full candidate set when a query straddles more than max_hits triangles.	16

Returns:

Type	Description
(candidate_indices_flat, offsets, hit_counts)	offsets is the exclusive prefix sum of per-query hit counts. Query k owns candidate_indices_flat[offsets[k] : offsets[k] + hit_counts[k]]. All three are empty when m == 0.

Raises:

Type	Description
ValueError	If query_lower and query_upper have different lengths, or max_hits < 1.
RuntimeError	If mesh, query_lower and query_upper are not all on one device.

Source code in triwarp/proximity.py

def query_mesh_aabb_with_offsets(
    mesh: wp.Mesh,
    query_lower: wp.array[wp.vec3],
    query_upper: wp.array[wp.vec3],
    *,
    max_hits: int = 16,
) -> tuple[wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.int32]]:
    """
    Low-level mesh AABB query with per-query axis-aligned bounds.

    For each query primitive ``k``, tests intersection of ``[query_lower[k],
    query_upper[k]]`` against every triangle in ``mesh`` via ``wp.mesh_query_aabb``.
    At most ``max_hits`` candidate face indices are recorded per query.

    Requires the default Warp mesh BVH backend; ``bvh_constructor="cubql"`` meshes
    do not support AABB queries.

    Parameters
    ----------
    mesh
        Target ``warp.Mesh`` built with the default BVH backend.
    query_lower
        Length-``m`` lower corners of the query boxes, on the target device.
    query_upper
        Length-``m`` upper corners of the query boxes.
    max_hits
        Maximum candidate faces recorded per query. Hits past this cap are dropped, so the
        result is a bounded sample rather than the full candidate set when a query straddles
        more than ``max_hits`` triangles.

    Returns
    -------
    candidate_indices_flat, offsets, hit_counts
        ``offsets`` is the exclusive prefix sum of per-query hit counts.
        Query ``k`` owns ``candidate_indices_flat[offsets[k] : offsets[k] + hit_counts[k]]``.
        All three are empty when ``m == 0``.

    Raises
    ------
    ValueError
        If ``query_lower`` and ``query_upper`` have different lengths, or ``max_hits < 1``.
    RuntimeError
        If ``mesh``, ``query_lower`` and ``query_upper`` are not all on one device.
    """
    require_same_device(mesh=mesh, query_lower=query_lower, query_upper=query_upper)
    device = query_lower.device
    m = int(query_lower.shape[0])
    if int(query_upper.shape[0]) != m:
        raise ValueError("query_lower and query_upper must have the same length")
    if max_hits < 1:
        raise ValueError("max_hits must be >= 1")

    if m == 0:
        return (
            wp.empty(0, dtype=wp.int32, device=device),
            wp.empty(0, dtype=wp.int32, device=device),
            wp.empty(0, dtype=wp.int32, device=device),
        )

    hit_counts = wp.empty(m, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.query_mesh_aabb_count,
        dim=m,
        inputs=[query_lower, query_upper, mesh.id, wp.int32(max_hits), hit_counts],
        device=device,
    )

    # One scan pass yields both the row starts and their total; an all-zero ``hit_counts`` scans to
    # all-zero offsets, which is exactly what the empty case wants to return.
    offsets, total_hits = tw.array.counts_to_offsets(hit_counts)
    if total_hits == 0:
        return wp.empty(0, dtype=wp.int32, device=device), offsets, hit_counts

    candidate_indices_flat = wp.empty(total_hits, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.query_mesh_aabb_neighbors,
        dim=m,
        inputs=[
            query_lower,
            query_upper,
            mesh.id,
            wp.int32(max_hits),
            offsets,
            candidate_indices_flat,
        ],
        device=device,
    )

    return candidate_indices_flat, offsets, hit_counts

containing_faces_2d

containing_faces_2d(
    vertices: array[vec2],
    faces: array[int32],
    points: array[vec2],
) -> wp.array[wp.int32]

For each 2D query point, the triangle of a planar triangulation that contains it.

Point location in the plane -- the primitive an inverse UV lookup needs, and the counterpart of triwarp.texture's forward direction: remap_attribute_from_uv samples an image at a UV coordinate, where this answers which triangle of a UV atlas a coordinate falls in, and so which surface point it corresponds to. Pair it with points_to_barycentric on the returned face to finish the inverse map.

The triangulation must not overlap itself -- a UV atlas, a Delaunay triangulation, or the xy projection of a height field. Where triangles do overlap, each query gets one of the containing faces and which one is unspecified.

Parameters:

Name	Type	Description	Default
vertices	array[vec2]	(n_vertices,) planar vertex positions as wp.vec2.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
points	array[vec2]	(m,) planar query positions as wp.vec2.	required

Returns:

Type	Description
array[int32]	Length m on vertices.device: the containing triangle's index, or -1 where the query lies outside the triangulation. A query exactly on a shared edge is inside both its triangles and which one is returned is not specified.

Raises:

Type	Description
RuntimeError	If vertices, faces and points are not all on one device.

Notes

Two stages: a closest-point query against the triangulation lifted to the z = 0 plane picks a candidate face -- sufficient because a point inside any triangle is at distance zero from it, so the closest triangle contains it whenever one does -- and a barycentric sign test on the query's own coordinates decides. That reuses the BVH wp.Mesh already builds, at the cost of one lifted wp.vec3 copy of vertices.

The second stage is not redundant: deciding on the query radius alone misclassifies some queries, because an in-plane point's closest-point distance is not exactly zero in float32 and no radius threshold separates every interior point from every exterior one. The barycentric test is much sharper, so the radius is only a search bound and there is no tolerance to tune.

Building that BVH happens on every call, so a caller locating several point sets against the same triangulation pays for it each time -- there is no prebuilt-index entry point.

See Also

closest_point_on_mesh points_to_barycentric remap_attribute_from_uv

Source code in triwarp/proximity.py

def containing_faces_2d(
    vertices: wp.array[wp.vec2], faces: wp.array[wp.int32], points: wp.array[wp.vec2]
) -> wp.array[wp.int32]:
    """
    For each 2D query point, the triangle of a planar triangulation that contains it.

    Point location in the plane -- the primitive an inverse UV lookup needs, and the counterpart of
    [`triwarp.texture`][]'s forward direction: ``remap_attribute_from_uv`` *samples an image* at a
    UV coordinate, where this answers which triangle of a UV atlas a coordinate falls in, and so
    which surface point it corresponds to. Pair it with
    [`points_to_barycentric`][triwarp.triangles.points_to_barycentric] on the returned face to
    finish the inverse map.

    The triangulation must not overlap itself -- a UV atlas, a Delaunay triangulation, or the
    ``xy`` projection of a height field. Where triangles do overlap, each query gets one of the
    containing faces and which one is unspecified.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` planar vertex positions as ``wp.vec2``.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    points
        ``(m,)`` planar query positions as ``wp.vec2``.

    Returns
    -------
    wp.array[wp.int32]
        Length ``m`` on ``vertices.device``: the containing triangle's index, or ``-1`` where the
        query lies outside the triangulation. A query exactly on a shared edge is inside *both* its
        triangles and which one is returned is not specified.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces`` and ``points`` are not all on one device.

    Notes
    -----
    Two stages: a closest-point query against the triangulation lifted to the ``z = 0`` plane picks
    a candidate face -- sufficient because a point inside any triangle is at distance zero from it,
    so the *closest* triangle contains it whenever one does -- and a barycentric sign test on the
    query's own coordinates decides. That reuses the BVH ``wp.Mesh`` already builds, at the cost of
    one lifted ``wp.vec3`` copy of ``vertices``.

    The second stage is not redundant: deciding on the query radius alone misclassifies some
    queries, because an in-plane point's closest-point distance is not exactly zero in
    ``float32`` and no radius threshold separates every interior point from every exterior one.
    The barycentric test is much sharper, so the radius is only a search bound and there is no
    tolerance to tune.

    Building that BVH happens on every call, so a caller locating several point sets against the
    same triangulation pays for it each time -- there is no prebuilt-index entry point.

    See Also
    --------
    [`closest_point_on_mesh`][triwarp.proximity.closest_point_on_mesh]
    [`points_to_barycentric`][triwarp.triangles.points_to_barycentric]
    [`remap_attribute_from_uv`][triwarp.texture.remap_attribute_from_uv]
    """
    require_same_device(vertices=vertices, faces=faces, points=points)
    device = vertices.device
    m = int(points.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if m == 0:
        return wp.empty(0, dtype=wp.int32, device=device)
    if n_faces == 0:
        return wp.full(m, -1, dtype=wp.int32, device=device)

    lifted = wp.empty(int(vertices.shape[0]), dtype=wp.vec3, device=device)
    wp.map(kernel_array.lift_vec2, vertices, wp.float32(0.0), out=lifted)
    # One readback, the same one `closest_point_on_mesh` pays and for the same reason: the search
    # radius has to be in the triangulation's own units and nothing else knows its scale.
    search_radius = _CONTAINMENT_SEARCH_SCALE * tw.bounds.enclosing_diagonal(lifted)

    require_nonempty_mesh(faces, "containing_faces_2d")
    # Both buffers are local to this call (``lifted`` is built above), so no copy is needed.
    mesh = wp.Mesh(points=lifted, indices=faces)
    out_face = wp.empty(m, dtype=wp.int32, device=device)
    wp.launch(
        kernel_proximity.face_containing_point_2d,
        dim=m,
        inputs=[
            mesh.id,
            vertices,
            faces,
            points,
            wp.float32(search_radius),
            _CONTAINMENT_BARYCENTRIC_EPS,
            out_face,
        ],
        device=device,
    )
    return out_face