solidify_server.py

"""
Make solid: any model rebuilt as one closed, printable surface, with as much
of its detail as a CPU can keep, and its colour baked onto the new mesh.

  seal        solid.py at 512 or 768 voxels across: the model filled into a
              grid and wrapped in one sealed surface (marching cubes on a
              lightly smoothed field), so holes, overlapping parts and edges
              shared by three faces are gone whatever the input was
  snap        every vertex near the original surface moved onto it (90% of
              the way), which brings back the sharp edges the voxels rounded
              off; only positions move, so the seal cannot break
  reduce      a quadric decimator down to a triangle count a slicer takes
              happily; anything it tangles is cut out and the hole filled,
              so the result stays closed; then snapped once more
  colour      when the model has a texture or vertex colours: new UVs
              (axis-facing charts, packed in rows) and base colour plus metallic/roughness baked from the
              original. No normal map: the mesh is dense enough to carry its
              own shape.
"""

import json
import time

import fast_simplification
import numpy as np
import point_cloud_utils as pcu
import trimesh
from scipy import ndimage
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from skimage import measure

from glb import write_glb
from remesh import load, rasterize, seam_fill, source_look, to_image
from solid import PAD

# Voxels along the longest side, and the triangles the result is reduced to.
# On two cores of a desktop CPU, a textured 200k-triangle generated model (the
# chest, the teapot) takes 13 to 15 s at standard and 24 to 30 s at fine,
# colour bake included, with at most 3.5 GB of memory. A cpu-basic Space
# (2 vCPU, 16 GB) is about two to three times slower: about 45 s and 90 s.
DETAIL = {
    "standard": {"resolution": 512, "triangles": 300_000},
    "fine": {"resolution": 768, "triangles": 600_000},
}
TEXTURE_SIZE = 2048
# How far a vertex may be from the original to be snapped onto it, in voxels,
# and how far of the way it goes (as the site's snapTo: neighbours that would
# land on the same original corner stay a hair apart, so an STL, which has no
# connectivity, does not pinch there when a slicer merges equal points).
# The sealed surface lies up to about a voxel outside the original, and parts
# widened to three voxels lie a voxel and a half out; 2.5 voxels reaches them
# all (on a race car, the mean distance to the original drops from 0.7 to
# 0.26 per mille of its size) and is still less, in model units, than the
# browser's 1.5 voxels at 256 across.
SNAP_REACH = 2.5
SNAP_SHARE = 0.9
# The least cosine between a vertex's normal and the original face it snaps to.
SNAP_FACING = 0.3
# A separate piece with fewer triangles than this share of the largest is a
# crumb of the input (a floating shard), not a part: it is dropped.
CRUMB = 0.01


def edge_table(faces):
    """Every half-edge as (undirected key, direction, face)."""
    f = np.asarray(faces, dtype=np.int64)
    a = f[:, [0, 1, 2]].ravel()
    b = f[:, [1, 2, 0]].ravel()
    lo, hi = np.minimum(a, b), np.maximum(a, b)
    n = int(f.max()) + 1 if len(f) else 1
    return lo * n + hi, a < b, np.repeat(np.arange(len(f)), 3), a, b


def edge_report(faces):
    """(open edges, tangled edges): edges with one face, and edges with more than two or two the same way round."""
    key, forward, _, _, _ = edge_table(faces)
    order = np.argsort(key, kind="stable")
    key, forward = key[order], forward[order]
    uniq, start, count = np.unique(key, return_index=True, return_counts=True)
    ways = np.add.reduceat(forward.astype(np.int64), start)
    open_edges = int((count == 1).sum())
    tangled = int(((count > 2) | ((count == 2) & (ways != 1))).sum())
    return open_edges, tangled


