ordito.repair

def make_solid(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    keep_largest: bool = True,
    join_components: bool = False,
    max_iter: int = 10,
    inner_iter: int = 3,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Turn a broken digitised surface into a single watertight solid.

    The composite this module's pieces exist to make possible: debris removed, holes closed,
    degeneracies collapsed and self-intersections cut out and refilled, alternating until nothing is
    left to fix. Every stage is a public function here or in [`ordito.holes`][ordito.holes]; what
    this adds is the order and the loop, which is where the difficulty actually is -- closing a hole
    can create a self-intersection, and cutting one out reopens a hole.

    The stages, in order:

    0. The connectivity repair that does not look like a stage at all until it is missing:
       [`remove_unreferenced_vertices`][ordito.repair.remove_unreferenced_vertices],
       [`make_winding_consistent`][ordito.repair.make_winding_consistent] and
       [`split_non_manifold_vertices`][ordito.repair.split_non_manifold_vertices]. Without it a
       mesh whose defect is a non-manifold *edge* comes back with the right Euler characteristic and
       one component and is still not watertight, because no later stage looks at edge manifoldness.
    1. ``keep_largest`` -> [`remove_small_components`][ordito.repair.remove_small_components], so
       the scan debris goes before anything expensive runs on it.
    2. ``join_components`` ->
       [`holes.join_closest_components`][ordito.holes.join_closest_components], for an input
       whose pieces are meant to be one surface rather than a largest piece plus rubbish. Mutually
       useful with ``keep_largest`` rather than exclusive: keep the big piece *and* weld what is
       left.
    3. [`holes.fill_min_weight`][ordito.holes.fill_min_weight] if any boundary remains.
    4. Up to ``max_iter`` rounds of
       [`remove_degenerate_faces`][ordito.repair.remove_degenerate_faces] and
       [`collapse_small_triangles`][ordito.repair.collapse_small_triangles], then
       [`fix_self_intersections`][ordito.repair.fix_self_intersections] at ``max_iter=inner_iter``.
       The loop exits as soon as a round changes nothing.
    5. A final fill if stage 4 reopened a boundary -- which it routinely does, since cutting an
       intersecting region out is what opens one. Nothing geometric runs after it, and that is not
       an omission: filling a 3-vertex rim produces one sliver, a degeneracy pass deletes the sliver
       and reopens the rim, and the two trade the same faces indefinitely.

    Under ``keep_largest`` the component filter runs **inside** stage 4 as well as at the top, and
    that is not belt and braces: cutting an intersecting band out can disconnect the surface, so the
    extra piece does not exist yet when stage 1 looks.

    !!! warning "It returns the best it managed, not a guarantee"
        There is no success flag, deliberately. Convergence is not guaranteed for any input -- a
        patch can intersect something, and its repair can open another hole -- so on a stubborn mesh
        this returns a *partly* repaired surface rather than looping harder or raising. Ask
        [`validation.is_watertight`][ordito.validation.is_watertight] if the answer matters.

        ``join_components`` is where that bites in practice: welding several shells leaves **one**
        rim spanning all of them, and if that rim is badly non-planar the minimum-weight patch
        across it self-intersects, so stage 4 cuts the patch out and undoes the join. Hemispherical
        bowls joined at coplanar rims converge to one solid; the same bowls tilted relative to each
        other come back as separate shells. Rims far from coplanar want
        [`holes.stitch_loops`][ordito.holes.stitch_loops] or a per-pair
        [`holes.bridge_edges`][ordito.holes.bridge_edges] followed by a targeted fill, not this.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    keep_largest
        Drop every connected component but the one with the most faces, first.
    join_components
        Bridge the remaining open components together instead of leaving them separate.
    max_iter
        Cap on the alternating degeneracy / self-intersection rounds.
    inner_iter
        Cap on the cut-and-refill passes inside each self-intersection repair.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_new_vertices,)`` positions of the repaired mesh, on ``vertices.device``.
    new_faces : wp.array[wp.int32]
        ``(3 * n_new_faces,)`` flat face buffer of the repaired mesh.

    Raises
    ------
    ValueError
        If ``max_iter`` or ``inner_iter`` is negative, or if a face references a vertex index
        ``vertices`` does not cover.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`validation.is_watertight`][ordito.validation.is_watertight]
        The question this does not answer for you.
    [`remove_small_components`][ordito.repair.remove_small_components]
    [`fix_self_intersections`][ordito.repair.fix_self_intersections]
    [`holes.fill_min_weight`][ordito.holes.fill_min_weight]
    [`make_normals_outward`][ordito.repair.make_normals_outward]
        What to run afterwards if the winding has to face outward as well as be consistent.

    Examples
    --------
    An icosphere with three faces missing and a speck of debris beside it comes back as one
    closed sphere:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.icosphere(subdivisions=2)
    debris_vertices, debris_faces = od.creation.tetrahedron()
    broken_vertices = wp.array(
        np.concatenate([vertices.numpy(), 0.1 * debris_vertices.numpy() + 3.0]), dtype=wp.vec3
    )
    broken_faces = wp.array(
        np.concatenate([faces.numpy()[9:], debris_faces.numpy() + vertices.shape[0]]),
        dtype=wp.int32,
    )
    print(od.validation.is_watertight(broken_vertices, broken_faces))  # False
    print(len(od.boundary.boundary_loops(broken_vertices, broken_faces)))  # 3
    solid_vertices, solid_faces = od.repair.make_solid(broken_vertices, broken_faces)
    print(od.validation.is_watertight(solid_vertices, solid_faces))  # True
    print(od.measures.euler_characteristic(solid_faces))  # 2
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    if max_iter < 0 or inner_iter < 0:
        raise ValueError(
            f"max_iter and inner_iter must be non-negative, got {max_iter}, {inner_iter}"
        )

    if faces.size == 0:
        return vertices, faces

    # This is the trust boundary for "a broken digitised surface": every stage below indexes
    # ``vertices[faces]`` without a bound check, the same way every other per-face kernel wrapper in
    # this package does, so an out-of-range index arriving here would otherwise reach the first one
    # silently rather than raising a Python exception.
    require_valid_faces(faces, vertices.size, "make_solid")

    # Stage 0. The reference does this inside its *loader*, which is why it is easy to leave out and
    # why leaving it out is dangerous: skipping it can leave a mesh at the right Euler
    # characteristic and one component while still **not watertight**, because nothing downstream
    # addresses a non-manifold edge.
    vertices, faces, _remap = remove_unreferenced_vertices(vertices, faces)
    # ``make_winding_consistent`` then ``split_non_manifold_vertices`` over one pairing of the
    # halfedges: the flip carries the mates over to the wound faces, which the split reads.
    vertices, faces, _source = _split_corners(vertices, *_wound_with_mates(faces, vertices.size))

    # Every stage below is a deterministic function of the mesh it is handed, and ``version`` counts
    # the changes made to the mesh, so a stage is skipped while the mesh is still one it is known
    # to leave alone: the output of the last run of an idempotent stage (degeneracy removal, the
    # collapse to a fixpoint, the largest-component filter), or a mesh the cut-and-refill found
    # nothing to cut in. Each stage reports a change exactly: the removers only drop faces and
    # compact the vertices those keep, so equal counts are an identical mesh; the fill only adds
    # faces; the cut-and-refill returns its own input when it cuts nothing.
    version = 0
    clean_after = {"degenerate": -1, "collapse": -1, "cut": -1, "components": -1}

    def counts() -> tuple[int, int]:
        return vertices.size, faces.size

    if keep_largest:
        vertices, faces = remove_small_components(vertices, faces, keep_largest=True)
        clean_after["components"] = version
    if join_components:
        faces = od.holes.join_closest_components(vertices, faces)
        version += 1

    # The fill finds the rims itself and returns the faces unchanged (a copy) when there is none, so
    # no separate boundary test runs first.
    before = counts()
    # A component filter has compacted the faces, so their indices are below the vertex count by
    # construction and the fill needs no reduction to bound them; the caller's own faces keep it,
    # as the guard against a negative index.
    faces = od.holes.fill_min_weight(
        vertices, faces, n_vertices=vertices.size if keep_largest else None
    )
    version += int(counts() != before)
    for _ in range(max_iter):
        n_faces_before = faces.size
        for name, remover in (
            ("degenerate", remove_degenerate_faces),
            ("collapse", collapse_small_triangles),
        ):
            if clean_after[name] != version:
                before = counts()
                vertices, faces = remover(vertices, faces)
                version += int(counts() != before)
                clean_after[name] = version
        # The collapse runs on the degeneracy pass's output and keeps only faces passing the same
        # test (its remap applies the predicate on the positions it keeps), so what it returns is
        # also a mesh that pass leaves alone.
        if clean_after["collapse"] == version:
            clean_after["degenerate"] = version
        if clean_after["cut"] != version:
            uncut = faces
            vertices, faces = _cut_and_refill(vertices, faces, 1, inner_iter)
            if faces is uncut:
                clean_after["cut"] = version
            else:
                version += 1
        # Cutting an intersecting band out can *disconnect* the surface, so the component filter has
        # to run again here and not only at the top: a single pass at the start cannot see a
        # component that did not exist yet.
        if keep_largest and clean_after["components"] != version:
            before = counts()
            vertices, faces = remove_small_components(vertices, faces, keep_largest=True)
            if counts() != before:
                # A subset of the faces on the same positions: what the three stages above found
                # nothing to do in still holds nothing for them (no degenerate face, none below the
                # collapse's threshold, which a smaller box only lowers, no crossing pair).
                for name in ("degenerate", "collapse", "cut"):
                    if clean_after[name] == version:
                        clean_after[name] = version + 1
                version += 1
            clean_after["components"] = version
        if faces.size == n_faces_before:
            break

    # Stage 4 opens a rim whenever it cuts an intersecting region out, so the last fill is not a
    # belt-and-braces repeat of stage 3 -- it is what makes the common case come back closed. And
    # nothing geometric may run *after* it: filling a small rim produces one sliver, degeneracy
    # removal deletes that sliver and reopens the rim, and the two then trade the same faces
    # forever -- so the fill goes last.
    before = counts()
    # The faces here are a compacting stage's output, or else the ones the first fill left: it
    # range-checked them whenever it filled a rim, and with none filled there is none here to fill.
    # Either way the bound is the vertex count.
    faces = od.holes.fill_min_weight(vertices, faces, n_vertices=vertices.size)
    version += int(counts() != before)
    if keep_largest and clean_after["components"] != version:
        vertices, faces = remove_small_components(vertices, faces, keep_largest=True)
    return vertices, faces

