from __future__ import annotations

from collections import defaultdict
from dataclasses import dataclass
from dataclasses import replace
import heapq
import math
import time

import numpy as np

from .mesh import Meshlet, ObjMesh, compute_triangle_centers, compute_triangle_normals, unique_valid
from .progress import Heartbeat, log, step


@dataclass
class MeshletBuildOptions:
    builder: str = "greedy"
    target_vertices: int = 30
    max_vertices: int = 127
    max_texcoords: int = 255
    max_normals: int = 255
    max_triangles: int = 127
    max_normal_angle_deg: float = 90.0
    radius_weight: float = 7.0
    normal_weight: float = 18.0
    vertex_weight: float = 2.0
    stop_score: float = 4.8
    shared_vertex_weight: float = 4.0
    test_pos_weight: float = 3.0
    compatible_edge_weight: float = 0.0
    incompatible_edge_penalty: float = 0.0
    compatible_merge_weight: float = 0.0
    lkh_exe: str | None = None
    visibility_merge_samples: int = 1024
    visibility_merge_size: int = 1024
    visibility_merge_margin_deg: float = -1.0
    visibility_merge_cone_weight: float = 45.0
    visibility_pack_mode: str = "none"
    visibility_pack_compression: float = 50.0
    fill_holes: bool = True
    merge_passes: int = 1
    smooth_passes: int = 0
    merge_max_normal_angle_deg: float = 48.0


def _edge_key(a: int, b: int) -> tuple[int, int]:
    return (a, b) if a < b else (b, a)


def build_triangle_adjacency(mesh: ObjMesh) -> list[set[int]]:
    edge_to_triangles: dict[tuple[int, int], list[int]] = defaultdict(list)
    for ti, tri in enumerate(mesh.triangles):
        v = [c.v for c in tri.corners]
        edge_to_triangles[_edge_key(v[0], v[1])].append(ti)
        edge_to_triangles[_edge_key(v[1], v[2])].append(ti)
        edge_to_triangles[_edge_key(v[2], v[0])].append(ti)

    adjacency = [set() for _ in mesh.triangles]
    for tris in edge_to_triangles.values():
        if len(tris) < 2:
            continue
        for a in tris:
            adjacency[a].update(b for b in tris if b != a)
    return adjacency


def build_edge_to_triangles(mesh: ObjMesh) -> dict[tuple[int, int], list[int]]:
    edge_to_triangles: dict[tuple[int, int], list[int]] = defaultdict(list)
    for ti, tri in enumerate(mesh.triangles):
        v = [c.v for c in tri.corners]
        edge_to_triangles[_edge_key(v[0], v[1])].append(ti)
        edge_to_triangles[_edge_key(v[1], v[2])].append(ti)
        edge_to_triangles[_edge_key(v[2], v[0])].append(ti)
    return edge_to_triangles


def _corner_edge_key(a, b) -> tuple:
    ka = (a.v, a.vt, a.vn)
    kb = (b.v, b.vt, b.vn)
    return (ka, kb) if ka < kb else (kb, ka)


def build_compatible_edge_to_triangles(mesh: ObjMesh) -> dict[tuple, list[int]]:
    edge_to_triangles: dict[tuple, list[int]] = defaultdict(list)
    for ti, tri in enumerate(mesh.triangles):
        c = tri.corners
        edge_to_triangles[_corner_edge_key(c[0], c[1])].append(ti)
        edge_to_triangles[_corner_edge_key(c[1], c[2])].append(ti)
        edge_to_triangles[_corner_edge_key(c[2], c[0])].append(ti)
    return edge_to_triangles


def _triangle_sets(mesh: ObjMesh, tri_index: int) -> tuple[set[int], set[int], set[int]]:
    tri = mesh.triangles[tri_index]
    return (
        {c.v for c in tri.corners},
        unique_valid(c.vt for c in tri.corners),
        unique_valid(c.vn for c in tri.corners),
    )


def _sphere_from_vertices(points: np.ndarray) -> tuple[np.ndarray, float]:
    bb_min = np.min(points, axis=0)
    bb_max = np.max(points, axis=0)
    center = (bb_min + bb_max) * 0.5
    radius = float(np.max(np.linalg.norm(points - center, axis=1))) if len(points) else 0.0
    return center, radius


def _normal_cone(normals: np.ndarray) -> tuple[np.ndarray, float]:
    if len(normals) == 0:
        return np.array([0.0, 0.0, 1.0], dtype=np.float64), 1.0
    if float(np.max(np.linalg.norm(normals, axis=1))) <= 1e-12:
        return np.array([0.0, 0.0, 1.0], dtype=np.float64), 0.0
    if len(normals) == 1:
        length = float(np.linalg.norm(normals[0]))
        if length <= 1e-12:
            return np.array([0.0, 0.0, 1.0], dtype=np.float64), 0.0
        return normals[0].copy() / length, 1.0

    seeds = [
        np.array([1.0, 0.0, 0.0]),
        np.array([-1.0, 0.0, 0.0]),
        np.array([0.0, 1.0, 0.0]),
        np.array([0.0, -1.0, 0.0]),
        np.array([0.0, 0.0, 1.0]),
        np.array([0.0, 0.0, -1.0]),
    ]
    mean_axis = np.sum(normals, axis=0)
    mean_len = float(np.linalg.norm(mean_axis))
    if mean_len > 1e-12:
        seeds.append(mean_axis / mean_len)
    seeds.extend(normals)

    best_axis = seeds[0]
    best_dot = -2.0
    for axis in seeds:
        axis_len = float(np.linalg.norm(axis))
        if axis_len <= 1e-12:
            continue
        axis = axis / axis_len
        min_dot = float(np.min(normals @ axis))
        if min_dot > best_dot:
            best_dot = min_dot
            best_axis = axis

    # Local search on the unit sphere. This is the same spirit as the old prototype:
    # maximize the minimum dot product between the axis and all triangle normals.
    directions = np.asarray(
        [
            [1.0, 0.0, 0.0],
            [-1.0, 0.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, -1.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.0, 0.0, -1.0],
            [0.0, 0.52573111, 0.85065081],
            [0.0, -0.52573111, 0.85065081],
            [0.52573111, 0.85065081, 0.0],
            [-0.52573111, 0.85065081, 0.0],
            [0.85065081, 0.0, 0.52573111],
            [-0.85065081, 0.0, 0.52573111],
        ],
        dtype=np.float64,
    )
    eps = 0.1
    while eps > 0.0005:
        improved = False
        base_axis = best_axis.copy()
        for direction in directions:
            axis = base_axis + eps * direction
            axis_len = float(np.linalg.norm(axis))
            if axis_len <= 1e-12:
                continue
            axis /= axis_len
            min_dot = float(np.min(normals @ axis))
            if min_dot > best_dot:
                best_dot = min_dot
                best_axis = axis
                improved = True
        if not improved:
            eps *= 0.5
    return best_axis, best_dot