def piece_labels(n_vertices, faces):
    """Which connected piece every face belongs to, and how many pieces there are."""
    f = np.asarray(faces, dtype=np.int64)
    rows = np.concatenate([f[:, 0], f[:, 1]])
    cols = np.concatenate([f[:, 1], f[:, 2]])
    graph = coo_matrix((np.ones(len(rows), dtype=np.int8), (rows, cols)), shape=(n_vertices, n_vertices))
    count, label = connected_components(graph, directed=False)
    return label[f[:, 0]], count


def drop_crumbs(v, faces, share=CRUMB):
    """Pieces much smaller than the largest removed; returns (faces, pieces kept)."""
    label, count = piece_labels(len(v), faces)
    if count <= 1:
        return faces, 1
    sizes = np.bincount(label)
    keep_piece = sizes >= max(12, share * sizes.max())
    return faces[keep_piece[label]], int(keep_piece.sum())


def compact(v, faces):
    """Vertices nothing uses dropped, faces renumbered."""
    used = np.unique(faces)
    remap = np.full(len(v), -1, dtype=np.int64)
    remap[used] = np.arange(len(used))
    return v[used], remap[faces]


def boundary_loops(faces):
    """The open edges of a surface as closed loops of vertices, in the direction the faces run."""
    key, _, _, a, b = edge_table(faces)
    uniq, inverse, count = np.unique(key, return_inverse=True, return_counts=True)
    lone = count[inverse] == 1
    starts = {}
    for s, e in zip(a[lone].tolist(), b[lone].tolist()):
        starts.setdefault(s, []).append(e)
    loops = []
    while starts:
        first = next(iter(starts))
        loop, at = [first], first
        while True:
            ends = starts.get(at)
            if not ends:
                break
            nxt = ends.pop()
            if not ends:
                del starts[at]
            if nxt == first:
                loops.append(loop)
                break
            if nxt in loop:
                # Two holes touching at a vertex: close the part walked so far.
                cut = loop.index(nxt)
                loops.append(loop[cut:])
                loop = loop[: cut + 1]
                at = nxt
                continue
            loop.append(nxt)
            at = nxt
    return [loop for loop in loops if len(loop) >= 3]


def mend(v, faces, rounds=6):
    """Faces on tangled edges cut out, crumbs dropped, holes filled with a fan.

    A decimator now and then folds two sheets onto one edge where a part is
    about one edge thin; the few faces around it go, and the small hole they
    leave is closed again. Returns (vertices, faces, closed).
    """
    v = np.asarray(v, dtype=np.float64)
    faces = np.asarray(faces, dtype=np.int64)
    for _ in range(rounds):
        key, forward, face_of, _, _ = edge_table(faces)
        order = np.argsort(key, kind="stable")
        skey, sfwd, sface = key[order], forward[order], face_of[order]
        uniq, start, count = np.unique(skey, return_index=True, return_counts=True)
        ways = np.add.reduceat(sfwd.astype(np.int64), start)
        bad_edge = (count > 2) | ((count == 2) & (ways != 1))
        bad = np.zeros(len(faces), dtype=bool)
        if bad_edge.any():
            bad[sface[np.repeat(bad_edge, count)]] = True
            faces = faces[~bad]
        faces, _ = drop_crumbs(v, faces)
        loops = boundary_loops(faces)
        if not loops and not bad.any():
            break
        extra_v, extra_f = [], []
        for loop in loops:
            ring = loop[::-1]
            if len(ring) == 3:
                extra_f.append(ring)
                continue
            centre = len(v) + len(extra_v)
            extra_v.append(v[ring].mean(axis=0))
            extra_f.extend([ring[k], ring[(k + 1) % len(ring)], centre] for k in range(len(ring)))
        if extra_v:
            v = np.concatenate([v, np.asarray(extra_v)])
        if extra_f:
            faces = np.concatenate([faces, np.asarray(extra_f, dtype=np.int64)])
    v, faces = compact(v, faces)
    open_edges, tangled = edge_report(faces)
    return v, faces, open_edges == 0 and tangled == 0


