visibility

triwarp.visibility

How far the surface is from a point.

Outward until nothing is hit (ambient_occlusion, volumetric_obscurance); inward until it is hit again (shape_diameter, thickness); or in every direction at once, which is the radius at which the surface first appears (max_tangent_sphere).

All five take (mesh: wp.Mesh, points, *, normals=None, ...) -- a prebuilt BVH and a set of positions to measure at -- which is what separates this module from triwarp.proximity, where the queries take raw (vertices, faces) and ask where the surface is rather than how far away.

The two outward fields integrate the same bundle of rays over the outward hemisphere at each point, weighted by Lambert's cosine law, and differ only in what a blocked ray costs: ambient occlusion charges a hit its full weight however far away it is, obscurance discounts it by exp(-tau * distance) so that only nearby geometry darkens a point. Ambient occlusion is the tau -> 0 limit, which is why they share a kernel. Neither is normalized against a scene, so both are comparable across meshes and resolutions, and on a convex closed surface no ray can return and every point reads exactly 0 -- the cheapest available sanity check on a result.

The three inward measures differ in how much evidence they take. thickness is a dispatcher over the other two: one inward ray, or one tangent sphere. shape_diameter fires a whole cone and takes an outlier-trimmed mean, which is what makes it stable enough for skeleton extraction and part segmentation. max_tangent_sphere uses no rays at all: it shrinks a sphere until nothing but the surface touches it, so it answers the question for a volume rather than along a direction, and it is the reason this module's name is a loose fit for one of its five members.

RayWeight module-attribute

RayWeight = Literal['cosine', 'uniform']

Weighting of a ray; see ambient_occlusion.

ambient_occlusion

ambient_occlusion(
    mesh: Mesh,
    points: array[vec3],
    *,
    normals: array[vec3] | None = None,
    n_rays: int = 64,
    weight: RayWeight = "cosine",
    max_t: float | None = None,
) -> wp.array[wp.float32]

Fraction of the outward hemisphere at each point that is blocked by the mesh itself.

A Fibonacci hemisphere lattice is rotated into each point's tangent frame and every direction is traced against the mesh; the result is the weighted share of directions that hit something. This is MeshLab's compute_scalar_ambient_occlusion and libigl's igl::ambient_occlusion, which differ from each other exactly in the weight argument below.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangle mesh with a built BVH (wp.Mesh). The occluder and the surface being shaded.	required
points	array[vec3]	(m,) positions to shade, normally the mesh's own vertices.	required
normals	array[vec3] | None	(m,) outward unit normals, defining which hemisphere is "outward" at each point. When None they are taken from the closest face of mesh (normals_at_closest_faces), which is right for points on the surface and meaningless for points far off it — pass them explicitly in that case. For a smooth result on the mesh's own vertices, pass vertex_normals at weighting="area" instead: face normals make the field piecewise constant across each vertex's ring.	None
n_rays	int	Directions per point. Error falls as 1 / sqrt(n_rays); MeshLab's default is 64, which is this one. Cost is exactly linear in it.	64
weight	RayWeight	How a direction is weighted in the integral: "cosine" (default) — by dot(direction, normal), which is the physically correct weighting for ambient irradiance and MeshLab's. Grazing directions barely matter. "uniform" — every direction counts once, which is libigl's convention and reads as a solid-angle fraction rather than an irradiance one.	'cosine'
max_t	float | None	Maximum ray length. Anything beyond it does not occlude. When None, the diagonal of the AABB enclosing the mesh and the query points, so nothing in the mesh is missed.	None

Returns:

Type	Description
array[float32]	(m,) occlusion in [0, 1] on points.device: 0 where the hemisphere is fully open and 1 where it is fully blocked. Exactly 0 everywhere on a convex closed mesh, which is the cheapest available sanity check on a result.

Raises:

Type	Description
ValueError	If n_rays < 1, weight is not one of the two names, or normals has a different length from points.
RuntimeError	If mesh, points and normals are not all on one device.

See Also

volumetric_obscurance triwarp.visibility.shape_diameter triwarp.sample.sample_fibonacci_hemisphere

Notes

MeshLab reports the complement of this (higher is more exposed) and leaves it unnormalized — its number is roughly (1 - occlusion) * n_rays / 4 for a cosine-weighted bundle over n_rays whole-sphere directions. Its direction set is also its own, so the two agree in distribution and in ranking rather than value by value.