def _normal_cone_fast(normals: np.ndarray) -> tuple[np.ndarray, float]:
    if len(normals) == 0:
        return np.array([0.0, 0.0, 1.0], dtype=np.float64), 1.0
    axis = np.sum(normals, axis=0)
    length = float(np.linalg.norm(axis))
    if length <= 1e-12:
        axis = normals[0].copy()
        length = float(np.linalg.norm(axis))
    axis = axis / length if length > 0.0 else np.array([0.0, 0.0, 1.0], dtype=np.float64)
    return axis, float(np.min(normals @ axis))


def _backface_cull_cone(normal_axis: np.ndarray, normal_min_dot: float) -> tuple[np.ndarray, float]:
    """Return the object-to-camera visibility cone derived from a normal cone.

    The normal cone half-angle is alpha=acos(normal_min_dot). A face can be visible
    when the object-to-camera direction lies within 90deg+alpha from the normal-cone
    axis. The runtime stores visibility cones, so the cosine is -sin(alpha).
    cull_cos > 1 marks an empty/disabled culling cone.
    """
    if normal_min_dot <= 0.0:
        return normal_axis, 2.0
    cull_cos = -float(np.sqrt(max(0.0, 1.0 - normal_min_dot * normal_min_dot)))
    return normal_axis, cull_cos


def _make_meshlet(mesh: ObjMesh, tri_indices: list[int], triangle_normals: np.ndarray) -> Meshlet:
    vertices: set[int] = set()
    texcoords: set[int] = set()
    normals: set[int] = set()
    for ti in tri_indices:
        v, vt, vn = _triangle_sets(mesh, ti)
        vertices.update(v)
        texcoords.update(vt)
        normals.update(vn)

    center, radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
    normal_axis, normal_min_dot = _normal_cone(triangle_normals[tri_indices])
    cull_axis, cull_cos = _backface_cull_cone(normal_axis, normal_min_dot)
    material = mesh.triangles[tri_indices[0]].material
    return Meshlet(tri_indices, material, vertices, texcoords, normals, center, radius, normal_axis, normal_min_dot, cull_axis, cull_cos)


def build_meshlets(mesh: ObjMesh, options: MeshletBuildOptions | None = None) -> list[Meshlet]:
    """Build connected, local, normal-coherent meshlets under meshlet index limits."""
    opt = options or MeshletBuildOptions()
    if opt.builder == "visibility_merge":
        return build_meshlets_visibility_merge(mesh, opt)
    if opt.builder == "hybrid":
        meshlets = build_meshlets_hybrid(mesh, opt)
        return refine_meshlets(mesh, meshlets, opt)
    if opt.builder != "greedy":
        raise ValueError(f"unknown meshlet builder: {opt.builder}")

    centers = compute_triangle_centers(mesh)
    normals = compute_triangle_normals(mesh)
    adjacency = build_triangle_adjacency(mesh)
    scale = mesh.scale_hint()

    remaining: set[int] = set(range(len(mesh.triangles)))
    meshlets: list[Meshlet] = []
    heartbeat = Heartbeat("greedy meshlets")

    max_normal_cost = 1.0 - float(np.cos(np.radians(opt.max_normal_angle_deg)))

    def candidate_score(
        current: list[int],
        vertices: set[int],
        texcoords: set[int],
        normal_indices: set[int],
        current_radius: float,
        current_center_mean: np.ndarray,
        tri_index: int,
    ) -> float | None:
        tri = mesh.triangles[tri_index]
        if tri.material != mesh.triangles[current[0]].material:
            return None

        v, vt, vn = _triangle_sets(mesh, tri_index)
        new_vertices = vertices | v
        if len(new_vertices) > opt.max_vertices:
            return None
        if len(texcoords | vt) > opt.max_texcoords:
            return None
        if len(normal_indices | vn) > opt.max_normals:
            return None
        if len(current) >= opt.max_triangles:
            return None

        new_center, new_radius = _sphere_from_vertices(mesh.vertices[list(new_vertices)])
        _ = new_center
        radius_cost = max(0.0, new_radius - current_radius) / scale

        normal_axis, normal_min_dot = _normal_cone_fast(normals[current + [tri_index]])
        _ = normal_axis
        normal_cost = 1.0 - normal_min_dot
        if normal_cost > max_normal_cost:
            return None

        new_vertex_cost = len(new_vertices) - len(vertices)
        spatial_hint = float(np.linalg.norm(centers[tri_index] - current_center_mean)) / scale

        return (
            opt.vertex_weight * new_vertex_cost
            + opt.radius_weight * radius_cost
            + opt.normal_weight * normal_cost
            + 0.5 * spatial_hint
        )

    while remaining:
        seed = min(remaining)
        current = [seed]
        remaining.remove(seed)
        vertices, texcoords, normals_set = _triangle_sets(mesh, seed)
        frontier = {n for n in adjacency[seed] if n in remaining}

        while frontier:
            _, current_radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
            current_center_mean = np.mean(centers[current], axis=0)
            scored = []
            for ti in sorted(frontier):
                score = candidate_score(current, vertices, texcoords, normals_set, current_radius, current_center_mean, ti)
                if score is not None:
                    scored.append((score, ti))
            if not scored:
                break
            score, chosen = min(scored, key=lambda x: (x[0], x[1]))
            chosen_vertices, chosen_texcoords, chosen_normals = _triangle_sets(mesh, chosen)
            if len(vertices) >= opt.target_vertices and score > opt.stop_score:
                break

            current.append(chosen)
            remaining.remove(chosen)
            vertices.update(chosen_vertices)
            texcoords.update(chosen_texcoords)
            normals_set.update(chosen_normals)
            frontier.discard(chosen)
            for n in adjacency[chosen]:
                if n in remaining:
                    frontier.add(n)

        meshlets.append(_make_meshlet(mesh, current, normals))
        if len(meshlets) % 250 == 0:
            log(f"greedy meshlets: {len(meshlets)} meshlets, {len(mesh.triangles) - len(remaining)}/{len(mesh.triangles)} triangles")
        else:
            heartbeat.ping(f"{len(meshlets)} meshlets, {len(mesh.triangles) - len(remaining)}/{len(mesh.triangles)} triangles")

    return refine_meshlets(mesh, meshlets, opt)


def sort_meshlets_by_material(meshlets: list[Meshlet]) -> list[Meshlet]:
    """Return meshlets grouped by material while preserving local order inside each group."""
    return [meshlet for _, meshlet in sorted(enumerate(meshlets), key=lambda item: (item[1].material, item[0]))]


def _cone_angle(cos_angle: float) -> float:
    return float(math.acos(max(-1.0, min(1.0, cos_angle))))


def _meshlet_with_visibility(
    meshlet: Meshlet,
    visibility_axis: np.ndarray,
    visibility_cos: float,
    visibility_samples: int = 0,
) -> Meshlet:
    return replace(
        meshlet,
        visibility_axis=visibility_axis,
        visibility_cos=visibility_cos,
        visibility_cull_axis=-visibility_axis,
        visibility_cull_cos=-visibility_cos,
        visibility_samples=visibility_samples,
    )


