#!/usr/bin/env python3
"""Compare upward mesh heights on shared XY edges between adjacent LOD rings."""
import array
import collections
import math
import struct
import sys


def read_scene(path):
    tiles = []
    with open(path, "rb") as f:
        if f.read(8) != b"VFMSH001":
            raise SystemExit(f"{path}: invalid mesh header")
        version, voxel_um, count = struct.unpack("<3I", f.read(12))
        camera = struct.unpack("<3d", f.read(24))
        yaw, pitch = struct.unpack("<2d", f.read(16))
        voxel_m = voxel_um / 1_000_000.0
        for _ in range(count):
            key = struct.unpack("<4i", f.read(16))
            surface_vertices, skirt_count, vertex_count, index_count, mask_count = struct.unpack(
                "<5I", f.read(20))
            positions = array.array("f")
            positions.frombytes(f.read(vertex_count * 12))
            normals = array.array("f")
            normals.frombytes(f.read(vertex_count * 12))
            indices = array.array("I")
            indices.frombytes(f.read(index_count * 4))
            masks = f.read(mask_count)
            if sys.byteorder != "little":
                positions.byteswap(); normals.byteswap(); indices.byteswap()
            triangles = []
            for i in range(0, index_count, 3):
                if masks and masks[i // 3] != 0:
                    continue
                a, b, c = indices[i:i + 3]
                if max(a, b, c) >= surface_vertices:
                    continue
                nz = (normals[a * 3 + 2] + normals[b * 3 + 2] + normals[c * 3 + 2]) / 3.0
                if nz < 0.65:
                    continue
                tri = tuple((positions[j * 3], positions[j * 3 + 1], positions[j * 3 + 2])
                            for j in (a, b, c))
                triangles.append(tri)
            if triangles:
                lod = key[3]
                extent = 32 * (1 << lod) * voxel_m
                tiles.append({"key": key, "lod": lod, "extent": extent,
                              "x0": key[0] * extent, "y0": key[1] * extent,
                              "x1": (key[0] + 1) * extent, "y1": (key[1] + 1) * extent,
                              "z0": key[2] * extent, "z1": (key[2] + 1) * extent,
                              "triangles": triangles})
    return camera, yaw, pitch, tiles


def height_at(tile, x, y):
    best = None
    for tri in tile["triangles"]:
        (x0, y0, z0), (x1, y1, z1), (x2, y2, z2) = tri
        den = (y1 - y2) * (x0 - x2) + (x2 - x1) * (y0 - y2)
        if abs(den) < 1.0e-10:
            continue
        a = ((y1 - y2) * (x - x2) + (x2 - x1) * (y - y2)) / den
        b = ((y2 - y0) * (x - x2) + (x0 - x2) * (y - y2)) / den
        c = 1.0 - a - b
        if min(a, b, c) >= -1.0e-5:
            z = a * z0 + b * z1 + c * z2
            if best is None or z > best:
                best = z
    return best


path = sys.argv[1]
camera, yaw, pitch, tiles = read_scene(path)
by_lod = collections.defaultdict(list)
for tile in tiles:
    by_lod[tile["lod"]].append(tile)
print(f"scene={path} camera_m={camera} yaw={yaw} pitch={pitch} nonempty_mesh_tiles={len(tiles)}")
print("lod_pair shared_edges samples_both median_abs_cm p95_abs_cm max_abs_cm median_delta_cm gaps")
for lod in range(8):
    fine = by_lod[lod]
    coarse = by_lod[lod + 1]
    deltas = []
    edge_rows = []
    edge_count = 0
    gap_count = 0
    # Shared XY edges occur only where tile extents touch and overlap on the other axis.
    for a in fine:
        for b in coarse:
            zlo, zhi = max(a["z0"], b["z0"]), min(a["z1"], b["z1"])
            if zhi - zlo <= 1.0e-4:
                continue
            edge = None
            if abs(a["x1"] - b["x0"]) < 1.0e-4 or abs(b["x1"] - a["x0"]) < 1.0e-4:
                lo, hi = max(a["y0"], b["y0"]), min(a["y1"], b["y1"])
                if hi - lo > 1.0e-4:
                    edge = ("x", a["x1"] if abs(a["x1"] - b["x0"]) < 1.0e-4 else a["x0"], lo, hi)
            if edge is None and (abs(a["y1"] - b["y0"]) < 1.0e-4 or abs(b["y1"] - a["y0"]) < 1.0e-4):
                lo, hi = max(a["x0"], b["x0"]), min(a["x1"], b["x1"])
                if hi - lo > 1.0e-4:
                    edge = ("y", a["y1"] if abs(a["y1"] - b["y0"]) < 1.0e-4 else a["y0"], lo, hi)
            if edge is None:
                continue
            edge_count += 1
            axis, fixed, lo, hi = edge
            epsilon = min(0.02, a["extent"] * 0.001, b["extent"] * 0.001)
            samples = 9
            matched = 0
            edge_deltas = []
            for i in range(samples):
                t = (i + 0.5) / samples
                varying = lo + (hi - lo) * t
                if axis == "x":
                    side_a = 1.0 if abs(a["x1"] - fixed) < 1.0e-4 else -1.0
                    side_b = 1.0 if abs(b["x1"] - fixed) < 1.0e-4 else -1.0
                    za = height_at(a, fixed - side_a * epsilon, varying)
                    zb = height_at(b, fixed - side_b * epsilon, varying)
                else:
                    side_a = 1.0 if abs(a["y1"] - fixed) < 1.0e-4 else -1.0
                    side_b = 1.0 if abs(b["y1"] - fixed) < 1.0e-4 else -1.0
                    za = height_at(a, varying, fixed - side_a * epsilon)
                    zb = height_at(b, varying, fixed - side_b * epsilon)
                if za is not None and not zlo - 1.0e-3 <= za <= zhi + 1.0e-3:
                    za = None
                if zb is not None and not zlo - 1.0e-3 <= zb <= zhi + 1.0e-3:
                    zb = None
                if za is not None and zb is not None:
                    delta = (za - zb) * 100.0
                    deltas.append(delta)
                    edge_deltas.append((delta, za, zb, varying))
                    matched += 1
            if matched == 0:
                gap_count += 1
            elif edge_deltas:
                edge_rows.append((max(abs(row[0]) for row in edge_deltas),
                                  a["key"], b["key"], edge, edge_deltas))
    if deltas:
        absolute = sorted(abs(x) for x in deltas)
        signed = sorted(deltas)
        p95 = absolute[math.floor(0.95 * (len(absolute) - 1))]
        median_abs = absolute[len(absolute) // 2]
        median_delta = signed[len(signed) // 2]
        print(f"{lod}->{lod+1} {edge_count} {len(deltas)} {median_abs:.3f} {p95:.3f} {max(absolute):.3f} {median_delta:.3f} {gap_count}")
        for worst, fine_key, coarse_key, edge, rows in sorted(edge_rows, reverse=True)[:2]:
            row = max(rows, key=lambda item: abs(item[0]))
            print(f"  worst_edge max_abs_cm={worst:.3f} fine={fine_key} coarse={coarse_key} "
                  f"axis={edge[0]} fixed_m={edge[1]:.3f} span_m={edge[2]:.3f}..{edge[3]:.3f} "
                  f"sample_v_m={row[3]:.3f} fine_z_m={row[1]:.3f} "
                  f"coarse_z_m={row[2]:.3f} delta_cm={row[0]:.3f}")
    else:
        print(f"{lod}->{lod+1} {edge_count} 0 n/a n/a n/a n/a {gap_count}")