Blame for vli/simplify.py

"""The simplification ladder: how far a mesh can be reduced and still read as its word.

Rungs run from the full asset down to an oriented box. The low rungs are the
ones that can be materialized (fewer faces = fewer papercraft pieces).
"""
import numpy as np
import trimesh

FACE_LEVELS = [20000, 5000, 1000, 250, 60]


def base_color(mesh):
    """Mean RGBA of a mesh's colours, or neutral grey."""
    vis = getattr(mesh, "visual", None)
    # an uncoloured mesh still reports trimesh's default dark grey, so ask whether colours are defined
    if isinstance(vis, trimesh.visual.ColorVisuals) and vis.defined:
        try:
            cols = vis.vertex_colors if vis.kind == "vertex" else vis.face_colors
            if cols is not None and len(cols):
                return np.asarray(cols, dtype=np.float64).mean(axis=0).astype(np.uint8)
        except Exception:
            pass
    return np.array([180, 180, 180, 255], dtype=np.uint8)


def _welded(mesh):
    """Geometry only, with coincident vertices merged, in fresh writable arrays.

    Asset files split vertices at UV and normal seams; the decimator can't collapse
    edges that aren't shared, so an unwelded mesh stalls far above its target.
    fast-simplification also rejects read-only arrays (cached views), hence the copies.
    """
    m = trimesh.Trimesh(np.array(mesh.vertices, dtype=np.float64), np.array(mesh.faces, dtype=np.int64))
    return trimesh.Trimesh(np.array(m.vertices, dtype=np.float64), np.array(m.faces, dtype=np.int64),
                           process=False)


def _clustered(mesh, cell):
    """Vertex clustering: snap vertices to a grid of `cell` metres and merge each cell."""
    v, f = np.asarray(mesh.vertices, dtype=np.float64), np.asarray(mesh.faces, dtype=np.int64)
    inv = np.unique(np.floor((v - v.min(axis=0)) / cell).astype(np.int64), axis=0, return_inverse=True)[1].ravel()
    nv = np.zeros((inv.max() + 1, 3))
    np.add.at(nv, inv, v)
    nv /= np.bincount(inv)[:, None]
    nf = inv[f]
    nf = nf[(nf[:, 0] != nf[:, 1]) & (nf[:, 1] != nf[:, 2]) & (nf[:, 0] != nf[:, 2])]
    return trimesh.Trimesh(nv, nf, process=False)


def _cluster_to(mesh, faces):
    """The finest grid that gets to `faces` or fewer (log-scale bisection on the cell size)."""
    ext = float(mesh.extents.max()) or 1.0
    lo, hi, best = ext / 2000, ext / 2, None
    for _ in range(20):
        cell = (lo * hi) ** 0.5
        m = _clustered(mesh, cell)
        if len(m.faces) > faces:
            lo = cell
        else:
            hi, best = cell, m
    return best


def decimate(mesh, faces):
    """Reduce toward `faces`, re-welding between passes while that still helps.

    Quadric decimation stalls on some CAD meshes (long sliver triangles: a SimReady
    lamp stopped at ~21k faces whatever the settings). Then vertex clustering takes
    over: cruder, but it always reduces and keeps thin parts that quadric drops.
    """
    color = base_color(mesh)
    out = _welded(mesh)
    while len(out.faces) > faces:
        nxt = _welded(out.simplify_quadric_decimation(face_count=faces))
        if len(nxt.faces) > 0.95 * len(out.faces):
            break
        out = nxt
    if len(out.faces) > faces:
        out = _cluster_to(out, faces) or out
    out.visual = trimesh.visual.ColorVisuals(out, face_colors=np.tile(color, (len(out.faces), 1)))
    return out


def _recolor(mesh, color):
    mesh.visual = trimesh.visual.ColorVisuals(mesh, face_colors=np.tile(color, (len(mesh.faces), 1)))
    return mesh


def _ray_hits(center, dirs, a, b):
    """Farthest hit of each 2D ray (center + t*dir) on segments a->b; nan where a ray misses."""
    e, w = b - a, a - center
    denom = dirs[:, None, 0] * e[None, :, 1] - dirs[:, None, 1] * e[None, :, 0]
    with np.errstate(divide="ignore", invalid="ignore"):
        t = (w[:, 0] * e[:, 1] - w[:, 1] * e[:, 0])[None] / denom
        s = (w[:, 0] * dirs[:, None, 1] - w[:, 1] * dirs[:, None, 0]) / denom
    far = np.where((np.abs(denom) > 1e-12) & (t >= 0) & (s >= 0) & (s <= 1), t, -np.inf).max(axis=1)
    return np.where(np.isfinite(far), far, np.nan)