def _visibility_cone_from_mask(views: np.ndarray, mask: np.ndarray, margin_deg: float, *, fast: bool) -> tuple[np.ndarray, float, int]:
    from .cones import direction_cone, expand_cone_cos

    visible_count = int(np.count_nonzero(mask))
    if visible_count == 0:
        return np.array([0.0, 0.0, 1.0], dtype=np.float64), -1.0, 0
    # The TGX helper renders along camera->object directions. Runtime meshlet cones
    # store object->camera visibility directions, hence the sign flip.
    visible_camera_dirs = -views[mask]
    if fast:
        axis, cos_angle = _normal_cone_fast(visible_camera_dirs)
    else:
        axis, cos_angle = direction_cone(visible_camera_dirs)
    return axis, expand_cone_cos(cos_angle, margin_deg), visible_count


def _apply_mask_visibility_cones(
    meshlets: list[Meshlet],
    views: np.ndarray,
    masks: list[np.ndarray],
    margin_deg: float,
    *,
    fast: bool,
) -> list[Meshlet]:
    updated: list[Meshlet] = []
    for meshlet, mask in zip(meshlets, masks):
        axis, cos_angle, visible_count = _visibility_cone_from_mask(views, mask, margin_deg, fast=fast)
        updated.append(_meshlet_with_visibility(meshlet, axis, cos_angle, visible_count))
    return updated


def _triangle_visibility_meshlets(mesh: ObjMesh, opt: MeshletBuildOptions, tri_normals: np.ndarray) -> tuple[list[Meshlet], np.ndarray, list[np.ndarray], float, float]:
    """Create one meshlet per triangle and keep their TGX visibility masks."""
    from .cones import auto_visibility_margin_deg
    from .quality import fibonacci_sphere
    from .visibility import compute_tgx_visibility

    meshlets = [_make_meshlet(mesh, [ti], tri_normals) for ti in range(len(mesh.triangles))]
    views = fibonacci_sphere(max(1, int(opt.visibility_merge_samples)))
    margin = float(opt.visibility_merge_margin_deg)
    if margin < 0.0:
        margin = auto_visibility_margin_deg(len(views))
    print("visibility-merge seed cones:")
    print(f"  triangle seeds: {len(meshlets)}")
    print(f"  views         : {len(views)}")
    print(f"  render size   : {opt.visibility_merge_size}x{opt.visibility_merge_size}")
    print(f"  margin        : {margin:.2f} deg" + (" (auto)" if opt.visibility_merge_margin_deg < 0.0 else ""))
    t0 = time.perf_counter()
    visibility = compute_tgx_visibility(
        mesh,
        meshlets,
        views,
        width=int(opt.visibility_merge_size),
        height=int(opt.visibility_merge_size),
    )
    tgx_time = time.perf_counter() - t0
    visible_counts = np.sum(visibility, axis=0)
    print(f"  visible views : min {int(np.min(visible_counts))}, avg {float(np.mean(visible_counts)):.1f}, max {int(np.max(visible_counts))}")
    print(f"  never visible : {int(np.sum(visible_counts == 0))}")
    print(f"  TGX time      : {tgx_time:.3f} s")
    masks = [visibility[:, ti].astype(bool, copy=True) for ti in range(visibility.shape[1])]
    t1 = time.perf_counter()
    meshlets = _apply_mask_visibility_cones(meshlets, views, masks, margin, fast=True)
    cone_time = time.perf_counter() - t1
    print(f"  seed cones    : {cone_time:.3f} s")
    return meshlets, views, masks, margin, tgx_time


def _visibility_pair_links(mesh: ObjMesh, a: Meshlet, b: Meshlet, compatible_edges: dict[tuple, list[int]], geometric_edges: dict[tuple[int, int], list[int]]) -> tuple[int, int]:
    b_tris = set(b.triangles)
    compatible = 0
    geometric = 0
    for ti in a.triangles:
        tri = mesh.triangles[ti]
        c = tri.corners
        for key in (_corner_edge_key(c[0], c[1]), _corner_edge_key(c[1], c[2]), _corner_edge_key(c[2], c[0])):
            if any(tj in b_tris for tj in compatible_edges.get(key, [])):
                compatible += 1
        v = [corner.v for corner in tri.corners]
        for key in (_edge_key(v[0], v[1]), _edge_key(v[1], v[2]), _edge_key(v[2], v[0])):
            if any(tj in b_tris for tj in geometric_edges.get(key, [])):
                geometric += 1
    return compatible, max(0, geometric - compatible)


def _visibility_merge_candidate(
    mesh: ObjMesh,
    a: Meshlet,
    b: Meshlet,
    mask_a: np.ndarray,
    mask_b: np.ndarray,
    nb_views: int,
    tri_normals: np.ndarray,
    opt: MeshletBuildOptions,
    scale: float,
    compatible_edges: dict[tuple, list[int]],
    geometric_edges: dict[tuple[int, int], list[int]],
) -> tuple[float, Meshlet, np.ndarray] | None:
    if a.material != b.material:
        return None
    tri_indices = a.triangles + b.triangles
    if len(tri_indices) > opt.max_triangles:
        return None
    vertices = a.vertices | b.vertices
    texcoords = a.texcoords | b.texcoords
    normals = a.normals | b.normals
    if len(vertices) > opt.max_vertices or len(texcoords) > opt.max_texcoords or len(normals) > opt.max_normals:
        return None

    normal_axis, normal_min_dot = _normal_cone_fast(tri_normals[tri_indices])
    normal_angle = _cone_angle(normal_min_dot)
    if math.degrees(normal_angle) > opt.max_normal_angle_deg:
        return None

    candidate_mask = np.logical_or(mask_a, mask_b)
    visible_count = int(np.count_nonzero(candidate_mask))
    overlap_count = int(np.count_nonzero(np.logical_and(mask_a, mask_b)))
    visible_growth = max(0, visible_count - max(a.visibility_samples, b.visibility_samples)) / max(1, nb_views)
    overlap_bonus = overlap_count / max(1, nb_views)

    old_radius = max(a.radius, b.radius)
    center, radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
    radius_growth = max(0.0, radius - old_radius) / scale
    normal_cost = 1.0 - normal_min_dot
    new_vertex_cost = len(vertices) - max(len(a.vertices), len(b.vertices))
    compatible_links, incompatible_links = _visibility_pair_links(mesh, a, b, compatible_edges, geometric_edges)
    shared_refs = len(a.vertices & b.vertices) + len(a.texcoords & b.texcoords) + len(a.normals & b.normals)

    score = (
        opt.vertex_weight * new_vertex_cost
        + opt.radius_weight * radius_growth
        + opt.normal_weight * normal_cost
        + opt.visibility_merge_cone_weight * visible_growth
        + opt.incompatible_edge_penalty * incompatible_links
        - opt.compatible_edge_weight * compatible_links
        - 0.10 * opt.visibility_merge_cone_weight * overlap_bonus
        - 0.05 * shared_refs
    )
    if len(vertices) >= opt.target_vertices and score > opt.stop_score:
        return None

    cull_axis, cull_cos = _backface_cull_cone(normal_axis, normal_min_dot)
    candidate = Meshlet(
        tri_indices,
        a.material,
        vertices,
        texcoords,
        normals,
        center,
        radius,
        normal_axis,
        normal_min_dot,
        cull_axis,
        cull_cos,
    )
    visibility_axis = a.visibility_axis if a.visibility_axis is not None else np.array([0.0, 0.0, 1.0], dtype=np.float64)
    candidate = _meshlet_with_visibility(candidate, visibility_axis, -1.0, visible_count)
    return score, candidate, candidate_mask


