ordito.metrics

def chamfer_points_to_points(
    x: wp.array[wp.vec3],
    y: wp.array[wp.vec3],
    *,
    point_reduction: _PointReduction | None = "mean",
    single_directional: bool = False,
    squared: bool = False,
) -> float | _UnreducedChamfer:
    """
    Chamfer distance between two point clouds.

    For each point in ``x`` the Euclidean distance to its nearest neighbor in ``y`` is
    accumulated (and symmetrically for ``y`` into ``x`` unless ``single_directional``).
    ``pytorch3d.loss.chamfer_distance`` is this call with ``squared=True``; its
    ``single_directional`` keyword maps onto this one exactly.

    Parameters
    ----------
    x
        ``(n,)`` query point cloud.
    y
        ``(m,)`` target point cloud.
    point_reduction
        How to reduce the per-point distances of each direction:
        ``"mean"`` (default), ``"sum"``, ``"max"``, or ``None`` to return the
        per-point distances unreduced. ``"max"`` yields a directed
        Hausdorff-like value.
    single_directional
        If ``True``, only the ``x -> y`` term is computed.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.

    Returns
    -------
    float or wp.array[wp.float32] or tuple
        The reduced Chamfer distance as a Python ``float`` when ``point_reduction``
        is not ``None``. Otherwise the per-point distances: a single
        ``(n,)`` array when ``single_directional``, else an ``((n,), (m,))`` tuple.
        Empty inputs yield ``0.0`` (reduced) or empty arrays.

    Raises
    ------
    RuntimeError
        If ``x`` and ``y`` are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"``, ``"sum"``, ``"max"`` or ``None``.

    See Also
    --------
    [`chamfer_points_to_mesh`][ordito.metrics.chamfer_points_to_mesh]
    [`chamfer_mesh_to_mesh`][ordito.metrics.chamfer_mesh_to_mesh]
    [`query_nearest`][ordito.neighbors.query_nearest]

    Examples
    --------
    The ``x -> y`` distances are ``[0, 1]`` and the ``y -> x`` ones ``[0, 3]``; the symmetric value
    is the sum of the two means:

    ```python
    import warp as wp

    import ordito as od

    x = wp.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], dtype=wp.vec3)
    y = wp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 3.0]], dtype=wp.vec3)
    print(od.metrics.chamfer_points_to_points(x, y))  # 2.0
    print(od.metrics.chamfer_points_to_points(x, y, single_directional=True))  # 0.5
    print(od.metrics.chamfer_points_to_points(x, y, squared=True))  # 5.0
    forward, backward = od.metrics.chamfer_points_to_points(x, y, point_reduction=None)
    print(forward.numpy().tolist(), backward.numpy().tolist())  # [0.0, 1.0] [0.0, 3.0]
    ```
    """
    require_same_device(x=x, y=y)
    _validate_point_reduction(point_reduction)
    distances = _distances_points_to_points(x, y, single_directional)
    if distances is None:
        return _empty_chamfer(point_reduction, single_directional, x.device)
    return _chamfer(distances, point_reduction, single_directional, squared)