def reduce(v, faces, target):
    """About `target` triangles by quadric decimation, closed and edge-manifold.

    If mending cannot close what the decimator did, a gentler reduction is
    tried, and at last the unreduced surface is kept: larger, but sealed.
    """
    for want in (target, int(target * 1.5)):
        if len(faces) <= want:
            break
        rv, rf = fast_simplification.simplify(v, faces, target_count=want, agg=7)
        rv, rf, closed = mend(rv, rf)
        if closed:
            return rv, rf, True
        print(f"Reduction to {want} could not be closed, trying gentler", flush=True)
    rv, rf, closed = mend(v, faces)
    return rv, rf, closed


def snap(v, f, source_v, source_f, source_normals, reach, share=SNAP_SHARE):
    """Vertices within `reach` of the original surface moved `share` of the way onto it.

    Only onto original faces turned roughly the same way (either way round,
    as generated models carry patches wound backwards): at a sharp edge, the
    nearest original face can be the side around the corner, and a vertex
    pulled there shows as a hair across the edge.
    """
    normals = trimesh.Trimesh(v, f, process=False).vertex_normals
    distance, face, bc = pcu.closest_points_on_mesh(np.asarray(v, dtype=np.float64), source_v, source_f)
    closest = np.einsum("nk,nkd->nd", bc, source_v[source_f[face]])
    facing = np.abs(np.einsum("ij,ij->i", source_normals[face], normals)) > SNAP_FACING
    move = ((distance < reach) & facing)[:, None]
    return np.where(move, v + share * (closest - v), v)


AXES = np.array([[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]], dtype=np.float64)
# A face may join a chart facing up to this far from its own normal (cosine),
# so charts do not break up where the normal hovers between two axes.
CHART_COS = 0.35
CHART_PADDING = 3


def charts(v, f):
    """Faces grouped into charts that each project flat onto one axis plane: (chart per face, axis per chart)."""
    a, b, c = v[f[:, 0]], v[f[:, 1]], v[f[:, 2]]
    n = np.cross(b - a, c - a)
    n /= np.linalg.norm(n, axis=1, keepdims=True) + 1e-20
    facing = n @ AXES.T
    label = facing.argmax(1)
    allowed = facing > CHART_COS
    adjacency = trimesh.graph.face_adjacency(f)
    # A few rounds of majority vote among the three neighbours, within what each face allows.
    for _ in range(4):
        votes = np.zeros((len(f), 6), dtype=np.int32)
        np.add.at(votes, (adjacency[:, 0], label[adjacency[:, 1]]), 1)
        np.add.at(votes, (adjacency[:, 1], label[adjacency[:, 0]]), 1)
        votes[np.arange(len(f)), label] += 1
        votes[~allowed] = -1
        label = votes.argmax(1)
    same = label[adjacency[:, 0]] == label[adjacency[:, 1]]
    pairs = adjacency[same]
    graph = coo_matrix((np.ones(len(pairs), dtype=np.int8), (pairs[:, 0], pairs[:, 1])), shape=(len(f), len(f)))
    _, chart = connected_components(graph, directed=False)
    axis = np.zeros(chart.max() + 1, dtype=np.int64)
    axis[chart] = label
    return chart, axis


def shelf_pack(widths, heights, size, padding):
    """Rectangles (in model units) packed in rows into a square of `size` texels.

    Returns (scale in texels per unit, x, y of each rectangle in texels). The
    largest scale that fits is found by bisection.
    """
    order = np.argsort(-heights, kind="stable")

    def place(scale):
        w = np.ceil(widths * scale) + padding
        h = np.ceil(heights * scale) + padding
        x = np.zeros(len(w))
        y = np.zeros(len(w))
        cx = cy = row = 0.0
        for i in order.tolist():
            if w[i] > size:
                return None
            if cx + w[i] > size:
                cy += row
                cx = row = 0.0
            x[i], y[i] = cx, cy
            cx += w[i]
            row = max(row, h[i])
        return (x, y) if cy + row <= size else None

    lo, hi = 0.0, size / max(float(np.sqrt((widths * heights).sum())), 1e-12) * 1.5
    best = None
    for _ in range(18):
        mid = 0.5 * (lo + hi)
        placed = place(mid)
        if placed is None:
            hi = mid
        else:
            lo, best = mid, placed
    if best is None:
        raise ValueError("The model has too many separate surfaces to lay out in one texture.")
    return lo, best[0], best[1]