Source code in triwarp/visibility.py

def ambient_occlusion(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    normals: wp.array[wp.vec3] | None = None,
    n_rays: int = 64,
    weight: RayWeight = "cosine",
    max_t: float | None = None,
) -> wp.array[wp.float32]:
    """
    Fraction of the outward hemisphere at each point that is blocked by the mesh itself.

    A Fibonacci hemisphere lattice is rotated into each point's tangent frame and every direction is
    traced against the mesh; the result is the weighted share of directions that hit something. This
    is MeshLab's ``compute_scalar_ambient_occlusion`` and libigl's ``igl::ambient_occlusion``, which
    differ from each other exactly in the ``weight`` argument below.

    Parameters
    ----------
    mesh
        Triangle mesh with a built BVH (``wp.Mesh``). The occluder *and* the surface being shaded.
    points
        ``(m,)`` positions to shade, normally the mesh's own vertices.
    normals
        ``(m,)`` outward unit normals, defining which hemisphere is "outward" at each point. When
        ``None`` they are taken from the closest face of ``mesh``
        ([`normals_at_closest_faces`][triwarp.proximity.normals_at_closest_faces]), which is right
        for points on the surface and meaningless for points far off it — pass them explicitly in
        that case. For a smooth result on the mesh's own vertices, pass
        [`vertex_normals`][triwarp.vertices.vertex_normals] at ``weighting="area"`` instead:
        face normals make the field piecewise constant across each vertex's ring.
    n_rays
        Directions per point. Error falls as ``1 / sqrt(n_rays)``; MeshLab's default is ``64``,
        which is this one. Cost is exactly linear in it.
    weight
        How a direction is weighted in the integral:

        - ``"cosine"`` (default) — by ``dot(direction, normal)``, which is the physically correct
          weighting for ambient irradiance and MeshLab's. Grazing directions barely matter.
        - ``"uniform"`` — every direction counts once, which is libigl's convention and reads as a
          solid-angle fraction rather than an irradiance one.
    max_t
        Maximum ray length. Anything beyond it does not occlude. When ``None``, the diagonal of the
        AABB enclosing the mesh and the query points, so nothing in the mesh is missed.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` occlusion in ``[0, 1]`` on ``points.device``: ``0`` where the hemisphere is fully
        open and ``1`` where it is fully blocked. **Exactly ``0`` everywhere on a convex closed
        mesh**, which is the cheapest available sanity check on a result.

    Raises
    ------
    ValueError
        If ``n_rays < 1``, ``weight`` is not one of the two names, or ``normals`` has a different
        length from ``points``.
    RuntimeError
        If ``mesh``, ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`volumetric_obscurance`][triwarp.visibility.volumetric_obscurance]
    [`triwarp.visibility.shape_diameter`][triwarp.visibility.shape_diameter]
    [`triwarp.sample.sample_fibonacci_hemisphere`][triwarp.sample.sample_fibonacci_hemisphere]

    Notes
    -----
    MeshLab reports the *complement* of this (higher is more exposed) and leaves it unnormalized —
    its number is roughly ``(1 - occlusion) * n_rays / 4`` for a cosine-weighted bundle over
    ``n_rays`` whole-sphere directions. Its direction set is also its own, so the two agree in
    distribution and in ranking rather than value by value.
    """
    require_same_device(mesh=mesh, points=points, normals=normals)
    return _occlusion_bundle(mesh, points, normals, n_rays, weight, 0.0, max_t, "ambient_occlusion")

volumetric_obscurance

volumetric_obscurance(
    mesh: Mesh,
    points: array[vec3],
    *,
    normals: array[vec3] | None = None,
    n_rays: int = 64,
    tau: float = 0.1,
    weight: RayWeight = "cosine",
    max_t: float | None = None,
) -> wp.array[wp.float32]

Distance-attenuated ambient occlusion: an occluder at range t counts exp(-tau * t).