def chamfer_points_to_mesh(
    points: wp.array[wp.vec3],
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    point_reduction: _PointReduction | None = "mean",
    single_directional: bool = False,
    squared: bool = False,
    mesh: wp.Mesh | None = None,
) -> float | _UnreducedChamfer:
    """
    Chamfer distance between a point cloud and a triangle mesh.

    The forward term is the distance from each point to the mesh surface
    (exact, via [`closest_point_on_mesh`][ordito.proximity.closest_point_on_mesh]);
    the backward term is the distance from each mesh vertex to its nearest
    point in the cloud.

    Parameters
    ----------
    points
        ``(n,)`` query point cloud.
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index array.
    point_reduction
        ``"mean"`` (default), ``"sum"``, ``"max"``, or ``None``. See
        [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points].
    single_directional
        If ``True``, only the ``points -> mesh surface`` term is computed.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the BVH build each
        call otherwise pays (worth it when one mesh is measured repeatedly). Not checked against
        ``vertices`` / ``faces``: the answer is for the geometry ``mesh`` holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    float or wp.array[wp.float32] or tuple
        Same layout as [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points].
        The forward array is ``(n,)`` and the backward array ``(n_vertices,)``, except when
        ``points``, ``vertices`` or ``faces`` is empty -- a mesh with no faces counts as empty here
        even when ``vertices`` is not, since there is no surface for the forward term to measure
        against -- in which case both arrays are length 0, matching
        [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points]'s own convention.

    Raises
    ------
    RuntimeError
        If ``points``, ``vertices``, ``faces`` and ``mesh`` are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"``, ``"sum"``, ``"max"`` or ``None``.

    See Also
    --------
    [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points]
    [`chamfer_mesh_to_mesh`][ordito.metrics.chamfer_mesh_to_mesh]
    [`closest_point_on_mesh`][ordito.proximity.closest_point_on_mesh]

    Examples
    --------
    A point ``1.0`` outside a unit cube and one at its centre, ``0.5`` from every face. The backward
    term adds the mean distance from the cube's corners to the cloud, ``sqrt(3) / 2``:

    ```python
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    points = wp.array([[1.5, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=wp.vec3)
    forward = od.metrics.chamfer_points_to_mesh(points, vertices, faces, single_directional=True)
    print(forward)  # 0.75
    print(f"{od.metrics.chamfer_points_to_mesh(points, vertices, faces):.4f}")  # 1.6160
    ```
    """
    require_same_device(points=points, vertices=vertices, faces=faces, mesh=mesh)
    _validate_point_reduction(point_reduction)
    distances = _distances_points_to_mesh(points, vertices, faces, single_directional, mesh)
    if distances is None:
        return _empty_chamfer(point_reduction, single_directional, points.device)
    return _chamfer(distances, point_reduction, single_directional, squared)

def chamfer_mesh_to_mesh(
    vertices_a: wp.array[wp.vec3],
    faces_a: wp.array[wp.int32],
    vertices_b: wp.array[wp.vec3],
    faces_b: wp.array[wp.int32],
    *,
    point_reduction: _PointReduction | None = "mean",
    single_directional: bool = False,
    squared: bool = False,
    mesh_a: wp.Mesh | None = None,
    mesh_b: wp.Mesh | None = None,
) -> float | _UnreducedChamfer:
    """
    Chamfer distance between two triangle meshes (vertex-to-surface).

    The forward term is the distance from each vertex of mesh ``A`` to the
    surface of mesh ``B``; the backward term is the reverse. This mirrors the
    vertex-based approach of libigl while accumulating (mean/sum) instead of maximizing.

    Parameters
    ----------
    vertices_a, faces_a
        ``(n_vertices_a,)`` vertices and ``(3 * n_faces_a,)`` flat faces of mesh ``A``.
    vertices_b, faces_b
        ``(n_vertices_b,)`` vertices and ``(3 * n_faces_b,)`` flat faces of mesh ``B``.
    point_reduction
        ``"mean"`` (default), ``"sum"``, ``"max"``, or ``None``. See
        [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points].
    single_directional
        If ``True``, only the ``A -> surface(B)`` term is computed.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.
    mesh_a, mesh_b
        ``wp.Mesh`` objects already built over mesh ``A`` and mesh ``B``, to spare the BVH build
        each call otherwise pays per mesh queried (worth it when a mesh is measured repeatedly).
        Not checked against the arrays: the answer is for the geometry each holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    float or wp.array[wp.float32] or tuple
        Same layout as [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points].
        The forward array is ``(n_vertices_a,)`` and the backward array ``(n_vertices_b,)``, except
        when any of ``vertices_a``, ``faces_a``, ``vertices_b`` or ``faces_b`` is empty -- a mesh
        with no faces counts as empty here even when its vertices are not -- in which case both
        arrays are length 0.

    Raises
    ------
    RuntimeError
        If ``vertices_a``, ``faces_a``, ``vertices_b``, ``faces_b``, ``mesh_a`` and ``mesh_b``
        are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"``, ``"sum"``, ``"max"`` or ``None``.

    See Also
    --------
    [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points]
    [`chamfer_points_to_mesh`][ordito.metrics.chamfer_points_to_mesh]
    [`hausdorff_mesh_to_mesh`][ordito.metrics.hausdorff_mesh_to_mesh]

    Examples
    --------
    A unit cube inside a cube of side 2: each inner corner is ``0.5`` from the outer surface, each
    outer corner ``sqrt(3) / 2`` from the inner one:

    ```python
    import ordito as od

    vertices_a, faces_a = od.creation.box()
    vertices_b, faces_b = od.creation.box(extents=(2.0, 2.0, 2.0))
    forward = od.metrics.chamfer_mesh_to_mesh(
        vertices_a, faces_a, vertices_b, faces_b, single_directional=True
    )
    print(forward)  # 0.5
    symmetric = od.metrics.chamfer_mesh_to_mesh(vertices_a, faces_a, vertices_b, faces_b)
    print(f"{symmetric:.4f}")  # 1.3660
    ```
    """
    require_same_device(
        vertices_a=vertices_a,
        faces_a=faces_a,
        vertices_b=vertices_b,
        faces_b=faces_b,
        mesh_a=mesh_a,
        mesh_b=mesh_b,
    )
    _validate_point_reduction(point_reduction)
    distances = _distances_mesh_to_mesh(
        vertices_a, faces_a, vertices_b, faces_b, single_directional, mesh_a, mesh_b
    )
    if distances is None:
        return _empty_chamfer(point_reduction, single_directional, vertices_a.device)
    return _chamfer(distances, point_reduction, single_directional, squared)