def remove_unreferenced_vertices(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], *, return_inverse: bool = False
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32]]
    | tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.int32]]
):
    """
    Remove vertices not referenced by any face and remap face indices.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    return_inverse
        If ``True``, also return ``inverse`` with ``new_vertices[inverse]`` sourcing
        ``vertices``.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_referenced,)`` referenced vertices only, compacted from index zero.
    new_faces : wp.array[wp.int32]
        ``(3 * n_faces,)`` face buffer with indices remapped into ``new_vertices``.
    remap : wp.array[wp.int32]
        ``(n_vertices,)`` old-to-new map (``-1`` when unreferenced).
    inverse : wp.array[wp.int32], optional
        ``(n_referenced,)`` source index of each new vertex. Present when
        ``return_inverse=True``.

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

    Examples
    --------
    A stray vertex no face uses is dropped and maps to ``-1``:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    padded = wp.array(np.vstack([vertices.numpy(), [5.0, 5.0, 5.0]]), dtype=wp.vec3)
    new_vertices, new_faces, remap = od.repair.remove_unreferenced_vertices(padded, faces)
    print(padded.shape, new_vertices.shape)  # (9,) (8,)
    print(remap.numpy().tolist())  # [0, 1, 2, 3, 4, 5, 6, 7, -1]
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    device = faces.device
    n_vertices = vertices.size

    # The referenced flags are marked as ``int32`` and scanned in place, and one pass reads the scan
    # to write all three maps and the remapped faces: the tail of the scan is the referenced count,
    # which sizes the outputs and is the one readback.
    remap = _launch.empty(n_vertices, dtype=wp.int32, device=device)
    if n_vertices == 0:
        new_vertices = _launch.empty(0, dtype=wp.vec3, device=device)
        inverse = _launch.empty(0, dtype=wp.int32, device=device)
        new_faces = od.array.remap_indices(faces, remap)
    else:
        inclusive = _launch.zeros(n_vertices, dtype=wp.int32, device=device)
        n_indices = faces.size
        if n_indices > 0:
            _launch.launch(
                kernel_repair.mark_referenced,
                dim=n_indices,
                inputs=[faces, inclusive],
                device=device,
            )
        _launch.array_scan(inclusive, out_array=inclusive, inclusive=True)
        n_referenced = int(read_scalar(inclusive))
        new_vertices = _launch.empty(n_referenced, dtype=wp.vec3, device=device)
        inverse = _launch.empty(n_referenced, dtype=wp.int32, device=device)
        new_faces = _launch.empty(n_indices, dtype=wp.int32, device=device)
        _launch.launch(
            kernel_repair.compact_referenced,
            dim=max(n_vertices, n_indices),
            inputs=[vertices, faces, inclusive, remap, new_vertices, inverse, new_faces],
            device=device,
        )

    if return_inverse:
        return new_vertices, new_faces, remap, inverse
    return new_vertices, new_faces, remap

def remove_duplicated_vertices(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], epsilon: float = 0.0
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.int32]]:
    """
    Merge duplicate vertex positions up to a coordinate tolerance and remap face indices.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    epsilon
        Uniqueness tolerance. Positive values snap coordinates to ``round(v / epsilon)``, so the
        tolerance is absolute and is the one you chose. ``0`` instead groups by a *relative*
        bucket about ``2.4e-4`` wide (see Notes) — pass an explicit ``epsilon`` unless that is
        what you want.

    Returns
    -------
    unique_vertices : wp.array[wp.vec3]
        ``(n_unique,)`` deduplicated vertex positions (first occurrence per equivalence class).
    unique_indices : wp.array[wp.int32]
        ``(n_unique,)`` original indices into ``vertices`` for each output row.
    inverse : wp.array[wp.int32]
        ``(n_vertices,)`` map from each input vertex to its slot in ``unique_vertices``.
    unique_faces : wp.array[wp.int32]
        ``(3 * n_faces,)`` face buffer with indices remapped into ``unique_vertices``.

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

    Notes
    -----
    Both tolerance modes quantize positions and group by cell, so both separate a pair straddling a
    cell boundary however close it is. That is worth knowing about ``epsilon=0`` in particular,
    which is **not** an exact-equality test despite requiring no tolerance: it buckets by the high
    bits of each coordinate's ``float32`` representation, giving a *relative* cell about ``2.4e-4``
    wide. So ``1.0`` and ``1.000244`` merge, while two adjacent ``float32`` values on either side of
    a bucket edge do not. Prefer an explicit ``epsilon`` whenever the tolerance matters; use ``0``
    only to collapse positions that are already bitwise equal, which it does reliably (including
    across ``+0.0`` / ``-0.0``).

    Duplicates that are known from construction rather than measured are better collapsed directly:
    see [`revolve`][ordito.creation.revolve], which derives them from its profile instead of
    hashing positions.

    See Also
    --------
    [`duplicate_vertex_inverse`][ordito.repair.duplicate_vertex_inverse]
    [`hash_vector_rows`][ordito.grouping.hash_vector_rows]
    [`ordito.seams.cut_along_edges`][ordito.seams.cut_along_edges]
        The inverse operation: welding coincident positions back together closes every cut it makes.

    Examples
    --------
    A cube stored as a triangle soup, three private vertices per face, is not watertight until its
    coincident corners are welded:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    soup_vertices = wp.array(vertices.numpy()[faces.numpy()], dtype=wp.vec3)
    soup_faces = wp.array(np.arange(36), dtype=wp.int32)
    print(od.validation.is_watertight(soup_vertices, soup_faces))  # False
    welded, _, inverse, welded_faces = od.repair.remove_duplicated_vertices(
        soup_vertices, soup_faces, epsilon=1e-6
    )
    print(soup_vertices.shape, welded.shape)  # (36,) (8,)
    print(od.validation.is_watertight(welded, welded_faces))  # True
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    # Each class's first occurrence comes with the inverse, read off the grouping's stable sort.
    inverse, unique_indices = _vertex_classes(vertices, epsilon)
    unique_vertices = od.array.gather(vertices, unique_indices)
    unique_faces = od.array.remap_indices(faces, inverse)
    return unique_vertices, unique_indices, inverse, unique_faces

