#!/usr/bin/env python3
"""Audit every LOD0 mesh triangle in a cavern corridor for clustered defects."""
from __future__ import annotations

import argparse
import array
import json
import math
from pathlib import Path
import struct
import sys


ROOT = Path(__file__).resolve().parents[4]
CELL_M = 1.0
CORRIDOR_MARGIN_M = 3.0
OPEN_GROUND_EXTENSION_M = 5.0
WELD_M = 0.0001


class Triangle:
    __slots__ = ("verts", "center", "area", "normal", "smooth", "quality",
                 "route_distance", "route_lateral", "route_s", "floor_z",
                 "in_region", "reference")

    def __init__(self, verts, area, normal, smooth, quality):
        self.verts = verts
        self.center = tuple(sum(v[i] for v in verts) / 3.0 for i in range(3))
        self.area = area
        self.normal = normal
        self.smooth = smooth
        self.quality = quality
        self.route_distance = 0.0
        self.route_lateral = 0.0
        self.route_s = 0.0
        self.floor_z = 0.0
        self.in_region = False
        self.reference = False


def normalize(value):
    length = math.sqrt(sum(item * item for item in value))
    return tuple(item / length for item in value) if length > 1.0e-12 else (0.0, 0.0, 0.0)


def dot(a, b):
    return sum(a[i] * b[i] for i in range(3))


def angle_degrees(a, b):
    return math.degrees(math.acos(max(-1.0, min(1.0, dot(a, b)))))


def quantile(values, fraction):
    if not values:
        return 0.0
    ordered = sorted(values)
    return ordered[min(len(ordered) - 1, max(0, int(math.ceil(fraction * len(ordered))) - 1))]