Iones et al.'s obscurance, and MeshLab's compute_scalar_by_volumetric_obscurance. It exists because binary ambient_occlusion treats a wall across the room like a crevice wall a millimetre away, which darkens the interior of any closed room uniformly and hides exactly the small-scale detail the field is usually wanted for. Attenuating by distance keeps the response local.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangle mesh with a built BVH (wp.Mesh).	required
points	array[vec3]	(m,) positions to shade.	required
normals	array[vec3] | None	(m,) outward unit normals; see ambient_occlusion for the default.	None
n_rays	int	Directions per point.	64
tau	float	Attenuation rate, in inverse length units of the mesh — so it is not scale-invariant, and a mesh scaled by k wants tau / k for the same result. MeshLab's default is 0.1, which suits a mesh of extent order 1. As tau -> 0 this becomes ambient_occlusion; as tau -> inf everything reads 0. Must be positive.	0.1
weight	RayWeight	Ray weighting, as in ambient_occlusion.	'cosine'
max_t	float | None	Maximum ray length.	None

Returns:

Type	Description
array[float32]	(m,) obscurance-weighted occlusion in [0, 1] on points.device, 0 where nothing nearby blocks. Exactly 0 on a convex closed mesh.

Raises:

Type	Description
ValueError	If tau <= 0, n_rays < 1, weight is unknown, or normals is the wrong length.
RuntimeError	If mesh, points and normals are not all on one device.

See Also

ambient_occlusion

Source code in triwarp/visibility.py

def volumetric_obscurance(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    normals: wp.array[wp.vec3] | None = None,
    n_rays: int = 64,
    tau: float = 0.1,
    weight: RayWeight = "cosine",
    max_t: float | None = None,
) -> wp.array[wp.float32]:
    """
    Distance-attenuated ambient occlusion: an occluder at range ``t`` counts ``exp(-tau * t)``.

    Iones et al.'s obscurance, and MeshLab's ``compute_scalar_by_volumetric_obscurance``. It exists
    because binary [`ambient_occlusion`][triwarp.visibility.ambient_occlusion] treats a wall across
    the room like a crevice wall a millimetre away, which darkens the interior of any closed room
    uniformly and hides exactly the small-scale detail the field is usually wanted for. Attenuating
    by distance keeps the response local.

    Parameters
    ----------
    mesh
        Triangle mesh with a built BVH (``wp.Mesh``).
    points
        ``(m,)`` positions to shade.
    normals
        ``(m,)`` outward unit normals; see
        [`ambient_occlusion`][triwarp.visibility.ambient_occlusion] for the default.
    n_rays
        Directions per point.
    tau
        Attenuation rate, in inverse length units of the mesh — so it is **not** scale-invariant,
        and a mesh scaled by ``k`` wants ``tau / k`` for the same result. MeshLab's default is
        ``0.1``, which suits a mesh of extent order 1. As ``tau -> 0`` this becomes
        [`ambient_occlusion`][triwarp.visibility.ambient_occlusion]; as ``tau -> inf`` everything
        reads ``0``. Must be positive.
    weight
        Ray weighting, as in [`ambient_occlusion`][triwarp.visibility.ambient_occlusion].
    max_t
        Maximum ray length.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` obscurance-weighted occlusion in ``[0, 1]`` on ``points.device``, ``0`` where
        nothing nearby blocks. Exactly ``0`` on a convex closed mesh.

    Raises
    ------
    ValueError
        If ``tau <= 0``, ``n_rays < 1``, ``weight`` is unknown, or ``normals`` is the wrong length.
    RuntimeError
        If ``mesh``, ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`ambient_occlusion`][triwarp.visibility.ambient_occlusion]
    """
    require_same_device(mesh=mesh, points=points, normals=normals)
    if tau <= 0.0:
        raise ValueError(f"tau must be positive, got {tau}; use ambient_occlusion for the limit")
    return _occlusion_bundle(
        mesh, points, normals, n_rays, weight, tau, max_t, "volumetric_obscurance"
    )

shape_diameter

shape_diameter(
    mesh: Mesh,
    points: array[vec3],
    *,
    normals: array[vec3] | None = None,
    n_rays: int = 64,
    cone_angle: float = math.pi / 3.0,
    trim: float = 1.0,
    max_t: float | None = None,
) -> wp.array[wp.float32]

Shape diameter function: local thickness of the volume, from an inward cone of rays.

Shapira et al.'s SDF and MeshLab's compute_scalar_by_shape_diameter_function_per_vertex. A cone of n_rays rays is fired into the volume about -normal, each is traced to the far side, and the result is the cosine-weighted mean of those distances after discarding the outliers — the rays that escaped through a nearby opening or crossed the entire model, which would otherwise dominate the average near any concavity.