def _visibility_pack_mode(opt: MeshletBuildOptions) -> str:
    mode = str(opt.visibility_pack_mode or "none").strip().lower().replace("-", "_")
    if mode in ("", "0", "false", "no", "none", "off", "disabled"):
        return "none"
    compression = max(0.0, min(100.0, float(opt.visibility_pack_compression)))
    if compression <= 0.0:
        return "none"
    if mode == "compact" and compression >= 100.0:
        return "ultra"
    if mode in ("compact", "ultra"):
        return mode
    raise ValueError(f"unknown visibility non-adjacent packing mode: {opt.visibility_pack_mode!r}")


def _visibility_pack_thresholds(opt: MeshletBuildOptions) -> tuple[float, float, float, float]:
    compression = max(0.0, min(100.0, float(opt.visibility_pack_compression)))
    if compression <= 50.0:
        t = compression / 50.0
        return (
            1.0 - 0.25 * t,  # bad-visible ratio: 1.00 -> 0.75
            0.08 * t,        # normal visible growth: 0.00 -> 0.08
            0.18 * t,        # relaxed visible growth: 0.00 -> 0.18
            0.15 * t,        # radius growth: 0.00 -> 0.15
        )
    t = (compression - 50.0) / 50.0
    return (
        0.75 * (1.0 - t),          # bad-visible ratio: 0.75 -> 0.00
        0.08 + (1.00 - 0.08) * t,  # normal visible growth: 0.08 -> 1.00
        0.18 + (1.00 - 0.18) * t,  # relaxed visible growth: 0.18 -> 1.00
        0.15 + (10.0 - 0.15) * t,  # radius growth: 0.15 -> effectively unbounded
    )


def _visibility_pack_candidate(
    mesh: ObjMesh,
    a: Meshlet,
    b: Meshlet,
    mask_a: np.ndarray,
    mask_b: np.ndarray,
    nb_views: int,
    tri_normals: np.ndarray,
    opt: MeshletBuildOptions,
    scale: float,
    mode: str,
) -> tuple[float, Meshlet, np.ndarray] | None:
    if a.material != b.material:
        return None

    vertices = a.vertices | b.vertices
    texcoords = a.texcoords | b.texcoords
    normals = a.normals | b.normals
    if len(vertices) > opt.max_vertices or len(texcoords) > opt.max_texcoords or len(normals) > opt.max_normals:
        return None

    candidate_mask = np.logical_or(mask_a, mask_b)
    visible_a = a.visibility_samples if a.visibility_samples > 0 else int(np.count_nonzero(mask_a))
    visible_b = b.visibility_samples if b.visibility_samples > 0 else int(np.count_nonzero(mask_b))
    visible_count = int(np.count_nonzero(candidate_mask))
    nb_views = max(1, nb_views)
    visible_growth = max(0, visible_count - max(visible_a, visible_b)) / nb_views
    union_visible_ratio = visible_count / nb_views
    min_visible_ratio = min(visible_a, visible_b) / nb_views

    old_radius = max(float(a.radius), float(b.radius), 1.0e-20)
    safe_scale = max(float(scale), 1.0e-20)
    center, radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
    radius_growth = max(0.0, float(radius) - old_radius) / safe_scale
    radius_ratio_growth = max(0.0, float(radius) / old_radius - 1.0)

    if mode == "compact":
        bad_visible_ratio, max_growth, bad_max_growth, max_radius_growth = _visibility_pack_thresholds(opt)
        already_bad = min_visible_ratio >= bad_visible_ratio
        low_visibility_cost = visible_growth <= max_growth
        close_enough = radius_growth <= max_radius_growth
        relaxed_visibility_cost = visible_growth <= bad_max_growth
        if not (low_visibility_cost or (already_bad and relaxed_visibility_cost) or (close_enough and relaxed_visibility_cost)):
            return None
    elif mode != "ultra":
        return None

    tri_indices = a.triangles + b.triangles
    normal_axis, normal_min_dot = _normal_cone_fast(tri_normals[tri_indices])
    cull_axis, cull_cos = _backface_cull_cone(normal_axis, normal_min_dot)
    candidate = Meshlet(
        tri_indices,
        a.material,
        vertices,
        texcoords,
        normals,
        center,
        radius,
        normal_axis,
        normal_min_dot,
        cull_axis,
        cull_cos,
    )
    visibility_axis = a.visibility_axis if a.visibility_axis is not None else np.array([0.0, 0.0, 1.0], dtype=np.float64)
    candidate = _meshlet_with_visibility(candidate, visibility_axis, -1.0, visible_count)

    fill = max(
        len(vertices) / max(1, opt.max_vertices),
        len(texcoords) / max(1, opt.max_texcoords),
        len(normals) / max(1, opt.max_normals),
    )
    tri_fill = min(1.0, len(tri_indices) / max(1, opt.max_triangles))
    shared_refs = len(a.vertices & b.vertices) + len(a.texcoords & b.texcoords) + len(a.normals & b.normals)
    new_ref_cost = (
        (len(vertices) - max(len(a.vertices), len(b.vertices)))
        + (len(texcoords) - max(len(a.texcoords), len(b.texcoords)))
        + (len(normals) - max(len(a.normals), len(b.normals)))
    )
    if mode == "ultra":
        score = -120.0 * fill - 12.0 * tri_fill + 4.0 * radius_growth + 0.5 * union_visible_ratio + 0.08 * new_ref_cost - 0.04 * shared_refs
    else:
        score = (
            120.0 * visible_growth
            + 18.0 * radius_growth
            + 4.0 * radius_ratio_growth
            + 3.0 * union_visible_ratio
            + 0.12 * new_ref_cost
            - 18.0 * fill
            - 4.0 * tri_fill
            - 5.0 * min_visible_ratio
            - 0.06 * shared_refs
        )
    return score, candidate, candidate_mask