def read_mesh(path: Path):
    triangles = []
    selected_tiles = 0
    lod_counts = {}
    with path.open("rb") as stream:
        prefix = stream.read(60)
        if len(prefix) != 60 or prefix[:8] != b"VFMSH001":
            raise ValueError(f"{path}: invalid VFMSH001 mesh header")
        version, voxel_um, tile_count = struct.unpack_from("<3I", prefix, 8)
        if version != 1 or not voxel_um or tile_count > 1_000_000:
            raise ValueError(f"{path}: unsupported mesh version or tile count")
        for _ in range(tile_count):
            header = stream.read(36)
            if len(header) != 36:
                raise ValueError(f"{path}: truncated tile header")
            key = struct.unpack_from("<4i", header, 0)
            surface_count, _surface_aux, vertex_count, index_count, mask_count = \
                struct.unpack_from("<5I", header, 16)
            lod_counts[key[3]] = lod_counts.get(key[3], 0) + 1
            position_bytes = stream.read(vertex_count * 12)
            normal_bytes = stream.read(vertex_count * 12)
            index_bytes = stream.read(index_count * 4)
            masks = stream.read(mask_count)
            if (len(position_bytes) != vertex_count * 12
                    or len(normal_bytes) != vertex_count * 12
                    or len(index_bytes) != index_count * 4
                    or len(masks) != mask_count):
                raise ValueError(f"{path}: truncated tile payload")
            if key[3] != 0:
                continue
            selected_tiles += 1
            positions = array.array("f")
            positions.frombytes(position_bytes)
            normals = array.array("f")
            normals.frombytes(normal_bytes)
            indices = array.array("I")
            indices.frombytes(index_bytes)
            if sys.byteorder != "little":
                positions.byteswap()
                normals.byteswap()
                indices.byteswap()
            for offset in range(0, index_count - 2, 3):
                if masks and masks[offset // 3]:
                    continue
                ids = indices[offset:offset + 3]
                if max(ids) >= surface_count:
                    continue
                verts = tuple(tuple(float(positions[index * 3 + axis]) for axis in range(3))
                              for index in ids)
                edge_a = tuple(verts[1][axis] - verts[0][axis] for axis in range(3))
                edge_b = tuple(verts[2][axis] - verts[0][axis] for axis in range(3))
                cross = (edge_a[1] * edge_b[2] - edge_a[2] * edge_b[1],
                         edge_a[2] * edge_b[0] - edge_a[0] * edge_b[2],
                         edge_a[0] * edge_b[1] - edge_a[1] * edge_b[0])
                cross_length = math.sqrt(dot(cross, cross))
                if cross_length <= 1.0e-12:
                    continue
                face_normal = tuple(value / cross_length for value in cross)
                area = cross_length * 0.5
                smooth_normal = normalize(tuple(
                    sum(float(normals[index * 3 + axis]) for index in ids) for axis in range(3)))
                if dot(face_normal, smooth_normal) < 0.0:
                    face_normal = tuple(-value for value in face_normal)
                edge_squares = 0.0
                for a, b in ((verts[0], verts[1]), (verts[1], verts[2]),
                             (verts[2], verts[0])):
                    edge_squares += sum((a[i] - b[i]) ** 2 for i in range(3))
                quality = 4.0 * math.sqrt(3.0) * area / max(edge_squares, 1.0e-15)
                triangles.append(Triangle(verts, area, face_normal, smooth_normal, quality))
        if stream.read(1):
            raise ValueError(f"{path}: trailing bytes after mesh tiles")
    if not triangles:
        raise ValueError(f"{path}: no LOD0 surface triangles")
    return triangles, voxel_um, selected_tiles, lod_counts


def bezier(points, t):
    one = 1.0 - t
    weights = (one ** 3, 3.0 * one * one * t, 3.0 * one * t * t, t ** 3)
    return tuple(sum(weights[i] * points[i][axis] for i in range(4)) for axis in range(3))


def make_route(layout):
    edges = layout.get("cavern_edges", [])
    rooms = layout.get("cavern_rooms", [])
    if len(edges) != 1 or len(rooms) != 1:
        raise ValueError("cavern audit expects one resolved room and one edge in its layout")
    edge = edges[0]
    points = [[value / 1_000_000.0 for value in row]
              for row in edge["control_points_um"]]
    bore_radius = edge["bore_radius_um"] / 1_000_000.0
    bore_height = edge["bore_height_um"] / 1_000_000.0
    samples = [bezier(points, index / 160.0) for index in range(161)]
    horizontal_s = [0.0]
    for a, b in zip(samples, samples[1:]):
        horizontal_s.append(horizontal_s[-1] + math.hypot(b[0] - a[0], b[1] - a[1]))
    first_dx = samples[1][0] - samples[0][0]
    first_dy = samples[1][1] - samples[0][1]
    first_len = max(1.0e-9, math.hypot(first_dx, first_dy))
    extension = (samples[0][0] - first_dx / first_len * OPEN_GROUND_EXTENSION_M,
                 samples[0][1] - first_dy / first_len * OPEN_GROUND_EXTENSION_M,
                 samples[0][2] - (samples[1][2] - samples[0][2]) / first_len
                    * OPEN_GROUND_EXTENSION_M)
    route_points = [extension, *samples]
    route_s = [-OPEN_GROUND_EXTENSION_M, *horizontal_s]
    segments = []
    for i, (a, b) in enumerate(zip(route_points, route_points[1:])):
        segments.append((a, b, route_s[i], route_s[i + 1]))
    room = rooms[0]
    room_center = tuple(value / 1_000_000.0 for value in room["center_um"])
    room_half = tuple(value / 1_000_000.0 for value in room["half_extents_um"])
    room_floor = room["floor_um"] / 1_000_000.0
    return {
        "segments": segments,
        "edge": edge,
        "bore_radius": bore_radius,
        "bore_height": bore_height,
        "room_center": room_center,
        "room_half": room_half,
        "room_floor": room_floor,
        "route_points": samples,
        "route_length": horizontal_s[-1],
    }


def route_nearest(point, route):
    best = None
    for a, b, s0, s1 in route["segments"]:
        dx, dy = b[0] - a[0], b[1] - a[1]
        denominator = dx * dx + dy * dy
        t = 0.0 if denominator <= 1.0e-12 else max(0.0, min(1.0,
            ((point[0] - a[0]) * dx + (point[1] - a[1]) * dy) / denominator))
        x, y = a[0] + dx * t, a[1] + dy * t
        distance = math.hypot(point[0] - x, point[1] - y)
        if best is None or distance < best[0]:
            segment_length = math.sqrt(max(denominator, 1.0e-12))
            lateral = (dx * (point[1] - y) - dy * (point[0] - x)) / segment_length
            best = (distance, s0 + (s1 - s0) * t,
                    a[2] + (b[2] - a[2]) * t, lateral)
    return best


def mark_region(triangles, route):
    room_c = route["room_center"]
    room_h = route["room_half"]
    radius = route["bore_radius"] + CORRIDOR_MARGIN_M
    for tri in triangles:
        x, y, z = tri.center
        distance, distance_s, floor_z, lateral = route_nearest(tri.center, route)
        tri.route_distance = distance
        tri.route_lateral = lateral
        tri.route_s = distance_s
        tri.floor_z = floor_z
        corridor = (-OPEN_GROUND_EXTENSION_M <= distance_s <= route["route_length"]
            and distance <= radius
            and floor_z - CORRIDOR_MARGIN_M <= z
            <= floor_z + route["bore_height"] + CORRIDOR_MARGIN_M)
        room_box = (abs(x - room_c[0]) <= room_h[0] + CORRIDOR_MARGIN_M
            and abs(y - room_c[1]) <= room_h[1] + CORRIDOR_MARGIN_M
            and route["room_floor"] - CORRIDOR_MARGIN_M <= z
            <= route["room_floor"] + 2.0 * room_h[2] + CORRIDOR_MARGIN_M)
        tri.in_region = corridor or room_box
        outside_bore = (distance > route["bore_radius"] + 0.45
            or abs(x - room_c[0]) > room_h[0] + 0.45
            or abs(y - room_c[1]) > room_h[1] + 0.45)
        # Keep the reference on native, ground-facing terrain. The outer mesh
        # also contains steep cavern walls and room edges, which are features.
        tri.reference = tri.in_region and outside_bore and tri.smooth[2] >= 0.45


def vertex_key(point):
    return tuple(int(round(value / WELD_M)) for value in point)


def make_edges(triangles, indices):
    edge_map = {}
    for index in indices:
        tri = triangles[index]
        points = tri.verts
        for a, b in ((0, 1), (1, 2), (2, 0)):
            ka, kb = vertex_key(points[a]), vertex_key(points[b])
            if ka == kb:
                continue
            key = (ka, kb) if ka < kb else (kb, ka)
            delta = tuple(points[a][axis] - points[b][axis] for axis in range(3))
            length = math.sqrt(dot(delta, delta))
            midpoint = tuple((points[a][axis] + points[b][axis]) * 0.5 for axis in range(3))
            edge_map.setdefault(key, []).append((index, length, midpoint))
    pairs = []
    for incident in edge_map.values():
        unique = {}
        for index, length, midpoint in incident:
            unique[index] = (length, midpoint)
        incident = list(unique.items())
        if len(incident) < 2:
            continue
        if len(incident) > 4:
            incident = incident[:4]
        for left in range(len(incident)):
            for right in range(left + 1, len(incident)):
                i, (length, midpoint) = incident[left]
                j, _ = incident[right]
                angle = max(angle_degrees(triangles[i].normal, triangles[j].normal),
                            0.75 * angle_degrees(triangles[i].smooth, triangles[j].smooth))
                pairs.append((i, j, angle, length, midpoint))
    return pairs


def cell_of(point):
    return tuple(math.floor(value / CELL_M) for value in point)


def score_surface(triangles, considered, pairs, edge_cutoff, rough_cutoff,
                  ref_cell_density):
    considered_set = set(considered)
    variation = {index: 0.0 for index in considered}
    max_jump = {index: 0.0 for index in considered}
    severe_count = {index: 0.0 for index in considered}
    for i, j, angle, length, _midpoint in pairs:
        if i not in considered_set or j not in considered_set:
            continue
        weighted = angle * length * 0.5
        variation[i] += weighted
        variation[j] += weighted
        max_jump[i] = max(max_jump[i], angle)
        max_jump[j] = max(max_jump[j], angle)
        if angle >= edge_cutoff:
            severe_count[i] += 0.5
            severe_count[j] += 0.5
    cells = {}
    for index in considered:
        tri = triangles[index]
        key = cell_of(tri.center)
        cell = cells.setdefault(key, {"area": 0.0, "variation": 0.0, "severe": 0.0,
            "sliver_area": 0.0, "slivers": 0, "triangles": 0, "max_jump": 0.0,
            "reference_area": 0.0, "reference_variation": 0.0,
            "reference_severe": 0.0, "reference_sliver_area": 0.0,
            "reference_slivers": 0})
        cell["area"] += tri.area
        cell["variation"] += variation[index]
        cell["severe"] += severe_count[index]
        cell["triangles"] += 1
        cell["max_jump"] = max(cell["max_jump"], max_jump[index])
        sliver = tri.quality < 0.07 and tri.area < 0.12
        if sliver:
            cell["sliver_area"] += tri.area
            cell["slivers"] += 1
        if tri.reference:
            cell["reference_area"] += tri.area
            cell["reference_variation"] += variation[index]
            cell["reference_severe"] += severe_count[index]
            if sliver:
                cell["reference_sliver_area"] += tri.area
                cell["reference_slivers"] += 1
    ref_rough = []
    ref_density = []
    ref_sliver_fraction = []
    for cell in cells.values():
        if cell["reference_area"] >= max(0.15, cell["area"] * 0.65):
            ref_rough.append(cell["reference_variation"] / max(cell["reference_area"], 1.0e-9))
            ref_density.append(cell["reference_severe"] / max(cell["reference_area"], 1.0e-9))
            ref_sliver_fraction.append(cell["reference_sliver_area"]
                / max(cell["reference_area"], 1.0e-9))
    reference_median = quantile(ref_rough, 0.5)
    reference_mad = quantile([abs(value - reference_median) for value in ref_rough], 0.5)
    cell_rough_limit = max(10.0, quantile(ref_rough, 0.995) * 1.5,
                           reference_median + 8.0 * reference_mad)
    edge_density_limit = max(1.5, quantile(ref_density, 0.995) * 1.7,
                             ref_cell_density * 1.7)
    sliver_limit = max(0.22, quantile(ref_sliver_fraction, 0.995) + 0.18)
    hot = set()
    for key, cell in cells.items():
        area = max(cell["area"], 1.0e-9)
        roughness = cell["variation"] / area
        edge_density = cell["severe"] / area
        sliver_fraction = cell["sliver_area"] / area
        rough_cluster = roughness > cell_rough_limit and cell["severe"] >= 2.0
        edge_cluster = cell["severe"] >= 4.0 and edge_density > edge_density_limit
        sliver_cluster = cell["slivers"] >= 5 and sliver_fraction > sliver_limit
        if rough_cluster or edge_cluster or sliver_cluster:
            hot.add(key)
        cell["roughness_deg_per_m"] = roughness
        cell["severe_edge_density_m2"] = edge_density
        cell["sliver_area_fraction"] = sliver_fraction
        cell["reasons"] = [name for name, active in (("normal-variation", rough_cluster),
            ("edge-density", edge_cluster), ("slivers", sliver_cluster)) if active]
    clusters = []
    flagged_cells = set()
    remaining = set(hot)
    while remaining:
        stack = [remaining.pop()]
        members = []
        while stack:
            current = stack.pop()
            members.append(current)
            x, y, z = current
            neighbors = [(x + dx, y + dy, z + dz)
                for dx in (-1, 0, 1) for dy in (-1, 0, 1) for dz in (-1, 0, 1)
                if dx or dy or dz]
            for neighbor in neighbors:
                if neighbor in remaining:
                    remaining.remove(neighbor)
                    stack.append(neighbor)
        cells_in_cluster = [cells[key] for key in members]
        cluster_area = sum(cell["area"] for cell in cells_in_cluster)
        cluster_severe = sum(cell["severe"] for cell in cells_in_cluster)
        cluster_slivers = sum(cell["slivers"] for cell in cells_in_cluster)
        accepted = len(members) >= 2 and cluster_area >= 0.25 and \
            (cluster_severe >= 8.0 or cluster_slivers >= 12)
        if accepted:
            flagged_cells.update(members)
            clusters.append({"cells": len(members), "area_m2": cluster_area,
                "severe_edges": cluster_severe, "sliver_triangles": cluster_slivers,
                "bounds_m": [[min(key[i] for key in members) * CELL_M,
                              (max(key[i] for key in members) + 1) * CELL_M]
                             for i in range(3)],
                "reasons": sorted({reason for cell in cells_in_cluster
                                   for reason in cell["reasons"]})})
    clusters.sort(key=lambda cluster: (cluster["severe_edges"], cluster["area_m2"]), reverse=True)
    flagged_triangles = [index for index in considered
                         if cell_of(triangles[index].center) in flagged_cells]
    return {"cells": cells, "hot_cells": hot, "flagged_cells": flagged_cells,
            "clusters": clusters, "hot_triangles": flagged_triangles,
            "thresholds": {"edge_angle_deg": edge_cutoff,
                "roughness_deg_per_m": cell_rough_limit,
                "severe_edge_density_m2": edge_density_limit,
                "sliver_area_fraction": sliver_limit,
                "reference_cell_roughness_p995_deg_per_m": quantile(ref_rough, 0.995),
                "reference_cell_roughness_median_deg_per_m": reference_median}}


def floor_profile_audit(triangles, considered, region_indices, route):
    """Measure the whole floor patch as area-weighted bins, not selected rails."""
    if route is None:
        return {"present": False, "detected": False, "triangles": 0, "bins": 0,
                "lateral_bands": 0, "ripple_amplitude_m": 0.0,
                "max_adjacent_step_m": 0.0}
    along_step = 0.25
    lateral_step = 0.5
    bins = {}
    floor_triangles = 0
    wall_triangles = 0
    wall_bins = set()
    for index in region_indices:
        tri = triangles[index]
        height = tri.center[2] - tri.floor_z
        if (0.0 <= tri.route_s <= 10.0
                and route["bore_radius"] - 0.9 <= tri.route_distance
                    <= route["bore_radius"] + 1.0
                and 0.45 <= height <= route["bore_height"] + 0.5
                and abs(tri.smooth[2]) <= 0.45):
            wall_triangles += 1
            wall_bins.add(int(tri.route_s))
    for index in considered:
        tri = triangles[index]
        if (tri.route_s < 0.0 or tri.route_s > min(10.0, route["route_length"])
                or tri.route_distance > min(1.9, route["bore_radius"] - 0.1)
                or abs(tri.center[2] - tri.floor_z) > 0.60
                or tri.smooth[2] < 0.55):
            continue
        floor_triangles += 1
        along_bin = int(math.floor(tri.route_s / along_step))
        lateral_bin = int(math.floor((tri.route_lateral + 2.0) / lateral_step))
        key = (lateral_bin, along_bin)
        row = bins.setdefault(key, [0.0, 0.0, 0])
        residual = tri.center[2] - tri.floor_z
        row[0] += residual * tri.area
        row[1] += tri.area
        row[2] += 1
    by_lateral = {}
    for (lateral_bin, along_bin), (weighted_residual, area, count) in bins.items():
        if area < 0.002:
            continue
        by_lateral.setdefault(lateral_bin, {})[along_bin] = weighted_residual / area
    amplitudes = []
    steps = []
    valid_bins = 0
    for values in by_lateral.values():
        ordered = sorted(values.items())
        run = []
        runs = []
        for item in ordered:
            if run and item[0] != run[-1][0] + 1:
                runs.append(run)
                run = []
            run.append(item)
        if run:
            runs.append(run)
        for segment in runs:
            if len(segment) < 8:
                continue
            valid_bins += len(segment)
            raw = [value for _, value in segment]
            highpass = []
            for position, value in enumerate(raw):
                nearby = raw[max(0, position - 4):min(len(raw), position + 5)]
                highpass.append(value - sum(nearby) / len(nearby))
            amplitudes.append(max(highpass) - min(highpass))
            steps.extend(abs(raw[i] - raw[i - 1]) for i in range(1, len(raw)))
    max_ripple = max(amplitudes, default=0.0) * 0.5
    max_step = max(steps, default=0.0)
    present = (floor_triangles >= 24 and valid_bins >= 8
        and wall_triangles >= 24 and len(wall_bins) >= 3)
    detected = present and (max_ripple > 0.015 or max_step > 0.050)
    return {"present": present, "detected": detected, "triangles": floor_triangles,
        "bins": valid_bins, "lateral_bands": len(by_lateral),
        "wall_triangles": wall_triangles, "wall_bins": len(wall_bins),
        "ripple_amplitude_m": max_ripple, "max_adjacent_step_m": max_step,
        "thresholds": {"ripple_amplitude_m": 0.015, "adjacent_step_m": 0.050}}


def wall_band_measure(pairs, route, edge_cutoff):
    """Locate high dihedrals around the bore-foot and the L47 0.25 m band edge."""
    if route is None:
        return {}
    bins = {}
    radius = route["bore_radius"]
    for _i, _j, angle, _length, midpoint in pairs:
        distance, route_s, floor_z, _lateral = route_nearest(midpoint, route)
        height = midpoint[2] - floor_z
        if (0.0 <= route_s <= 10.0 and abs(distance - radius) <= 0.75
                and -0.25 <= height <= 0.875):
            index = int(math.floor((height + 0.25) / 0.125))
            entry = bins.setdefault(index, [])
            entry.append(angle)
    result = {}
    for index, angles in sorted(bins.items()):
        center_height = (index + 0.5) * 0.125 - 0.25
        result[f"{center_height:.3f}"] = {
            "edges": len(angles),
            "over_reference_cutoff": sum(value >= edge_cutoff for value in angles),
            "p95_dihedral_deg": quantile(angles, 0.95),
            "max_dihedral_deg": max(angles, default=0.0),
        }
    return result


def radial_wall_measure(pairs, route, edge_cutoff):
    """Measure dihedrals across the hard bore-radius gate at the entrance foot."""
    if route is None:
        return {}
    bins = {}
    radius = route["bore_radius"]
    for _i, _j, angle, _length, midpoint in pairs:
        distance, route_s, floor_z, _lateral = route_nearest(midpoint, route)
        height = midpoint[2] - floor_z
        if (-5.0 <= route_s <= 3.0 and -1.0 <= height <= 0.75
                and abs(distance - radius) <= 1.0):
            index = int(math.floor((distance - (radius - 1.0)) / 0.125))
            bins.setdefault(index, []).append(angle)
    return {f"{radius - 1.0 + (index + 0.5) * 0.125:.3f}": {
        "edges": len(angles),
        "over_reference_cutoff": sum(value >= edge_cutoff for value in angles),
        "p95_dihedral_deg": quantile(angles, 0.95),
        "max_dihedral_deg": max(angles, default=0.0),
    } for index, angles in sorted(bins.items())}


def svg_for(path, triangles, route, hot_triangles, label):
    selected = [triangles[index] for index in hot_triangles]
    all_region = [tri for tri in triangles if tri.in_region]
    points = [point for tri in all_region for point in tri.verts]
    if not points:
        return
    xmin, xmax = min(p[0] for p in points), max(p[0] for p in points)
    ymin, ymax = min(p[1] for p in points), max(p[1] for p in points)
    width, height, margin = 1200.0, 900.0, 35.0
    sx = (width - 2 * margin) / max(1.0e-6, xmax - xmin)
    sy = (height - 2 * margin) / max(1.0e-6, ymax - ymin)
    scale = min(sx, sy)
    draw_w, draw_h = (xmax - xmin) * scale, (ymax - ymin) * scale
    ox, oy = (width - draw_w) * 0.5, (height - draw_h) * 0.5
    def project(point):
        return ox + (point[0] - xmin) * scale, height - oy - (point[1] - ymin) * scale
    lines = [f'<svg xmlns="http://www.w3.org/2000/svg" width="{int(width)}" height="{int(height)}" viewBox="0 0 {int(width)} {int(height)}">',
        '<rect width="100%" height="100%" fill="#fbfbfa"/>',
        f'<text x="28" y="28" font-family="sans-serif" font-size="19" fill="#222">{label}: whole-surface defect clusters (XY plan)</text>',
        '<g fill="#c8ced2" fill-opacity="0.24" stroke="#8b959b" stroke-opacity="0.12" stroke-width="0.3">']
    for tri in all_region:
        coords = [project(p) for p in tri.verts]
        lines.append('<polygon points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in coords)
                     + '"/>')
    lines.append('</g><g fill="#d7191c" fill-opacity="0.78" stroke="#7f0000" stroke-width="0.4">')
    for tri in selected:
        coords = [project(p) for p in tri.verts]
        lines.append('<polygon points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in coords)
                     + '"/>')
    route_coords = [project(p) for p in route["route_points"]]
    lines.append('</g><polyline points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in route_coords)
                 + '" fill="none" stroke="#1769aa" stroke-width="2"/>')
    lines.append(f'<text x="28" y="{int(height - 24)}" font-family="sans-serif" font-size="15" fill="#333">Blue: intended route · Red: flagged triangle cells · Gray: remaining ROI triangles · flagged triangles: {len(selected)}</text>')
    lines.append('</svg>')
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")


def audit_case(label, mesh_path, layout_path, expect, focus, svg_dir):
    triangles, voxel_um, tile_count, lod_counts = read_mesh(mesh_path)
    route = None
    if layout_path != Path("-"):
        layout = json.loads(layout_path.read_text(encoding="utf-8"))
        route = make_route(layout)
        mark_region(triangles, route)
        region_indices = [i for i, tri in enumerate(triangles) if tri.in_region]
    else:
        region_indices = list(range(len(triangles)))
        for tri in triangles:
            tri.in_region = True
            tri.reference = True
    if len(region_indices) < 24:
        raise ValueError(f"{label}: only {len(region_indices)} triangles intersect the ROI")
    reference_indices = [i for i in region_indices if triangles[i].reference]
    if len(reference_indices) < 24:
        reference_indices = region_indices
        for index in region_indices:
            triangles[index].reference = True
    edge_pairs = make_edges(triangles, region_indices)
    reference_edges = [edge[2] for edge in edge_pairs
        if triangles[edge[0]].reference and triangles[edge[1]].reference]
    edge_cutoff = max(35.0, quantile(reference_edges, 0.995))
    reference_set = set(reference_indices)
    ref_roughness = {index: 0.0 for index in reference_indices}
    for i, j, angle, length, _midpoint in edge_pairs:
        if i in reference_set and j in reference_set:
            ref_roughness[i] += angle * length * 0.5
            ref_roughness[j] += angle * length * 0.5
    ref_values = [ref_roughness[index] / max(triangles[index].area, 1.0e-9)
                  for index in reference_indices]
    rough_cutoff = 10.0
    ref_density = 0.0
    ref_area = sum(triangles[index].area for index in reference_indices)
    ref_severe = 0.0
    for i, j, angle, _length, _midpoint in edge_pairs:
        if i in reference_set and j in reference_set and angle >= edge_cutoff:
            ref_severe += 1.0
    ref_density = ref_severe / max(ref_area, 1.0e-9)
    if focus == "floor":
        considered = [i for i in region_indices if route is not None
            and triangles[i].route_s >= 1.5
            and triangles[i].route_s <= 10.0
            and triangles[i].route_distance <= 1.25
            and abs(triangles[i].center[2] - triangles[i].floor_z) <= 0.24
            and triangles[i].smooth[2] >= 0.55]
    else:
        considered = region_indices
    if len(considered) < 24:
        raise ValueError(f"{label}: focus={focus} has only {len(considered)} triangles")
    scores = score_surface(triangles, considered, edge_pairs, edge_cutoff,
                           rough_cutoff, ref_density)
    floor_profile = floor_profile_audit(triangles, considered, region_indices, route)
    wall_band = wall_band_measure(edge_pairs, route, edge_cutoff)
    radial_wall = radial_wall_measure(edge_pairs, route, edge_cutoff)
    found = bool(scores["clusters"]) or floor_profile["detected"]
    passed = expect == "measure" or (expect == "fail" and found) or (expect == "pass" and not found)
    if svg_dir and route is not None:
        svg_for(svg_dir / f"flagged-regions-{label}.svg", triangles, route,
                scores["hot_triangles"], label)
    cluster_summary = scores["clusters"][:8]
    result = {
        "label": label,
        "mesh": str(mesh_path),
        "layout": None if layout_path == Path("-") else str(layout_path),
        "focus": focus,
        "expected": expect,
        "observed": "fail" if found else "pass",
        "status": "PASS" if passed else "FAIL",
        "voxel_um": voxel_um,
        "lod0_tiles": tile_count,
        "lod_tile_counts": lod_counts,
        "triangles_loaded_lod0": len(triangles),
        "triangles_in_roi": len(region_indices),
        "triangles_audited": len(considered),
        "adjacent_pairs": len(edge_pairs),
        "reference_triangles": len(reference_indices),
        "reference_edge_p995_deg": quantile(reference_edges, 0.995),
        "reference_p995_normal_variation_deg_per_m": quantile(ref_values, 0.995),
        "thresholds": scores["thresholds"],
        "flagged_cells": len(scores["hot_cells"]),
        "flagged_triangles": len(scores["hot_triangles"]),
        "floor_profile": floor_profile,
        "wall_band_height_histogram": wall_band,
        "radial_wall_histogram": radial_wall,
        "clusters": cluster_summary,
    }
    print(f"[WholeSurfaceAudit] label={label} focus={focus} lod0_tiles={tile_count} "
        f"triangles_roi={len(region_indices)} audited={len(considered)} "
        f"adjacent_pairs={len(edge_pairs)} reference={len(reference_indices)} "
        f"edge_cutoff_deg={edge_cutoff:.2f} "
        f"roughness_cutoff_deg_m={scores['thresholds']['roughness_deg_per_m']:.2f} "
        f"hot_cells={len(scores['hot_cells'])} flagged_triangles={len(scores['hot_triangles'])} "
        f"clusters={len(scores['clusters'])} observed={result['observed'].upper()} "
        f"expected={expect.upper()} status={result['status']}")
    print(f"[WholeSurfaceFloorProfile] label={label} triangles={floor_profile['triangles']} "
        f"bins={floor_profile['bins']} lateral_bands={floor_profile['lateral_bands']} "
        f"wall_triangles={floor_profile.get('wall_triangles', 0)} "
        f"wall_bins={floor_profile.get('wall_bins', 0)} "
        f"ripple_amplitude_m={floor_profile['ripple_amplitude_m']:.4f} "
        f"max_adjacent_step_m={floor_profile['max_adjacent_step_m']:.4f} "
        f"observed={'FAIL' if floor_profile['detected'] else ('PASS' if floor_profile['present'] else 'NOT_PRESENT')}")
    if wall_band:
        print(f"[WholeSurfaceWallBandMeasure] label={label} voxel_m={voxel_um / 1_000_000.0:.3f} "
            f"recarve_edge_height_m=0.250 radial_band_m={route['bore_radius']:.3f} "
            f"height_histogram={json.dumps(wall_band, sort_keys=True, separators=(',', ':'))}")
        print(f"[WholeSurfaceRadialWallMeasure] label={label} route_s_m=-5..3 "
            f"height_above_floor_m=-1..0.75 bore_radius_m={route['bore_radius']:.3f} "
            f"radial_histogram={json.dumps(radial_wall, sort_keys=True, separators=(',', ':'))}")
    for index, cluster in enumerate(cluster_summary):
        print(f"[WholeSurfaceCluster] label={label} index={index} cells={cluster['cells']} "
            f"area_m2={cluster['area_m2']:.3f} severe_edges={cluster['severe_edges']:.1f} "
            f"slivers={cluster['sliver_triangles']} bounds_m={cluster['bounds_m']} "
            f"reasons={','.join(cluster['reasons'])}")
    hot_cells = scores["cells"]
    for key in sorted(scores["hot_cells"], key=lambda item: (
            hot_cells[item]["severe"], hot_cells[item]["variation"]), reverse=True)[:12]:
        cell = hot_cells[key]
        local = route_nearest(tuple((axis + 0.5) * CELL_M for axis in key), route) if route else None
        print(f"[WholeSurfaceHotCell] label={label} cell_m=({key[0]},{key[1]},{key[2]}) "
            f"route_s_m={local[1] if local else 0.0:.2f} radial_m={local[0] if local else 0.0:.2f} "
            f"height_above_floor_m={(key[2] + 0.5 - local[2]) if local else 0.0:.2f} "
            f"area_m2={cell['area']:.3f} roughness_deg_m={cell['roughness_deg_per_m']:.2f} "
            f"severe_edges={cell['severe']:.1f} edge_density_m2={cell['severe_edge_density_m2']:.2f} "
            f"slivers={cell['slivers']} max_jump_deg={cell['max_jump']:.2f} "
            f"reasons={','.join(cell['reasons'])}")
    return result


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--case", nargs=5, action="append", required=True,
        metavar=("LABEL", "MESH", "LAYOUT_OR_DASH", "EXPECT", "FOCUS"))
    parser.add_argument("--json", type=Path, required=True)
    parser.add_argument("--svg-dir", type=Path)
    args = parser.parse_args()
    results = []
    for label, mesh, layout, expect, focus in args.case:
        if expect not in ("pass", "fail", "measure") or focus not in ("whole", "floor"):
            parser.error("EXPECT must be pass/fail/measure and FOCUS must be whole/floor")
        results.append(audit_case(label, Path(mesh), Path(layout), expect, focus, args.svg_dir))
    args.json.parent.mkdir(parents=True, exist_ok=True)
    args.json.write_text(json.dumps({"audit": "whole-surface-l48-v1", "cases": results},
                                    indent=2) + "\n", encoding="utf-8")
    failures = [result for result in results if result["status"] != "PASS"]
    print(f"[WholeSurfaceAuditSummary] cases={len(results)} failed={len(failures)} "
          f"json={args.json} status={'FAIL' if failures else 'PASS'}")
    return 1 if failures else 0


if __name__ == "__main__":
    raise SystemExit(main())