def chamfer_points_to_points_loss(
    x: wp.array[wp.vec3],
    y: wp.array[wp.vec3],
    *,
    tape: wp.Tape | None = None,
    point_reduction: _DiffReduction = "mean",
    single_directional: bool = False,
    squared: bool = False,
) -> wp.array[wp.float32]:
    """
    Differentiable Chamfer loss between two point clouds.

    Squared-distance analogue of
    [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points] that returns
    a length-1 ``wp.float32`` device array instead of a host ``float``, so the loss can be
    back-propagated to ``x`` and ``y`` through a caller-owned ``wp.Tape``.

    The nearest-neighbor assignment (via
    [`query_nearest`][ordito.neighbors.query_nearest]) is computed
    outside ``tape`` and held constant during the backward pass, matching pytorch3d's
    ``chamfer_distance`` gradient.

    Parameters
    ----------
    x
        ``(n,)`` query point cloud. Set ``requires_grad=True`` for gradients.
    y
        ``(m,)`` target point cloud. Set ``requires_grad=True`` for gradients.
    tape
        Caller-owned ``wp.Tape`` into which the differentiable kernels are recorded. When
        ``None`` the loss is still computed but no operations are taped (no gradient).
    point_reduction
        ``"mean"`` (default) or ``"sum"``; reduces the per-point distances of each
        direction. ``"max"`` and ``None`` are not supported (see
        [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points]).
    single_directional
        If ``True``, only the ``x -> y`` term is accumulated.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.

    Returns
    -------
    wp.array[wp.float32]
        ``(1,)`` device array holding the Chamfer loss (``requires_grad=True``). Empty
        inputs yield a length-1 zero array.

    Raises
    ------
    RuntimeError
        If ``x`` and ``y`` are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"`` or ``"sum"``.

    See Also
    --------
    [`chamfer_points_to_points`][ordito.metrics.chamfer_points_to_points]
    [`chamfer_points_to_mesh_loss`][ordito.metrics.chamfer_points_to_mesh_loss]
    [`chamfer_mesh_to_mesh_loss`][ordito.metrics.chamfer_mesh_to_mesh_loss]

    Examples
    --------
    The point ``(1, 0, 0)`` is 1 from its nearest neighbour at the origin, so the mean loss is
    ``0.5`` and its gradient pushes it towards the origin with weight ``1 / 2``:

    ```python
    import warp as wp

    import ordito as od

    x = wp.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], dtype=wp.vec3, requires_grad=True)
    y = wp.array([[0.0, 0.0, 0.0], [1.0, 0.0, 2.0]], dtype=wp.vec3)
    tape = wp.Tape()
    loss = od.metrics.chamfer_points_to_points_loss(x, y, tape=tape, single_directional=True)
    tape.backward(loss=loss)
    print(loss.numpy())  # [0.5]
    print(x.grad.numpy()[1])  # [0.5 0.  0. ]
    ```
    """
    require_same_device(x=x, y=y)
    _validate_diff_reduction(point_reduction)
    device = x.device
    n = x.size
    m = y.size
    loss = _zero_loss(device)
    if n == 0 or m == 0:
        return loss

    # Non-differentiable nearest-neighbor indices (computed outside the tape). The backward search
    # is chosen from the forward one's own distances -- see `_backward_nearest`.
    nearest_xy, distance_xy = _nearest(y, x)
    terms: list[_ChamferTerm] = [_NnTerm(x, y, nearest_xy, _reduction_scale(point_reduction, n))]
    if not single_directional:
        nearest_yx = _backward_nearest(x, y, distance_xy, nearest_xy)[0][0]
        terms.append(_NnTerm(y, x, nearest_yx, _reduction_scale(point_reduction, m)))

    _accumulate_chamfer_terms(tape, terms, loss, squared=squared)
    return loss