def duplicate_vertex_inverse(vertices: wp.array[wp.vec3], epsilon: float) -> wp.array[wp.int32]:
    """
    Map each vertex to the slot of its coincident-vertex equivalence class.

    The inverse map produced by welding vertices at ``epsilon`` tolerance, without also
    computing the deduplicated vertex/face buffers — useful for remapping per-vertex
    attributes (colors, UVs, ...) to match a [`remove_duplicated_vertices`]
    [ordito.repair.remove_duplicated_vertices] call made with the same ``epsilon``.

    This is the shared equivalence map the dedup remaps by, not a remover -- which is why it lives
    here rather than in [`ordito.grouping`][ordito.grouping] beside
    [`first_occurrence_indices`][ordito.grouping.first_occurrence_indices], whose shape it has. Its
    key is a *tolerance-quantized position*, and choosing that tolerance is mesh-repair policy;
    ``grouping`` is deliberately dtype-generic and geometry-free, so an ``epsilon`` there would be
    the first crack in that contract.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    epsilon
        Uniqueness tolerance, with the same meaning as in
        [`remove_duplicated_vertices`][ordito.repair.remove_duplicated_vertices]: positive values
        snap coordinates to ``round(v / epsilon)``, while ``0`` groups by a *relative* bucket about
        ``2.4e-4`` wide rather than testing for equality.

    Returns
    -------
    wp.array[wp.int32]
        ``(n_vertices,)`` map from each input vertex to its slot in the deduplicated set.

    See Also
    --------
    [`remove_duplicated_vertices`][ordito.repair.remove_duplicated_vertices]
        The remover this map is the shared half of.
    [`hash_vector_rows`][ordito.grouping.hash_vector_rows]
    [`grouping.first_occurrence_indices`][ordito.grouping.first_occurrence_indices]
        The geometry-free counterpart: the same class-to-representative reduction over any key.

    Examples
    --------
    The 36 corners of a cube's triangle soup fall into 8 classes:

    ```python
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    soup_vertices = wp.array(vertices.numpy()[faces.numpy()], dtype=wp.vec3)
    inverse = od.repair.duplicate_vertex_inverse(soup_vertices, 1e-6)
    print(inverse.shape, int(inverse.numpy().max()) + 1)  # (36,) 8
    ```
    """
    return _vertex_classes(vertices, epsilon)[0]