def unwrap(v, f, size, padding=CHART_PADDING):
    """UVs for a dense mesh, fast: (vertex map, faces, uvs with v up).

    xatlas takes minutes on a few hundred thousand triangles. A mesh this
    dense carries its shape in the geometry, so the texture only needs even
    colour, not tidy seams: faces are grouped into charts that each face one
    of the six axis directions, every chart is projected straight onto its
    plane (no stretching beyond the slope of its faces), and the charts are
    packed in rows at one texel density.
    """
    chart, axis = charts(v, f)
    # Every corner gets its chart's copy of its vertex.
    corner_key = chart[np.repeat(np.arange(len(f)), 3)] * len(v) + f.ravel()
    keys, new_index = np.unique(corner_key, return_inverse=True)
    vmap = keys % len(v)
    vchart = keys // len(v)
    faces = new_index.reshape(-1, 3).astype(np.int64)
    # The two coordinates across each chart's axis.
    dims = np.array([[1, 2], [2, 1], [2, 0], [0, 2], [0, 1], [1, 0]])[axis[vchart]]
    p = v[vmap]
    uv = np.stack([p[np.arange(len(p)), dims[:, 0]], p[np.arange(len(p)), dims[:, 1]]], axis=1)
    count = axis.size
    lo = np.full((count, 2), np.inf)
    hi = np.full((count, 2), -np.inf)
    np.minimum.at(lo, vchart, uv)
    np.maximum.at(hi, vchart, uv)
    extent = hi - lo
    scale, x, y = shelf_pack(extent[:, 0], extent[:, 1], size, padding)
    px = x[vchart] + padding / 2 + (uv[:, 0] - lo[vchart, 0]) * scale
    py = y[vchart] + padding / 2 + (uv[:, 1] - lo[vchart, 1]) * scale
    uvs = np.stack([px / size, 1.0 - py / size], axis=1)
    return vmap, faces, uvs


def bake_colour(positions, uvs, faces, high, size):
    """Base colour (and metallic/roughness when the original has it) baked per texel from the nearest point of the original.

    remesh.bake without the normal map, which a mesh this dense does not need.
    """
    rows, cols, tri, bary = rasterize(uvs, faces, size)
    point = np.einsum("nk,nkd->nd", bary, positions[faces[tri]])
    _, face, bc = pcu.closest_points_on_mesh(point, np.asarray(high.vertices, dtype=np.float64), np.asarray(high.faces, dtype=np.int32))
    lookup, has_mr = source_look(high)
    color, mr = lookup(face, bc)
    mask = np.zeros((size, size), dtype=bool)
    mask[rows, cols] = True
    near, from_y, from_x = seam_fill(mask)

    def image(values, background):
        img = np.tile(np.asarray(background, dtype=np.float32), (size, size, 1))
        img[rows, cols] = values
        img[near] = img[from_y, from_x]
        return img

    maps = {"color": image(color, [0, 0, 0, 1])}
    if has_mr:
        maps["metallicRoughness"] = image(mr, [0, 1, 0])
    return maps, float(mask.mean())


# Points closer than this (in the model's units) are one point to a reader:
# the site welds on this grid, and slicers weld nearly equal points too.
WELD = 1e-5