def chamfer_points_to_mesh_loss(
    points: wp.array[wp.vec3],
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    tape: wp.Tape | None = None,
    point_reduction: _DiffReduction = "mean",
    single_directional: bool = False,
    squared: bool = False,
    mesh: wp.Mesh | None = None,
) -> wp.array[wp.float32]:
    """
    Differentiable Chamfer loss between a point cloud and a triangle mesh.

    Squared-distance analogue of
    [`chamfer_points_to_mesh`][ordito.metrics.chamfer_points_to_mesh]. The forward term
    is the exact point-to-surface distance (closest triangle held constant during
    backprop); the backward term is the distance from each mesh vertex to its
    nearest point in the cloud. Gradients flow to ``points`` **and** ``vertices``.

    Parameters
    ----------
    points
        ``(n,)`` query point cloud.
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index array.
    tape
        Caller-owned ``wp.Tape`` for the backward pass, or ``None`` to skip taping.
    point_reduction
        ``"mean"`` (default) or ``"sum"``.
    single_directional
        If ``True``, only the ``points -> mesh surface`` term is accumulated.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the BVH build each
        call otherwise pays (worth it when one mesh is measured repeatedly). Not checked against
        ``vertices`` / ``faces``: the answer is for the geometry ``mesh`` holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    wp.array[wp.float32]
        ``(1,)`` device array holding the Chamfer loss (``requires_grad=True``).

    Raises
    ------
    RuntimeError
        If ``points``, ``vertices``, ``faces`` and ``mesh`` are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"`` or ``"sum"``.

    See Also
    --------
    [`chamfer_points_to_mesh`][ordito.metrics.chamfer_points_to_mesh]
    [`chamfer_points_to_points_loss`][ordito.metrics.chamfer_points_to_points_loss]
    [`chamfer_mesh_to_mesh_loss`][ordito.metrics.chamfer_mesh_to_mesh_loss]

    Examples
    --------
    A point ``1.0`` outside the ``x = 0.5`` face of a unit cube; the gradient of its distance points
    away from the face:

    ```python
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    points = wp.array([[1.5, 0.1, 0.2]], dtype=wp.vec3, requires_grad=True)
    tape = wp.Tape()
    loss = od.metrics.chamfer_points_to_mesh_loss(
        points, vertices, faces, tape=tape, single_directional=True
    )
    tape.backward(loss=loss)
    print(loss.numpy())  # [1.]
    print(points.grad.numpy()[0])  # [1. 0. 0.]
    ```
    """
    require_same_device(points=points, vertices=vertices, faces=faces, mesh=mesh)
    _validate_diff_reduction(point_reduction)
    device = points.device
    n = points.size
    v = vertices.size
    n_faces = faces.size // 3
    loss = _zero_loss(device)
    if n == 0 or v == 0 or n_faces == 0:
        return loss

    # Non-differentiable closest-face and nearest-neighbor assignment (outside the tape). The
    # backward search is chosen from the forward one's own distances -- see `_backward_nearest`.
    distance_forward, face_id = _closest_on_mesh(points, vertices, faces, mesh)
    terms: list[_ChamferTerm] = [
        _SurfaceTerm(points, vertices, faces, face_id, _reduction_scale(point_reduction, n))
    ]
    if not single_directional:
        nearest_vp = _backward_nearest(points, vertices, distance_forward, face_id, faces)[0][0]
        terms.append(_NnTerm(vertices, points, nearest_vp, _reduction_scale(point_reduction, v)))

    _accumulate_chamfer_terms(tape, terms, loss, squared=squared)
    return loss