This is the many-ray generalization of thickness(method="ray"): at n_rays=1 with a vanishing cone_angle the bundle collapses to the inward normal and the two agree to float32. The extra rays are what make it stable — a single ray through a thin sliver of geometry reads a thickness the neighbourhood does not have — and the reason it is the quantity used for skeleton extraction and part segmentation rather than the one-ray version.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangle mesh with a built BVH (wp.Mesh). Should be closed: on an open surface the rays that find nothing on the far side are simply absent from the mean.	required
points	array[vec3]	(m,) surface positions to measure at, normally the mesh's own vertices.	required
normals	array[vec3] | None	(m,) outward unit normals; the cone opens along -normals. When None they are taken from the closest face of mesh. Pass vertex_normals at weighting="area" for a smooth field over a mesh's own vertices.	None
n_rays	int	Rays per point. MeshLab's default is 64, which is this one. Note that the single ray of n_rays=1 is the centroid of the cone's Fibonacci lattice rather than its axis, so it only coincides with the inward normal as cone_angle goes to zero.	64
cone_angle	float	Half-angle of the cone in radians, measured from the inward normal. The default pi / 3 (60 degrees) is Shapira's 120-degree cone. Must be in (0, pi / 2] — beyond that the cone reaches around to the outside of the surface and the distances stop meaning thickness.	pi / 3.0
trim	float	Keep only rays whose distance is within trim standard deviations of the mean before averaging. 1.0 (the default) is Shapira's rule; a large value keeps everything and turns this into a plain weighted mean. Must be non-negative.	1.0
max_t	float | None	Maximum ray length. When None, the diagonal of the AABB enclosing the mesh and the query points.	None

Returns:

Type	Description
array[float32]	(m,) diameters in the mesh's own length units on points.device. inf at a point where no ray in the cone found the far side at all.

Raises:

Type	Description
ValueError	If n_rays < 1, cone_angle is outside (0, pi / 2], trim < 0, or normals has a different length from points.
RuntimeError	If mesh, points and normals are not all on one device.

See Also

thickness max_tangent_sphere triwarp.visibility.ambient_occlusion triwarp.sample.sample_fibonacci_cone

Notes

MeshLab's cone_amplitude parameter is a no-op in the 2025.07 build — its output is byte-identical at 90 and 120 degrees — and its trimming rule is not the one documented in the paper, so its values differ from these by a roughly constant factor on a given mesh. Compare against it by rank rather than by value; the exactly-checkable statements are the reduction to thickness at n_rays=1 and the analytic 2 R on a sphere.

Source code in triwarp/visibility.py