def as_slicer_sees_it(v, faces):
    """The mesh with nearly equal points merged, as the site and slicers read it, and anything that tangles mended.

    Snapping can bring two vertices that sat a hair apart onto the same
    point. Connected, they are fine; in an STL, which has no connectivity,
    the reader merges them and two sheets share an edge. Merging here first,
    on the same grid the site uses, means the file is checked as it will be read.
    """
    keys = np.round(np.asarray(v, dtype=np.float64) / WELD).astype(np.int64)
    _, first, inverse = np.unique(keys, axis=0, return_index=True, return_inverse=True)
    points = np.asarray(v, dtype=np.float64)[first]
    f = inverse.ravel()[faces]
    f = f[(f[:, 0] != f[:, 1]) & (f[:, 1] != f[:, 2]) & (f[:, 0] != f[:, 2])]
    v = points
    open_edges, tangled = edge_report(f)
    if open_edges or tangled:
        v, f, _ = mend(v, f)
    else:
        v, f = compact(v, f)
    return v, f


def has_look(mesh):
    """True when the model carries colour worth baking: a texture or vertex colours."""
    visual = mesh.visual
    if getattr(visual, "kind", None) == "vertex":
        return True
    material = getattr(visual, "material", None)
    if material is None or getattr(visual, "uv", None) is None:
        return False
    image = getattr(material, "baseColorTexture", None)
    if image is None:
        image = getattr(material, "image", None)
    return image is not None


def _shifted(grid, combine):
    """`grid` combined with its six face neighbours (the cross ndimage uses by default), outside the grid counting as empty."""
    out = grid.copy()
    for axis in range(3):
        for side in (1, -1):
            src = [slice(None)] * 3
            dst = [slice(None)] * 3
            src[axis] = slice(0, -1) if side == 1 else slice(1, None)
            dst[axis] = slice(1, None) if side == 1 else slice(0, -1)
            combine(out[tuple(dst)], grid[tuple(src)], out=out[tuple(dst)])
    if combine is np.logical_and:
        # Erosion: a voxel on the grid's own border has an empty neighbour.
        for axis in range(3):
            edge = [slice(None)] * 3
            for end in (0, -1):
                edge[axis] = end
                out[tuple(edge)] = False
    return out


def dilate(grid):
    return _shifted(grid, np.logical_or)


def erode(grid):
    return _shifted(grid, np.logical_and)