def chamfer_mesh_to_mesh_loss(
    vertices_a: wp.array[wp.vec3],
    faces_a: wp.array[wp.int32],
    vertices_b: wp.array[wp.vec3],
    faces_b: wp.array[wp.int32],
    *,
    tape: wp.Tape | None = None,
    point_reduction: _DiffReduction = "mean",
    single_directional: bool = False,
    squared: bool = False,
    mesh_a: wp.Mesh | None = None,
    mesh_b: wp.Mesh | None = None,
) -> wp.array[wp.float32]:
    """
    Differentiable Chamfer loss between two triangle meshes (vertex-to-surface).

    Squared-distance analogue of
    [`chamfer_mesh_to_mesh`][ordito.metrics.chamfer_mesh_to_mesh]. The forward term is
    the distance from each vertex of mesh ``A`` to the surface of mesh ``B`` (with
    the closest triangle held constant during backprop); the backward term is the reverse.
    Gradients flow to both meshes' vertices.

    Parameters
    ----------
    vertices_a, faces_a
        ``(n_vertices_a,)`` vertices and ``(3 * n_faces_a,)`` flat faces of mesh ``A``.
    vertices_b, faces_b
        ``(n_vertices_b,)`` vertices and ``(3 * n_faces_b,)`` flat faces of mesh ``B``.
    tape
        Caller-owned ``wp.Tape`` for the backward pass, or ``None`` to skip taping.
    point_reduction
        ``"mean"`` (default) or ``"sum"``.
    single_directional
        If ``True``, only the ``A -> surface(B)`` term is accumulated.
    squared
        Square each distance before reducing, as ``pytorch3d`` does. Off by default: the result
        is then in squared units, and its ``"mean"`` is not a mean distance, which is easy to
        misread.
    mesh_a, mesh_b
        ``wp.Mesh`` objects already built over mesh ``A`` and mesh ``B``, to spare the BVH build
        each call otherwise pays per mesh queried (worth it when a mesh is measured repeatedly).
        Not checked against the arrays: the answer is for the geometry each holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    wp.array[wp.float32]
        ``(1,)`` device array holding the Chamfer loss (``requires_grad=True``).

    Raises
    ------
    RuntimeError
        If ``vertices_a``, ``faces_a``, ``vertices_b``, ``faces_b``, ``mesh_a`` and ``mesh_b``
        are not all on one device.
    ValueError
        If ``point_reduction`` is not ``"mean"`` or ``"sum"``.

    See Also
    --------
    [`chamfer_mesh_to_mesh`][ordito.metrics.chamfer_mesh_to_mesh]
    [`chamfer_points_to_points_loss`][ordito.metrics.chamfer_points_to_points_loss]
    [`chamfer_points_to_mesh_loss`][ordito.metrics.chamfer_points_to_mesh_loss]

    Examples
    --------
    A unit cube under a slab whose bottom face is the plane ``z = 1``: the cube's top corners are
    ``0.5`` from it and the bottom ones ``1.5``, a mean of 1, and every corner is pulled up with
    weight ``1 / 8``:

    ```python
    import warp as wp

    import ordito as od

    vertices_a, faces_a = od.creation.box()
    vertices_b, faces_b = od.creation.box(bounds=[[-2.0, -2.0, 1.0], [2.0, 2.0, 2.0]])
    vertices_a.requires_grad = True
    tape = wp.Tape()
    loss = od.metrics.chamfer_mesh_to_mesh_loss(
        vertices_a, faces_a, vertices_b, faces_b, tape=tape, single_directional=True
    )
    tape.backward(loss=loss)
    print(loss.numpy())  # [1.]
    print(set(vertices_a.grad.numpy()[:, 2].tolist()))  # {-0.125}
    ```
    """
    require_same_device(
        vertices_a=vertices_a,
        faces_a=faces_a,
        vertices_b=vertices_b,
        faces_b=faces_b,
        mesh_a=mesh_a,
        mesh_b=mesh_b,
    )
    _validate_diff_reduction(point_reduction)
    device = vertices_a.device
    va = vertices_a.size
    vb = vertices_b.size
    n_faces_a = faces_a.size // 3
    n_faces_b = faces_b.size // 3
    loss = _zero_loss(device)
    if va == 0 or vb == 0 or n_faces_a == 0 or n_faces_b == 0:
        return loss

    # Non-differentiable closest-face assignment for each direction (outside the tape).
    _distances, face_ids = _closest_both_ways(
        vertices_a, faces_a, vertices_b, faces_b, single_directional, mesh_a, mesh_b
    )
    face_id_ab = odt.as_dense(face_ids[:va])
    terms: list[_ChamferTerm] = [
        _SurfaceTerm(
            vertices_a, vertices_b, faces_b, face_id_ab, _reduction_scale(point_reduction, va)
        )
    ]
    if not single_directional:
        face_id_ba = odt.as_dense(face_ids[va:])
        terms.append(
            _SurfaceTerm(
                vertices_b, vertices_a, faces_a, face_id_ba, _reduction_scale(point_reduction, vb)
            )
        )

    _accumulate_chamfer_terms(tape, terms, loss, squared=squared)
    return loss

