remesh

triwarp.remesh

Changing a mesh's triangulation: subdividing it, coarsening it, and improving its triangle shapes.

Four families, in decreasing order of how much they rearrange:

• Remeshing. isotropic_remesh runs the Botsch-Kobbelt split / collapse / flip / smooth / reproject loop until every edge is near a target length. It is the only entry point here that does all four of the others' jobs at once.
• Decimation. quadric_decimate collapses edges in quadric-error order to a target face count; cluster_decimate instead welds each voxel of a uniform grid to a single vertex, which is far cheaper and far blunter.
• Edge flipping, which moves no vertex and changes no vertex count: flip_to_delaunay toward the Delaunay criterion, flip_by_objective toward a triangle-shape or flatness objective, and intrinsic_delaunay toward the intrinsic Delaunay triangulation, which flips the connectivity a Laplacian sees without touching the embedding.
• Subdivision, which only ever adds: subdivide (one-to-four splits), subdivide_loop (Loop's approximating scheme), subdivide_to_size and subdivide_region_to_size (until every edge, or every edge of a face region, is under a target length).

isotropic_remesh

isotropic_remesh(
    vertices: array[vec3],
    faces: array[int32],
    *,
    target_length: float | array[float32] | None = None,
    iterations: int = 10,
    feature_angle: float = 30.0,
    split: bool = True,
    collapse: bool = True,
    swap: bool = True,
    smooth: bool = True,
    reproject: bool = True,
    max_deviation: float | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Isotropic explicit remeshing (Botsch-Kobbelt split / collapse / flip / smooth / reproject).

GPU port of the classic incremental isotropic remesher (PyMeshLab's meshing_isotropic_explicit_remeshing, vcglib IsotropicRemeshing): each iteration drives all edge lengths toward target_length by (1) splitting every edge longer than 4/3 * target_length at its midpoint (crack-free, reusing subdivide_to_size); (2) collapsing every edge shorter than 4/5 * target_length (a parallel primitive with full 1-ring locking and a manifold link-condition guard); (3) flipping interior edges toward the ideal vertex valence (6 interior, 4 boundary); (4) tangentially smoothing free vertices (area-equalizing Laplacian projected onto the tangent plane); and (5) reprojecting free vertices back onto the original surface. Feature and boundary structure is preserved: each vertex is classified FREE / CREASE (on a boundary loop or a dihedral crease sharper than feature_angle) / CORNER (feature junction or endpoint), corners are frozen, crease vertices only move along their feature, and feature edges are never flipped. Passing max_deviation relaxes that last guarantee on purpose — see its own entry below.

Inputs are cloned and never mutated.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions on the target device.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer. Whenever a reference wp.Mesh is needed — under reproject, under max_deviation, or with an array target_length, any one of which is enough — it aliases vertices and faces rather than copying them; do not mutate them for the duration of the call.	required
target_length	float | array[float32] | None	Desired edge length. A scalar is the uniform target, defaulting to 1 % of the bounding-box diagonal. A (n_vertices,) wp.float32 array is an adaptive sizing field over the input vertices: every stage then reads its own local target, so the result is fine where the field is small and coarse where it is large. This is the general form of PyMeshLab's adaptive flag — rather than deriving the field from curvature internally, the caller supplies it, which also covers a painted field, a distance-to-feature field, or a field carried from another mesh. Build a curvature-driven one from triwarp.curvature, or an interpolated one from interpolate_from_points.	None
iterations	int	Number of full remeshing passes.	10
feature_angle	float	Dihedral angle in degrees above which an interior edge is treated as a sharp feature (protected from flipping and collapsing across).	30.0
split	bool	Enable/disable each stage of the per-iteration pipeline.	True
collapse	bool	Enable/disable each stage of the per-iteration pipeline.	True
swap	bool	Enable/disable each stage of the per-iteration pipeline.	True
smooth	bool	Enable/disable each stage of the per-iteration pipeline.	True
reproject	bool	Enable/disable each stage of the per-iteration pipeline.	True
max_deviation	float | None	Bound on how far the result may move off the input surface, in model units. At the end of every iteration each vertex further than this from the input surface is pulled straight back toward its own closest point until it is exactly this far, bounding the result's one-sided Hausdorff distance to the input. None (the default) leaves fidelity to reproject alone, as before. This is PyMeshLab's checksurfdist / maxsurfdist pair as a single optional bound. See the Notes for how tightly the bound actually holds. The clamp overrides feature preservation, deliberately, and it is the one thing here that does. It runs over every vertex with no regard for its FREE / CREASE / CORNER classification, because a crease or a corner is exactly the kind of vertex that drifts and that reproject refuses to touch. So a crease or corner vertex that has strayed past the bound is pulled back toward its closest point on the whole input surface, which need not lie on its own feature curve. It moves no further than the bound requires, and a vertex already inside the band is untouched — but if the feature curves must be honoured exactly, leave this None.	None

Returns:

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

Raises:

Type	Description
ValueError	If target_length is non-positive, a sizing field does not have one entry per vertex or holds a non-positive value, or max_deviation is non-positive.
RuntimeError	If vertices, faces and target_length are not all on one device.

See Also

subdivide_to_size split_edges flip_to_delaunay

Notes

Hysteresis (split above 4/3 t, collapse below 4/5 t) keeps split and collapse from fighting. The collapse primitive guarantees manifoldness through the link condition but has no normal-flip guard in this version, relying on the reprojection step to keep free vertices on the original surface.

An adaptive field is re-sampled from the input surface before each stage that reads it, with transfer_onto_vertices, rather than transported through the split / collapse operations. The field is a property of the input geometry, so re-sampling keeps it exact under an arbitrary sequence of operations where transport would accumulate error; the cost is one closest-point query per vertex per stage (two per iteration with both split and collapse on). Inside a single stage the field is transported, because there the correspondence is known exactly: a split midpoint takes the mean of the endpoints it splits, and a collapse compacts the bands alongside the vertices.

A constant field is not quite the scalar path: the two agree on the face buffer exactly and on positions to within float rounding. The gap is float rounding in the threshold alone — the array path forms 4/3 * t per vertex in float32 where the scalar path forms it in Python float64 and narrows once. Pass a scalar when the target is uniform; it is also one closest-point query per stage cheaper.

max_deviation is a positional bound applied per iteration, not a per-operation rejection test: an individual collapse or flip is never vetoed for moving the surface too far, it is the accumulated vertex position that is corrected afterwards. A mesh whose edges must never sweep past the bound mid-iteration needs the flip stage's own gate as well (flip_to_delaunay takes one).

How tightly the bound holds is worth stating, because it is set by wp.mesh_query_point_no_sign rather than by this function. Against that query — the one the clamp is implemented with, and the one reproject has always used — the result is within the bound almost exactly. Against an independent closest-point query, the bound controls deviation proportionally, but at a bound near Warp's own query accuracy it becomes approximate rather than hard, because the two queries can disagree by a small absolute amount and iterating the clamp does not converge further (Warp's answer is a fixed point). Ask for a bound comfortably above the scale of a single-precision closest-point query, or scale the model up.

The remaining limitation, stated because the parameter that used to advertise it is gone: PyMeshLab's selectedonly has no equivalent here, and a region-restricted refinement is subdivide_region_to_size rather than a mode of this function.

Source code in triwarp/remesh.py

def isotropic_remesh(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    target_length: float | wp.array[wp.float32] | None = None,
    iterations: int = 10,
    feature_angle: float = 30.0,
    split: bool = True,
    collapse: bool = True,
    swap: bool = True,
    smooth: bool = True,
    reproject: bool = True,
    max_deviation: float | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Isotropic explicit remeshing (Botsch-Kobbelt split / collapse / flip / smooth / reproject).

    GPU port of the classic incremental isotropic remesher (PyMeshLab's
    ``meshing_isotropic_explicit_remeshing``, vcglib ``IsotropicRemeshing``): each iteration drives
    all edge lengths toward ``target_length`` by (1) **splitting** every edge longer than
    ``4/3 * target_length`` at its midpoint (crack-free, reusing
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size]); (2) **collapsing** every edge shorter
    than ``4/5 * target_length`` (a parallel primitive with full 1-ring locking and a manifold
    link-condition guard); (3) **flipping** interior edges toward the ideal vertex valence (6
    interior, 4 boundary); (4) **tangentially smoothing** free vertices (area-equalizing Laplacian
    projected onto the tangent plane); and (5) **reprojecting** free vertices back onto the original
    surface. Feature and boundary structure is preserved: each vertex is classified FREE / CREASE
    (on a boundary loop or a dihedral crease sharper than ``feature_angle``) / CORNER (feature
    junction or endpoint), corners are frozen, crease vertices only move along their feature, and
    feature edges are never flipped. Passing ``max_deviation`` relaxes that last guarantee on
    purpose — see its own entry below.

    Inputs are cloned and never mutated.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions on the target device.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer. Whenever a reference
        ``wp.Mesh`` is needed — under ``reproject``, under ``max_deviation``, or with an array
        ``target_length``, any one of which is enough — it aliases ``vertices`` and ``faces``
        rather than copying them; do not mutate them for the duration of the call.
    target_length
        Desired edge length. A **scalar** is the uniform target, defaulting to ``1 %`` of the
        bounding-box diagonal. A ``(n_vertices,)`` ``wp.float32`` array is an **adaptive sizing
        field** over the *input* vertices: every stage then reads its own local target, so the
        result is fine where the field is small and coarse where it is large. This is the general
        form of PyMeshLab's ``adaptive`` flag — rather than deriving the field from curvature
        internally, the caller supplies it, which also covers a painted field, a distance-to-feature
        field, or a field carried from another mesh. Build a curvature-driven one from
        [`triwarp.curvature`][triwarp.curvature], or an interpolated one from
        [`interpolate_from_points`][triwarp.interpolation.interpolate_from_points].
    iterations
        Number of full remeshing passes.
    feature_angle
        Dihedral angle in **degrees** above which an interior edge is treated as a sharp feature
        (protected from flipping and collapsing across).
    split, collapse, swap, smooth, reproject
        Enable/disable each stage of the per-iteration pipeline.
    max_deviation
        Bound on how far the result may move off the input surface, in model units. At the end of
        every iteration each vertex further than this from the input surface is pulled straight back
        toward its own closest point until it is exactly this far, bounding the result's one-sided
        Hausdorff distance to the input. ``None`` (the default) leaves fidelity to ``reproject``
        alone, as before. This is PyMeshLab's ``checksurfdist`` / ``maxsurfdist`` pair as a single
        optional bound. See the Notes for how tightly the bound actually holds.

        **The clamp overrides feature preservation, deliberately, and it is the one thing here
        that does.** It runs over *every* vertex with no regard for its FREE / CREASE / CORNER
        classification, because a crease or a corner is exactly the kind of vertex that drifts and
        that ``reproject`` refuses to touch. So a crease or corner vertex that has strayed past the
        bound is pulled back toward its closest point on the whole input surface, which need not
        lie on its own feature curve. It moves no further than the bound requires, and a vertex
        already inside the band is untouched — but if the feature curves must be honoured exactly,
        leave this ``None``.

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

    Raises
    ------
    ValueError
        If ``target_length`` is non-positive, a sizing field does not have one entry per vertex or
        holds a non-positive value, or ``max_deviation`` is non-positive.
    RuntimeError
        If ``vertices``, ``faces`` and ``target_length`` are not all on one device.

    See Also
    --------
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size]
    [`split_edges`][triwarp.remesh.split_edges]
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]

    Notes
    -----
    Hysteresis (split above ``4/3 t``, collapse below ``4/5 t``) keeps split and collapse from
    fighting. The collapse primitive guarantees manifoldness through the link condition but has no
    normal-flip guard in this version, relying on the reprojection step to keep free vertices on the
    original surface.

    An adaptive field is **re-sampled from the input surface** before each stage that reads it, with
    [`transfer_onto_vertices`][triwarp.interpolation.transfer_onto_vertices], rather than
    transported through the split / collapse operations. The field is a property of the input
    geometry, so re-sampling keeps it exact under an arbitrary sequence of operations where
    transport would accumulate error; the cost is one closest-point query per vertex per stage (two
    per iteration with both ``split`` and ``collapse`` on). Inside a single stage the field *is*
    transported, because there the correspondence is known exactly: a split midpoint takes the mean
    of the endpoints it splits, and a collapse compacts the bands alongside the vertices.

    A **constant** field is not quite the scalar path: the two agree on the face buffer exactly and
    on positions to within float rounding. The gap is float rounding in the threshold alone — the
    array path forms ``4/3 * t`` per vertex in ``float32`` where the scalar path forms it in Python
    ``float64`` and narrows once. Pass a scalar when the target is uniform; it is also one
    closest-point query per stage cheaper.

    ``max_deviation`` is a **positional bound applied per iteration**, not a per-operation rejection
    test: an individual collapse or flip is never vetoed for moving the surface too far, it is the
    accumulated vertex position that is corrected afterwards. A mesh whose *edges* must never sweep
    past the bound mid-iteration needs the flip stage's own gate as well
    ([`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay] takes one).

    How tightly the bound holds is worth stating, because it is set by
    ``wp.mesh_query_point_no_sign`` rather than by this function. Against that query — the one the
    clamp is implemented with, and the one ``reproject`` has always used — the result is within the
    bound almost exactly. Against an independent closest-point query, the bound controls deviation
    proportionally, but at a bound near Warp's own query accuracy it becomes approximate rather
    than hard, because the two queries can disagree by a small absolute amount and iterating the
    clamp does not converge further (Warp's answer is a fixed point). Ask for a bound comfortably
    above the scale of a single-precision closest-point query, or scale the model up.

    The remaining limitation, stated because the parameter that used to advertise it is gone:
    PyMeshLab's ``selectedonly`` has no equivalent here, and a region-restricted refinement is
    [`subdivide_region_to_size`][triwarp.remesh.subdivide_region_to_size] rather than a mode of this
    function.
    """
    require_same_device(vertices=vertices, faces=faces, target_length=target_length)
    n_faces = int(faces.shape[0]) // 3
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    current_vertices = wp.clone(vertices)
    current_faces = wp.clone(faces)
    if n_faces == 0 or iterations <= 0:
        return current_vertices, current_faces

    diag = tw.bounds.enclosing_diagonal(vertices)
    sizing_input: wp.array[wp.float32] | None = None
    if target_length is None or isinstance(target_length, int | float):
        target = float(target_length) if target_length is not None else 0.01 * diag
        if target <= 0.0:
            raise ValueError(f"isotropic_remesh requires target_length > 0, got {target}.")
    else:
        field = target_length
        if int(field.shape[0]) != n_vertices:
            raise ValueError(
                f"isotropic_remesh requires one target_length per vertex ({n_vertices}), "
                f"got {int(field.shape[0])}."
            )
        # One readback, on a buffer the caller just built: a non-positive entry makes the split
        # stage diverge (every edge over-long), so it is worth catching here rather than at
        # ``max_iter``.
        smallest = float(tw.reduce.min(field))
        if smallest <= 0.0:
            raise ValueError(
                f"isotropic_remesh requires a positive target_length everywhere, got {smallest}."
            )
        target = float(tw.reduce.mean(field))
        sizing_input = field
    if max_deviation is not None and max_deviation <= 0.0:
        raise ValueError(f"isotropic_remesh requires max_deviation > 0, got {max_deviation}.")
    feature = wp.float32(math.radians(feature_angle))

    # Original surface, built once, for reprojecting free vertices and for the deviation bound
    # (never a 0-triangle mesh).
    original_mesh = None
    if reproject or max_deviation is not None or sizing_input is not None:
        require_nonempty_mesh(faces, "isotropic_remesh")
        # The mesh aliases the caller's buffers and is discarded here, so it needs no copy: the
        # loop below rebinds ``current_vertices`` / ``current_faces`` and never writes ``vertices``.
        original_mesh = wp.Mesh(points=vertices, indices=faces)
    query_radius = max(diag, 1.0)
    clamp_kernel = None
    if max_deviation is not None:
        # Hoisted out of the loop: the generated kernel is cached, but the per-call Python is not.
        clamp_kernel = wp.map(
            kernel_remesh.clamp_to_surface_band,
            current_vertices,
            wp.uint64(0),
            wp.float32(0.0),
            wp.float32(0.0),
            out=wp.empty_like(current_vertices),
            return_kernel=True,
        )

    for _ in range(iterations):
        # The sizing field is re-sampled from the *input* surface immediately before each stage that
        # reads it, because the stage before it changed the vertex set. Two closest-point passes per
        # iteration is the price of keeping the field exact; see this function's Notes.
        if split:
            # Only the *high* band gates a split, and with no sizing field it is one constant at
            # every vertex -- which ``subdivide_to_size`` takes as a scalar. So the uniform path
            # builds neither the two per-vertex band buffers nor the ``low`` band that only the
            # collapse stage below reads.
            split_limit: float | wp.array[wp.float32] = 4.0 / 3.0 * target
            if sizing_input is not None:
                split_limit = _length_bands(
                    _sizing_at(current_vertices, vertices, faces, sizing_input, query_radius),
                    target,
                    int(current_vertices.shape[0]),
                    device,
                )[1]
            current_vertices, current_faces = subdivide_to_size(
                current_vertices, current_faces, split_limit, max_iter=20
            )
        if collapse:
            low, high = _length_bands(
                _sizing_at(current_vertices, vertices, faces, sizing_input, query_radius),
                target,
                int(current_vertices.shape[0]),
                device,
            )
            current_vertices, current_faces = _collapse_pass(
                current_vertices, current_faces, low, high, feature
            )
        if int(current_faces.shape[0]) == 0:
            break
        if swap:
            _valence_flip_pass(current_vertices, current_faces, feature)
        if smooth or reproject:
            codes, _boundary = _classify(current_vertices, current_faces, feature)
            if smooth:
                current_vertices = _smooth_pass(current_vertices, current_faces, codes)
            if reproject and original_mesh is not None:
                current_vertices = _reproject_pass(
                    current_vertices, codes, original_mesh, query_radius
                )
        if clamp_kernel is not None and original_mesh is not None:
            bounded = wp.empty_like(current_vertices)
            wp.launch(
                clamp_kernel,
                dim=int(current_vertices.shape[0]),
                inputs=[
                    current_vertices,
                    wp.uint64(original_mesh.id),
                    wp.float32(max_deviation),
                    wp.float32(query_radius),
                ],
                outputs=[bounded],
                device=device,
            )
            current_vertices = bounded

    return current_vertices, current_faces

cluster_decimate

cluster_decimate(
    vertices: array[vec3],
    faces: array[int32],
    *,
    voxel_size: float | None = None,
    contraction: Literal["average", "closest"] = "average",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Decimate by snapping vertices to a uniform voxel grid and welding each cell to one vertex.

The one decimation scheme that is naturally parallel: there is no priority queue and no sequential dependence anywhere, so the whole thing is a handful of kernel launches whatever the mesh size. Every vertex is binned into a cell of width voxel_size, each occupied cell becomes a single output vertex, faces are remapped onto those, and the faces that collapsed (two or three corners in the same cell) or duplicated are dropped.

Ports MeshLab's meshing_decimation_clustering and Open3D's simplify_vertex_clustering; the grid is anchored half a cell below the bounding box, which is Open3D's convention, so both libraries produce the same cell assignment for the same voxel_size. The vertex output alone is also MeshLab's generate_sampling_clustered_vertex.

This does not preserve topology

Two sheets of the surface that pass within voxel_size of each other get welded together, and a thin feature narrower than a cell disappears. That is the point of the algorithm — it is a resampling, not a simplification — but it means the result can be non-manifold even when the input is not. Use isotropic_remesh when the topology matters and quadric_decimate when a face budget does.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions on the target device.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
voxel_size	float | None	Cell width. Defaults to 1 % of the bounding-box diagonal, matching MeshLab's threshold default of 1 %. Larger cells decimate harder.	None
contraction	Literal['average', 'closest']	How each cell picks its output position: "average" (default) — the mean of the cell's vertices, which is Open3D's SimplificationContraction.Average. Smooths slightly and cannot land off the input's convex hull. "closest" — the input vertex nearest the cell centre, which is MeshLab's 'Closest to center' sampling. Keeps every output vertex on the input surface, so it is the right choice when the positions must stay exact (ties break to the lowest index).	'average'

Returns:

Name	Type	Description
vertices	array[vec3]	One position per occupied cell that still carries a face, compacted from index zero.
faces	array[int32]	Flat 3 * n_faces triangle index buffer, free of collapsed and duplicate faces.

Raises:

Type	Description
ValueError	If voxel_size is not positive, or contraction is not one of the two names.
RuntimeError	If vertices and faces are not all on one device.

See Also

quadric_decimate The other way to simplify: a face budget and a quadric error metric, rather than a voxel size. isotropic_remesh triwarp.repair.remove_duplicated_vertices triwarp.grouping.unique_faces triwarp.voxels.voxel_down_sample

Notes

Cells whose every face collapsed are dropped from the output, where Open3D keeps them as unreferenced vertices. So the face counts agree exactly and the vertex counts can differ by the number of such cells — usually zero, and never in a way that changes the surface.

The binning goes through triwarp.voxels.cell_indices, but the dedup deliberately stays on triwarp.grouping.unique_rows rather than moving onto a NanoVDB grid: a grid numbers its clusters leaf-major, so preserving today's vertex order would need a restoring sort that gives back most of any gain, and it is not worth changing the public output convention for what remains.

Source code in triwarp/remesh.py

def cluster_decimate(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    voxel_size: float | None = None,
    contraction: Literal["average", "closest"] = "average",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Decimate by snapping vertices to a uniform voxel grid and welding each cell to one vertex.

    The one decimation scheme that is *naturally* parallel: there is no priority queue and no
    sequential dependence anywhere, so the whole thing is a handful of kernel launches whatever the
    mesh size. Every vertex is binned into a cell of width ``voxel_size``, each occupied cell
    becomes a single output vertex, faces are remapped onto those, and the faces that collapsed
    (two or three corners in the same cell) or duplicated are dropped.

    Ports MeshLab's ``meshing_decimation_clustering`` and Open3D's ``simplify_vertex_clustering``;
    the grid is anchored half a cell below the bounding box, which is Open3D's convention, so both
    libraries produce the same cell assignment for the same ``voxel_size``. The vertex output alone
    is also MeshLab's ``generate_sampling_clustered_vertex``.

    !!! warning "This does not preserve topology"
        Two sheets of the surface that pass within ``voxel_size`` of each other get welded
        together, and a thin feature narrower than a cell disappears. That is the *point* of the
        algorithm — it is a resampling, not a simplification — but it means the result can be
        non-manifold even when the input is not. Use
        [`isotropic_remesh`][triwarp.remesh.isotropic_remesh] when the topology matters and
        [`quadric_decimate`][triwarp.remesh.quadric_decimate] when a face budget does.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions on the target device.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    voxel_size
        Cell width. Defaults to ``1 %`` of the bounding-box diagonal, matching MeshLab's
        ``threshold`` default of ``1 %``. Larger cells decimate harder.
    contraction
        How each cell picks its output position:

        - ``"average"`` (default) — the mean of the cell's vertices, which is Open3D's
          ``SimplificationContraction.Average``. Smooths slightly and cannot land off the input's
          convex hull.
        - ``"closest"`` — the input vertex nearest the cell centre, which is MeshLab's
          ``'Closest to center'`` sampling. Keeps every output vertex *on* the input surface, so it
          is the right choice when the positions must stay exact (ties break to the lowest index).

    Returns
    -------
    vertices : wp.array[wp.vec3]
        One position per occupied cell that still carries a face, compacted from index zero.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer, free of collapsed and duplicate faces.

    Raises
    ------
    ValueError
        If ``voxel_size`` is not positive, or ``contraction`` is not one of the two names.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`quadric_decimate`][triwarp.remesh.quadric_decimate]
        The other way to simplify: a face budget and a quadric error metric, rather than a voxel
        size.
    [`isotropic_remesh`][triwarp.remesh.isotropic_remesh]
    [`triwarp.repair.remove_duplicated_vertices`][triwarp.repair.remove_duplicated_vertices]
    [`triwarp.grouping.unique_faces`][triwarp.grouping.unique_faces]
    [`triwarp.voxels.voxel_down_sample`][triwarp.voxels.voxel_down_sample]

    Notes
    -----
    Cells whose every face collapsed are dropped from the output, where Open3D keeps them as
    unreferenced vertices. So the face counts agree exactly and the vertex counts can differ by the
    number of such cells — usually zero, and never in a way that changes the surface.

    The binning goes through [`triwarp.voxels.cell_indices`][triwarp.voxels.cell_indices], but the
    *dedup* deliberately stays on [`triwarp.grouping.unique_rows`][triwarp.grouping.unique_rows]
    rather than moving onto a NanoVDB grid: a grid numbers its clusters leaf-major, so preserving
    today's vertex order would need a restoring sort that gives back most of any gain, and it is not
    worth changing the public output convention for what remains.
    """
    require_same_device(vertices=vertices, faces=faces)
    if contraction not in ("average", "closest"):
        raise ValueError(f"contraction must be 'average' or 'closest', got {contraction!r}")

    device = vertices.device
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if n_vertices == 0 or n_faces == 0:
        return wp.clone(vertices), wp.clone(faces)

    # One definition of "the default voxel grid for these points", shared with ``triwarp.voxels``
    # so the two modules cannot drift on the cell size or on Open3D's half-cell anchor.
    voxel_size, origin = tw.voxels.resolve_voxel_grid(
        vertices, voxel_size, caller="cluster_decimate"
    )
    cells = tw.voxels.cell_indices(vertices, voxel_size, origin=origin)
    unique_cells, labels = tw.grouping.unique_rows(cells, return_inverse=True)
    # The unique rows *are* the clusters, so their count is the answer -- a device reduction over
    # ``labels`` plus its readback would re-derive a number ``unique_rows`` has already paid for.
    n_clusters = int(unique_cells.shape[0])

    cluster_vertices = _cluster_positions(
        vertices, labels, n_clusters, origin, voxel_size, contraction
    )
    remapped = tw.array.remap_indices(faces, labels)
    keep_mask = wp.empty(n_faces, dtype=wp.bool, device=device)
    wp.launch(
        kernel_remesh.faces_with_distinct_indices,
        dim=n_faces,
        inputs=[remapped, keep_mask],
        device=device,
    )
    kept_vertices, kept_faces = tw.selection.submesh_from_face_mask(
        cluster_vertices, remapped, keep_mask
    )
    # Welding can map two distinct input faces onto the same triple, which would leave a duplicated
    # face rather than a manifold one, so the dedup is part of the algorithm rather than polish.
    out_vertices, out_faces, _remap = tw.repair.remove_unreferenced_vertices(
        kept_vertices, tw.grouping.unique_faces(kept_faces)
    )
    return out_vertices, out_faces

quadric_decimate

quadric_decimate(
    vertices: array[vec3],
    faces: array[int32],
    *,
    target_faces: int | None = ...,
    target_ratio: float | None = ...,
    feature_angle: float = ...,
    max_iter: int = ...,
    return_index: Literal[False] = False,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

quadric_decimate(
    vertices: array[vec3],
    faces: array[int32],
    *,
    target_faces: int | None = ...,
    target_ratio: float | None = ...,
    feature_angle: float = ...,
    max_iter: int = ...,
    return_index: Literal[True],
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.int32],
    wp.array[wp.int32],
    wp.array[wp.int32],
]

quadric_decimate(
    vertices: array[vec3],
    faces: array[int32],
    *,
    target_faces: int | None = None,
    target_ratio: float | None = None,
    feature_angle: float = 30.0,
    max_iter: int = 100,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[
        wp.array[wp.vec3],
        wp.array[wp.int32],
        wp.array[wp.int32],
        wp.array[wp.int32],
    ]
)

Simplify to a target face count by quadric-error edge collapses (Garland-Heckbert).

The decimation to reach for when the requirement is a face budget rather than an edge length: isotropic_remesh targets a length and cluster_decimate a voxel size, and neither lets a caller ask for "this mesh at 10 % of its triangles". This does, and it is the method every comparable library exposes for the purpose (MeshLab's meshing_decimation_quadric_edge_collapse, igl.decimate, Open3D's simplify_quadric_decimation).

Each vertex accumulates the area-weighted plane quadrics of its incident faces; the cost of collapsing an edge is the residual of the summed quadric at its own minimizer, which is also where the surviving vertex is placed. Cheap collapses are the ones that barely move the surface, so the flat regions go first and the features last — the property that makes this the standard method.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions on the target device. Never mutated.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
target_faces	int | None	Desired face count. Mutually exclusive with target_ratio; exactly one must be given. A value at or above the input count collapses nothing, but still returns an independent copy with the same compaction and provenance every other target gets — see Returns.	None
target_ratio	float | None	Desired face count as a fraction of the input's, so 0.1 is MeshLab's usual "10 %". Must be in (0, 1].	None
feature_angle	float	Dihedral angle in degrees above which an edge is a feature. Feature and boundary structure is preserved exactly as in isotropic_remesh: corners are frozen, a crease vertex only collapses along its own feature, and a crease is never dragged off it.	30.0
max_iter	int	Cap on collapse passes. Each pass commits a conflict-free independent set, so a large reduction needs many; the loop also stops early once the target is met or a pass commits nothing.	100
return_index	bool	If True, also return the two provenance maps below, which is how a per-vertex or per-face attribute survives the decimation.	False

Returns:

Name	Type	Description
vertices	array[vec3]	Simplified vertex positions on vertices.device, compacted from index zero.
faces	array[int32]	Flat 3 * n_faces triangle index buffer. The count is a best effort at target_faces and is not bounded by it: a mesh whose remaining edges all fail the link condition or the normal-flip guard stops above the target — see Notes. Read the returned count rather than assuming it.
vertex_index	array[int32]	Only when return_index is True: length n_input_vertices, the output vertex each input vertex ended up in, or -1 for an input vertex that survives in no output face (one that was already unreferenced). Many-to-one, since that is what a collapse is, so it is the direction a scatter or a segmented reduction wants.
face_index	array[int32]	Only when return_index is True: length n_output_faces, the input face each output face came from -- the same output-to-input direction split_edges uses, so a per-face attribute follows through tw.array.gather. A collapse only deletes faces and never creates one, so every output face has exactly one source.

Raises:

Type	Description
ValueError	If neither or both of target_faces / target_ratio is given, target_faces is negative, or target_ratio is outside (0, 1].
RuntimeError	If vertices and faces are not all on one device.

See Also

cluster_decimate isotropic_remesh triwarp.repair.collapse_small_triangles

Notes

This is a batched-parallel greedy method, not the textbook serial one, and the difference is visible in the output. Textbook QEM pops one edge at a time from a global priority queue, which is inherently sequential. Here each pass scores every edge, ranks the candidates by cost, and commits the cheapest independent set of them — two collapses may commit together only if their closed 1-rings are disjoint. So the sequence of collapses differs from a serial run's and the resulting triangulation is not the same mesh, even though both are driven by the same metric. Compare the two by deviation from the input rather than by equality.

In exchange the quality is competitive with the sequential method despite the different collapse order: committing an independent set spreads the error over the surface where draining a priority queue concentrates it, and a max-norm error measure rewards that.

A pass commits several independent sets against one scoring, not one. A single hashed-key round takes on the order of m / 50 of the candidates, because each winner locks the closed 1-rings of both its endpoints, and rebuilding the geometry between rounds is comparatively expensive. So the pass retires only the candidates the previous round's commits actually invalidated — those whose closed 1-rings touch a collapsed neighbourhood, for which the cached quadric, cost, target position, link condition and normal-flip verdict are the only things that went stale — and runs another round until one finds nothing new.

That round loop runs entirely on device, as one wp.capture_while graph — see _run_collapse_rounds — which removes the host readback that would otherwise happen once per round.

The per-pass rebuild cost is dominated by the number of wrapper calls it issues rather than by the mesh size, which is why it shares its edge grouping (an _EdgeIncidence) with _classify instead of letting each re-derive it, and why _DecimationBuffers goes further still and replays the whole rebuild as one captured graph with fixed-width buffers and live sizes carried in a device array. Read that class before changing anything here.

return_index costs nothing when it is off and next to nothing when it is on. The two provenance maps are folded per pass by two launches and one copy, and the branch that adds them is evaluated when the pass is issued, so the captured graph a CUDA run replays does not even contain it.

Four consequences to plan around:

• The target is reached exactly whenever it is reachable, and it is feature_angle that decides whether it is. A pass is budgeted at half the remaining surplus (an interior collapse removes two faces), shared across its rounds, and the loop stops early when a pass can commit nothing. What stops it is almost always the feature rule rather than the link condition or the normal-flip guard: a surface's own dihedral angles grow as it is coarsened, so past some face count every edge of a smooth mesh is sharper than feature_angle, every vertex becomes a frozen corner, and no collapse is legal at any max_iter. That floor is the parameter working, not a limitation to route around — it is the same rule that keeps a cylinder's rim and a box's creases intact. Raising feature_angle lowers it; at 180 degrees nothing is a feature and the target is reached. Check the returned face count if it matters.
• Every collapse is also checked against a normal-flip guard: an incident face whose normal would turn by more than ~78 degrees vetoes it. That is what keeps the output free of the inverted, self-intersecting triangles an unguarded quadric method produces at high reduction ratios. It is not usually what stops a decimation short, and is deliberately not exposed as a keyword — see COLLAPSE_MIN_NORMAL_DOT in kernels/remesh.py, which records the veto census this claim rests on.
• The independent set is chosen under a hashed lock key rather than by cost rank. That looks like a detail and is not: on a structured mesh both the edge index and the quadric cost are spatially monotone fields, and a monotone key has one local minimum, so either of those keys would commit only a single collapse per pass. See scramble_index in kernels/remesh.py.
• The output is not bit-reproducible on a mesh with tied costs, and never was. The vertex-face incidence CSR is built by an atomic counting scatter, so a row's order varies run to run; where two candidate edges tie on cost, which one the sort keeps varies with it. On a mesh with many tied costs this can move the two-sided Hausdorff distance noticeably between otherwise identical runs, so treat the max-norm as a band, not a value — the mean deviation is far more stable. Compare a change to this function on the mean, or on many repeats.

Source code in triwarp/remesh.py

def quadric_decimate(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    target_faces: int | None = None,
    target_ratio: float | None = None,
    feature_angle: float = 30.0,
    max_iter: int = 100,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.int32]]
):
    """
    Simplify to a target face count by quadric-error edge collapses (Garland-Heckbert).

    The decimation to reach for when the requirement is a **face budget** rather than an edge
    length: [`isotropic_remesh`][triwarp.remesh.isotropic_remesh] targets a length and
    [`cluster_decimate`][triwarp.remesh.cluster_decimate] a voxel size, and neither lets a caller
    ask for "this mesh at 10 % of its triangles". This does, and it is the method every comparable
    library exposes for the purpose (MeshLab's ``meshing_decimation_quadric_edge_collapse``,
    ``igl.decimate``, Open3D's ``simplify_quadric_decimation``).

    Each vertex accumulates the area-weighted plane quadrics of its incident faces; the cost of
    collapsing an edge is the residual of the summed quadric at its own minimizer, which is also
    where the surviving vertex is placed. Cheap collapses are the ones that barely move the surface,
    so the flat regions go first and the features last — the property that makes this the standard
    method.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions on the target device. Never mutated.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    target_faces
        Desired face count. Mutually exclusive with ``target_ratio``; exactly one must be given.
        A value at or above the input count collapses nothing, but still returns an independent
        copy with the same compaction and provenance every other target gets — see Returns.
    target_ratio
        Desired face count as a fraction of the input's, so ``0.1`` is MeshLab's usual "10 %". Must
        be in ``(0, 1]``.
    feature_angle
        Dihedral angle in **degrees** above which an edge is a feature. Feature and boundary
        structure is preserved exactly as in
        [`isotropic_remesh`][triwarp.remesh.isotropic_remesh]: corners are frozen, a crease vertex
        only collapses along its own feature, and a crease is never dragged off it.
    max_iter
        Cap on collapse passes. Each pass commits a conflict-free independent set, so a large
        reduction needs many; the loop also stops early once the target is met or a pass commits
        nothing.
    return_index
        If ``True``, also return the two provenance maps below, which is how a per-vertex or
        per-face attribute survives the decimation.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        Simplified vertex positions on ``vertices.device``, compacted from index zero.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer. The count is a best effort at ``target_faces``
        and is **not** bounded by it: a mesh whose remaining edges all fail the link condition or
        the normal-flip guard stops above the target — see Notes. Read the returned count rather
        than assuming it.
    vertex_index : wp.array[wp.int32]
        Only when ``return_index`` is ``True``: length ``n_input_vertices``, the **output** vertex
        each input vertex ended up in, or ``-1`` for an input vertex that survives in no output face
        (one that was already unreferenced). Many-to-one, since that is what a collapse is, so it is
        the direction a scatter or a segmented reduction wants.
    face_index : wp.array[wp.int32]
        Only when ``return_index`` is ``True``: length ``n_output_faces``, the **input** face each
        output face came from -- the same output-to-input direction
        [`split_edges`][triwarp.remesh.split_edges] uses, so a per-face attribute follows through
        ``tw.array.gather``. A collapse only deletes faces and never creates one, so every output
        face has exactly one source.

    Raises
    ------
    ValueError
        If neither or both of ``target_faces`` / ``target_ratio`` is given, ``target_faces`` is
        negative, or ``target_ratio`` is outside ``(0, 1]``.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`cluster_decimate`][triwarp.remesh.cluster_decimate]
    [`isotropic_remesh`][triwarp.remesh.isotropic_remesh]
    [`triwarp.repair.collapse_small_triangles`][triwarp.repair.collapse_small_triangles]

    Notes
    -----
    **This is a batched-parallel greedy method, not the textbook serial one, and the difference is
    visible in the output.** Textbook QEM pops one edge at a time from a global priority queue,
    which is inherently sequential. Here each pass scores every edge, ranks the candidates by cost,
    and commits the cheapest *independent set* of them — two collapses may commit together only if
    their closed 1-rings are disjoint. So the sequence of collapses differs from a serial run's and
    the resulting triangulation is not the same mesh, even though both are driven by the same
    metric. Compare the two by deviation from the input rather than by equality.

    In exchange the *quality* is competitive with the sequential method despite the different
    collapse order: committing an independent set spreads the error over the surface where draining
    a priority queue concentrates it, and a max-norm error measure rewards that.

    A pass commits **several independent sets against one scoring**, not one. A single hashed-key
    round takes on the order of ``m / 50`` of the candidates, because each winner locks the closed
    1-rings of both its endpoints, and rebuilding the geometry between rounds is comparatively
    expensive. So the pass retires only the candidates the previous round's commits actually
    invalidated — those whose closed 1-rings touch a collapsed neighbourhood, for which the cached
    quadric, cost, target position, link condition and normal-flip verdict are the only things that
    went stale — and runs another round until one finds nothing new.

    That round loop runs **entirely on device**, as one ``wp.capture_while`` graph — see
    ``_run_collapse_rounds`` — which removes the host readback that would otherwise happen once per
    round.

    The per-pass rebuild cost is dominated by the number of wrapper calls it issues rather than by
    the mesh size, which is why it shares its edge grouping (an ``_EdgeIncidence``) with
    ``_classify`` instead of letting each re-derive it, and why ``_DecimationBuffers`` goes further
    still and replays the whole rebuild as one captured graph with fixed-width buffers and live
    sizes carried in a device array. Read that class before changing anything here.

    **``return_index`` costs nothing when it is off and next to nothing when it is on.** The two
    provenance maps are folded per pass by two launches and one copy, and the branch that adds them
    is evaluated when the pass is *issued*, so the captured graph a CUDA run replays does not even
    contain it.

    Four consequences to plan around:

    - **The target is reached exactly whenever it is reachable, and it is ``feature_angle`` that
      decides whether it is.** A pass is budgeted at half the remaining surplus (an interior
      collapse removes two faces), shared across its rounds, and the loop stops early when a pass
      can commit nothing. What stops it is almost always the *feature* rule rather than the link
      condition or the normal-flip guard: a surface's own dihedral angles grow as it is coarsened,
      so past some face count every edge of a smooth mesh is sharper than ``feature_angle``, every
      vertex becomes a frozen corner, and no collapse is legal at any ``max_iter``. **That floor is
      the parameter working, not a limitation to route around** — it is the same rule that keeps a
      cylinder's rim and a box's creases intact. Raising ``feature_angle`` lowers it; at 180 degrees
      nothing is a feature and the target is reached. Check the returned face count if it matters.
    - Every collapse is also checked against a **normal-flip guard**: an incident face whose normal
      would turn by more than ~78 degrees vetoes it. That is what keeps the output free of the
      inverted, self-intersecting triangles an unguarded quadric method produces at high reduction
      ratios. It is **not** usually what stops a decimation short, and is deliberately not exposed
      as a keyword — see ``COLLAPSE_MIN_NORMAL_DOT`` in ``kernels/remesh.py``, which records the
      veto census this claim rests on.
    - The independent set is chosen under a **hashed** lock key rather than by cost rank. That looks
      like a detail and is not: on a structured mesh both the edge index and the quadric cost are
      spatially monotone fields, and a monotone key has one local minimum, so either of those keys
      would commit only a single collapse per pass. See ``scramble_index`` in
      ``kernels/remesh.py``.
    - **The output is not bit-reproducible on a mesh with tied costs, and never was.** The
      vertex-face incidence CSR is built by an atomic counting scatter, so a row's order varies run
      to run; where two candidate edges tie on cost, which one the sort keeps varies with it. On a
      mesh with many tied costs this can move the two-sided Hausdorff distance noticeably between
      otherwise identical runs, so **treat the max-norm as a band, not a value** — the mean
      deviation is far more stable. Compare a change to this function on the mean, or on many
      repeats.
    """
    require_same_device(vertices=vertices, faces=faces)
    n_faces = int(faces.shape[0]) // 3
    target = _resolve_decimation_target(target_faces, target_ratio, n_faces)

    if n_faces == 0 or target >= n_faces:
        # Nothing to collapse, but the *output* contract still holds: vertices compacted from index
        # zero, and ``vertex_index`` reporting -1 for an input vertex no output face references.
        # Returning the buffer verbatim with an identity map would make the shape of the answer
        # depend on whether the target happened to clear the input's face count -- a caller sweeping
        # a ratio would see an already-unreferenced vertex appear and disappear across that
        # boundary. This is the same compaction ``_DecimationBuffers`` runs at the end of every
        # pass, reached through the shared helper so the two cannot drift apart.
        kept_vertices, kept_faces, remap = tw.repair.remove_unreferenced_vertices(vertices, faces)
        if not return_index:
            return kept_vertices, kept_faces
        return kept_vertices, kept_faces, remap, tw.array.arange(n_faces, device=faces.device)

    buffers = _DecimationBuffers(
        vertices, faces, target, wp.float32(math.radians(feature_angle)), track_index=return_index
    )
    for _ in range(max_iter):
        if not buffers.run_pass():
            break
    out_vertices, out_faces, face_source = buffers.result()
    if not return_index:
        return out_vertices, out_faces
    return out_vertices, out_faces, buffers.vertex_index, face_source

flip_to_delaunay

flip_to_delaunay(
    vertices: array[vec3],
    faces: array[int32],
    region: array[bool] | None = None,
    max_angle_change: float | None = None,
    max_deviation: float | None = None,
    critical_aspect_ratio: float = 1000.0,
    max_iter: int = 100,
) -> wp.array[wp.int32]

Improve triangle quality by flipping interior edges toward the Delaunay criterion.

For every interior edge whose two incident faces are both in region, the shared diagonal is flipped when doing so satisfies the local Delone (empty-circumcircle) test — subject to an optional dihedral-angle-change gate and a surface-deviation gate, so the flips never distort the surface. Rim edges (with a face outside the region, or on the mesh boundary) are never flipped. Vertices, face count and region membership are unchanged; only the triangulation of the region is rewritten.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
region	array[bool] | None	Optional length-n_faces wp.bool mask; only edges interior to the True faces are flippable. None treats the whole mesh as flippable.	None
max_angle_change	float | None	Maximum dihedral-angle change (radians) a flip may introduce. None disables the gate.	None
max_deviation	float | None	Maximum surface deviation a flip may introduce. None disables the gate.	None
critical_aspect_ratio	float	Triangle aspect ratio above which the dihedral-angle gate is lifted, so degenerate triangles can still be repaired.	1000.0
max_iter	int	Maximum number of parallel flip passes.	100

Returns:

Type	Description
array[int32]	Flat face buffer with the region re-triangulated, on faces.device (a copy; the input is not modified).

Raises:

Type	Description
ValueError	If region is given and is not length n_faces.
RuntimeError	If vertices, faces and region are not all on one device.

See Also

flip_by_objective The same flip engine with a shape or flatness predicate in front of it instead of the Delone test. subdivide_region_to_size delaunay_triangulation face_adjacency

Source code in triwarp/remesh.py

def flip_to_delaunay(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    region: wp.array[wp.bool] | None = None,
    max_angle_change: float | None = None,
    max_deviation: float | None = None,
    critical_aspect_ratio: float = 1000.0,
    max_iter: int = 100,
) -> wp.array[wp.int32]:
    """
    Improve triangle quality by flipping interior edges toward the Delaunay criterion.

    For every interior edge whose two incident faces are both in ``region``, the shared diagonal
    is flipped when doing so satisfies the local Delone (empty-circumcircle) test — subject to an
    optional dihedral-angle-change gate and a surface-deviation gate, so the flips never distort
    the surface. Rim edges (with a face
    outside the region, or on the mesh boundary) are never flipped. Vertices, face count and
    region membership are unchanged; only the triangulation of the region is rewritten.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    region
        Optional length-``n_faces`` ``wp.bool`` mask; only edges interior to the ``True`` faces
        are flippable. ``None`` treats the whole mesh as flippable.
    max_angle_change
        Maximum dihedral-angle change (radians) a flip may introduce. ``None`` disables the gate.
    max_deviation
        Maximum surface deviation a flip may introduce. ``None`` disables the gate.
    critical_aspect_ratio
        Triangle aspect ratio above which the dihedral-angle gate is lifted, so degenerate
        triangles can still be repaired.
    max_iter
        Maximum number of parallel flip passes.

    Returns
    -------
    wp.array[wp.int32]
        Flat face buffer with the region re-triangulated, on ``faces.device`` (a copy; the
        input is not modified).

    Raises
    ------
    ValueError
        If ``region`` is given and is not length ``n_faces``.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`flip_by_objective`][triwarp.remesh.flip_by_objective]
        The same flip engine with a shape or flatness predicate in front of it instead of the
        Delone test.
    [`subdivide_region_to_size`][triwarp.remesh.subdivide_region_to_size]
    [`delaunay_triangulation`][triwarp.reconstruction.delaunay_triangulation]
    [`face_adjacency`][triwarp.adjacency.face_adjacency]
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    device = faces.device
    setup = _flip_setup(faces, region)
    if setup is None:
        return wp.clone(faces)
    out_faces, n_vertices, region_flags = setup
    mac, mdsq, car = _flip_gates(max_angle_change, max_deviation, critical_aspect_ratio)

    def launch(adjacency, adjacency_edges, unshared, sorted_keys, key_base, out_flip, out_quad):
        wp.launch(
            kernel_remesh.delone_flip_candidates,
            dim=int(adjacency.shape[0]),
            inputs=[
                vertices,
                out_faces,
                adjacency,
                adjacency_edges,
                unshared,
                region_flags,
                sorted_keys,
                key_base,
                mac,
                mdsq,
                car,
                out_flip,
                out_quad,
            ],
            device=device,
        )

    _flip_interior_edges(out_faces, n_vertices, launch, max_iter)
    return out_faces

flip_by_objective

flip_by_objective(
    vertices: array[vec3],
    faces: array[int32],
    *,
    objective: Literal[
        "planarity", "curvature", "t_vertex"
    ] = "planarity",
    region: array[bool] | None = None,
    planar_angle: float = 1.0,
    metric: Literal[
        "radius_ratio", "area_max_side", "mean_ratio"
    ] = "area_max_side",
    aspect_threshold: float = 40.0,
    max_iter: int = 100,
) -> wp.array[wp.int32]

Flip interior edges to optimize triangle shape or surface flatness, instead of the Delone test.

Runs the same parallel flip engine as flip_to_delaunay — independent-set selection over the face adjacency, iterated until no edge is a candidate — with a different predicate at the front. That is the whole port: the machinery was already there, and these two objectives are what MeshLab's meshing_edge_flip_by_planar_optimization and meshing_edge_flip_by_curvature_optimization put in front of it.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions. Never modified — only the triangulation changes.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
objective	Literal['planarity', 'curvature', 't_vertex']	Which predicate to flip on: "planarity" (default) — shape only, surface preserved. A quad is eligible when its two triangles meet within planar_angle of flat, and its diagonal is flipped when doing so raises the metric score of the worse of the two triangles. Because the quad is near-planar to begin with, the rewrite is a retriangulation and not a deformation. This is the one to reach for after any operation that leaves thin triangles across a flat region. "curvature" — flatness, surface changed. The diagonal is flipped whenever the other one bends less, i.e. whenever the dihedral angle across it is smaller. This does move the surface (it chooses between two interpolations of the same four points) and is what makes a coarse triangulation of a curved shape follow its principal directions. planar_angle and metric are ignored. "t_vertex" — slivers only. A quad is eligible only when one of its two triangles is a sliver: its aspect_ratio (circumradius over twice the inradius) exceeds aspect_threshold; the diagonal is then flipped if that improves the worse of the two. A T-vertex — a vertex sitting in the interior of a neighbouring edge — is exactly what produces such a sliver, which is why this is the repair for one; see flip_t_vertices for the wrapper that says so. planar_angle and metric are ignored.	'planarity'
region	array[bool] | None	Optional length-n_faces wp.bool mask; only edges interior to the True faces are flippable. None treats the whole mesh as flippable.	None
planar_angle	float	Planarity tolerance in degrees for objective="planarity": a quad whose dihedral exceeds it is left alone. MeshLab's pthreshold, whose default of 1 is this one. Must be in [0, 180].	1.0
metric	Literal['radius_ratio', 'area_max_side', 'mean_ratio']	Which face_quality measure objective="planarity" maximizes. Only the three larger-is-better shape measures are accepted — "aspect_ratio" runs the other way and "area" is not a shape measure at all. MeshLab's planartype, whose default 'area/max side' is this one.	'area_max_side'
aspect_threshold	float	Sliver threshold for objective="t_vertex": only a quad whose worse triangle has an aspect_ratio above this is eligible. MeshLab's meshing_remove_t_vertices threshold, whose default of 40 is this one. Must be positive.	40.0
max_iter	int	Maximum number of parallel flip passes. Each pass commits a conflict-free independent set, so a mesh needing many local rewrites needs several.	100

Returns:

Type	Description
array[int32]	Flat face buffer with the region re-triangulated, on faces.device (a copy; the input is not modified).

Raises:

Type	Description
ValueError	If objective or metric is unknown, planar_angle is outside [0, 180], or region has the wrong length.
RuntimeError	If vertices, faces and region are not all on one device.

See Also

flip_to_delaunay triwarp.triangles.face_quality isotropic_remesh

Notes

A quad whose two diagonals score equally — every quad of a regular grid — would flip back and forth forever, one pass each way, so a flip must beat the incumbent by a relative 1e-6 rather than merely tie it. That margin is what makes max_iter a safety net rather than the normal stopping condition.

Source code in triwarp/remesh.py

def flip_by_objective(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    objective: Literal["planarity", "curvature", "t_vertex"] = "planarity",
    region: wp.array[wp.bool] | None = None,
    planar_angle: float = 1.0,
    metric: Literal["radius_ratio", "area_max_side", "mean_ratio"] = "area_max_side",
    aspect_threshold: float = 40.0,
    max_iter: int = 100,
) -> wp.array[wp.int32]:
    """
    Flip interior edges to optimize triangle shape or surface flatness, instead of the Delone test.

    Runs the same parallel flip engine as [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay] —
    independent-set selection over the face adjacency, iterated until no edge is a candidate — with
    a different predicate at the front. That is the whole port: the machinery was already there, and
    these two objectives are what MeshLab's ``meshing_edge_flip_by_planar_optimization`` and
    ``meshing_edge_flip_by_curvature_optimization`` put in front of it.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions. Never modified — only the triangulation changes.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    objective
        Which predicate to flip on:

        - ``"planarity"`` (default) — **shape only, surface preserved.** A quad is eligible when its
          two triangles meet within ``planar_angle`` of flat, and its diagonal is flipped when doing
          so raises the ``metric`` score of the *worse* of the two triangles. Because the quad is
          near-planar to begin with, the rewrite is a retriangulation and not a deformation. This is
          the one to reach for after any operation that leaves thin triangles across a flat region.
        - ``"curvature"`` — **flatness, surface changed.** The diagonal is flipped whenever the
          other one bends less, i.e. whenever the dihedral angle across it is smaller. This *does*
          move the surface (it chooses between two interpolations of the same four points) and is
          what makes a coarse triangulation of a curved shape follow its principal directions.
          ``planar_angle`` and ``metric`` are ignored.
        - ``"t_vertex"`` — **slivers only.** A quad is eligible only when one of its two triangles
          is a sliver: its ``aspect_ratio`` (circumradius over twice the inradius) exceeds
          ``aspect_threshold``; the diagonal is then flipped if that improves the worse of the two.
          A T-vertex — a vertex sitting in the interior of a neighbouring edge — is exactly what
          produces such a sliver, which is why this is the repair for one; see
          [`flip_t_vertices`][triwarp.repair.flip_t_vertices] for the wrapper that says so.
          ``planar_angle`` and ``metric`` are ignored.
    region
        Optional length-``n_faces`` ``wp.bool`` mask; only edges interior to the ``True`` faces are
        flippable. ``None`` treats the whole mesh as flippable.
    planar_angle
        Planarity tolerance in **degrees** for ``objective="planarity"``: a quad whose dihedral
        exceeds it is left alone. MeshLab's ``pthreshold``, whose default of ``1`` is this one. Must
        be in ``[0, 180]``.
    metric
        Which [`face_quality`][triwarp.triangles.face_quality] measure ``objective="planarity"``
        maximizes. Only the three larger-is-better shape measures are accepted — ``"aspect_ratio"``
        runs the other way and ``"area"`` is not a shape measure at all. MeshLab's ``planartype``,
        whose default ``'area/max side'`` is this one.
    aspect_threshold
        Sliver threshold for ``objective="t_vertex"``: only a quad whose worse triangle has an
        ``aspect_ratio`` above this is eligible. MeshLab's ``meshing_remove_t_vertices`` threshold,
        whose default of ``40`` is this one. Must be positive.
    max_iter
        Maximum number of parallel flip passes. Each pass commits a conflict-free independent set,
        so a mesh needing many local rewrites needs several.

    Returns
    -------
    wp.array[wp.int32]
        Flat face buffer with the region re-triangulated, on ``faces.device`` (a copy; the input is
        not modified).

    Raises
    ------
    ValueError
        If ``objective`` or ``metric`` is unknown, ``planar_angle`` is outside ``[0, 180]``, or
        ``region`` has the wrong length.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]
    [`triwarp.triangles.face_quality`][triwarp.triangles.face_quality]
    [`isotropic_remesh`][triwarp.remesh.isotropic_remesh]

    Notes
    -----
    A quad whose two diagonals score *equally* — every quad of a regular grid — would flip back and
    forth forever, one pass each way, so a flip must beat the incumbent by a relative ``1e-6``
    rather than merely tie it. That margin is what makes ``max_iter`` a safety net rather than the
    normal stopping condition.
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    if objective not in ("planarity", "curvature", "t_vertex"):
        raise ValueError(
            f"objective must be 'planarity', 'curvature' or 't_vertex', got {objective!r}"
        )
    if aspect_threshold <= 0.0:
        raise ValueError(f"aspect_threshold must be positive, got {aspect_threshold}")
    if metric not in _OBJECTIVE_QUALITY_METRICS:
        raise ValueError(
            f"metric must be one of {list(_OBJECTIVE_QUALITY_METRICS)}, got {metric!r}"
        )
    if not 0.0 <= planar_angle <= 180.0:
        raise ValueError(f"planar_angle must be in [0, 180] degrees, got {planar_angle}")

    device = faces.device
    setup = _flip_setup(faces, region)
    if setup is None:
        return wp.clone(faces)
    out_faces, n_vertices, region_flags = setup

    objective_flag = {
        "planarity": kernel_remesh.OBJECTIVE_PLANARITY,
        "curvature": kernel_remesh.OBJECTIVE_CURVATURE,
        "t_vertex": kernel_remesh.OBJECTIVE_T_VERTEX,
    }[objective]
    metric_flag = tw.triangles._QUALITY_METRICS[metric]
    # The gate is on the dihedral's cosine so the kernel needs no inverse trigonometry.
    planar_cos = wp.float32(math.cos(math.radians(planar_angle)))

    def launch(adjacency, adjacency_edges, unshared, sorted_keys, key_base, out_flip, out_quad):
        wp.launch(
            kernel_remesh.objective_flip_candidates,
            dim=int(adjacency.shape[0]),
            inputs=[
                vertices,
                out_faces,
                adjacency,
                adjacency_edges,
                unshared,
                region_flags,
                sorted_keys,
                key_base,
                objective_flag,
                metric_flag,
                planar_cos,
                wp.float32(aspect_threshold),
                out_flip,
                out_quad,
            ],
            device=device,
        )

    _flip_interior_edges(out_faces, n_vertices, launch, max_iter)
    return out_faces

intrinsic_delaunay

intrinsic_delaunay(
    vertices: array[vec3],
    faces: array[int32],
    epsilon: float = TOLERANCE_MOLLIFY,
    max_iter: int = 100,
) -> tuple[wp.array[wp.int32], twt.Array2dFloat32, int]

Retriangulate to the intrinsic Delaunay triangulation, without moving a vertex.

Flips edges whose two opposite angles sum past pi — exactly the edges whose cotangent weight is negative — until none is left. The flip is intrinsic: the new edge is not a straight line in space but the geodesic across the two triangles, and its length comes from unfolding them into a plane and measuring the other diagonal. The surface, its vertices and its metric are all untouched; only which pairs of vertices count as connected changes, so every operator built from the result is a better-behaved operator for the same geometry.

Its practical effect is that the cotangent weights all become non-negative, which is what a Laplacian needs to satisfy a maximum principle: no spurious extrema, no negative diffusion, far better conditioned solves on a badly-shaped mesh.

Flips run in parallel rounds, each committing a conflict-free independent set (the same engine behind flip_to_delaunay) with the edge-length table carried alongside the connectivity.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer. Not modified.	required
epsilon	float	Mollification margin applied before flipping, relative to the mean edge length: a degenerate triangle has no well-defined angles to test.	TOLERANCE_MOLLIFY
max_iter	int	Cap on the number of parallel flip rounds.	100

Returns:

Name	Type	Description
intrinsic_faces	array[int32]	Length-3 * n_faces connectivity of the intrinsic triangulation, over the same vertices.
edge_lengths	Array2dFloat32	(n_faces, 3) intrinsic edge lengths for those faces, column e opposite corner e.
n_flips	int	How many edges were flipped. Zero means the input was already intrinsically Delaunay.

Raises:

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

Notes

Unlike the other flip passes in this module, this one can create a second edge between two vertices some other edge already connects — intrinsically that is a different geodesic, not a duplicate, and the flip loop tracks connectivity through an incrementally-maintained halfedge twin table rather than a vertex-pair key, so the two never collide. The one thing that still cannot be flipped away is a negative weight on a boundary edge: the Delaunay two-opposite- angles condition has nothing to compare a boundary edge's one incident angle against, so no flip of any kind addresses it. A caller that needs the maximum principle should check the weights it got (cotmatrix_entries_intrinsic) rather than inferring them from convergence.

See Also

robust_laplacian mollify_intrinsic flip_to_delaunay

Source code in triwarp/remesh.py

def intrinsic_delaunay(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    epsilon: float = TOLERANCE_MOLLIFY,
    max_iter: int = 100,
) -> tuple[wp.array[wp.int32], twt.Array2dFloat32, int]:
    """
    Retriangulate to the intrinsic Delaunay triangulation, without moving a vertex.

    Flips edges whose two opposite angles sum past ``pi`` — exactly the edges whose cotangent weight
    is negative — until none is left. The flip is *intrinsic*: the new edge is not a straight line
    in space but the geodesic across the two triangles, and its length comes from unfolding them
    into a plane and measuring the other diagonal. The surface, its vertices and its metric are all
    untouched; only which pairs of vertices count as connected changes, so every operator built from
    the result is a better-behaved operator for the *same* geometry.

    Its practical effect is that the cotangent weights all become non-negative, which is what a
    Laplacian needs to satisfy a maximum principle: no spurious extrema, no negative diffusion, far
    better conditioned solves on a badly-shaped mesh.

    Flips run in parallel rounds, each committing a conflict-free independent set (the same engine
    behind [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]) with the edge-length table carried
    alongside the connectivity.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer. Not modified.
    epsilon
        Mollification margin applied before flipping, relative to the mean edge length: a degenerate
        triangle has no well-defined angles to test.
    max_iter
        Cap on the number of parallel flip rounds.

    Returns
    -------
    intrinsic_faces : wp.array[wp.int32]
        Length-``3 * n_faces`` connectivity of the intrinsic triangulation, over the same vertices.
    edge_lengths : twt.Array2dFloat32
        ``(n_faces, 3)`` intrinsic edge lengths for those faces, column ``e`` opposite corner ``e``.
    n_flips : int
        How many edges were flipped. Zero means the input was already intrinsically Delaunay.

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

    Notes
    -----
    Unlike the other flip passes in this module, this one *can* create a second edge between two
    vertices some other edge already connects — intrinsically that is a different geodesic, not a
    duplicate, and the flip loop tracks connectivity through an incrementally-maintained halfedge
    twin table rather than a vertex-pair key, so the two never collide. The one thing that still
    cannot be flipped away is a negative weight on a *boundary* edge: the Delaunay two-opposite-
    angles condition has nothing to compare a boundary edge's one incident angle against, so no
    flip of any kind addresses it. A caller that needs the maximum principle should check the
    weights it got ([`cotmatrix_entries_intrinsic`][triwarp.laplacian.cotmatrix_entries_intrinsic])
    rather than inferring them from convergence.

    See Also
    --------
    [`robust_laplacian`][triwarp.laplacian.robust_laplacian]
    [`mollify_intrinsic`][triwarp.laplacian.mollify_intrinsic]
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    lengths, _ = mollify_intrinsic(vertices, faces, epsilon=epsilon)
    intrinsic_faces = wp.clone(faces)
    if n_faces == 0:
        return intrinsic_faces, lengths, 0

    n_half = n_faces * 3
    # The one point in this loop where a vertex-pair-keyed adjacency table is trustworthy: the
    # caller's input is still a simplicial complex, so ``_FlipTopology``'s ordinary rebuild can
    # derive it, once. Every later round instead maintains ``twin`` in place, because a flip can
    # make the vertex-pair key ambiguous (see ``kernel_remesh.build_intrinsic_twins``).
    initial = _FlipTopology(intrinsic_faces, n_vertices)
    m0 = initial.rebuild()
    twin = wp.full(n_half, -1, dtype=wp.int32, device=device)
    if m0 > 0:
        wp.launch(
            kernel_remesh.build_intrinsic_twins,
            dim=m0,
            inputs=[intrinsic_faces, initial.adjacency, initial.unshared, twin],
            device=device,
        )
    del initial  # its vertex-pair-keyed tables cannot represent what a flip may do from here on

    flip = wp.empty(n_half, dtype=wp.bool, device=device)
    quad = twt.empty_2d((n_half, 4), wp.int32, device=device)
    new_length = wp.empty(n_half, dtype=wp.float32, device=device)
    neighbors = twt.empty_2d((n_half, 4), wp.int32, device=device)
    face_claim = wp.empty(n_faces, dtype=wp.int32, device=device)
    remap = wp.empty(n_half, dtype=wp.int32, device=device)
    no_remap = wp.zeros(n_half, dtype=wp.bool, device=device)
    count = wp.zeros(1, dtype=wp.int32, device=device)

    total = 0
    for _ in range(max_iter):
        wp.launch(
            kernel_remesh.intrinsic_delaunay_candidates,
            dim=n_half,
            inputs=[intrinsic_faces, lengths, twin, flip, quad, new_length, neighbors],
            device=device,
        )
        # Independent-set selection over just the flipping pair -- ``claim_intrinsic_flips``'s own
        # docstring says why that is enough here, unlike the vertex-pair-keyed flip loops.
        face_claim.fill_(INT32_MAX)
        wp.launch(
            kernel_remesh.claim_intrinsic_flips,
            dim=n_half,
            inputs=[flip, twin, face_claim],
            device=device,
        )
        count.zero_()
        remap.fill_(INT32_MAX)
        no_remap.zero_()
        wp.launch(
            kernel_remesh.commit_intrinsic_flips,
            dim=n_half,
            inputs=[
                flip,
                quad,
                neighbors,
                face_claim,
                new_length,
                intrinsic_faces,
                lengths,
                twin,
                remap,
                no_remap,
                count,
            ],
            device=device,
        )
        wp.launch(
            kernel_remesh.fixup_twin_remap,
            dim=n_half,
            inputs=[remap, no_remap, twin],
            device=device,
        )
        n = int(read_scalar(count, 0))
        total += n
        if n == 0:
            break
    return intrinsic_faces, lengths, total

subdivide

subdivide(
    vertices: array[vec3], faces: array[int32]
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Subdivide a mesh by splitting every face into four triangles.

Each triangle is split by placing a new vertex at the midpoint of each edge. The four child triangles share these midpoints and preserve the original winding order, matching trimesh.remesh.subdivide exactly.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions on the target device.	required
faces	array[int32]	Length-3 * n_faces flat triangle index buffer.	required

Returns:

Type	Description
tuple[array[vec3], array[int32]]	(new_vertices, new_faces) on vertices.device.

Raises:

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

See Also

trimesh.remesh.subdivide

Source code in triwarp/remesh.py

def subdivide(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Subdivide a mesh by splitting every face into four triangles.

    Each triangle is split by placing a new vertex at the midpoint of each
    edge. The four child triangles share these midpoints and preserve the
    original winding order, matching [`trimesh.remesh.subdivide`][] exactly.

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

    Returns
    -------
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
        ``(new_vertices, new_faces)`` on ``vertices.device``.

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

    See Also
    --------
    [`trimesh.remesh.subdivide`][]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_vertices = int(vertices.shape[0])

    n_faces = int(faces.shape[0]) // 3
    if n_faces == 0:
        # Cloned, not aliased: every other entry point in this module returns independent buffers,
        # and a caller that mutates a "subdivided" mesh must not reach back into its own input.
        return wp.clone(vertices), wp.clone(faces)

    unique_edges, inverse = tw.edges.edges_unique(faces, n_vertices=n_vertices)
    n_unique = int(unique_edges.shape[0])

    # Compute midpoint vertex for each unique edge
    out_midpoints = wp.empty(n_unique, dtype=wp.vec3, device=device)
    wp.launch(
        kernel_remesh.compute_midpoints,
        dim=n_unique,
        inputs=[vertices, unique_edges, out_midpoints],
        device=device,
    )

    new_vertices, _ = tw.array.pack_1d_arrays([vertices, out_midpoints])
    return new_vertices, _split_faces_four(faces, inverse, n_vertices)

subdivide_loop

subdivide_loop(
    vertices: array[vec3],
    faces: array[int32],
    *,
    return_operator: Literal[False] = False,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

subdivide_loop(
    vertices: array[vec3],
    faces: array[int32],
    *,
    return_operator: Literal[True],
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.int32],
    wps.BsrMatrix[wp.float32],
]

subdivide_loop(
    vertices: array[vec3],
    faces: array[int32],
    *,
    return_operator: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[
        wp.array[wp.vec3],
        wp.array[wp.int32],
        wps.BsrMatrix[wp.float32],
    ]
)

Subdivide a mesh with one pass of Loop subdivision.

Same 1-to-4 split as subdivide -- identical face table, identical index layout -- but the positions are the smooth Loop stencils rather than midpoints, so the surface is approximated instead of interpolated: original vertices move, and repeated application converges to a C² limit surface (C¹ at irregular vertices).

• Odd (edge) vertices, one per unique edge: 3/8 on each endpoint and 1/8 on each of the two vertices opposite the edge. A boundary edge takes the midpoint instead.
• Even (original) vertices, interior: (1 - n * beta) * v + beta * sum(ring) with Warren's beta -- 3/16 at valence 3, 3/(8n) above -- which is the variant igl.loop uses, rather than Loop's original trigonometric weight.
• Even vertices on a boundary: 3/4 * v plus 1/8 of each of the two neighbours along the boundary, with interior neighbours excluded, so a boundary curve subdivides identically from either side of a seam.

The mesh should be edge-manifold. Where it is not, the stencils are not defined and the fallbacks are conservative rather than arbitrary: an edge with three or more incident faces takes the midpoint rule, and a vertex where one or three-plus boundary edges meet -- or one with no edges at all -- keeps its position.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions on the target device.	required
faces	array[int32]	Length-3 * n_faces flat triangle index buffer.	required
return_operator	bool	If True, also return the sparse interpolation operator P this pass applies, so that new_vertices == P @ vertices and any per-vertex attribute can be carried through the subdivision by the same matrix (see triwarp.interpolation.transfer_through_operator).	False

Returns:

Name	Type	Description
new_vertices	array[vec3]	Positions on vertices.device: the n_vertices relocated originals first and then one vertex per unique edge, so the leading n_vertices rows are the input vertex set displaced -- unlike subdivide, where that prefix is unchanged.
new_faces	array[int32]	Flat 3 * 4 * n_faces triangle index buffer.
operator	BsrMatrix	Only when return_operator is True: the (n_vertices + n_edges, n_vertices) float32 matrix of Loop weights, one row per output vertex in the same layout as new_vertices. Every row sums to 1.

Raises:

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

Notes

Why an operator rather than an index map. Every other topology edit here reports provenance as one int32 per output element, because each output comes from exactly one input. A Loop vertex does not: an odd vertex is an affine combination of four inputs and an even vertex of its whole 1-ring, so no index map can express it, and an attribute cannot otherwise be carried through this function at all. The operator is the honest form, it is the standard prolongation object (kernels/algorithms/multigrid.py builds one for a different purpose), and it costs nothing unless asked for.

The operator is assembled from the same three grids and through the same two weight functions (kernels/remesh.loop_odd_weights / loop_even_weights) that the position kernels use, so the two cannot drift onto different surfaces. It is not used to compute the positions -- those stay two direct kernels, since a bsr_from_triplets build plus a bsr_mv would make every caller pay for the matrix.

Four launches over three grids -- faces, unique edges, vertices -- plus the shared topology. Neither stencil needs an ordered 1-ring: the valence and the ring sum come from an atomic pass over the unique edges, which is the deduplicated neighbour count igl::loop reads off a sorted adjacency list, and the two boundary neighbours are found as the ones joined by boundary edges rather than as the ends of that list.

For several passes, call this repeatedly -- that is what igl.loop's number_of_subdivs does internally, and each pass multiplies the face count by four.

See Also

subdivide subdivide_to_size igl.loop

Source code in triwarp/remesh.py

def subdivide_loop(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], *, return_operator: bool = False
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[wp.array[wp.vec3], wp.array[wp.int32], wps.BsrMatrix[wp.float32]]
):
    """
    Subdivide a mesh with one pass of Loop subdivision.

    Same 1-to-4 split as [`subdivide`][triwarp.remesh.subdivide] -- identical face table, identical
    index layout -- but the positions are the smooth Loop stencils rather than midpoints, so the
    surface is *approximated* instead of interpolated: original vertices move, and repeated
    application converges to a C² limit surface (C¹ at irregular vertices).

    - **Odd (edge) vertices**, one per unique edge: ``3/8`` on each endpoint and ``1/8`` on each of
      the two vertices opposite the edge. A boundary edge takes the midpoint instead.
    - **Even (original) vertices**, interior: ``(1 - n * beta) * v + beta * sum(ring)`` with
      **Warren's** ``beta`` -- ``3/16`` at valence 3, ``3/(8n)`` above -- which is the variant
      ``igl.loop`` uses, rather than Loop's original trigonometric weight.
    - **Even vertices on a boundary**: ``3/4 * v`` plus ``1/8`` of each of the two neighbours along
      the boundary, with interior neighbours excluded, so a boundary curve subdivides identically
      from either side of a seam.

    The mesh should be edge-manifold. Where it is not, the stencils are not defined and the
    fallbacks are conservative rather than arbitrary: an edge with three or more incident faces
    takes the midpoint rule, and a vertex where one or three-plus boundary edges meet -- or one with
    no edges at all -- keeps its position.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions on the target device.
    faces
        Length-``3 * n_faces`` flat triangle index buffer.
    return_operator
        If ``True``, also return the sparse interpolation operator ``P`` this pass applies, so that
        ``new_vertices == P @ vertices`` and **any** per-vertex attribute can be carried through the
        subdivision by the same matrix (see
        [`triwarp.interpolation.transfer_through_operator`][triwarp.interpolation.transfer_through_operator]).

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        Positions on ``vertices.device``: the ``n_vertices`` relocated originals first and then one
        vertex per unique edge, so the leading ``n_vertices`` rows are the input vertex set
        *displaced* -- unlike ``subdivide``, where that prefix is unchanged.
    new_faces : wp.array[wp.int32]
        Flat ``3 * 4 * n_faces`` triangle index buffer.
    operator : warp.sparse.BsrMatrix
        Only when ``return_operator`` is ``True``: the ``(n_vertices + n_edges, n_vertices)``
        ``float32`` matrix of Loop weights, one row per output vertex in the same layout as
        ``new_vertices``. Every row sums to 1.

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

    Notes
    -----
    **Why an operator rather than an index map.** Every other topology edit here reports provenance
    as one ``int32`` per output element, because each output comes from exactly one input. A Loop
    vertex does not: an odd vertex is an affine combination of four inputs and an even vertex of its
    whole 1-ring, so no index map can express it, and an attribute cannot otherwise be carried
    through this function at all. The operator is the honest form, it is the standard prolongation
    object (``kernels/algorithms/multigrid.py`` builds one for a different purpose), and it costs
    nothing unless asked for.

    The operator is assembled from the same three grids and through the same two weight functions
    (``kernels/remesh.loop_odd_weights`` / ``loop_even_weights``) that the position kernels use, so
    the two cannot drift onto different surfaces. It is *not* used to compute the positions -- those
    stay two direct kernels, since a ``bsr_from_triplets`` build plus a ``bsr_mv`` would make every
    caller pay for the matrix.

    Four launches over three grids -- faces, unique edges, vertices -- plus the shared topology.
    Neither stencil needs an ordered 1-ring: the valence and the ring sum come from an atomic pass
    over the *unique* edges, which is the deduplicated neighbour count ``igl::loop`` reads off a
    sorted adjacency list, and the two boundary neighbours are found as the ones joined by boundary
    edges rather than as the ends of that list.

    For several passes, call this repeatedly -- that is what ``igl.loop``'s ``number_of_subdivs``
    does internally, and each pass multiplies the face count by four.

    See Also
    --------
    [`subdivide`][triwarp.remesh.subdivide]
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size]
    ``igl.loop``
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if n_faces == 0:
        # No faces means no edges and no relocation, so the pass is the identity -- but it
        # returns independent buffers all the same, as every other entry point here does.
        if return_operator:
            return (
                wp.clone(vertices),
                wp.clone(faces),
                wps.bsr_identity(n_vertices, wp.float32, device=device),
            )
        return wp.clone(vertices), wp.clone(faces)

    unique_edges, inverse = tw.edges.edges_unique(faces, n_vertices=n_vertices)
    n_unique = int(unique_edges.shape[0])

    # How many faces each edge carries, and the sum of the vertices opposite it.
    edge_opposite_sum = wp.zeros(n_unique, dtype=wp.vec3, device=device)
    edge_face_count = wp.zeros(n_unique, dtype=wp.int32, device=device)
    wp.launch(
        kernel_remesh.loop_edge_opposites,
        dim=n_faces,
        inputs=[vertices, faces, inverse, edge_opposite_sum, edge_face_count],
        device=device,
    )

    valence = wp.zeros(n_vertices, dtype=wp.int32, device=device)
    ring_sum = wp.zeros(n_vertices, dtype=wp.vec3, device=device)
    boundary_count = wp.zeros(n_vertices, dtype=wp.int32, device=device)
    boundary_sum = wp.zeros(n_vertices, dtype=wp.vec3, device=device)
    wp.launch(
        kernel_remesh.loop_vertex_rings,
        dim=n_unique,
        inputs=[
            vertices,
            unique_edges,
            edge_face_count,
            valence,
            ring_sum,
            boundary_count,
            boundary_sum,
        ],
        device=device,
    )

    # One buffer sized for its final use, written through two views: the relocated originals in the
    # prefix and the new edge vertices after them, which is the index layout `_split_faces_four`
    # assumes and the one ``igl.loop`` returns.
    new_vertices = wp.empty(n_vertices + n_unique, dtype=wp.vec3, device=device)
    wp.launch(
        kernel_remesh.loop_even_positions,
        dim=n_vertices,
        inputs=[vertices, valence, ring_sum, boundary_count, boundary_sum],
        outputs=[new_vertices[:n_vertices]],
        device=device,
    )
    wp.launch(
        kernel_remesh.loop_odd_positions,
        dim=n_unique,
        inputs=[vertices, unique_edges, edge_opposite_sum, edge_face_count],
        outputs=[new_vertices[n_vertices:]],
        device=device,
    )
    new_faces = _split_faces_four(faces, inverse, n_vertices)
    if not return_operator:
        return new_vertices, new_faces
    return (
        new_vertices,
        new_faces,
        _loop_operator(
            faces, unique_edges, inverse, edge_face_count, valence, boundary_count, n_vertices
        ),
    )

subdivide_to_size

subdivide_to_size(
    vertices: array[vec3],
    faces: array[int32],
    max_edge: float | array[float32],
    max_iter: int = 10,
    return_index: Literal[False] = False,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

subdivide_to_size(
    vertices: array[vec3],
    faces: array[int32],
    max_edge: float | array[float32],
    max_iter: int = 10,
    *,
    return_index: Literal[True],
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.int32],
    wp.array[wp.int32],
]

subdivide_to_size(
    vertices: array[vec3],
    faces: array[int32],
    max_edge: float | array[float32],
    max_iter: int = 10,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[
        wp.array[wp.vec3],
        wp.array[wp.int32],
        wp.array[wp.int32],
    ]
)

Subdivide a mesh until every edge is at most max_edge long.

Every edge longer than max_edge is bisected at a single shared midpoint, so the two faces on either side stay in sync and a watertight input stays watertight — no T-junctions (cracks) are introduced. Faces already small enough are left untouched. Each pass splits every over-long edge once and re-triangulates the incident faces with per-face templates (1, 2, or 3 split edges; the 2-split quad is cut along its shorter diagonal), iterating until no edge exceeds the threshold, matching trimesh.remesh.subdivide_to_size.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions on the target device.	required
faces	array[int32]	Length-3 * n_faces flat triangle index buffer.	required
max_edge	float | array[float32]	Maximum length of any edge in the result. A scalar gives the uniform target; a (n_vertices,) wp.float32 array is a per-vertex sizing field, and an edge's own target is then the mean of its two endpoints', so the refinement is fine where the field is small and coarse where it is large. The field is extended to each inserted midpoint as the mean of the endpoints it splits, so no resampling is needed between passes and a field that satisfies the target cannot be driven past it by a later pass.	required
max_iter	int	Maximum number of subdivision passes. A ValueError is raised if the mesh still has an over-long edge after this many passes.	10
return_index	bool	If True, also return the source face index of each output face.	False

Returns:

Name	Type	Description
new_vertices	array[vec3]	Refined vertex positions on vertices.device (original vertices first, then the inserted edge midpoints).
new_faces	array[int32]	Flat buffer of the refined faces.
index	array[int32]	Only returned when return_index is True: length n_out_faces, the index of the original face each output face was refined from.

Raises:

Type	Description
ValueError	If any edge is still longer than max_edge after max_iter passes.
RuntimeError	If vertices, faces and max_edge are not all on one device.

See Also

subdivide trimesh.remesh.subdivide_to_size

Source code in triwarp/remesh.py

def subdivide_to_size(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    max_edge: float | wp.array[wp.float32],
    max_iter: int = 10,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32]]
):
    """
    Subdivide a mesh until every edge is at most ``max_edge`` long.

    Every edge longer than ``max_edge`` is bisected at a single shared midpoint,
    so the two faces on either side stay in sync and a watertight input stays
    watertight — no T-junctions (cracks) are introduced. Faces already small
    enough are left untouched. Each pass splits every over-long edge once and
    re-triangulates the incident faces with per-face templates (1, 2, or 3 split
    edges; the 2-split quad is cut along its shorter diagonal), iterating until
    no edge exceeds the threshold, matching
    [`trimesh.remesh.subdivide_to_size`][].

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions on the target device.
    faces
        Length-``3 * n_faces`` flat triangle index buffer.
    max_edge
        Maximum length of any edge in the result. A **scalar** gives the uniform target; a
        ``(n_vertices,)`` ``wp.float32`` array is a per-vertex **sizing field**, and an edge's own
        target is then the mean of its two endpoints', so the refinement is fine where the field is
        small and coarse where it is large. The field is *extended* to each inserted midpoint as the
        mean of the endpoints it splits, so no resampling is needed between passes and a field that
        satisfies the target cannot be driven past it by a later pass.
    max_iter
        Maximum number of subdivision passes. A ``ValueError`` is raised if the
        mesh still has an over-long edge after this many passes.
    return_index
        If ``True``, also return the source face index of each output face.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        Refined vertex positions on ``vertices.device`` (original vertices first,
        then the inserted edge midpoints).
    new_faces : wp.array[wp.int32]
        Flat buffer of the refined faces.
    index : wp.array[wp.int32]
        Only returned when ``return_index`` is ``True``: length ``n_out_faces``,
        the index of the original face each output face was refined from.

    Raises
    ------
    ValueError
        If any edge is still longer than ``max_edge`` after ``max_iter`` passes.
    RuntimeError
        If ``vertices``, ``faces`` and ``max_edge`` are not all on one device.

    See Also
    --------
    [`subdivide`][triwarp.remesh.subdivide]
    [`trimesh.remesh.subdivide_to_size`][]
    """
    require_same_device(vertices=vertices, faces=faces, max_edge=max_edge)
    device = vertices.device
    sizing = max_edge if isinstance(max_edge, wp.array) else None
    max_edge_f = wp.float32(0.0) if sizing is not None else wp.float32(max_edge)

    current_vertices = vertices
    current_faces = faces
    n_faces = int(faces.shape[0]) // 3
    index = tw.array.arange(n_faces, device=device)

    if n_faces == 0:
        if return_index:
            return wp.clone(vertices), wp.clone(faces), index
        return wp.clone(vertices), wp.clone(faces)

    for i in range(max_iter + 1):
        n_vertices = int(current_vertices.shape[0])

        unique_edges, inverse = tw.edges.edges_unique(current_faces, n_vertices=n_vertices)
        m = int(unique_edges.shape[0])
        lengths = tw.edges.edges_unique_length(
            current_vertices, current_faces, unique_edges=unique_edges
        )

        # Flag the edges that are longer than the target length.
        long_mask = wp.empty(m, dtype=wp.bool, device=device)
        if sizing is None:
            wp.map(kernel_array.greater, lengths, max_edge_f, out=long_mask)
        else:
            wp.launch(
                kernel_remesh.mark_edges_over_sizing_field,
                dim=m,
                inputs=[unique_edges, lengths, sizing, long_mask],
                device=device,
            )

        # A sizing field must grow with the vertex buffer, and it is *extended* rather than
        # re-sampled: a midpoint's target is the mean of the endpoints it splits, which is the same
        # value the edge was tested against, so a run of passes cannot drift the field. Computed
        # before the split because it reads the pre-split edge rows.
        next_sizing = (
            None if sizing is None else _extend_sizing_field(sizing, unique_edges, long_mask)
        )

        # ``index`` rides through the split rather than being gathered afterwards, and the "did
        # anything split" test reads the resulting face count rather than reducing the mask: a face
        # count strictly grows when a mask is non-empty and is unchanged when it is empty, so this
        # is exact and costs nothing.
        new_vertices, new_faces, new_index = split_edges(
            current_vertices,
            current_faces,
            long_mask,
            unique_edges=unique_edges,
            inverse=inverse,
            index=index,
            return_index=True,
        )
        # Every edge is short enough: we are done.
        if int(new_faces.shape[0]) == int(current_faces.shape[0]):
            break
        # Ran out of passes with over-long edges still present.
        if i >= max_iter:
            raise ValueError("max_iter exceeded!")
        current_vertices, current_faces, index, sizing = (
            new_vertices,
            new_faces,
            new_index,
            next_sizing,
        )

    # Nothing ever split, so ``current_*`` are still the caller's own buffers. Every entry point in
    # this module returns independent ones -- ``subdivide`` says so in as many words -- and a caller
    # that mutates a "subdivided" mesh must not reach back into its own input. Cloning here rather
    # than up front keeps the path that *did* split free, since ``split_edges`` already handed back
    # fresh buffers there.
    if current_vertices is vertices:
        current_vertices, current_faces = wp.clone(vertices), wp.clone(faces)
    if return_index:
        return current_vertices, current_faces, index
    return current_vertices, current_faces

subdivide_region_to_size

subdivide_region_to_size(
    vertices: array[vec3],
    faces: array[int32],
    region: array[bool],
    max_edge: float,
    max_iter: int = 10,
    max_splits: int | None = None,
    delaunay: bool = True,
    max_angle_change: float | None = math.pi / 6.0,
    max_deviation: float | None = None,
) -> tuple[
    wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]
]

Subdivide only a face region until its edges are at most max_edge long.

The region-restricted form of subdivide_to_size, and the refinement stage of the smooth-patch pipeline: every edge with at least one incident region face and length greater than max_edge is bisected, the incident faces are re-triangulated crack-free (the subdivide_to_size 1/2/3-split templates, so faces outside the region that touch a split edge stay watertight), and — unless disabled — a parallel Delaunay edge-flip pass (flip_to_delaunay) improves the region triangulation after each pass. New vertices are appended after the originals, so the caller derives the new-vertex set as the index range [len(vertices), len(new_vertices)).

Where a sequential implementation would drive this from a longest-edge-first priority queue, splitting is done in parallel passes; max_splits is honoured as a soft budget by keeping only the longest eligible edges of the pass that would exceed it.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
region	array[bool]	Length-n_faces wp.bool mask; only edges touching a True face are refined.	required
max_edge	float	Target maximum edge length inside the region.	required
max_iter	int	Maximum number of subdivision passes.	10
max_splits	int | None	Optional soft cap on the total number of edge splits. None keeps splitting until convergence and raises if max_iter is exhausted first.	None
delaunay	bool	When True (default), interleave and finish with the Delaunay flip pass.	True
max_angle_change	float | None	Dihedral-angle-change gate (radians) for the flip pass (default 30°). None disables the gate.	pi / 6.0
max_deviation	float | None	Surface-deviation gate for the flip pass. None disables it.	None

Returns:

Name	Type	Description
new_vertices	array[vec3]	Original vertices followed by the inserted midpoints, on vertices.device.
new_faces	array[int32]	Flat buffer of the refined faces.
new_region	array[bool]	Length n_out_faces region mask; child faces inherit their parent's membership.

Raises:

Type	Description
ValueError	If region length does not match the face count, or if over-long region edges remain after max_iter passes and max_splits is None.
RuntimeError	If vertices, faces and region are not all on one device.

See Also

subdivide_to_size flip_to_delaunay fill_smooth

Source code in triwarp/remesh.py

def subdivide_region_to_size(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    region: wp.array[wp.bool],
    max_edge: float,
    max_iter: int = 10,
    max_splits: int | None = None,
    delaunay: bool = True,
    max_angle_change: float | None = math.pi / 6.0,
    max_deviation: float | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]]:
    """
    Subdivide only a face region until its edges are at most ``max_edge`` long.

    The region-restricted form of [`subdivide_to_size`][triwarp.remesh.subdivide_to_size], and the
    refinement stage of the smooth-patch pipeline: every edge with at least one incident region
    face and length
    greater than ``max_edge`` is bisected, the incident faces are re-triangulated crack-free
    (the [`subdivide_to_size`][triwarp.remesh.subdivide_to_size] 1/2/3-split templates, so faces
    outside the region that touch a split edge stay watertight), and — unless disabled — a
    parallel Delaunay edge-flip pass ([`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay])
    improves the region triangulation after each pass. New vertices are appended after the
    originals, so the caller derives the new-vertex set as the index range
    ``[len(vertices), len(new_vertices))``.

    Where a sequential implementation would drive this from a longest-edge-first priority queue,
    splitting is done in parallel passes; ``max_splits`` is honoured as a soft budget by keeping
    only the longest eligible edges of the pass that would exceed it.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    region
        Length-``n_faces`` ``wp.bool`` mask; only edges touching a ``True`` face are refined.
    max_edge
        Target maximum edge length inside the region.
    max_iter
        Maximum number of subdivision passes.
    max_splits
        Optional soft cap on the total number of edge splits. ``None`` keeps
        splitting until convergence and raises if ``max_iter`` is exhausted first.
    delaunay
        When ``True`` (default), interleave and finish with the Delaunay flip pass.
    max_angle_change
        Dihedral-angle-change gate (radians) for the flip pass (default 30°). ``None`` disables
        the gate.
    max_deviation
        Surface-deviation gate for the flip pass. ``None`` disables it.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        Original vertices followed by the inserted midpoints, on ``vertices.device``.
    new_faces : wp.array[wp.int32]
        Flat buffer of the refined faces.
    new_region : wp.array[wp.bool]
        Length ``n_out_faces`` region mask; child faces inherit their parent's membership.

    Raises
    ------
    ValueError
        If ``region`` length does not match the face count, or if over-long region edges remain
        after ``max_iter`` passes and ``max_splits`` is ``None``.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size]
    [`flip_to_delaunay`][triwarp.remesh.flip_to_delaunay]
    [`fill_smooth`][triwarp.holes.fill_smooth]
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    device = vertices.device
    n_faces = int(faces.shape[0]) // 3
    if int(region.shape[0]) != n_faces:
        raise ValueError(f"region must have length n_faces={n_faces}, got {int(region.shape[0])}")
    if n_faces == 0:
        return wp.clone(vertices), wp.clone(faces), wp.clone(region)

    max_edge_f = wp.float32(max_edge)
    current_vertices = vertices
    current_faces = faces
    region_flags = tw.array.astype(region, wp.int32)
    splits_done = 0

    for i in range(max_iter + 1):
        n_faces = int(current_faces.shape[0]) // 3
        n_vertices = int(current_vertices.shape[0])

        unique_edges, inverse = tw.edges.edges_unique(current_faces, n_vertices=n_vertices)
        m = int(unique_edges.shape[0])
        lengths = tw.edges.edges_unique_length(
            current_vertices, current_faces, unique_edges=unique_edges
        )

        edge_in_region = wp.zeros(m, dtype=wp.bool, device=device)
        wp.launch(
            kernel_remesh.mark_region_edges,
            dim=3 * n_faces,
            inputs=[region_flags, inverse, edge_in_region],
            device=device,
        )
        long_mask = wp.empty(m, dtype=wp.bool, device=device)
        wp.map(kernel_remesh.long_region_edge, lengths, max_edge_f, edge_in_region, out=long_mask)

        # Only the *count* is wanted here -- for the stopping test, the budget and the running
        # total. ``split_edges`` derives the per-edge vertex slots from ``long_mask`` itself.
        n_long = int(tw.reduce.sum(tw.array.astype(long_mask, wp.int32)))

        if n_long == 0:
            break
        if i >= max_iter:
            if max_splits is None:
                raise ValueError("max_iter exceeded!")
            break

        if max_splits is not None:
            remaining = max_splits - splits_done
            if remaining <= 0:
                break
            if n_long > remaining:
                # It keeps exactly ``remaining`` of the flagged edges, by construction, so the
                # new count needs no second reduction.
                long_mask = _keep_longest_edges(long_mask, lengths, remaining, m, device)
                n_long = remaining

        # The crack-free split itself is [`split_edges`][triwarp.remesh.split_edges], which is
        # exactly what this loop contributes nothing new to: all this function decides is *which*
        # edges (long, and inside the region) and what rides along (``region_flags``, carried
        # through ``index`` so the grown region comes back resolved onto the new faces).
        current_vertices, current_faces, region_flags = split_edges(
            current_vertices,
            current_faces,
            long_mask,
            unique_edges=unique_edges,
            inverse=inverse,
            index=region_flags,
            return_index=True,
        )
        splits_done += n_long

        if delaunay:
            _flip_region_faces(
                current_vertices, current_faces, region_flags, max_angle_change, max_deviation, 8
            )

    # Nothing in the region needed splitting, so ``current_*`` are still the caller's own buffers;
    # see ``subdivide_to_size``'s tail for why that has to be broken here. ``new_region`` is always
    # freshly allocated by ``astype``, so only the two mesh buffers are at stake. This has to run
    # *above* the closing flip pass, not below it: ``_flip_region_faces`` rewrites its face buffer
    # in place, so cloning afterwards would hand back a copy of an already-mutated input.
    if current_vertices is vertices:
        current_vertices, current_faces = wp.clone(vertices), wp.clone(faces)

    if delaunay:
        _flip_region_faces(
            current_vertices, current_faces, region_flags, max_angle_change, max_deviation, 50
        )

    new_region = tw.array.astype(region_flags, wp.bool)
    return current_vertices, current_faces, new_region