def seal_grid(vertices, faces, resolution):
    """solid.solid_grid followed by solid.thicken(grid, 1.0), with the same result, faster.

    At 768 voxels across the grid has a few hundred million cells; ndimage's
    general morphology takes seconds per pass there, while a pass with the
    six-neighbour cross (all these steps use) is six shifted array operations.
    The wall is also sampled in single precision, in batches.
    """
    lo, hi = vertices.min(0), vertices.max(0)
    h = (hi - lo).max() / resolution
    shape = np.ceil((hi - lo) / h).astype(int) + 2 * PAD + 1

    a, b, c = vertices[faces[:, 0]], vertices[faces[:, 1]], vertices[faces[:, 2]]
    longest = np.max(np.stack([np.linalg.norm(b - a, axis=1), np.linalg.norm(c - a, axis=1), np.linalg.norm(c - b, axis=1)]), axis=0)
    n = np.maximum(1, np.ceil(longest / (0.5 * h))).astype(int)
    # Positions in voxel units, so a sample is floored straight into its cell.
    a = ((a - lo) / h + PAD).astype(np.float32)
    ab = ((b - lo) / h + PAD).astype(np.float32) - a
    ac = ((c - lo) / h + PAD).astype(np.float32) - a
    wall = np.zeros(shape, dtype=bool)
    flat = wall.reshape(-1)
    strides = np.array([shape[1] * shape[2], shape[2], 1], dtype=np.int64)
    for steps in np.unique(n):
        sel = np.flatnonzero(n == steps)
        i, j = np.meshgrid(np.arange(steps + 1), np.arange(steps + 1), indexing="ij")
        keep = i + j <= steps
        wi = (i[keep] / steps).astype(np.float32)[None, :, None]
        wj = (j[keep] / steps).astype(np.float32)[None, :, None]
        batch = max(1, 4_000_000 // wi.shape[1])
        for start in range(0, len(sel), batch):
            part = sel[start : start + batch]
            pts = a[part][:, None] + ab[part][:, None] * wi + ac[part][:, None] * wj
            idx = pts.reshape(-1, 3).astype(np.int64)
            flat[idx @ strides] = True

    thick = dilate(wall)
    labels, _ = ndimage.label(~thick)
    outside = labels == labels[0, 0, 0]
    del labels, thick
    solid = ~outside
    solid &= ~(dilate(outside) & ~wall)
    del outside, wall
    # thicken(solid, 1.0): ball(1) is the same cross.
    thin = solid & ~dilate(erode(solid))
    if thin.any():
        solid |= dilate(thin)
    return solid, h, lo


def sealed_surface(solid, h, lo, block=64, halo=4):
    """solid.grid_surface with the same result, faster on a large grid.

    Marching cubes only looks at cells near the solid's boundary (a band three
    voxels wide each side, which catches every crossing of the smoothed
    field), and the field is only smoothed in the blocks of the grid that
    band touches; elsewhere it is 0 or 1 as the voxels are.
    """
    band = solid & ~erode(solid)
    band = dilate(dilate(dilate(band)))
    reach = dilate(band)
    field = solid.astype(np.float32)
    shape = solid.shape
    blocks = [(shape[k] + block - 1) // block for k in range(3)]
    padded = np.pad(reach, [(0, blocks[k] * block - shape[k]) for k in range(3)])
    active = padded.reshape(blocks[0], block, blocks[1], block, blocks[2], block).any(axis=(1, 3, 5))
    del padded, reach
    for corner in np.argwhere(active):
        start = corner * block
        stop = np.minimum(start + block, shape)
        gstart = np.maximum(start - halo, 0)
        gstop = np.minimum(stop + halo, shape)
        sub = solid[gstart[0] : gstop[0], gstart[1] : gstop[1], gstart[2] : gstop[2]].astype(np.float32)
        smooth = ndimage.gaussian_filter(sub, sigma=1.0, mode="constant")
        a, b = start - gstart, stop - gstart
        field[start[0] : stop[0], start[1] : stop[1], start[2] : stop[2]] = smooth[a[0] : b[0], a[1] : b[1], a[2] : b[2]]
    verts, tris, _, _ = measure.marching_cubes(field, level=0.5, allow_degenerate=False, mask=band)
    del field, band
    # Cell i covers [lo + (i - PAD) h, lo + (i - PAD + 1) h): its value sits at
    # the cell's centre. (grid_surface places it at the cell's corner, half a
    # voxel off along every axis; the remesher bakes its look afterwards, so
    # it does not show there, but here every vertex counts.)
    verts = (verts - PAD + 0.5) * h + lo
    # Faces turned so normals point out, checked by volume (as grid_surface).
    tris = tris[:, [0, 2, 1]]
    a, b, c = verts[tris[:, 0]], verts[tris[:, 1]], verts[tris[:, 2]]
    if np.einsum("ij,ij->i", a, np.cross(b, c)).sum() < 0:
        tris = tris[:, [0, 2, 1]]
    return verts.astype(np.float64), tris.astype(np.int64)


def solid_mesh(high, detail="standard", step=None, timings=None):
    """The sealed, reduced and snapped surface of a loaded model: (vertices, faces, info)."""
    step = step or (lambda fraction, text: None)
    timings = {} if timings is None else timings
    clock = [time.perf_counter()]

    def lap(name):
        now = time.perf_counter()
        timings[name] = round(now - clock[0], 1)
        clock[0] = now

    size = DETAIL[detail]
    source_v = np.asarray(high.vertices, dtype=np.float64)
    source_f = np.asarray(high.faces, dtype=np.int32)

    step(0.1, "Sealing the model in a voxel grid")
    grid, h, lo = seal_grid(source_v, source_f, size["resolution"])
    lap("seal")
    step(0.3, "Wrapping it in one surface")
    v, f = sealed_surface(grid, h, lo)
    del grid
    f, _ = drop_crumbs(v, f)
    v, f = compact(v, f)
    sealed = len(f)
    lap("surface")

    # Snapped before the reduction too, so the decimator works on the sharp
    # shape and keeps its edges, rather than on the rounded voxel skin
    # (measured on the chest: mean distance to the original 0.03 per mille of
    # its size this way, against 0.11 snapping only after, and no streaks
    # from long thin triangles pulled onto it).
    step(0.4, "Snapping onto the original surface")
    source_normals = np.asarray(high.face_normals, dtype=np.float64)
    v = snap(v, f, source_v, source_f, source_normals, SNAP_REACH * h)
    lap("snap")

    step(0.5, "Reducing to a printable size")
    v, f, closed = reduce(v, f, size["triangles"])
    lap("reduce")

    step(0.6, "Snapping again for the sharp edges")
    v = snap(v, f, source_v, source_f, source_normals, SNAP_REACH * h)
    v, f = as_slicer_sees_it(v, f)
    lap("snap_again")

    _, pieces = piece_labels(len(v), f)
    open_edges, tangled = edge_report(f)
    info = {
        "detail": detail,
        "resolution": size["resolution"],
        "voxel": float(h),
        "sealed_triangles": int(sealed),
        "triangles": int(len(f)),
        "vertices": int(len(v)),
        "manifold": bool(open_edges == 0 and tangled == 0),
        "pieces": int(pieces),
    }
    return v, f, info


def solidify(path, glb_out, stl_out, detail="standard", texture_size=TEXTURE_SIZE, progress=None):
    """Make solid, from a model file to a GLB (textured when the input was) and an STL. Returns the details."""
    if detail not in DETAIL:
        raise ValueError(f"Detail must be one of {', '.join(DETAIL)}.")
    step = progress or (lambda fraction, text: None)
    timings = {}
    clock = time.perf_counter()
    started = clock

    def lap(name):
        nonlocal clock
        now = time.perf_counter()
        timings[name] = round(now - clock, 1)
        clock = now

    step(0.02, "Reading the model")
    high = load(path)
    if len(high.faces) == 0:
        raise ValueError("The file has no triangles in it.")
    lap("load")

    v, f, info = solid_mesh(high, detail, step, timings)
    clock = time.perf_counter()

    mesh = trimesh.Trimesh(v, f, process=False)
    normals = np.asarray(mesh.vertex_normals, dtype=np.float64)
    mesh.export(stl_out, file_type="stl")

    textured = has_look(high)
    if textured:
        step(0.7, "Unwrapping UVs")
        vmap, faces, uvs = unwrap(v, f, texture_size)
        positions, vnormals = v[vmap], normals[vmap]
        lap("uv")

        step(0.8, "Baking the colour from the original")
        maps, coverage = bake_colour(positions, uvs, faces, high, texture_size)
        lap("bake")

        step(0.95, "Writing the files")
        write_glb(
            glb_out, positions, vnormals, None, uvs, faces,
            {
                "color": to_image(maps["color"], "RGBA"),
                "metallicRoughness": to_image(maps["metallicRoughness"], "RGB") if "metallicRoughness" in maps else None,
                "normal": None,
            },
            generator="3dvalley.com make solid",
        )
        info["uv_coverage"] = round(float(coverage), 3)
    else:
        step(0.95, "Writing the files")
        mesh.export(glb_out, file_type="glb")
    lap("write")

    info["textured"] = bool(textured)
    info["timings"] = timings
    info["seconds"] = round(time.perf_counter() - started, 1)
    return info


if __name__ == "__main__":
    import sys

    src, dst = sys.argv[1], sys.argv[2]
    level = sys.argv[3] if len(sys.argv) > 3 else "standard"
    print(json.dumps(solidify(src, dst, dst.replace(".glb", ".stl"), level), indent=1))