def shape_diameter(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    normals: wp.array[wp.vec3] | None = None,
    n_rays: int = 64,
    cone_angle: float = math.pi / 3.0,
    trim: float = 1.0,
    max_t: float | None = None,
) -> wp.array[wp.float32]:
    """
    Shape diameter function: local thickness of the volume, from an inward cone of rays.

    Shapira et al.'s SDF and MeshLab's ``compute_scalar_by_shape_diameter_function_per_vertex``. A
    cone of ``n_rays`` rays is fired *into* the volume about ``-normal``, each is traced to the far
    side, and the result is the cosine-weighted mean of those distances **after discarding the
    outliers** — the rays that escaped through a nearby opening or crossed the entire model, which
    would otherwise dominate the average near any concavity.

    This is the many-ray generalization of
    [`thickness(method="ray")`][triwarp.visibility.thickness]: at ``n_rays=1`` with a vanishing
    ``cone_angle`` the bundle collapses to the inward normal and the two agree to float32. The extra
    rays are what make it stable — a single ray through a thin sliver of geometry reads a thickness
    the neighbourhood does not have — and the reason it is the quantity used for skeleton extraction
    and part segmentation rather than the one-ray version.

    Parameters
    ----------
    mesh
        Triangle mesh with a built BVH (``wp.Mesh``). Should be closed: on an open surface the rays
        that find nothing on the far side are simply absent from the mean.
    points
        ``(m,)`` surface positions to measure at, normally the mesh's own vertices.
    normals
        ``(m,)`` **outward** unit normals; the cone opens along ``-normals``. When ``None`` they are
        taken from the closest face of ``mesh``. Pass
        [`vertex_normals`][triwarp.vertices.vertex_normals] at ``weighting="area"`` for a smooth
        field over a mesh's own vertices.
    n_rays
        Rays per point. MeshLab's default is ``64``, which is this one. Note that the single ray of
        ``n_rays=1`` is the *centroid* of the cone's Fibonacci lattice rather than its axis, so it
        only coincides with the inward normal as ``cone_angle`` goes to zero.
    cone_angle
        Half-angle of the cone in **radians**, measured from the inward normal. The default
        ``pi / 3`` (60 degrees) is Shapira's 120-degree cone. Must be in ``(0, pi / 2]`` — beyond
        that the cone reaches around to the outside of the surface and the distances stop meaning
        thickness.
    trim
        Keep only rays whose distance is within ``trim`` standard deviations of the mean before
        averaging. ``1.0`` (the default) is Shapira's rule; a large value keeps everything and turns
        this into a plain weighted mean. Must be non-negative.
    max_t
        Maximum ray length. When ``None``, the diagonal of the AABB enclosing the mesh and the query
        points.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` diameters in the mesh's own length units on ``points.device``. ``inf`` at a point
        where no ray in the cone found the far side at all.

    Raises
    ------
    ValueError
        If ``n_rays < 1``, ``cone_angle`` is outside ``(0, pi / 2]``, ``trim < 0``, or ``normals``
        has a different length from ``points``.
    RuntimeError
        If ``mesh``, ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`thickness`][triwarp.visibility.thickness]
    [`max_tangent_sphere`][triwarp.visibility.max_tangent_sphere]
    [`triwarp.visibility.ambient_occlusion`][triwarp.visibility.ambient_occlusion]
    [`triwarp.sample.sample_fibonacci_cone`][triwarp.sample.sample_fibonacci_cone]

    Notes
    -----
    MeshLab's ``cone_amplitude`` parameter is a **no-op** in the 2025.07 build — its output is
    byte-identical at ``90`` and ``120`` degrees — and its trimming rule is not the one documented
    in the paper, so its values differ from these by a roughly constant factor on a given mesh.
    Compare against it by rank rather than by value; the exactly-checkable statements are the
    reduction to [`thickness`][triwarp.visibility.thickness] at ``n_rays=1`` and the analytic
    ``2 R`` on a sphere.
    """
    require_same_device(mesh=mesh, points=points, normals=normals)
    if n_rays < 1:
        raise ValueError(f"shape_diameter requires n_rays >= 1, got {n_rays}")
    if not 0.0 < cone_angle <= math.pi / 2.0:
        raise ValueError(f"cone_angle must be in (0, pi / 2] radians, got {cone_angle}")
    if trim < 0.0:
        raise ValueError(f"trim must be non-negative, got {trim}")

    device = points.device
    m = int(points.shape[0])
    # See `_occlusion_bundle`'s identical comment: resolved before the empty-input check so a
    # mismatched `normals` length is caught even when there is nothing to measure.
    normals, diagonal = _resolve_normals_and_radius(mesh, points, normals, "shape_diameter")
    if m == 0:
        return wp.empty(0, dtype=wp.float32, device=device)

    directions = tw.sample.sample_fibonacci_cone(n_rays, cone_angle, device=device)
    # Distances are kept so the trimming pass can revisit them against a mean the first pass had not
    # finished computing; re-tracing instead would double the only expensive part of the kernel.
    scratch = twt.empty_2d((m, n_rays), wp.float32, device=device)
    out_diameter = wp.empty(m, dtype=wp.float32, device=device)
    wp.launch_tiled(
        kernel_visibility.shape_diameter,
        dim=(m,),
        inputs=[
            mesh.id,
            points,
            normals,
            directions,
            wp.float32(max_t if max_t is not None else diagonal),
            wp.float32(_SURFACE_OFFSET * max(diagonal, 1e-12)),
            wp.float32(trim),
            scratch,
            out_diameter,
        ],
        block_dim=kernel_visibility.BUNDLE_BLOCK,
        device=device,
    )
    return out_diameter

thickness

thickness(
    mesh: Mesh,
    points: array[vec3],
    *,
    exterior: bool = False,
    normals: array[vec3] | None = None,
    method: Literal["max_sphere", "ray"] = "max_sphere",
) -> wp.array[wp.float32]

Local thickness of the volume at each point, by one inward ray or one tangent sphere.