def _visibility_nonadjacent_pack_pass(
    mesh: ObjMesh,
    meshlets: list[Meshlet],
    masks: list[np.ndarray],
    views: np.ndarray,
    tri_normals: np.ndarray,
    opt: MeshletBuildOptions,
    scale: float,
) -> tuple[list[Meshlet], list[np.ndarray], int, float, float, int]:
    mode = _visibility_pack_mode(opt)
    if mode == "none" or len(meshlets) < 2:
        return meshlets, masks, 0, 0.0, 0.0, 0

    t_total = time.perf_counter()
    print(f"visibility-pack: mode={mode}, start {len(meshlets)} meshlets")
    alive: set[int] = set(range(len(meshlets)))
    versions = [0 for _ in meshlets]
    heap: list[tuple[float, int, int, int, int, int, Meshlet, np.ndarray]] = []
    serial = 0
    candidate_time = 0.0
    candidate_count = 0

    def push_pair(i: int, j: int) -> None:
        nonlocal serial, candidate_time, candidate_count
        if i == j or i not in alive or j not in alive:
            return
        if i > j:
            i, j = j, i
        if meshlets[i].material != meshlets[j].material:
            return
        t_candidate = time.perf_counter()
        result = _visibility_pack_candidate(
            mesh,
            meshlets[i],
            meshlets[j],
            masks[i],
            masks[j],
            len(views),
            tri_normals,
            opt,
            scale,
            mode,
        )
        candidate_time += time.perf_counter() - t_candidate
        candidate_count += 1
        if result is None:
            return
        score, candidate, candidate_mask = result
        heapq.heappush(heap, (score, versions[i], versions[j], serial, i, j, candidate, candidate_mask))
        serial += 1

    by_material: dict[str, list[int]] = defaultdict(list)
    for idx, meshlet in enumerate(meshlets):
        by_material[meshlet.material].append(idx)
    for group in by_material.values():
        for pos, i in enumerate(group):
            for j in group[pos + 1 :]:
                push_pair(i, j)

    merges = 0
    while heap:
        score, version_i, version_j, _, i, j, candidate, candidate_mask = heapq.heappop(heap)
        if i not in alive or j not in alive:
            continue
        if versions[i] != version_i or versions[j] != version_j:
            continue

        keep, kill = (i, j) if len(meshlets[i].triangles) >= len(meshlets[j].triangles) else (j, i)
        if keep != i:
            t_candidate = time.perf_counter()
            result = _visibility_pack_candidate(
                mesh,
                meshlets[keep],
                meshlets[kill],
                masks[keep],
                masks[kill],
                len(views),
                tri_normals,
                opt,
                scale,
                mode,
            )
            candidate_time += time.perf_counter() - t_candidate
            candidate_count += 1
            if result is None:
                continue
            _, candidate, candidate_mask = result

        meshlets[keep] = candidate
        masks[keep] = candidate_mask
        versions[keep] += 1
        alive.remove(kill)
        versions[kill] += 1
        merges += 1
        if (merges % 100) == 0:
            print(f"visibility-pack: {merges} merges, {len(alive)} meshlets, best score {score:.3f}")

        for n in list(alive):
            if n != keep and meshlets[n].material == candidate.material:
                push_pair(keep, n)

    ordered_alive = sorted(alive, key=lambda idx: min(meshlets[idx].triangles))
    packed_meshlets = [meshlets[i] for i in ordered_alive]
    packed_masks = [masks[i] for i in ordered_alive]
    total_time = time.perf_counter() - t_total
    print(f"visibility-pack: {merges} merges, {len(packed_meshlets)} meshlets")
    return packed_meshlets, packed_masks, merges, total_time, candidate_time, candidate_count


def build_meshlets_visibility_merge(mesh: ObjMesh, options: MeshletBuildOptions | None = None) -> list[Meshlet]:
    """Build meshlets by merging triangle visibility masks computed by the TGX helper.

    This experimental builder starts with one triangle per meshlet, renders a TGX visibility
    matrix for those seeds, then greedily merges adjacent meshlets. Candidate visibility is
    computed as the exact discrete union of the triangle visibility masks.
    """
    t_total = time.perf_counter()
    opt = options or MeshletBuildOptions(builder="visibility_merge")
    pack_mode = _visibility_pack_mode(opt)
    tri_normals = compute_triangle_normals(mesh)
    meshlets, views, masks, visibility_margin, tgx_time = _triangle_visibility_meshlets(mesh, opt, tri_normals)
    scale = mesh.scale_hint()
    tri_adj = build_triangle_adjacency(mesh)
    geometric_edges = build_edge_to_triangles(mesh)
    compatible_edges = build_compatible_edge_to_triangles(mesh)

    alive: set[int] = set(range(len(meshlets)))
    versions = [0 for _ in meshlets]
    tri_to_meshlet = {ti: ti for ti in range(len(mesh.triangles))}
    adjacency: list[set[int]] = [set() for _ in meshlets]
    for ti, neighbors in enumerate(tri_adj):
        for tj in neighbors:
            if mesh.triangles[ti].material == mesh.triangles[tj].material:
                adjacency[ti].add(tj)

    heap: list[tuple[float, int, int, int, int, int, Meshlet, np.ndarray]] = []
    serial = 0
    candidate_time = 0.0
    candidate_count = 0

    def push_pair(i: int, j: int) -> None:
        nonlocal serial, candidate_time, candidate_count
        if i == j or i not in alive or j not in alive:
            return
        if i > j:
            i, j = j, i
        t_candidate = time.perf_counter()
        result = _visibility_merge_candidate(
            mesh,
            meshlets[i],
            meshlets[j],
            masks[i],
            masks[j],
            len(views),
            tri_normals,
            opt,
            scale,
            compatible_edges,
            geometric_edges,
        )
        candidate_time += time.perf_counter() - t_candidate
        candidate_count += 1
        if result is None:
            return
        score, candidate, candidate_mask = result
        heapq.heappush(heap, (score, versions[i], versions[j], serial, i, j, candidate, candidate_mask))
        serial += 1

    for i in range(len(meshlets)):
        for j in adjacency[i]:
            if i < j:
                push_pair(i, j)

    t_merge = time.perf_counter()
    merges = 0
    while heap:
        score, version_i, version_j, _, i, j, candidate, candidate_mask = heapq.heappop(heap)
        if i not in alive or j not in alive:
            continue
        if versions[i] != version_i or versions[j] != version_j:
            continue

        keep, kill = (i, j) if len(meshlets[i].triangles) >= len(meshlets[j].triangles) else (j, i)
        if keep != i:
            t_candidate = time.perf_counter()
            candidate = _visibility_merge_candidate(
                mesh,
                meshlets[keep],
                meshlets[kill],
                masks[keep],
                masks[kill],
                len(views),
                tri_normals,
                opt,
                scale,
                compatible_edges,
                geometric_edges,
            )
            candidate_time += time.perf_counter() - t_candidate
            candidate_count += 1
            if candidate is None:
                continue
            _, candidate_meshlet, candidate_mask = candidate
        else:
            candidate_meshlet = candidate

        old_neighbors = (adjacency[keep] | adjacency[kill]) - {keep, kill}
        meshlets[keep] = candidate_meshlet
        masks[keep] = candidate_mask
        versions[keep] += 1
        alive.remove(kill)
        versions[kill] += 1
        adjacency[keep] = {n for n in old_neighbors if n in alive}
        adjacency[kill].clear()
        for n in list(old_neighbors):
            adjacency[n].discard(kill)
            adjacency[n].discard(keep)
            if n in alive:
                adjacency[n].add(keep)
        for ti in candidate_meshlet.triangles:
            tri_to_meshlet[ti] = keep
        merges += 1
        if (merges % 250) == 0:
            print(f"visibility-merge: {merges} merges, {len(alive)} meshlets, best score {score:.3f}")
        for n in adjacency[keep]:
            push_pair(keep, n)

    merge_time = time.perf_counter() - t_merge
    ordered_alive = sorted(alive, key=lambda idx: min(meshlets[idx].triangles))
    result = [meshlets[i] for i in ordered_alive]
    result_masks = [masks[i] for i in ordered_alive]
    pack_merges = 0
    pack_time = 0.0
    pack_candidate_time = 0.0
    pack_candidate_count = 0
    if pack_mode != "none":
        result, result_masks, pack_merges, pack_time, pack_candidate_time, pack_candidate_count = _visibility_nonadjacent_pack_pass(
            mesh,
            result,
            result_masks,
            views,
            tri_normals,
            opt,
            scale,
        )
    t_final = time.perf_counter()
    result = _apply_mask_visibility_cones(result, views, result_masks, visibility_margin, fast=False)
    final_cone_time = time.perf_counter() - t_final
    if pack_mode == "none":
        print(f"visibility-merge: {merges} merges, {len(result)} meshlets")
    else:
        print(f"visibility-merge: {merges} adjacent merges, {pack_merges} non-adjacent merges, {len(result)} meshlets")
    total_time = time.perf_counter() - t_total
    print("visibility-merge timings:")
    print(f"  TGX seed masks : {tgx_time:.3f} s")
    print(f"  merge loop     : {merge_time:.3f} s")
    print(f"  candidates     : {candidate_time:.3f} s for {candidate_count}")
    if pack_mode != "none":
        print(f"  pack loop      : {pack_time:.3f} s")
        print(f"  pack candidates: {pack_candidate_time:.3f} s for {pack_candidate_count}")
    print(f"  final cones    : {final_cone_time:.3f} s")
    print(f"  total builder  : {total_time:.3f} s")
    return result