def resolve_duplicated_faces(
    faces: wp.array[wp.int32], *, n_vertices: int | None = None
) -> tuple[wp.array[wp.int32], wp.array[wp.int32]]:
    """
    Resolve duplicated triangles by orientation-aware cancellation rules.

    For each geometric duplicate:
    - equal positive and negative counts: remove all copies;
    - one extra positive copy: keep one positively oriented face;
    - one extra negative copy: keep one negatively oriented face;
    - otherwise raise ``ValueError`` when counts are not orientable.

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    n_vertices
        Optional exclusive bound on the vertex indices, typically the vertex count. Trusted, not
        checked: it spares the reduction that otherwise range-checks ``faces`` and infers the
        bound, and its host readback.

    Returns
    -------
    resolved_faces : wp.array[wp.int32]
        ``(3 * n_kept,)`` flat buffer of kept faces, in libigl's order: one per surviving group,
        the groups ascending lexicographically by their sorted corner indices.
    kept_indices : wp.array[wp.int32]
        ``(n_kept,)`` original face indices into the input ``faces`` buffer, in the same order.

    Raises
    ------
    ValueError
        If a duplicate group's signed count is not orientable, i.e. its positive and negative
        copies differ by more than one, or if ``n_vertices`` is ``None`` and a face index is
        negative.

    Examples
    --------
    A face and its reversed copy cancel, opening a hole; a third copy settles the pair in its own
    orientation:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    rows = faces.numpy().reshape(-1, 3)
    cancelled = wp.array(np.concatenate([rows, rows[:1, ::-1]]).ravel(), dtype=wp.int32)
    resolved, kept = od.repair.resolve_duplicated_faces(cancelled)
    print(resolved.shape[0] // 3, sorted(kept.numpy().tolist())[:3])  # 11 [1, 2, 3]
    tripled = wp.array(np.concatenate([rows, rows[:1, ::-1], rows[:1]]).ravel(), dtype=wp.int32)
    resolved, kept = od.repair.resolve_duplicated_faces(tripled)
    print(resolved.shape[0] // 3)  # 12
    ```
    """
    n_faces = faces.size // 3
    device = faces.device
    if n_faces == 0:
        empty = _launch.empty(0, dtype=wp.int32, device=device)
        return empty, empty

    if n_vertices is None:
        # The range check and the key radix ``grouping.unique_faces`` takes without a bound.
        min_index, max_value = od.reduce.minmax(faces)
        if min_index < 0:
            raise ValueError(f"faces must be non-negative, got a minimum of {min_index}")
        n_vertices = max_value + 1
    # Orientation-free keys ordered as the sorted rows are lexicographically, with the identity
    # payload, sorted, and then each run of equal keys -- one duplicate group, its faces in
    # ascending index -- decided by the thread at its start: no unique table, inverse or per-group
    # statistics are written between the sort and the decision. Three sorted indices below
    # ``n_vertices`` pack below ``n_vertices ** 3``.
    keys = _launch.empty(2 * n_faces, dtype=wp.uint64, device=device)
    order = _launch.empty(2 * n_faces, dtype=wp.int32, device=device)
    _launch.launch(
        kernel_repair.pack_lexicographic_face_keys_and_order,
        dim=n_faces,
        inputs=[faces, wp.uint64(n_vertices), keys, order],
        device=device,
    )
    end_bit = max(1, (min(n_vertices**3, 1 << 64) - 1).bit_length())
    _launch.radix_sort_pairs(keys, order, count=n_faces, end_bit=end_bit)
    # The keep flags at the sorted positions with the error word right behind them, so one readback
    # takes the kept count and the error together; the decided face of each run start beside them,
    # in the payload's upper half, the sort's scratch until now.
    flags = _launch.zeros(n_faces + 1, dtype=wp.int32, device=device)
    keep = _launch.view(order, n_faces, 2 * n_faces)
    _launch.launch(
        kernel_repair.resolve_duplicate_runs,
        dim=n_faces,
        inputs=[faces, keys, order, n_faces],
        outputs=[keep, flags],
        device=device,
    )
    # Scanned in place, a kept group's flag becomes its output row plus one, in ascending key order
    # (the reference's emission order); the last entry is the kept total.
    ranks = odt.as_dense(flags[:n_faces])
    _launch.array_scan(ranks, out_array=ranks, inclusive=True)
    n_kept, error_word = read_values(flags, n_faces - 1, 2)
    if error_word != 0:
        _raise_duplicate_group(faces, keys, order, n_faces - error_word)

    resolved = _launch.empty(3 * n_kept, dtype=wp.int32, device=device)
    kept_wp = _launch.empty(n_kept, dtype=wp.int32, device=device)
    if n_kept > 0:
        _launch.launch(
            kernel_repair.compact_kept_duplicate_runs,
            dim=n_faces,
            inputs=[faces, keep, ranks],
            outputs=[resolved, kept_wp],
            device=device,
        )
    return resolved, kept_wp