A dispatcher over the module's other two inward measures, and the cheapest of the three: it takes a single piece of evidence per point where shape_diameter fires a whole cone and trims the outliers. method="max_sphere" returns twice the radius of max_tangent_sphere, which answers the question for a volume rather than along a direction; method="ray" returns longest_ray along -normals (or +normals with exterior=True), which is one ray and therefore reads whatever thin sliver of geometry it happens to cross.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangle mesh with a built BVH (wp.Mesh).	required
points	array[vec3]	(m,) surface positions to measure at.	required
exterior	bool	When True measure outward (the reach) instead of inward (the thickness).	False
normals	array[vec3] | None	(m,) outward unit normals. When None they are taken from the closest face of mesh; see ambient_occlusion.	None
method	Literal['max_sphere', 'ray']	"max_sphere" (default) or "ray"; see the summary for the difference.	'max_sphere'

Returns:

Type	Description
array[float32]	(m,) thickness values in the mesh's own length units on points.device. inf where the measure is unbounded (no far side was found).

Raises:

Type	Description
ValueError	If method is neither "max_sphere" nor "ray", or if normals has a different length from points.
RuntimeError	If mesh, points and normals are not all on one device.

See Also

max_tangent_sphere shape_diameter The stable many-ray generalization of method="ray". longest_ray

Source code in triwarp/visibility.py

def thickness(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    exterior: bool = False,
    normals: wp.array[wp.vec3] | None = None,
    method: Literal["max_sphere", "ray"] = "max_sphere",
) -> wp.array[wp.float32]:
    """
    Local thickness of the volume at each point, by one inward ray or one tangent sphere.

    A dispatcher over the module's other two inward measures, and the cheapest of the three: it
    takes a single piece of evidence per point where
    [`shape_diameter`][triwarp.visibility.shape_diameter] fires a whole cone and trims the outliers.
    ``method="max_sphere"`` returns twice the radius of
    [`max_tangent_sphere`][triwarp.visibility.max_tangent_sphere], which answers the question for a
    *volume* rather than along a direction; ``method="ray"`` returns
    [`longest_ray`][triwarp.ray.longest_ray] along ``-normals`` (or ``+normals`` with
    ``exterior=True``), which is one ray and therefore reads whatever thin sliver of geometry it
    happens to cross.

    Parameters
    ----------
    mesh
        Triangle mesh with a built BVH (``wp.Mesh``).
    points
        ``(m,)`` surface positions to measure at.
    exterior
        When ``True`` measure outward (the reach) instead of inward (the thickness).
    normals
        ``(m,)`` **outward** unit normals. When ``None`` they are taken from the closest face of
        ``mesh``; see [`ambient_occlusion`][triwarp.visibility.ambient_occlusion].
    method
        ``"max_sphere"`` (default) or ``"ray"``; see the summary for the difference.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` thickness values in the mesh's own length units on ``points.device``. ``inf``
        where the measure is unbounded (no far side was found).

    Raises
    ------
    ValueError
        If ``method`` is neither ``"max_sphere"`` nor ``"ray"``, or if ``normals`` has a different
        length from ``points``.
    RuntimeError
        If ``mesh``, ``points`` and ``normals`` are not all on one device.

    See Also
    --------
    [`max_tangent_sphere`][triwarp.visibility.max_tangent_sphere]
    [`shape_diameter`][triwarp.visibility.shape_diameter]
        The stable many-ray generalization of ``method="ray"``.
    [`longest_ray`][triwarp.ray.longest_ray]
    """
    require_same_device(mesh=mesh, points=points, normals=normals)
    if method not in _THICKNESS_METHODS:
        raise ValueError(f"method must be one of {sorted(_THICKNESS_METHODS)}, got {method!r}")

    if method == "max_sphere":
        _centers, radii = max_tangent_sphere(mesh, points, inwards=not exterior, normals=normals)
        wp.map(wp.mul, radii, wp.float32(2.0), out=radii)
        return radii

    normals, max_t = _resolve_normals_and_radius(mesh, points, normals, "thickness")
    ray_dirs = normals
    if not exterior:
        ray_dirs = wp.empty(int(points.shape[0]), dtype=wp.vec3, device=points.device)
        wp.map(wp.neg, normals, out=ray_dirs)
    return tw.ray.longest_ray(mesh, points, ray_dirs, max_t=max_t)

max_tangent_sphere

max_tangent_sphere(
    mesh: Mesh,
    points: array[vec3],
    *,
    inwards: bool = True,
    normals: array[vec3] | None = None,
    threshold: float = 1e-06,
    max_iter: int = 100,
) -> tuple[wp.array[wp.vec3], wp.array[wp.float32]]