def _meshlets_share_compatible_edge(mesh: ObjMesh, a: Meshlet, b: Meshlet) -> bool:
    edge_to_tris = build_compatible_edge_to_triangles(mesh)
    b_tris = set(b.triangles)
    for ti in a.triangles:
        tri = mesh.triangles[ti]
        c = tri.corners
        for key in (_corner_edge_key(c[0], c[1]), _corner_edge_key(c[1], c[2]), _corner_edge_key(c[2], c[0])):
            if any(tj in b_tris for tj in edge_to_tris.get(key, [])):
                return True
    return False


def build_meshlets_hybrid(mesh: ObjMesh, options: MeshletBuildOptions | None = None) -> list[Meshlet]:
    """Build meshlets with a Meshlete-inspired border growth strategy.

    The original greedy builder scores all frontier triangles similarly. This variant is closer
    to Meshlete's generator: triangles are seeded in spatial order along the mesh major axis, then
    each meshlet grows from its active border, strongly preferring triangles sharing two or three
    vertices with the current meshlet. TGX-specific constraints keep normal cones and local index
    counts under control.
    """
    opt = options or MeshletBuildOptions(builder="hybrid")
    centers = compute_triangle_centers(mesh)
    tri_normals = compute_triangle_normals(mesh)
    adjacency = build_triangle_adjacency(mesh)
    edge_to_triangles = build_edge_to_triangles(mesh)
    compatible_edge_to_triangles = build_compatible_edge_to_triangles(mesh)
    scale = mesh.scale_hint()

    bb_min, bb_max = mesh.bounding_box()
    axis_id = int(np.argmax(bb_max - bb_min))
    axis = np.zeros(3, dtype=np.float64)
    axis[axis_id] = 1.0

    tri_order: list[tuple[float, int]] = []
    for ti in range(len(mesh.triangles)):
        pts = mesh.triangle_vertices(ti)
        tri_order.append((float(np.min(pts @ axis)), ti))
    tri_order.sort()

    remaining: set[int] = set(range(len(mesh.triangles)))
    meshlets: list[Meshlet] = []
    max_normal_cost = 1.0 - float(np.cos(np.radians(opt.max_normal_angle_deg)))
    heartbeat = Heartbeat("hybrid meshlets")

    def triangle_edges(ti: int) -> list[tuple[int, int]]:
        v = [c.v for c in mesh.triangles[ti].corners]
        return [_edge_key(v[0], v[1]), _edge_key(v[1], v[2]), _edge_key(v[2], v[0])]

    def geometric_link_count(ti: int, current_set: set[int]) -> int:
        count = 0
        for key in triangle_edges(ti):
            if any(tj in current_set for tj in edge_to_triangles.get(key, [])):
                count += 1
        return count

    def compatible_link_count(ti: int, current_set: set[int]) -> int:
        tri = mesh.triangles[ti]
        c = tri.corners
        count = 0
        for key in (_corner_edge_key(c[0], c[1]), _corner_edge_key(c[1], c[2]), _corner_edge_key(c[2], c[0])):
            if any(tj in current_set for tj in compatible_edge_to_triangles.get(key, [])):
                count += 1
        return count

    def next_seed() -> int:
        while tri_order:
            _, ti = tri_order.pop(0)
            if ti in remaining:
                return ti
        return min(remaining)

    def score_candidate(
        current: list[int],
        vertices: set[int],
        texcoords: set[int],
        normal_indices: set[int],
        current_set: set[int],
        test_pos: np.ndarray,
        old_radius: float,
        ti: int,
    ) -> tuple[float, set[int], set[int], set[int], int] | None:
        tri = mesh.triangles[ti]
        if tri.material != mesh.triangles[current[0]].material:
            return None

        v, vt, vn = _triangle_sets(mesh, ti)
        shared = len(v & vertices)
        new_vertices = vertices | v
        new_texcoords = texcoords | vt
        new_normals = normal_indices | vn
        if len(new_vertices) > opt.max_vertices:
            return None
        if len(new_texcoords) > opt.max_texcoords:
            return None
        if len(new_normals) > opt.max_normals:
            return None
        if len(current) >= opt.max_triangles:
            return None

        normal_axis, normal_min_dot = _normal_cone_fast(tri_normals[current + [ti]])
        _ = normal_axis
        normal_cost = 1.0 - normal_min_dot
        if normal_cost > max_normal_cost:
            return None

        _, new_radius = _sphere_from_vertices(mesh.vertices[list(new_vertices)])
        radius_cost = max(0.0, new_radius - old_radius) / scale
        test_cost = min(float(np.linalg.norm(mesh.vertices[vi] - test_pos)) for vi in v) / scale
        new_vertex_cost = len(new_vertices) - len(vertices)
        centroid_cost = float(np.linalg.norm(centers[ti] - np.mean(centers[current], axis=0))) / scale
        compatible_links = compatible_link_count(ti, current_set)
        incompatible_links = max(0, geometric_link_count(ti, current_set) - compatible_links)
        score = (
            opt.vertex_weight * new_vertex_cost
            + opt.radius_weight * radius_cost
            + opt.normal_weight * normal_cost
            + opt.test_pos_weight * test_cost
            + 0.75 * centroid_cost
            + opt.incompatible_edge_penalty * incompatible_links
            - opt.shared_vertex_weight * shared
            - opt.compatible_edge_weight * compatible_links
        )
        return score, new_vertices, new_texcoords, new_normals, shared + compatible_links

    while remaining:
        seed = next_seed()
        current = [seed]
        remaining.remove(seed)
        vertices, texcoords, normals_set = _triangle_sets(mesh, seed)
        current_set = {seed}
        test_pos = centers[seed].copy()
        border_edges = set(triangle_edges(seed))

        while border_edges:
            _, old_radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
            candidates: dict[int, tuple[float, set[int], set[int], set[int], int]] = {}
            dead_edges = []
            for edge in sorted(border_edges):
                found = False
                for ti in edge_to_triangles.get(edge, []):
                    if ti not in remaining:
                        continue
                    found = True
                    scored = score_candidate(current, vertices, texcoords, normals_set, current_set, test_pos, old_radius, ti)
                    if scored is not None and (ti not in candidates or scored[0] < candidates[ti][0]):
                        candidates[ti] = scored
                if not found:
                    dead_edges.append(edge)
            for edge in dead_edges:
                border_edges.discard(edge)
            if not candidates:
                break

            chosen, data = min(candidates.items(), key=lambda item: (item[1][0], -item[1][4], item[0]))
            score, new_vertices, new_texcoords, new_normals, shared = data
            if len(vertices) >= opt.target_vertices and score > opt.stop_score and shared < 3:
                break

            current.append(chosen)
            remaining.remove(chosen)
            current_set.add(chosen)
            vertices, texcoords, normals_set = new_vertices, new_texcoords, new_normals
            alpha = 1.0 / (len(current) + 1.0)
            test_pos = (1.0 - alpha) * test_pos + alpha * centers[chosen]
            for edge in triangle_edges(chosen):
                border_edges.add(edge)

        if opt.fill_holes:
            # Cheap Meshlete-like cleanup: add adjacent unassigned triangles whose vertices are
            # already local. This usually improves stripification without increasing vertex count.
            changed = True
            while changed and len(current) < opt.max_triangles:
                changed = False
                for base in list(current):
                    for ti in sorted(adjacency[base]):
                        if ti not in remaining:
                            continue
                        if opt.compatible_edge_weight > 0.0 and compatible_link_count(ti, current_set) <= 0:
                            continue
                        v, vt, vn = _triangle_sets(mesh, ti)
                        if not v.issubset(vertices):
                            continue
                        if len(texcoords | vt) > opt.max_texcoords or len(normals_set | vn) > opt.max_normals:
                            continue
                        normal_axis, normal_min_dot = _normal_cone_fast(tri_normals[current + [ti]])
                        _ = normal_axis
                        if 1.0 - normal_min_dot > max_normal_cost:
                            continue
                        current.append(ti)
                        remaining.remove(ti)
                        current_set.add(ti)
                        texcoords.update(vt)
                        normals_set.update(vn)
                        changed = True
                        if len(current) >= opt.max_triangles:
                            break
                    if len(current) >= opt.max_triangles:
                        break

        meshlets.append(_make_meshlet(mesh, current, tri_normals))
        if len(meshlets) % 250 == 0:
            log(f"hybrid meshlets: {len(meshlets)} meshlets, {len(mesh.triangles) - len(remaining)}/{len(mesh.triangles)} triangles")
        else:
            heartbeat.ping(f"{len(meshlets)} meshlets, {len(mesh.triangles) - len(remaining)}/{len(mesh.triangles)} triangles")

    return meshlets