def remove_degenerate_faces(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Drop degenerate (zero-area) triangles and reindex, keeping vertex positions unchanged.

    Mirrors ``trimesh.Trimesh.nondegenerate_faces`` + ``update_faces``: a face is degenerate when
    two of its vertices coincide or its three vertices are collinear, detected by
    [`face_nondegenerate_mask`][ordito.triangles.face_nondegenerate_mask] (both triangle altitudes
    exceed the merge tolerance). Surviving faces are unchanged; vertices left unreferenced after the
    drop are removed by the reindexing in
    [`submesh_from_face_mask`][ordito.selection.submesh_from_face_mask].

    Unlike [`collapse_small_triangles`][ordito.repair.collapse_small_triangles], no vertices are
    merged and no edges are collapsed: this only removes faces already degenerate in the input.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_new_vertices,)`` vertices still referenced by a non-degenerate face, compacted from
        index zero.
    new_faces : wp.array[wp.int32]
        ``(3 * n_new_faces,)`` flat buffer of the non-degenerate faces, remapped into
        ``new_vertices``.

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

    See Also
    --------
    [`collapse_small_triangles`][ordito.repair.collapse_small_triangles]
    [`face_nondegenerate_mask`][ordito.triangles.face_nondegenerate_mask]
    [`trimesh.triangles.nondegenerate`][]

    Examples
    --------
    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    with_sliver = wp.array(np.concatenate([faces.numpy(), [0, 1, 1]]), dtype=wp.int32)
    new_vertices, new_faces = od.repair.remove_degenerate_faces(vertices, with_sliver)
    print(with_sliver.shape[0] // 3, new_faces.shape[0] // 3)  # 13 12
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    n_faces = faces.size // 3
    if n_faces == 0:
        return _launch.clone(vertices), _launch.clone(faces)

    # The degeneracy test marks the survivors for the compaction in the pass that decides them.
    keep_mask = _launch.empty(n_faces, dtype=wp.bool, device=faces.device)
    ranks = _launch.zeros(n_faces + vertices.size, dtype=wp.int32, device=faces.device)
    _launch.launch(
        kernel_repair.mark_nondegenerate_faces,
        dim=n_faces,
        inputs=[vertices, faces, keep_mask, ranks],
        device=faces.device,
    )
    new_vertices, new_faces, _vertex_index, _ranks = _compact_marked_submesh(
        vertices, faces, keep_mask, ranks, keep_masked=True, with_index=False
    )
    return new_vertices, new_faces

def remove_non_manifold_faces(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], max_iter: int = 3
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Remove faces touching a non-manifold (>2-incident) edge, leaving the mesh edge-manifold.

    Keeps only faces whose three edges are each used by at most two faces
    ([`edge_manifold_mask`][ordito.validation.edge_manifold_mask]). Dropping a face only lowers the
    share count of its own edges, and every edge used three or more times loses all its faces, so
    one pass leaves every edge used at most twice: a further pass would keep every face.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    max_iter
        Maximum number of removal passes. One pass already reaches the fixpoint, so every value
        from ``1`` gives the same answer; ``0`` removes nothing.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_new_vertices,)`` vertices still referenced after non-manifold faces are dropped,
        compacted from index zero; ``vertices`` itself, not a copy, when no face is dropped.
    new_faces : wp.array[wp.int32]
        ``(3 * n_new_faces,)`` flat buffer of the surviving (edge-manifold) faces; ``faces``
        itself, not a copy, when no face is dropped.

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

    See Also
    --------
    [`remove_degenerate_faces`][ordito.repair.remove_degenerate_faces]
    [`remove_degenerate_and_non_manifold_faces`][ordito.repair.remove_degenerate_and_non_manifold_faces]
    [`edge_manifold_mask`][ordito.validation.edge_manifold_mask]

    Examples
    --------
    A fin glued to one edge of a cube gives that edge three faces; all three go:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    a, b = faces.numpy()[:2]
    fin_vertices = wp.array(np.vstack([vertices.numpy(), [0.0, 0.0, 2.0]]), dtype=wp.vec3)
    fin_faces = wp.array(np.concatenate([faces.numpy(), [a, b, 8]]), dtype=wp.int32)
    print(od.validation.is_edge_manifold(fin_faces))  # False
    new_vertices, new_faces = od.repair.remove_non_manifold_faces(fin_vertices, fin_faces)
    print(new_faces.shape[0] // 3, od.validation.is_edge_manifold(new_faces))  # 10 True
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    # The radix of the edge keys; any bound above every index groups the edges identically, so the
    # vertex count serves every pass over a subset of the faces, as it does for
    # ``remove_degenerate_and_non_manifold_faces``, with no reduction and readback to infer one.
    survivors, keep = _edge_manifold_survivors(faces, None, vertices.size, max_iter)
    if keep is None:
        return vertices, faces
    return od.selection.submesh_from_face_mask(vertices, survivors, keep)

def remove_degenerate_and_non_manifold_faces(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], max_iter: int = 3
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Drop degenerate faces, then faces on a non-manifold edge, compacting the vertices once.

    The same answer as [`remove_degenerate_faces`][ordito.repair.remove_degenerate_faces]
    followed by [`remove_non_manifold_faces`][ordito.repair.remove_non_manifold_faces], byte for
    byte. The manifold test reads only face indices and every compaction preserves the order of
    both the faces and the vertices, so the manifold pass runs over the input's own vertex
    numbering and the referenced vertices are compacted once, at the end.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer, every index in
        ``[0, n_vertices)``.
    max_iter
        Maximum number of non-manifold removal passes, as in
        [`remove_non_manifold_faces`][ordito.repair.remove_non_manifold_faces]: every value from
        ``1`` gives the same answer.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_new_vertices,)`` vertices still referenced by a surviving face, compacted from index
        zero.
    new_faces : wp.array[wp.int32]
        ``(3 * n_new_faces,)`` flat buffer of the surviving faces, remapped into
        ``new_vertices``.

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

    See Also
    --------
    [`remove_degenerate_faces`][ordito.repair.remove_degenerate_faces]
    [`remove_non_manifold_faces`][ordito.repair.remove_non_manifold_faces]

    Examples
    --------
    A cube with a fin on one edge and a zero-area face:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.box()
    a, b = faces.numpy()[:2]
    broken_vertices = wp.array(np.vstack([vertices.numpy(), [0.0, 0.0, 2.0]]), dtype=wp.vec3)
    broken_faces = wp.array(np.concatenate([faces.numpy(), [a, b, 8, 0, 1, 1]]), dtype=wp.int32)
    new_vertices, new_faces = od.repair.remove_degenerate_and_non_manifold_faces(
        broken_vertices, broken_faces
    )
    print(broken_faces.shape[0] // 3, new_faces.shape[0] // 3)  # 14 10
    print(od.validation.is_edge_manifold(new_faces))  # True
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    if faces.size == 0:
        return _launch.clone(vertices), _launch.clone(faces)
    # ``triangles.face_nondegenerate_mask`` with the flags its survivors are compacted from, in the
    # one pass.
    n_faces = faces.size // 3
    keep = _launch.empty(n_faces, dtype=wp.bool, device=faces.device)
    keep_flags = _launch.empty(n_faces, dtype=wp.int32, device=faces.device)
    _launch.launch(
        kernel_repair.nondegenerate_face_flags,
        dim=n_faces,
        inputs=[vertices, faces],
        outputs=[keep, keep_flags],
        device=faces.device,
    )
    survivors, keep = _edge_manifold_survivors(faces, (keep, keep_flags), vertices.size, max_iter)
    assert keep is not None
    return od.selection.submesh_from_face_mask(vertices, survivors, keep)

def remove_small_components(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    keep_largest: bool = False,
    min_faces: int | None = None,
    min_area: float | None = None,
    min_diameter: float | None = None,
    labels: wp.array[wp.int32] | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Drop face-connected components that are too small, by one of four measures of "small".

    Debris -- a stray shell from a scan, a sliver left by a boolean, a component a decimation
    stranded -- is what every repair pipeline removes first, and until now
    [`combine.split`][ordito.combine.split] handed back *every* component and left the caller to
    write the argmax. This is that step, and it stays on the device: the per-component statistic is
    accumulated by one scatter and thresholded by one kernel, so nothing is read back and no
    component is ever materialized as its own mesh.

    Exactly one criterion may be given, because they are four different questions rather than four
    spellings of one and combining them would hide which one rejected a component:

    - ``keep_largest`` keeps the single component with the **most faces** and drops every other, so
      the result always has exactly one component. Face count, not area or diameter: a small dense
      shell outranks a large coarse one. A tie goes to the component whose lowest face index is
      smallest, which makes the choice deterministic rather than dependent on thread order.
    - ``min_faces`` keeps components with at least that many faces -- scale-free, and the measure
      that tracks *how much data* a component carries.
    - ``min_area`` keeps components whose summed triangle area reaches it.
    - ``min_diameter`` keeps components whose axis-aligned bounding-box diagonal reaches it, which
      is the measure that survives a component being a thin sheet of many tiny triangles.

    All three ``min_*`` bounds are **inclusive**, matching the reference implementations: a
    component of exactly the threshold face count, or of exactly the threshold diagonal, survives.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    keep_largest
        Keep only the component with the most faces.
    min_faces
        Minimum face count for a component to be kept.
    min_area
        Minimum summed triangle area for a component to be kept. Under
        [`is_deterministic`][ordito.config.is_deterministic] each component's area is summed
        exactly, so a component at the threshold is kept or dropped the same way on every call.
    min_diameter
        Minimum bounding-box diagonal for a component to be kept.
    labels
        ``(n_faces,)`` face component labels of ``faces``, each a face id naming its component, as
        [`face_connected_component_labels`][ordito.adjacency.face_connected_component_labels]
        returns them
        ([`Trimesh.face_connected_component_labels`][ordito.mesh.Trimesh.face_connected_component_labels]),
        to spare the labelling. Trusted, not checked.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        ``(n_new_vertices,)`` vertices of the surviving components, compacted from index zero.
    new_faces : wp.array[wp.int32]
        ``(3 * n_new_faces,)`` flat buffer of the surviving faces, in their input order.

    Raises
    ------
    ValueError
        If no criterion is given, or if more than one is.
    RuntimeError
        If ``vertices``, ``faces`` and ``labels`` are not all on one device.

    See Also
    --------
    [`combine.split`][ordito.combine.split]
        The whole decomposition, when every component is wanted rather than a subset.
    [`holes.join_closest_components`][ordito.holes.join_closest_components]
        The other answer to a multi-component mesh -- weld the pieces together instead of
        discarding them.
    [`adjacency.face_connected_component_labels`][ordito.adjacency.face_connected_component_labels]
        The labelling this thresholds.
    [`remove_non_manifold_faces`][ordito.repair.remove_non_manifold_faces]

    Notes
    -----
    A component's label is a representative *face* index rather than a dense ``0..k-1`` id, so the
    per-component statistic is an ``n_faces``-long array of which only ``k`` slots are ever written.
    That is deliberate: densifying the labels first would cost a sort and a search per face to save
    an allocation, and it is the allocation that is cheap.

    Examples
    --------
    A small tetrahedron floating beside a sphere is dropped, whether by keeping the largest
    component or by a face-count threshold:

    ```python
    import numpy as np
    import warp as wp

    import ordito as od

    vertices, faces = od.creation.icosphere(subdivisions=2)
    debris_vertices, debris_faces = od.creation.tetrahedron()
    noisy_vertices = wp.array(
        np.concatenate([vertices.numpy(), 0.1 * debris_vertices.numpy() + 3.0]), dtype=wp.vec3
    )
    noisy_faces = wp.array(
        np.concatenate([faces.numpy(), debris_faces.numpy() + vertices.shape[0]]), dtype=wp.int32
    )
    print(noisy_vertices.shape, noisy_faces.shape[0] // 3)  # (166,) 324
    kept_vertices, kept_faces = od.repair.remove_small_components(
        noisy_vertices, noisy_faces, keep_largest=True
    )
    print(kept_vertices.shape, kept_faces.shape[0] // 3)  # (162,) 320
    kept_vertices, kept_faces = od.repair.remove_small_components(
        noisy_vertices, noisy_faces, min_faces=10
    )
    print(kept_vertices.shape, kept_faces.shape[0] // 3)  # (162,) 320
    ```
    """
    require_same_device(vertices=vertices, faces=faces, labels=labels)
    given = (keep_largest, min_faces is not None, min_area is not None, min_diameter is not None)
    if sum(given) != 1:
        raise ValueError(
            "pass exactly one of keep_largest, min_faces, min_area or min_diameter, "
            f"got {sum(given)}"
        )

    device = vertices.device
    n_faces = faces.size // 3
    if n_faces == 0:
        return vertices, faces

    counts_faces = min_area is None and min_diameter is None
    counts = None
    if labels is None and counts_faces:
        # The face-count criteria histogram the labels in the flatten that writes them.
        parents = _face_component_forest(faces, vertices.size, None)
        labels = _launch.empty(n_faces, dtype=wp.int32, device=device)
        counts = _launch.zeros(n_faces, dtype=wp.int32, device=device)
        _launch.launch(
            kernel_repair.flatten_and_count_labels,
            dim=n_faces,
            inputs=[parents, labels, counts],
            device=device,
        )
    elif labels is None:
        labels = od.adjacency.face_connected_component_labels(faces, n_vertices=vertices.size)

    # Each criterion accumulates one per-component table; the others stay ``None`` and are never
    # read. The table is then thresholded inside the submesh extraction's own marking pass, so no
    # separate keep-mask pass runs.
    statistic = corners = best = None
    statistic_sums = _launch.empty_group_sum_pass()
    if min_area is not None:
        criterion = kernel_repair.COMPONENT_MIN_AREA
        statistic, statistic_sums = _component_measures(vertices, faces, labels, "area")
    elif min_diameter is not None:
        criterion = kernel_repair.COMPONENT_MIN_DIAMETER
        # ``+inf`` in all six slots seeds both ends at once: the packing stores the upper corner
        # negated so every update is a ``wp.atomic_min``. Each face then decodes its own
        # component's box as it reads it.
        corners = _launch.full(6 * n_faces, value=math.inf, dtype=wp.float32, device=device)
        _launch.launch(
            kernel_scatter.scatter_group_bounds,
            dim=n_faces,
            inputs=[vertices, faces, labels, corners],
            device=device,
        )
    else:
        if counts is None:
            counts = _launch.zeros(n_faces, dtype=wp.int32, device=device)
            _launch.launch(
                kernel_scatter.count_occurrences,
                dim=n_faces,
                inputs=[labels, counts],
                device=device,
            )
        if min_faces is not None:
            criterion = kernel_repair.COMPONENT_MIN_FACES
        else:
            criterion = kernel_repair.COMPONENT_LARGEST
            # ``-1`` is below every packed key, so the reduction needs no separate seeding pass and
            # the winning label never reaches the host -- the marking pass recomputes its key.
            best = _launch.full(1, -1, dtype=wp.int64, device=device)
            _launch.launch(
                kernel_repair.reduce_largest_group,
                dim=n_faces,
                inputs=[counts, best],
                device=device,
            )

    # ``wp.empty``: the marking pass writes every face's entry.
    keep = _launch.empty(n_faces, dtype=wp.bool, device=device)
    ranks = _launch.zeros(n_faces + vertices.size, dtype=wp.int32, device=device)
    _launch.launch(
        kernel_repair.mark_kept_components,
        dim=n_faces,
        inputs=[
            faces,
            labels,
            criterion,
            counts,
            statistic,
            statistic_sums,
            corners,
            best,
            int(min_faces) if min_faces is not None else 0,
            float(min_area if min_area is not None else (min_diameter or 0.0)),
            keep,
            ranks,
        ],
        device=device,
    )
    new_vertices, new_faces, _vertex_index, _ranks = _compact_marked_submesh(
        vertices, faces, keep, ranks, keep_masked=True, with_index=False
    )
    return new_vertices, new_faces