refine_region_to_density

refine_region_to_density(
    vertices: array[vec3],
    faces: array[int32],
    region: array[bool],
    *,
    max_iter: int = 10,
    alpha: float = math.sqrt(2.0),
    delaunay: bool = True,
    max_angle_change: float | None = math.radians(30.0),
    max_deviation: float | None = None,
) -> tuple[
    wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]
]

Refine a face region until its sampling matches the surrounding mesh's, not a target length.

The density-driven sibling of subdivide_region_to_size, and the difference is the criterion rather than the mechanism. That one bisects every region edge longer than a single global max_edge; this one gives each vertex a scale attribute -- the average length of the edges incident to it -- and splits a region triangle only while its own scale is coarse relative to its corners'. On a uniformly sampled neighbourhood the two agree; on a graded one they do not, because a single length cannot be right at both ends of the grading.

The rule is Liepa's (see Notes). Writing sigma(v) for the scale attribute, c for a triangle's centroid and sigma(c) for the mean of its three corners' attributes, the triangle is split at c when

``alpha * |c - v_m| > sigma(c)``  and  ``alpha * sigma(c) > sigma(v_m)``

holds for every corner m. The first clause refines; the second is what makes the process terminate at the surrounding sampling rather than at a tolerance. A pass that splits nothing ends the loop.