Find the center and radius of the sphere tangent to the mesh at each point.

Implements the shrinking-sphere algorithm (Inui et al. 2016): iteratively finds the largest sphere tangent to the mesh at points with no non-tangential intersections.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Warp mesh (BVH built by caller).	required
points	array[vec3]	(m,) surface points as wp.vec3.	required
inwards	bool	If True, sphere grows inward (into the mesh interior). If False, grows outward.	True
normals	array[vec3] | None	(m,) unit surface normals at points. If None, computed from the closest triangle. A caller-supplied array is normalized defensively (the default never needs it), since a non-unit normal would silently scale the reported radius.	None
threshold	float	Convergence threshold as a fraction of the mesh's own bounding-box diagonal, not the (possibly larger) box enclosing the query points too -- a query far outside the mesh should not loosen how tightly the sphere converges.	1e-06
max_iter	int	Maximum number of shrink iterations.	100

Returns:

Type	Description
centers	(m,) sphere center positions as wp.vec3.
radii	(m,) sphere radii as float32. inf when the sphere is unbounded.

Raises:

Type	Description
ValueError	If normals has a different length from points.
RuntimeError	If mesh, points and normals are not all on one device.

Source code in triwarp/visibility.py

def max_tangent_sphere(
    mesh: wp.Mesh,
    points: wp.array[wp.vec3],
    *,
    inwards: bool = True,
    normals: wp.array[wp.vec3] | None = None,
    threshold: float = 1e-6,
    max_iter: int = 100,
) -> tuple[wp.array[wp.vec3], wp.array[wp.float32]]:
    """
    Find the center and radius of the sphere tangent to the mesh at each point.

    Implements the shrinking-sphere algorithm (Inui et al. 2016): iteratively
    finds the largest sphere tangent to the mesh at ``points`` with no
    non-tangential intersections.

    Parameters
    ----------
    mesh
        Warp mesh (BVH built by caller).
    points
        ``(m,)`` surface points as ``wp.vec3``.
    inwards
        If ``True``, sphere grows inward (into the mesh interior). If ``False``,
        grows outward.
    normals
        ``(m,)`` unit surface normals at ``points``. If ``None``, computed from
        the closest triangle. A caller-supplied array is normalized defensively (the default
        never needs it), since a non-unit normal would silently scale the reported radius.
    threshold
        Convergence threshold as a fraction of the **mesh's own** bounding-box diagonal, not the
        (possibly larger) box enclosing the query points too -- a query far outside the mesh
        should not loosen how tightly the sphere converges.
    max_iter
        Maximum number of shrink iterations.

    Returns
    -------
    centers
        ``(m,)`` sphere center positions as ``wp.vec3``.
    radii
        ``(m,)`` sphere radii as ``float32``. ``inf`` when the sphere is
        unbounded.

    Raises
    ------
    ValueError
        If ``normals`` has a different length from ``points``.
    RuntimeError
        If ``mesh``, ``points`` and ``normals`` are not all on one device.
    """
    require_same_device(mesh=mesh, points=points, normals=normals)
    device = points.device
    m = int(points.shape[0])
    # Resolved (and, for a caller-supplied array, normalized) before the empty-input early return
    # and before the AABB reductions below -- a mismatched or non-unit ``normals`` is a property of
    # ``normals`` alone, and there is no reason to pay for two reductions first only to reject it.
    normals = _resolve_normals(mesh, points, normals, "max_tangent_sphere")
    if m == 0:
        return (
            wp.empty(0, dtype=wp.vec3, device=device),
            wp.empty(0, dtype=wp.float32, device=device),
        )
    # Every ray-bundle measure in this module normalizes defensively per-thread (see
    # `hemisphere_frame`); this path shrinks a sphere along `normals` directly with no such guard,
    # so a caller-supplied non-unit normal silently scales `sphere_center`'s radius away from what
    # `step_sphere_shrink`'s own convergence test expects. One `wp.map` fixes it for the whole
    # iterative loop; `normals_at_closest_faces`'s own output is already unit, so this is a no-op
    # there, but cheap enough not to special-case.
    unit_normals = wp.empty(m, dtype=wp.vec3, device=device)
    wp.map(wp.normalize, normals, out=unit_normals)
    normals = unit_normals

    # One reduction of ``mesh.points``, not two: ``max_t`` needs the box enclosing the mesh *and*
    # the queries, while the convergence threshold is a fraction of the mesh's own diagonal. Taking
    # the mesh corners once and deriving both saves an ``aabb`` pass and its host sync.
    mesh_lower, mesh_upper = tw.bounds.aabb(mesh.points)
    query_lower, query_upper = tw.bounds.aabb(points)
    union_lower, union_upper = tw.bounds.aabb_union(
        mesh_lower, mesh_upper, query_lower, query_upper
    )
    max_t = float(wp.length(union_upper - union_lower))
    mesh_diagonal = float(wp.length(mesh_upper - mesh_lower))

    ray_dirs = normals
    if inwards:
        ray_dirs = wp.empty(m, dtype=wp.vec3, device=device)
        wp.map(wp.neg, normals, out=ray_dirs)

    distances = tw.ray.longest_ray(mesh, points, ray_dirs, max_t=max_t)

    n_verts = int(mesh.points.shape[0])
    radii = wp.empty(m, dtype=wp.float32, device=device)
    not_converged = wp.empty(m, dtype=wp.bool, device=device)
    needs_support = wp.empty(m, dtype=wp.bool, device=device)
    wp.map(
        kernel_visibility.init_sphere_radii_finite,
        distances,
        out=[radii, not_converged, needs_support],
    )
    # Escaped rays (typically exterior/reach queries) need the support point of the vertex
    # cloud in the ray direction. Compact them first — interior queries usually leave the
    # subset empty — then run one grid-stride packed-argmax pass over the vertices for just
    # that subset instead of a serial all-vertices loop per query thread.
    support_indices = tw.array.flatnonzero(needs_support)
    k = int(support_indices.shape[0])
    if k > 0:
        n_vert_slices = max(1, (n_verts + ITEMS_PER_QUERY_SLICE - 1) // ITEMS_PER_QUERY_SLICE)
        packed_support = wp.zeros(k, dtype=wp.uint64, device=device)
        wp.launch(
            kernel_visibility.support_argmax_sliced,
            dim=(k, n_vert_slices),
            inputs=[
                mesh.points,
                wp.int32(n_verts),
                wp.int32(n_vert_slices),
                ray_dirs,
                support_indices,
                packed_support,
            ],
            device=device,
        )
        wp.launch(
            kernel_visibility.init_sphere_radii_support,
            dim=k,
            inputs=[
                mesh.points,
                points,
                ray_dirs,
                support_indices,
                packed_support,
                radii,
                not_converged,
            ],
            device=device,
        )

    centers = wp.empty(m, dtype=wp.vec3, device=device)
    wp.map(kernel_visibility.sphere_center, points, ray_dirs, radii, out=centers)

    convergence_threshold = wp.float32(threshold * mesh_diagonal)

    # All per-iteration buffers are preallocated once and ping-ponged (the step kernel writes
    # every lane, passing converged state through). The convergence count is checked every
    # iteration on purpose: an extra iteration runs a full BVH closest-point pass, far more
    # expensive than the 8-byte readback the check costs.
    n_pts_wp = wp.empty(m, dtype=wp.vec3, device=device)
    n_dists_wp = wp.empty(m, dtype=wp.float32, device=device)
    n_face_wp = wp.empty(m, dtype=wp.int32, device=device)
    new_radii = wp.empty(m, dtype=wp.float32, device=device)
    new_centers = wp.empty(m, dtype=wp.vec3, device=device)
    new_nc = wp.empty(m, dtype=wp.bool, device=device)

    for _ in range(max_iter):
        if tw.reduce.sum(not_converged) == 0:
            break

        wp.launch(
            kernel_proximity.closest_point_on_mesh,
            dim=m,
            inputs=[mesh.id, centers, wp.float32(max_t), n_pts_wp, n_dists_wp, n_face_wp],
            device=device,
        )
        wp.launch(
            kernel_visibility.step_sphere_shrink,
            dim=m,
            inputs=[
                points,
                ray_dirs,
                n_pts_wp,
                n_dists_wp,
                centers,
                radii,
                convergence_threshold,
                not_converged,
                new_radii,
                new_centers,
                new_nc,
            ],
            device=device,
        )
        radii, new_radii = new_radii, radii
        centers, new_centers = new_centers, centers
        not_converged, new_nc = new_nc, not_converged

    return centers, radii