def hausdorff_points_to_points(
    x: wp.array[wp.vec3], y: wp.array[wp.vec3], *, single_directional: bool = False
) -> float:
    """
    Hausdorff distance between two point clouds.

    The directed distance ``d(x, y) = max_i min_j ||x_i - y_j||`` is computed (and
    symmetrized with ``d(y, x)`` unless ``single_directional``). Equivalent to
    symmetrizing [`scipy.spatial.distance.directed_hausdorff`][].

    Parameters
    ----------
    x
        ``(n,)`` query point cloud.
    y
        ``(m,)`` target point cloud.
    single_directional
        If ``True``, return only the directed distance ``d(x, y)``.

    Returns
    -------
    float
        The (symmetric or directed) Hausdorff distance. ``0.0`` for empty inputs.

    Raises
    ------
    RuntimeError
        If ``x`` and ``y`` are not all on one device.

    See Also
    --------
    [`hausdorff_points_to_mesh`][ordito.metrics.hausdorff_points_to_mesh]
    [`hausdorff_mesh_to_mesh`][ordito.metrics.hausdorff_mesh_to_mesh]
    [`scipy.spatial.distance.directed_hausdorff`][]

    Examples
    --------
    Every point of ``x`` is within 1 of ``y``, but ``(3, 0, 0)`` is 2 from ``x``:

    ```python
    import warp as wp

    import ordito as od

    x = wp.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], dtype=wp.vec3)
    y = wp.array([[0.0, 0.0, 0.0], [3.0, 0.0, 0.0]], dtype=wp.vec3)
    print(od.metrics.hausdorff_points_to_points(x, y))  # 2.0
    print(od.metrics.hausdorff_points_to_points(x, y, single_directional=True))  # 1.0
    ```
    """
    require_same_device(x=x, y=y)
    distances = _distances_points_to_points(x, y, single_directional)
    if distances is None:
        return 0.0
    return _hausdorff(distances, single_directional)

def hausdorff_points_to_mesh(
    points: wp.array[wp.vec3],
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    single_directional: bool = False,
    mesh: wp.Mesh | None = None,
) -> float:
    """
    Hausdorff distance between a point cloud and a triangle mesh.

    The forward distance is the maximum over points of the distance to the mesh
    surface; the backward distance is the maximum over mesh vertices of the distance
    to the nearest point in the cloud.

    Parameters
    ----------
    points
        ``(n,)`` query point cloud.
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index array.
    single_directional
        If ``True``, return only the directed ``points -> mesh surface`` distance.
    mesh
        A ``wp.Mesh`` already built over ``vertices`` and ``faces``, to spare the BVH build each
        call otherwise pays (worth it when one mesh is measured repeatedly). Not checked against
        ``vertices`` / ``faces``: the answer is for the geometry ``mesh`` holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    float
        The (symmetric or directed) Hausdorff distance. ``0.0`` for empty inputs.

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

    See Also
    --------
    [`hausdorff_points_to_points`][ordito.metrics.hausdorff_points_to_points]
    [`hausdorff_mesh_to_mesh`][ordito.metrics.hausdorff_mesh_to_mesh]

    Examples
    --------
    A point ``1.0`` outside a unit cube and one at its centre. The farthest cube corner from the
    cloud is ``sqrt(3) / 2`` from the centre point, so the symmetric value is still 1:

    ```python
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    points = wp.array([[1.5, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=wp.vec3)
    directed = od.metrics.hausdorff_points_to_mesh(points, vertices, faces, single_directional=True)
    print(directed)  # 1.0
    print(f"{od.metrics.hausdorff_points_to_mesh(points, vertices, faces):.4f}")  # 1.0000
    ```
    """
    require_same_device(points=points, vertices=vertices, faces=faces, mesh=mesh)
    distances = _distances_points_to_mesh(points, vertices, faces, single_directional, mesh)
    if distances is None:
        return 0.0
    return _hausdorff(distances, single_directional)