def refine_meshlets(mesh: ObjMesh, meshlets: list[Meshlet], options: MeshletBuildOptions) -> list[Meshlet]:
    """Optional post-process inspired by the old prototype.

    Merging removes tiny adjacent meshlets when the combined meshlet still satisfies the
    local-index and normal-coherence constraints. Smoothing then moves boundary triangles to
    neighboring meshlets if it improves local connectivity while preserving constraints.
    """
    tri_normals = compute_triangle_normals(mesh)
    for pass_index in range(max(0, options.merge_passes)):
        before = len(meshlets)
        with step("meshlet merge pass", f"{pass_index + 1}/{max(0, options.merge_passes)}, {before} meshlets"):
            meshlets, changed = _merge_meshlet_pass(mesh, meshlets, tri_normals, options)
        log(f"meshlet merge pass {pass_index + 1}: {before} -> {len(meshlets)} meshlets")
        if not changed:
            break
    for pass_index in range(max(0, options.smooth_passes)):
        with step("meshlet smooth pass", f"{pass_index + 1}/{max(0, options.smooth_passes)}, {len(meshlets)} meshlets"):
            meshlets, changed = _smooth_meshlet_pass(mesh, meshlets, tri_normals, options)
        log(f"meshlet smooth pass {pass_index + 1}: {'changed' if changed else 'no change'}")
        if not changed:
            break
    return meshlets


def _meshlet_from_tris_if_valid(
    mesh: ObjMesh,
    tri_indices: list[int],
    tri_normals: np.ndarray,
    options: MeshletBuildOptions,
    max_normal_angle_deg: float | None = None,
) -> Meshlet | None:
    if not tri_indices:
        return None
    if len(tri_indices) > options.max_triangles:
        return None
    material = mesh.triangles[tri_indices[0]].material
    vertices: set[int] = set()
    texcoords: set[int] = set()
    normals: set[int] = set()
    for ti in tri_indices:
        if mesh.triangles[ti].material != material:
            return None
        v, vt, vn = _triangle_sets(mesh, ti)
        vertices.update(v)
        texcoords.update(vt)
        normals.update(vn)
    if len(vertices) > options.max_vertices or len(texcoords) > options.max_texcoords or len(normals) > options.max_normals:
        return None
    normal_axis, normal_min_dot = _normal_cone(tri_normals[tri_indices])
    angle = float(np.degrees(np.arccos(np.clip(normal_min_dot, -1.0, 1.0))))
    if angle > (max_normal_angle_deg if max_normal_angle_deg is not None else options.max_normal_angle_deg):
        return None
    center, radius = _sphere_from_vertices(mesh.vertices[list(vertices)])
    cull_axis, cull_cos = _backface_cull_cone(normal_axis, normal_min_dot)
    return Meshlet(tri_indices.copy(), material, vertices, texcoords, normals, center, radius, normal_axis, normal_min_dot, cull_axis, cull_cos)


def _meshlet_adjacency(mesh: ObjMesh, meshlets: list[Meshlet], tri_to_meshlet: dict[int, int]) -> list[set[int]]:
    tri_adj = build_triangle_adjacency(mesh)
    adj = [set() for _ in meshlets]
    for mi, meshlet in enumerate(meshlets):
        for ti in meshlet.triangles:
            for tj in tri_adj[ti]:
                mj = tri_to_meshlet.get(tj)
                if mj is not None and mj != mi:
                    adj[mi].add(mj)
    return adj


def _compatible_meshlet_links(mesh: ObjMesh, a: Meshlet, b: Meshlet) -> int:
    edge_to_triangles = build_compatible_edge_to_triangles(mesh)
    b_tris = set(b.triangles)
    links: set[tuple[int, int]] = set()
    for ti in a.triangles:
        tri = mesh.triangles[ti]
        c = tri.corners
        for key in (_corner_edge_key(c[0], c[1]), _corner_edge_key(c[1], c[2]), _corner_edge_key(c[2], c[0])):
            for tj in edge_to_triangles.get(key, []):
                if tj in b_tris:
                    links.add((min(ti, tj), max(ti, tj)))
    return len(links)