The split is a 1 -> 3 centroid split, which is what makes this cheap in parallel: the new vertex is interior to the triangle and no edge is divided, so there is nothing to agree with the neighbours about and no crack-free template is needed -- unlike edge bisection, which is why subdivide_region_to_size carries the 1/2/3 split families. One pass is two prefix scans and one kernel.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
region	array[bool]	Length-n_faces wp.bool mask; only True faces are refined. Child faces inherit their parent's membership.	required
max_iter	int	Cap on refinement passes. Unlike subdivide_region_to_size this does not raise when the cap is reached: the criterion is a density match rather than a hard bound, so stopping early leaves a coarser patch and not a wrong one.	10
alpha	float	The criterion's constant, sqrt(2) in the paper. Larger refines further, smaller stops sooner; it is the only real tuning knob here.	sqrt(2.0)
delaunay	bool	When True (default), run the parallel Delone edge-flip pass over the region after each split pass, which is the relaxation step the paper pairs with the criterion.	True
max_angle_change	float | None	Dihedral-angle-change gate for that flip pass (default 30 degrees). None disables the gate.	radians(30.0)
max_deviation	float | None	Surface-deviation gate for the flip pass. None disables it.	None

Returns:

Name	Type	Description
new_vertices	array[vec3]	Positions with the centroid vertices appended after the originals, so the caller derives the new-vertex set as the index range [len(vertices), len(new_vertices)).
new_faces	array[int32]	Flat buffer of the refined faces.
new_region	array[bool]	Length n_out_faces region mask.