def hausdorff_mesh_to_mesh(
    vertices_a: wp.array[wp.vec3],
    faces_a: wp.array[wp.int32],
    vertices_b: wp.array[wp.vec3],
    faces_b: wp.array[wp.int32],
    *,
    single_directional: bool = False,
    mesh_a: wp.Mesh | None = None,
    mesh_b: wp.Mesh | None = None,
) -> float:
    """
    Hausdorff distance between two triangle meshes.

    Direct port of libigl's ``igl::hausdorff`` (vertex-based overload): the distance
    from every vertex of one mesh to the surface of the other is computed in both
    directions and the overall maximum is returned. This equals
    ``sqrt(max(d_ba, d_ab))`` in libigl, where ``d_*`` are squared distances.

    Parameters
    ----------
    vertices_a, faces_a
        ``(n_vertices_a,)`` vertices and ``(3 * n_faces_a,)`` flat faces of mesh ``A``.
    vertices_b, faces_b
        ``(n_vertices_b,)`` vertices and ``(3 * n_faces_b,)`` flat faces of mesh ``B``.
    single_directional
        If ``True``, return only the directed ``A -> surface(B)`` distance.
    mesh_a, mesh_b
        ``wp.Mesh`` objects already built over mesh ``A`` and mesh ``B``, to spare the BVH build
        each call otherwise pays per mesh queried (worth it when a mesh is measured repeatedly).
        Not checked against the arrays: the answer is for the geometry each holds.
        [`Trimesh.warp_mesh`][ordito.mesh.Trimesh.warp_mesh] is the cached one to pass.

    Returns
    -------
    float
        The (symmetric or directed) Hausdorff distance. ``0.0`` for empty inputs.

    Raises
    ------
    RuntimeError
        If ``vertices_a``, ``faces_a``, ``vertices_b``, ``faces_b``, ``mesh_a`` and ``mesh_b``
        are not all on one device.

    See Also
    --------
    [`chamfer_mesh_to_mesh`][ordito.metrics.chamfer_mesh_to_mesh]
    [`hausdorff_points_to_mesh`][ordito.metrics.hausdorff_points_to_mesh]

    Examples
    --------
    A unit cube inside a cube of side 2: the inner corners are ``0.5`` from the outer surface, the
    outer corners ``sqrt(3) / 2`` from the inner one:

    ```python
    import ordito as od

    vertices_a, faces_a = od.creation.box()
    vertices_b, faces_b = od.creation.box(extents=(2.0, 2.0, 2.0))
    directed = od.metrics.hausdorff_mesh_to_mesh(
        vertices_a, faces_a, vertices_b, faces_b, single_directional=True
    )
    print(directed)  # 0.5
    symmetric = od.metrics.hausdorff_mesh_to_mesh(vertices_a, faces_a, vertices_b, faces_b)
    print(f"{symmetric:.4f}")  # 0.8660
    ```
    """
    require_same_device(
        vertices_a=vertices_a,
        faces_a=faces_a,
        vertices_b=vertices_b,
        faces_b=faces_b,
        mesh_a=mesh_a,
        mesh_b=mesh_b,
    )
    distances = _distances_mesh_to_mesh(
        vertices_a, faces_a, vertices_b, faces_b, single_directional, mesh_a, mesh_b
    )
    if distances is None:
        return 0.0
    return _hausdorff(distances, single_directional)