def _tri_to_meshlet(meshlets: list[Meshlet]) -> dict[int, int]:
    return {ti: mi for mi, meshlet in enumerate(meshlets) for ti in meshlet.triangles}


def _is_connected_subset(adjacency: list[set[int]], triangles: list[int]) -> bool:
    if len(triangles) <= 1:
        return True
    allowed = set(triangles)
    stack = [triangles[0]]
    seen = {triangles[0]}
    while stack:
        ti = stack.pop()
        for tj in adjacency[ti]:
            if tj in allowed and tj not in seen:
                seen.add(tj)
                stack.append(tj)
    return len(seen) == len(allowed)


def _merge_meshlet_pass(
    mesh: ObjMesh,
    meshlets: list[Meshlet],
    tri_normals: np.ndarray,
    options: MeshletBuildOptions,
) -> tuple[list[Meshlet], bool]:
    tri_to_meshlet = _tri_to_meshlet(meshlets)
    mesh_adj = _meshlet_adjacency(mesh, meshlets, tri_to_meshlet)
    alive = [True] * len(meshlets)
    changed = False
    heartbeat = Heartbeat("meshlet merge pass")

    ordered_meshlets = sorted(range(len(meshlets)), key=lambda k: len(meshlets[k].triangles))
    for processed, mi in enumerate(ordered_meshlets, start=1):
        if not alive[mi]:
            heartbeat.ping(f"{processed}/{len(ordered_meshlets)} scanned")
            continue
        best = None
        for mj in sorted(mesh_adj[mi], key=lambda k: len(meshlets[k].triangles)):
            if not alive[mj] or mi == mj:
                continue
            merged = meshlets[mj].triangles + meshlets[mi].triangles
            candidate = _meshlet_from_tris_if_valid(mesh, merged, tri_normals, options, options.merge_max_normal_angle_deg)
            if candidate is None:
                continue
            # Prefer compact merges with fewer vertices per triangle.
            compatible_links = _compatible_meshlet_links(mesh, meshlets[mi], meshlets[mj])
            score = candidate.radius / mesh.scale_hint() + 0.01 * len(candidate.vertices) - options.compatible_merge_weight * compatible_links
            if best is None or score < best[0]:
                best = (score, mj, candidate)
        if best is None:
            continue
        _, mj, candidate = best
        meshlets[mj] = candidate
        alive[mi] = False
        changed = True
        for ti in candidate.triangles:
            tri_to_meshlet[ti] = mj
        heartbeat.ping(f"{processed}/{len(ordered_meshlets)} scanned, merged {sum(1 for item in alive if not item)}")

    if not changed:
        return meshlets, False
    return [m for i, m in enumerate(meshlets) if alive[i]], True


def _smooth_meshlet_pass(
    mesh: ObjMesh,
    meshlets: list[Meshlet],
    tri_normals: np.ndarray,
    options: MeshletBuildOptions,
) -> tuple[list[Meshlet], bool]:
    tri_adj = build_triangle_adjacency(mesh)
    tri_to_meshlet = _tri_to_meshlet(meshlets)
    scale = mesh.scale_hint()
    moved = False

    def local_degree(ti: int, mi: int) -> int:
        return sum(1 for tj in tri_adj[ti] if tri_to_meshlet.get(tj) == mi)

    candidates = sorted(tri_to_meshlet, key=lambda ti: (local_degree(ti, tri_to_meshlet[ti]), len(meshlets[tri_to_meshlet[ti]].triangles)))
    heartbeat = Heartbeat("meshlet smooth pass")
    for processed, ti in enumerate(candidates, start=1):
        src = tri_to_meshlet.get(ti)
        if src is None or len(meshlets[src].triangles) <= 1:
            heartbeat.ping(f"{processed}/{len(candidates)} triangles scanned")
            continue
        neighbor_meshlets = sorted({tri_to_meshlet[tj] for tj in tri_adj[ti] if tri_to_meshlet.get(tj) != src})
        if not neighbor_meshlets:
            heartbeat.ping(f"{processed}/{len(candidates)} triangles scanned")
            continue

        src_without = [t for t in meshlets[src].triangles if t != ti]
        if not _is_connected_subset(tri_adj, src_without):
            heartbeat.ping(f"{processed}/{len(candidates)} triangles scanned")
            continue

        old_src = meshlets[src]
        best = None
        for dst in neighbor_meshlets:
            dst_with = meshlets[dst].triangles + [ti]
            new_dst = _meshlet_from_tris_if_valid(mesh, dst_with, tri_normals, options)
            if new_dst is None:
                continue
            new_src = _meshlet_from_tris_if_valid(mesh, src_without, tri_normals, options)
            if new_src is None:
                continue
            old_degree = local_degree(ti, src)
            new_degree = local_degree(ti, dst)
            if new_degree < old_degree:
                continue
            old_score = old_src.radius + meshlets[dst].radius
            new_score = new_src.radius + new_dst.radius
            improvement = (new_degree - old_degree) + 0.25 * ((old_score - new_score) / scale)
            if improvement <= 0.0:
                continue
            if best is None or improvement > best[0]:
                best = (improvement, dst, new_src, new_dst)
        if best is None:
            heartbeat.ping(f"{processed}/{len(candidates)} triangles scanned")
            continue
        _, dst, new_src, new_dst = best
        meshlets[src] = new_src
        meshlets[dst] = new_dst
        tri_to_meshlet[ti] = dst
        moved = True
        heartbeat.ping(f"{processed}/{len(candidates)} triangles scanned")

    if not moved:
        return meshlets, False
    return [m for m in meshlets if len(m.triangles) > 0], True


def meshlet_stats(meshlets: list[Meshlet]) -> dict[str, float]:
    if not meshlets:
        return {}
    tri_counts = np.asarray([len(m.triangles) for m in meshlets], dtype=np.float64)
    vertex_counts = np.asarray([len(m.vertices) for m in meshlets], dtype=np.float64)
    radii = np.asarray([m.radius for m in meshlets], dtype=np.float64)
    angles = np.asarray([m.normal_angle_deg() for m in meshlets], dtype=np.float64)
    cull_angles = np.asarray([m.cull_angle_deg() for m in meshlets], dtype=np.float64)
    cull_empty = np.asarray([m.cull_cos > 1.0 for m in meshlets], dtype=bool)
    return {
        "meshlets": float(len(meshlets)),
        "triangles_min": float(np.min(tri_counts)),
        "triangles_avg": float(np.mean(tri_counts)),
        "triangles_max": float(np.max(tri_counts)),
        "vertices_min": float(np.min(vertex_counts)),
        "vertices_avg": float(np.mean(vertex_counts)),
        "vertices_max": float(np.max(vertex_counts)),
        "radius_avg": float(np.mean(radii)),
        "normal_angle_avg": float(np.mean(angles)),
        "normal_angle_max": float(np.max(angles)),
        "cull_angle_avg": float(np.mean(cull_angles)),
        "cull_angle_max": float(np.max(cull_angles)),
        "cull_angle_min": float(np.min(cull_angles)),
        "cull_empty": float(np.sum(cull_empty)),
    }