def rings(mesh, count=8, samples=16):
    """Planar circumferences: `count` horizontal sections, each sampled at `samples` angles.

    Each ring is the section's outer outline seen from its own centre: the farthest
    crossing along `samples` evenly spaced rays. Neighbouring rings join into bands,
    so every band is one paper strip and every vertex is (ring, angle). Cruder than
    decimation, but each number can be checked by hand. Returns (mesh, profile).
    """
    lo, hi = mesh.bounds[:, 1]
    inset = 0.01 * (hi - lo)
    ang = np.linspace(0, 2 * np.pi, samples, endpoint=False)
    dirs = np.c_[np.cos(ang), np.sin(ang)]  # in the (x, z) plane
    profile = []
    for y in np.linspace(lo + inset, hi - inset, count):
        seg = trimesh.intersections.mesh_plane(mesh, [0, 1, 0], [0, y, 0])
        if len(seg) == 0:
            continue  # a gap between parts at this height
        a, b = seg[:, 0][:, [0, 2]], seg[:, 1][:, [0, 2]]
        pts = np.r_[a, b]
        center = (pts.min(axis=0) + pts.max(axis=0)) / 2
        r = _ray_hits(center, dirs, a, b)
        hit = ~np.isnan(r)
        if not hit.any():
            continue
        r = np.interp(ang, ang[hit], r[hit], period=2 * np.pi)  # fill rays that found nothing
        profile.append({"y": round(float(y), 3), "center": [round(float(c), 3) for c in center],
                        "radii": [round(float(x), 3) for x in r]})
    if len(profile) < 2:
        return None, profile
    verts = [[p["center"][0] + rr * dx, p["y"], p["center"][1] + rr * dz]
             for p in profile for rr, (dx, dz) in zip(p["radii"], dirs)]
    n, faces = len(profile), []
    for i in range(n - 1):
        for j in range(samples):
            v00, v01 = i * samples + j, i * samples + (j + 1) % samples
            faces += [[v00, v00 + samples, v01 + samples], [v00, v01 + samples, v01]]
    for k, p in ((0, profile[0]), (n - 1, profile[-1])):  # flat caps: fan to the ring centre
        c = len(verts)
        verts.append([p["center"][0], p["y"], p["center"][1]])
        faces += [[c, k * samples + (j + 1) % samples, k * samples + j] for j in range(samples)]
    out = trimesh.Trimesh(np.array(verts), np.array(faces), process=False)
    if out.volume < 0:
        out.invert()
    return out, profile


def _collapsed(rung, original, prev):
    """A rung that shrank in space, not just in faces: under half the original's size,
    or under half the previous decimated rung's area. The original is left out of the
    area test, since welding away double-sided faces can halve area legitimately.
    """
    if rung.extents.max() < 0.5 * original.extents.max():
        return True
    return prev is not None and rung.area < 0.5 * prev.area


def ladder(mesh, ring_count=8, ring_samples=16):
    color = base_color(mesh)
    rungs = [("original", mesh)]
    for n in FACE_LEVELS:
        prev = len(rungs[-1][1].faces)
        if n >= 0.8 * prev:
            continue
        rung = decimate(mesh, n)
        if len(rung.faces) > 0.8 * prev:
            break  # decimation stalled: lower levels would repeat this mesh
        if _collapsed(rung, mesh, rungs[-1][1] if len(rungs) > 1 else None):
            break  # a degenerate blob would pass for a rung; lower levels only get worse
        rungs.append((f"{len(rung.faces)} faces", rung))
    ringed, profile = rings(mesh, ring_count, ring_samples)
    if ringed is not None:
        ringed.metadata["rings"] = profile
        rungs.append((f"rings {len(profile)}x{ring_samples}", _recolor(ringed, color)))
    rungs.append(("convex hull", _recolor(mesh.convex_hull, color)))
    box = mesh.bounding_box_oriented
    rungs.append(("bounding box", _recolor(trimesh.Trimesh(box.vertices, box.faces), color)))
    return rungs


def stats(mesh):
    return {"faces": int(len(mesh.faces)), "vertices": int(len(mesh.vertices)),
            "watertight": bool(mesh.is_watertight), "convex": bool(mesh.is_convex),
            "surface_area_m2": round(float(mesh.area), 4)}