#!/usr/bin/env python3
"""Fail on contour terraces, stepped profiles, or insufficient Gradient amplitude."""
import argparse
import array
import math
import resource
import struct
import sys
from collections import defaultdict


def read_scene(path, wanted_lod, height_min, height_max, bounds):
    triangles = []
    selected_tiles = 0
    x_min, x_max, y_min, y_max = bounds
    with open(path, "rb") as stream:
        if stream.read(8) != b"VFMSH001":
            raise SystemExit(f"{path}: invalid mesh scene header")
        version, voxel_um, tile_count = struct.unpack("<3I", stream.read(12))
        stream.read(40)
        if version != 1 or tile_count > 100000:
            raise SystemExit(f"{path}: unsupported mesh scene")
        for _ in range(tile_count):
            key = struct.unpack("<4i", stream.read(16))
            surface_count, _, vertex_count, index_count, mask_count = struct.unpack(
                "<5I", stream.read(20))
            tile_extent = 32.0 * (1 << key[3]) * voxel_um / 1_000_000.0
            tile_x0, tile_y0 = key[0] * tile_extent, key[1] * tile_extent
            tile_z0 = key[2] * tile_extent
            if (wanted_lod is not None and key[3] != wanted_lod) or \
                    tile_x0 > x_max or tile_x0 + tile_extent < x_min or \
                    tile_y0 > y_max or tile_y0 + tile_extent < y_min or \
                    tile_z0 > height_max or tile_z0 + tile_extent < height_min:
                stream.seek(vertex_count * 24 + index_count * 4 + mask_count, 1)
                continue
            selected_tiles += 1
            positions = array.array("f")
            positions.frombytes(stream.read(vertex_count * 12))
            normals = array.array("f")
            normals.frombytes(stream.read(vertex_count * 12))
            indices = array.array("I")
            indices.frombytes(stream.read(index_count * 4))
            masks = stream.read(mask_count)
            if sys.byteorder != "little":
                positions.byteswap()
                normals.byteswap()
                indices.byteswap()
            for offset in range(0, index_count, 3):
                if masks and masks[offset // 3]:
                    continue
                ids = indices[offset:offset + 3]
                if max(ids) >= surface_count:
                    continue
                verts = []
                for index in ids:
                    p, n = index * 3, index * 3
                    verts.append(((positions[p], positions[p + 1], positions[p + 2]),
                                  (normals[n], normals[n + 1], normals[n + 2])))
                average_nz = sum(vertex[1][2] for vertex in verts) / 3.0
                z_values = [vertex[0][2] for vertex in verts]
                # Select the requested surface branch by its face height. An
                # overlap or centroid test admits a giant coarse triangle that
                # merely clips this window and can belong to the neighboring face.
                mean_z = sum(z_values) / 3.0
                vertex_margin = voxel_um / 1_000_000.0 * 0.5
                if (mean_z < height_min or mean_z > height_max or
                        min(z_values) < height_min - vertex_margin or
                        max(z_values) > height_max + vertex_margin):
                    continue
                if (max(vertex[0][0] for vertex in verts) < x_min or
                        min(vertex[0][0] for vertex in verts) > x_max or
                        max(vertex[0][1] for vertex in verts) < y_min or
                        min(vertex[0][1] for vertex in verts) > y_max):
                    continue
                triangles.append((key[3], average_nz, tuple(verts)))
    return voxel_um, triangles, selected_tiles


def build_buckets(triangles, upward, bounds, cell_size=2.0):
    buckets = defaultdict(list)
    x_min, x_max, y_min, y_max = bounds
    for triangle in triangles:
        _, average_nz, verts = triangle
        if (upward and average_nz < 0.35) or (not upward and average_nz > -0.35):
            continue
        xs = [v[0][0] for v in verts]
        ys = [v[0][1] for v in verts]
        x0 = max(math.floor(min(xs) / cell_size), math.floor(x_min / cell_size))
        x1 = min(math.floor(max(xs) / cell_size), math.floor(x_max / cell_size))
        y0 = max(math.floor(min(ys) / cell_size), math.floor(y_min / cell_size))
        y1 = min(math.floor(max(ys) / cell_size), math.floor(y_max / cell_size))
        for iy in range(y0, y1 + 1):
            for ix in range(x0, x1 + 1):
                buckets[(ix, iy)].append(verts)
    return buckets


def sample_surface(buckets, x, y, upward, cell_size, preferred_height=None):
    found = None
    for verts in buckets.get((math.floor(x / cell_size), math.floor(y / cell_size)), ()):
        (a, na), (b, nb), (c, nc) = verts
        denominator = (b[1] - c[1]) * (a[0] - c[0]) + (c[0] - b[0]) * (a[1] - c[1])
        if abs(denominator) < 1.0e-12:
            continue
        u = ((b[1] - c[1]) * (x - c[0]) + (c[0] - b[0]) * (y - c[1])) / denominator
        v = ((c[1] - a[1]) * (x - c[0]) + (a[0] - c[0]) * (y - c[1])) / denominator
        w = 1.0 - u - v
        if min(u, v, w) < -1.0e-6:
            continue
        z = u * a[2] + v * b[2] + w * c[2]
        normal = tuple(u * na[axis] + v * nb[axis] + w * nc[axis] for axis in range(3))
        length = math.sqrt(sum(value * value for value in normal))
        if length < 1.0e-9:
            continue
        normal = tuple(value / length for value in normal)
        if (normal[2] > 0) != upward:
            continue
        if found is None or (abs(z - preferred_height) < abs(found[0] - preferred_height)
                             if preferred_height is not None else z > found[0]):
            found = (z, normal)
    return found


def make_lines(cx, cy, radius, step):
    lines = []
    for index in range(8):
        angle = index * math.tau / 8.0
        count = round(radius / step)
        points = [(cx + i * step * math.cos(angle), cy + i * step * math.sin(angle))
                  for i in range(count + 1)]
        lines.append((f"radial-{index}", points))
    count = round(2.0 * radius / step)
    t_values = [-radius + i * step for i in range(count + 1)]
    for name, points in (
        ("straight-x", [(cx + t, cy) for t in t_values]),
        ("straight-y", [(cx, cy + t) for t in t_values]),
        ("straight-diagonal", [(cx + t / math.sqrt(2.0), cy + t / math.sqrt(2.0))
                               for t in t_values]),
    ):
        lines.append((name, points))
    return lines


def audit_profile(points, buckets, upward, plateau_min_m, plateau_epsilon_m,
                  step_threshold_m, cell_size, preferred_height=None):
    samples = [sample_surface(buckets, x, y, upward, cell_size, preferred_height)
               for x, y in points]
    valid = [(index, value) for index, value in enumerate(samples) if value is not None]
    if len(valid) < max(8, len(points) * 0.75):
        return {"coverage": len(valid) / len(points), "plateaus": [], "max_angle": math.inf,
                "height_range": 0.0, "valid": False}
    heights = [value[0] for _, value in valid]
    plateau_runs = []
    run_start = valid[0][0]
    run_start_height = valid[0][1][0]
    sample_spacing = math.hypot(points[1][0] - points[0][0],
                                points[1][1] - points[0][1])
    for cursor in range(1, len(valid)):
        prior_index, prior = valid[cursor - 1]
        index, current = valid[cursor]
        if index != prior_index + 1:
            run_start = index
            run_start_height = current[0]
            continue
        if abs(current[0] - run_start_height) <= plateau_epsilon_m:
            continue
        run_length = (prior_index - run_start + 1) * sample_spacing
        if run_length >= plateau_min_m and abs(current[0] - prior[0]) >= step_threshold_m:
            plateau_runs.append((run_length, abs(current[0] - prior[0]), run_start,
                                 prior_index, index, prior[0], current[0]))
        run_start = index
        run_start_height = current[0]
    max_angle = 0.0
    for cursor in range(1, len(valid)):
        prior_index, prior = valid[cursor - 1]
        index, current = valid[cursor]
        if index != prior_index + 1:
            continue
        dot = max(-1.0, min(1.0, sum(a * b for a, b in zip(prior[1], current[1]))))
        max_angle = max(max_angle, math.degrees(math.acos(dot)))
    return {"coverage": len(valid) / len(points), "plateaus": plateau_runs,
            "max_angle": max_angle, "height_range": max(heights) - min(heights),
            "height_min": min(heights), "height_max": max(heights),
            "valid": True}


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("scene")
    parser.add_argument("--surface", choices=("floor", "ceiling"), default="floor")
    parser.add_argument("--center", nargs=2, type=float, default=(64.0, -32.0))
    parser.add_argument("--radius", type=float, default=16.0)
    parser.add_argument("--lod", required=True,
                        help="audit one LOD (0..8) or the streamed composite (`all`)")
    parser.add_argument("--step", type=float, default=0.125)
    parser.add_argument("--height-min-m", type=float, default=-130.0)
    parser.add_argument("--height-max-m", type=float, default=-110.0)
    parser.add_argument("--plateau-min-m", type=float, default=0.5)
    parser.add_argument("--plateau-epsilon-m", type=float, default=0.0005)
    parser.add_argument("--step-threshold-m", type=float, default=0.1)
    parser.add_argument("--normal-jump-max-deg", type=float, default=20.0)
    parser.add_argument("--amplitude-min-m", type=float, default=3.0)
    args = parser.parse_args()
    try:
        wanted_lod = None if args.lod == "all" else int(args.lod)
    except ValueError:
        parser.error("--lod must be 0..8 or all")
    if wanted_lod is not None and not 0 <= wanted_lod <= 8:
        parser.error("--lod must be 0..8 or all")
    bounds = (args.center[0] - args.radius, args.center[0] + args.radius,
              args.center[1] - args.radius, args.center[1] + args.radius)
    voxel_um, triangles, selected_tiles = read_scene(args.scene, wanted_lod,
        args.height_min_m, args.height_max_m, bounds)
    upward = args.surface == "floor"
    max_lod = max((triangle[0] for triangle in triangles), default=0)
    cell_size = max(2.0, 16.0 * voxel_um / 1_000_000.0 * (1 << max_lod))
    buckets = build_buckets(triangles, upward, bounds, cell_size)
    results = []
    # The requested face is the outer boundary of its authored band: highest Z
    # for a floor, lowest Z for a ceiling. Selecting the band midpoint can jump
    # between two overlapping downward/upward sheets at a cave edge.
    preferred_height = args.height_max_m if upward else args.height_min_m
    for name, points in make_lines(*args.center, args.radius, args.step):
        results.append((name, points, audit_profile(points, buckets, upward,
                         args.plateau_min_m, args.plateau_epsilon_m,
                         args.step_threshold_m, cell_size, preferred_height)))
    valid = [result for _, _, result in results if result["valid"]]
    if not valid:
        raise SystemExit("no sampled surface profiles; check scene, center, orientation and LOD")
    amplitude = max(result["height_max"] for result in valid) - min(
        result["height_min"] for result in valid)
    sampled_min = min(result["height_min"] for result in valid)
    sampled_max = max(result["height_max"] for result in valid)
    longest_plateau = max((run[0] for result in valid for run in result["plateaus"]),
                          default=0.0)
    largest_step = max((run[1] for result in valid for run in result["plateaus"]),
                       default=0.0)
    max_normal_jump = max(result["max_angle"] for result in valid)
    min_coverage = min(result["coverage"] for _, _, result in results)
    passed = (amplitude >= args.amplitude_min_m and longest_plateau < args.plateau_min_m
              and max_normal_jump <= args.normal_jump_max_deg and min_coverage >= 0.75)
    print(f"[RingTerraceAudit] scene={args.scene} surface={args.surface} lod={args.lod} "
          f"max_lod={max_lod} "
          f"height_window_m={args.height_min_m:.3f}..{args.height_max_m:.3f} "
          f"voxel_um={voxel_um} triangles={len(triangles)} lines={len(results)} "
          f"selected_tiles={selected_tiles} bucket_cell_m={cell_size:.3f} "
          f"coverage_min={min_coverage:.3f} sampled_height_m={sampled_min:.6f}.."
          f"{sampled_max:.6f} amplitude_range_m={amplitude:.6f} "
          f"longest_flat_run_m={longest_plateau:.6f} following_step_m={largest_step:.6f} "
          f"max_normal_delta_deg={max_normal_jump:.6f} status={'PASS' if passed else 'FAIL'}")
    for name, points, result in results:
        print(f"[RingTerraceProfile] line={name} coverage={result['coverage']:.3f} "
              f"range_m={result['height_range']:.6f} flat_step_runs={len(result['plateaus'])} "
              f"max_normal_delta_deg={result['max_angle']:.6f}")
        for run in result["plateaus"]:
            print(f"[RingTerraceViolation] line={name} flat_start_m=({points[run[2]][0]:.3f},"
                  f"{points[run[2]][1]:.3f}) flat_run_m={run[0]:.6f} "
                  f"step_m={run[1]:.6f} from_height_m={run[5]:.6f} "
                  f"to_height_m={run[6]:.6f}")
    print(f"[RingTerraceRSS] peak_rss_kib={resource.getrusage(resource.RUSAGE_SELF).ru_maxrss}")
    raise SystemExit(0 if passed else 1)


if __name__ == "__main__":
    main()