Raises:

Type	Description
ValueError	If region length does not match the face count.
RuntimeError	If vertices, faces and region are not all on one device.

See Also

subdivide_region_to_size The same operation driven by a target edge length instead. smoothing.refine_and_smooth_region Where the two criteria are selected between, and what the hole fillers reach through. holes.fill_smooth

Notes

The criterion is section 3 of P. Liepa, "Filling holes in meshes", Eurographics/ACM SIGGRAPH Symposium on Geometry Processing (2003).

The scale attribute is computed once, from the mesh as given, and only extended as vertices are added -- a centroid inherits the mean of its parents'. That is the point of it: it carries the surrounding sampling inward across the patch instead of being re-measured from the increasingly fine triangles it is producing, which would never converge.

Source code in triwarp/remesh.py

def refine_region_to_density(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    region: wp.array[wp.bool],
    *,
    max_iter: int = 10,
    alpha: float = math.sqrt(2.0),
    delaunay: bool = True,
    max_angle_change: float | None = math.radians(30.0),
    max_deviation: float | None = None,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]]:
    """
    Refine a face region until its sampling matches the surrounding mesh's, not a target length.

    The density-driven sibling of
    [`subdivide_region_to_size`][triwarp.remesh.subdivide_region_to_size], and the difference is the
    criterion rather than the mechanism. That one bisects every region edge longer than a single
    global ``max_edge``; this one gives each vertex a **scale attribute** -- the average length of
    the edges incident to it -- and splits a region triangle only while its own scale is coarse
    relative to its corners'. On a uniformly sampled neighbourhood the two agree; on a *graded* one
    they do not, because a single length cannot be right at both ends of the grading.

    The rule is Liepa's (see Notes). Writing ``sigma(v)`` for the scale attribute, ``c`` for a
    triangle's centroid and ``sigma(c)`` for the mean of its three corners' attributes, the triangle
    is split at ``c`` when

        ``alpha * |c - v_m| > sigma(c)``  and  ``alpha * sigma(c) > sigma(v_m)``

    holds for every corner ``m``. The first clause refines; the second is what makes the process
    *terminate at the surrounding sampling* rather than at a tolerance. A pass that splits nothing
    ends the loop.

    The split is a **1 -> 3 centroid split**, which is what makes this cheap in parallel: the new
    vertex is interior to the triangle and no edge is divided, so there is nothing to agree with the
    neighbours about and no crack-free template is needed -- unlike edge bisection, which is why
    ``subdivide_region_to_size`` carries the 1/2/3 split families. One pass is two prefix scans and
    one kernel.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    region
        Length-``n_faces`` ``wp.bool`` mask; only ``True`` faces are refined. Child faces inherit
        their parent's membership.
    max_iter
        Cap on refinement passes. Unlike ``subdivide_region_to_size`` this does **not** raise when
        the cap is reached: the criterion is a density match rather than a hard bound, so stopping
        early leaves a coarser patch and not a wrong one.
    alpha
        The criterion's constant, ``sqrt(2)`` in the paper. Larger refines further, smaller stops
        sooner; it is the only real tuning knob here.
    delaunay
        When ``True`` (default), run the parallel Delone edge-flip pass over the region after each
        split pass, which is the relaxation step the paper pairs with the criterion.
    max_angle_change
        Dihedral-angle-change gate for that flip pass (default 30 degrees). ``None`` disables the
        gate.
    max_deviation
        Surface-deviation gate for the flip pass. ``None`` disables it.

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        Positions with the centroid vertices appended after the originals, so the caller derives the
        new-vertex set as the index range ``[len(vertices), len(new_vertices))``.
    new_faces : wp.array[wp.int32]
        Flat buffer of the refined faces.
    new_region : wp.array[wp.bool]
        Length ``n_out_faces`` region mask.

    Raises
    ------
    ValueError
        If ``region`` length does not match the face count.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`subdivide_region_to_size`][triwarp.remesh.subdivide_region_to_size]
        The same operation driven by a target edge length instead.
    [`smoothing.refine_and_smooth_region`][triwarp.smoothing.refine_and_smooth_region]
        Where the two criteria are selected between, and what the hole fillers reach through.
    [`holes.fill_smooth`][triwarp.holes.fill_smooth]

    Notes
    -----
    The criterion is section 3 of P. Liepa, *"Filling holes in meshes"*, Eurographics/ACM SIGGRAPH
    Symposium on Geometry Processing (2003).

    The scale attribute is computed **once**, from the mesh as given, and only *extended* as
    vertices are added -- a centroid inherits the mean of its parents'. That is the point of it: it
    carries the surrounding sampling inward across the patch instead of being re-measured from the
    increasingly fine triangles it is producing, which would never converge.
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    device = vertices.device
    n_faces = int(faces.shape[0]) // 3
    if int(region.shape[0]) != n_faces:
        raise ValueError(f"region must have length n_faces={n_faces}, got {int(region.shape[0])}")
    if n_faces == 0:
        return wp.clone(vertices), wp.clone(faces), wp.clone(region)

    alpha_f = wp.float32(alpha)
    current_vertices = vertices
    current_faces = faces
    current_region = region
    scale = _vertex_scale_attribute(vertices, faces, region)

    for _ in range(max_iter):
        n_faces = int(current_faces.shape[0]) // 3
        split = wp.empty(n_faces, dtype=wp.int32, device=device)
        wp.launch(
            kernel_remesh.mark_density_splits,
            dim=n_faces,
            inputs=[current_vertices, current_faces, current_region, scale, alpha_f, split],
            device=device,
        )
        split_offsets, n_split = tw.array.counts_to_offsets(split)
        if n_split == 0:
            break

        counts = wp.empty(n_faces, dtype=wp.int32, device=device)
        wp.map(kernel_remesh.face_split_count, split, out=counts)
        face_offsets, n_out_faces = tw.array.counts_to_offsets(counts)

        positions = wp.empty(n_split, dtype=wp.vec3, device=device)
        new_scale = wp.empty(n_split, dtype=wp.float32, device=device)
        out_faces = wp.empty(3 * n_out_faces, dtype=wp.int32, device=device)
        out_region = wp.empty(n_out_faces, dtype=wp.bool, device=device)
        wp.launch(
            kernel_remesh.emit_density_splits,
            dim=n_faces,
            inputs=[
                current_vertices,
                current_faces,
                current_region,
                scale,
                split,
                split_offsets,
                face_offsets,
                wp.int32(int(current_vertices.shape[0])),
                positions,
                new_scale,
                out_faces,
                out_region,
            ],
            device=device,
        )
        current_vertices, _ = tw.array.pack_1d_arrays([current_vertices, positions])
        scale, _ = tw.array.pack_1d_arrays([scale, new_scale])
        current_faces = out_faces
        current_region = out_region

        if delaunay:
            _flip_region_faces(
                current_vertices,
                current_faces,
                tw.array.astype(current_region, wp.int32),
                max_angle_change,
                max_deviation,
                8,
            )

    # No face was dense enough to split, so all three are still the caller's own buffers -- this one
    # aliases ``region`` outright, with no ``astype`` in between. See ``subdivide_to_size``'s tail.
    if current_vertices is vertices:
        current_vertices, current_faces, current_region = (
            wp.clone(vertices),
            wp.clone(faces),
            wp.clone(region),
        )
    return current_vertices, current_faces, current_region

split_edges

split_edges(
    vertices: array[vec3],
    faces: array[int32],
    split_mask: array[bool],
    split_positions: array[vec3] | None = None,
    *,
    unique_edges: Array2dInt32 | None = None,
    inverse: array[int32] | None = None,
    index: array[int32] | None = None,
    return_index: Literal[False] = False,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

split_edges(
    vertices: array[vec3],
    faces: array[int32],
    split_mask: array[bool],
    split_positions: array[vec3] | None = None,
    *,
    unique_edges: Array2dInt32 | None = None,
    inverse: array[int32] | None = None,
    index: array[int32] | None = None,
    return_index: Literal[True],
) -> tuple[
    wp.array[wp.vec3],
    wp.array[wp.int32],
    wp.array[wp.int32],
]

split_edges(
    vertices: array[vec3],
    faces: array[int32],
    split_mask: array[bool],
    split_positions: array[vec3] | None = None,
    *,
    unique_edges: Array2dInt32 | None = None,
    inverse: array[int32] | None = None,
    index: array[int32] | None = None,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[
        wp.array[wp.vec3],
        wp.array[wp.int32],
        wp.array[wp.int32],
    ]
)

Split a chosen set of edges in one crack-free pass, inserting one vertex per edge.

The primitive the whole subdivide_* family is built from, exposed because the choice of edges and the position of the new vertex are the only things that differ between its members: subdivide_to_size iterates this with a length test and midpoints, and triwarp.intersection.split_mesh_with_plane calls it once with the edges a plane crosses and the crossing points. A caller with a different criterion — a curvature threshold, a paint selection, an isovalue — needs no new machinery.

Crack-free means the new vertex of an edge is inserted once and both incident faces reference it, so a watertight input stays watertight and no T-junction is introduced. Each face is re-triangulated by how many of its three edges were split: 1 gives two triangles, 2 gives three (the quad cut along its shorter diagonal), 3 gives the regular 1-to-4 split, and 0 passes through unchanged.

Every return is freshly allocated and independently owned, on the path where nothing was flagged and the mesh comes back unchanged as much as on the splitting one -- including index, which is a copy of the caller's rather than the array they passed in.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions on the target device. Never mutated.	required
faces	array[int32]	Length-3 * n_faces wp.int32 flat triangle index buffer.	required
split_mask	array[bool]	(n_edges,) wp.bool mask over the unique undirected edges in edges_unique order, True for each edge to split.	required
split_positions	array[vec3] | None	Where to put each new vertex, as a (n_split,) wp.vec3 array indexed by the exclusive scan of split_mask — that is, in ascending unique-edge order among the flagged edges, which is where a kernel writing out[offsets[e]] naturally puts them. None uses each edge's midpoint.	None
unique_edges	Array2dInt32 | None	The edges_unique pair for faces, when the caller has already built it to compute split_mask. Both must be given together; either being None rebuilds them.	None
inverse	Array2dInt32 | None	The edges_unique pair for faces, when the caller has already built it to compute split_mask. Both must be given together; either being None rebuilds them.	None
index	array[int32] | None	(n_faces,) wp.int32 per-face values to carry through the split: each output face receives the value of the input face it came from. None means the identity, so return_index then reports provenance into faces. Passing the previous pass's index is how an iterated caller composes provenance without a gather per pass.	None
return_index	bool	If True, also return index resolved onto the output faces (provenance into faces when index is None).	False

Returns:

Name	Type	Description
new_vertices	array[vec3]	Original vertices followed by the n_split inserted ones, in ascending edge order.
new_faces	array[int32]	Flat 3 * m triangle index buffer for the refined mesh.
index	array[int32]	Only when return_index is True: length m, the index into the input faces of the face each output face came from.

Raises:

Type	Description
ValueError	If split_mask does not have one entry per unique edge, split_positions does not have one entry per flagged edge, or index does not have one entry per face.
RuntimeError	If vertices, faces, split_mask, split_positions, unique_edges, inverse and index are not all on one device.

See Also

subdivide_to_size subdivide triwarp.intersection.split_mesh_with_plane triwarp.edges.edges_unique

Examples:

Splitting every edge is the regular 1-to-4 subdivision, so the face count quadruples:

unique_edges, inverse = tw.edges.edges_unique(f)
every_edge = wp.full(int(unique_edges.shape[0]), True, dtype=wp.bool, device=f.device)
fine_v, fine_f = tw.remesh.split_edges(
    v, f, every_edge, unique_edges=unique_edges, inverse=inverse
)
print(int(fine_f.shape[0]) // 3 == 4 * (int(f.shape[0]) // 3))

Source code in triwarp/remesh.py

def split_edges(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    split_mask: wp.array[wp.bool],
    split_positions: wp.array[wp.vec3] | None = None,
    *,
    unique_edges: twt.Array2dInt32 | None = None,
    inverse: wp.array[wp.int32] | None = None,
    index: wp.array[wp.int32] | None = None,
    return_index: bool = False,
) -> (
    tuple[wp.array[wp.vec3], wp.array[wp.int32]]
    | tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.int32]]
):
    """
    Split a chosen set of edges in one crack-free pass, inserting one vertex per edge.

    The primitive the whole ``subdivide_*`` family is built from, exposed because the *choice* of
    edges and the *position* of the new vertex are the only things that differ between its members:
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size] iterates this with a length test and
    midpoints, and
    [`triwarp.intersection.split_mesh_with_plane`][triwarp.intersection.split_mesh_with_plane] calls
    it once with the edges a plane crosses and the crossing points. A caller with a different
    criterion — a curvature threshold, a paint selection, an isovalue — needs no new machinery.

    Crack-free means the new vertex of an edge is inserted **once** and both incident faces
    reference it, so a watertight input stays watertight and no T-junction is introduced. Each face
    is re-triangulated by how many of its three edges were split: 1 gives two triangles, 2 gives
    three (the quad cut along its shorter diagonal), 3 gives the regular 1-to-4 split, and 0 passes
    through unchanged.

    Every return is freshly allocated and independently owned, on the path where nothing was
    flagged and the mesh comes back unchanged as much as on the splitting one -- including
    ``index``, which is a copy of the caller's rather than the array they passed in.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions on the target device. Never mutated.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    split_mask
        ``(n_edges,)`` ``wp.bool`` mask over the **unique undirected edges** in
        [`edges_unique`][triwarp.edges.edges_unique] order, ``True`` for each edge to split.
    split_positions
        Where to put each new vertex, as a ``(n_split,)`` ``wp.vec3`` array indexed by the
        **exclusive scan of** ``split_mask`` — that is, in ascending unique-edge order among the
        flagged edges, which is where a kernel writing ``out[offsets[e]]`` naturally puts them.
        ``None`` uses each edge's midpoint.
    unique_edges, inverse
        The [`edges_unique`][triwarp.edges.edges_unique] pair for ``faces``, when the caller has
        already built it to compute ``split_mask``. Both must be given together; either being
        ``None`` rebuilds them.
    index
        ``(n_faces,)`` ``wp.int32`` per-face values to carry through the split: each output face
        receives the value of the input face it came from. ``None`` means the identity, so
        ``return_index`` then reports provenance into ``faces``. Passing the *previous* pass's index
        is how an iterated caller composes provenance without a gather per pass.
    return_index
        If ``True``, also return ``index`` resolved onto the output faces (provenance into ``faces``
        when ``index`` is ``None``).

    Returns
    -------
    new_vertices : wp.array[wp.vec3]
        Original vertices followed by the ``n_split`` inserted ones, in ascending edge order.
    new_faces : wp.array[wp.int32]
        Flat ``3 * m`` triangle index buffer for the refined mesh.
    index : wp.array[wp.int32]
        Only when ``return_index`` is ``True``: length ``m``, the index into the **input** ``faces``
        of the face each output face came from.

    Raises
    ------
    ValueError
        If ``split_mask`` does not have one entry per unique edge, ``split_positions`` does not have
        one entry per flagged edge, or ``index`` does not have one entry per face.
    RuntimeError
        If ``vertices``, ``faces``, ``split_mask``, ``split_positions``, ``unique_edges``,
        ``inverse`` and ``index`` are not all on one device.

    See Also
    --------
    [`subdivide_to_size`][triwarp.remesh.subdivide_to_size]
    [`subdivide`][triwarp.remesh.subdivide]
    [`triwarp.intersection.split_mesh_with_plane`][triwarp.intersection.split_mesh_with_plane]
    [`triwarp.edges.edges_unique`][triwarp.edges.edges_unique]

    Examples
    --------
    Splitting *every* edge is the regular 1-to-4 subdivision, so the face count quadruples:

    ```python
    unique_edges, inverse = tw.edges.edges_unique(f)
    every_edge = wp.full(int(unique_edges.shape[0]), True, dtype=wp.bool, device=f.device)
    fine_v, fine_f = tw.remesh.split_edges(
        v, f, every_edge, unique_edges=unique_edges, inverse=inverse
    )
    print(int(fine_f.shape[0]) // 3 == 4 * (int(f.shape[0]) // 3))
    ```
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        split_mask=split_mask,
        split_positions=split_positions,
        unique_edges=unique_edges,
        inverse=inverse,
        index=index,
    )
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    # Bound in one expression rather than an ``if`` that reassigns the parameters, so the optional
    # annotations narrow for the type checker without an ``assert``.
    edges, corner_edge = (
        (unique_edges, inverse)
        if unique_edges is not None and inverse is not None
        else tw.edges.edges_unique(faces, n_vertices=n_vertices)
    )
    n_edges = int(edges.shape[0])
    if int(split_mask.shape[0]) != n_edges:
        raise ValueError(
            f"split_mask must have one entry per unique edge ({n_edges}), "
            f"got {int(split_mask.shape[0])}."
        )

    carried = index if index is not None else tw.array.arange(n_faces, device=device)
    if int(carried.shape[0]) != n_faces:
        raise ValueError(
            f"index must have one entry per face ({n_faces}), got {int(carried.shape[0])}."
        )

    # The exclusive scan both counts the split edges and assigns each one its new vertex slot, which
    # is the indexing ``split_positions`` is documented against.
    offsets, n_split = tw.array.counts_to_offsets(tw.array.astype(split_mask, wp.int32))
    if n_split == 0 or n_faces == 0:
        if return_index:
            # ``carried`` is still the caller's own ``index`` buffer when one was supplied, so it
            # is copied for the same reason ``vertices`` and ``faces`` are: every return of this
            # function is independently owned on the no-split path exactly as on the splitting one.
            # When ``index`` was ``None`` the ``arange`` above already allocated it here.
            owned = wp.clone(carried) if carried is index else carried
            return wp.clone(vertices), wp.clone(faces), owned
        return wp.clone(vertices), wp.clone(faces)

    if split_positions is None:
        new_points = wp.empty(n_split, dtype=wp.vec3, device=device)
        wp.launch(
            kernel_remesh.fill_edge_midpoints,
            dim=n_edges,
            inputs=[vertices, edges, split_mask, offsets, new_points],
            device=device,
        )
    else:
        if int(split_positions.shape[0]) != n_split:
            raise ValueError(
                f"split_positions must have one entry per flagged edge ({n_split}), "
                f"got {int(split_positions.shape[0])}."
            )
        new_points = split_positions

    # Appended before face emission so the new indices resolve against one buffer.
    new_vertices, _ = tw.array.pack_1d_arrays([vertices, new_points])

    new_index = wp.empty(n_edges, dtype=wp.int32, device=device)
    wp.launch(
        kernel_remesh.build_midpoint_index,
        dim=n_edges,
        inputs=[split_mask, offsets, wp.int32(n_vertices), new_index],
        device=device,
    )
    face_new = tw.array.gather(new_index, corner_edge).reshape((n_faces, 3))

    # Emit up to four triangles per face into fixed slots, then compact.
    out_faces = twt.empty_2d((n_faces * 4, 3), wp.int32, device=device)
    out_valid = wp.empty(n_faces * 4, dtype=wp.bool, device=device)
    out_slot_index = wp.empty(n_faces * 4, dtype=wp.int32, device=device)
    wp.launch(
        kernel_remesh.emit_size_faces,
        dim=n_faces,
        inputs=[faces, face_new, new_vertices, carried, out_faces, out_valid, out_slot_index],
        device=device,
    )

    kept = tw.array.flatnonzero(out_valid)
    new_faces = tw.array.gather(out_faces, kept).reshape(-1)
    if return_index:
        return new_vertices, new_faces, tw.array.gather(out_slot_index, kept)
    return new_vertices, new_faces