#!/usr/bin/env python3
"""Audit LOD0 triangles in the resolved monolith and platform surface ROIs."""
from __future__ import annotations

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


def solve3(matrix, vector):
    rows = [list(matrix[i]) + [vector[i]] for i in range(3)]
    for column in range(3):
        pivot = max(range(column, 3), key=lambda row: abs(rows[row][column]))
        if abs(rows[pivot][column]) < 1.0e-10:
            return None
        rows[column], rows[pivot] = rows[pivot], rows[column]
        scale = rows[column][column]
        rows[column] = [value / scale for value in rows[column]]
        for row in range(3):
            if row == column:
                continue
            scale = rows[row][column]
            rows[row] = [rows[row][i] - scale * rows[column][i] for i in range(4)]
    return [rows[row][3] for row in range(3)]


def fit_radius(points, center_x, center_z, center_y, window=1.05, y_window=1.3):
    selected = sorted({(round(x, 4), round(z, 4)) for x, y, z in points
        if abs(y - center_y) <= y_window
        and abs(x - center_x) <= window and abs(z - center_z) <= window})
    if len(selected) < 8:
        return None, len(selected)
    sums = [0.0] * 9
    rhs = [0.0] * 3
    for x, z in selected:
        row = (x, z, 1.0)
        target = -(x * x + z * z)
        for i in range(3):
            rhs[i] += row[i] * target
            for j in range(3):
                sums[i * 3 + j] += row[i] * row[j]
    solution = solve3([sums[0:3], sums[3:6], [sums[6], sums[7], sums[8]]], rhs)
    if solution is None:
        return None, len(selected)
    a, b, c = solution
    radius_squared = (a * a + b * b) * 0.25 - c
    if radius_squared <= 0.0:
        return None, len(selected)
    return math.sqrt(radius_squared), len(selected)


def read_mesh(path: Path):
    triangles = []
    vertices_all = []
    selected_tiles = 0
    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 header")
        version, voxel_um, tile_count = struct.unpack_from("<3I", prefix, 8)
        if version != 1 or voxel_um <= 0 or tile_count > 1_000_000:
            raise ValueError(f"{path}: invalid mesh version or size")
        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, _aux_count, vertex_count, index_count, mask_count = \
                struct.unpack_from("<5I", header, 16)
            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)
            indices = array.array("I")
            indices.frombytes(index_bytes)
            if sys.byteorder != "little":
                positions.byteswap()
                indices.byteswap()
            points = [tuple(float(positions[index * 3 + axis])
                            for axis in range(3)) for index in range(vertex_count)]
            vertices_all.extend(points)
            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(points[index] 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])
                area_twice = math.sqrt(sum(value * value for value in cross))
                if area_twice <= 1.0e-10:
                    continue
                center = tuple(sum(vertex[axis] for vertex in verts) / 3.0
                               for axis in range(3))
                triangles.append((center, area_twice * 0.5, verts))
        if stream.read(1):
            raise ValueError(f"{path}: trailing mesh data")
    return triangles, vertices_all, selected_tiles


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("meshes", type=Path, nargs="+")
    parser.add_argument("layout", type=Path)
    parser.add_argument("--view", required=True, choices=("combined",))
    args = parser.parse_args()
    layout = json.loads(args.layout.read_text(encoding="utf-8"))
    motifs = layout.get("stamps", [])
    monolith = next(row for row in motifs if row["kind"] == 5)
    platform = next(row for row in motifs if row["kind"] == 6)
    mesh_rows = [read_mesh(path) for path in args.meshes]
    triangles = [triangle for row in mesh_rows for triangle in row[0]]
    vertices = [vertex for row in mesh_rows for vertex in row[1]]
    tiles = sum(row[2] for row in mesh_rows)
    if not triangles:
        raise ValueError("meshes have no LOD0 surface triangles")

    def metres(point):
        return tuple(value / 1_000_000.0 for value in point)

    mono = metres(monolith["position_um"])
    floor = monolith["floor_um"] / 1_000_000.0
    radius = monolith["half_extents_um"][0] / 1_000_000.0
    half_height = monolith["half_extents_um"][2] / 1_000_000.0
    mono_round = monolith["rounding_radius_um"] / 1_000_000.0
    plat = metres(platform["position_um"])
    half = tuple(value / 1_000_000.0 for value in platform["half_extents_um"])
    plat_round = platform["rounding_radius_um"] / 1_000_000.0
    top = floor + 2.0 * half_height

    rois = {
        "monolith_foot": lambda p: math.hypot(p[0] - mono[0], p[1] - mono[1]) <= radius + 1.4
            and floor - 1.0 <= p[2] <= floor + 2.1,
        "monolith_cap": lambda p: math.hypot(p[0] - mono[0], p[1] - mono[1]) <= radius + 1.6
            and top - 2.0 <= p[2] <= top + 1.1,
        "platform_edges": lambda p: abs(p[0] - plat[0]) <= half[0] + 1.2
            and abs(p[1] - plat[1]) <= half[1] + 1.2
            and floor - 0.8 <= p[2] <= floor + 2.4 * half[2] + 0.8,
    }
    roi_counts = {name: sum(1 for center, _area, _verts in triangles if predicate(center))
                  for name, predicate in rois.items()}
    bad_area = sum(1 for _center, area, _verts in triangles if not math.isfinite(area) or area <= 0.0)
    foot_radius, foot_points = fit_radius(vertices, mono[0] + radius - mono_round,
        floor + mono_round, mono[1], window=1.1)
    cap_radius, cap_points = fit_radius(vertices, mono[0] + radius - mono_round,
        top - mono_round, mono[1], window=1.1)
    platform_radius, platform_points = fit_radius(vertices,
        plat[0] + half[0] - plat_round, floor + 2.0 * half[2] - plat_round,
        plat[1], window=1.1)
    radii = (foot_radius, cap_radius, platform_radius)
    fit_counts = (foot_points, cap_points, platform_points)
    def surface_pass(tile_count, degenerate_count, rois, measured_radii):
        covered = all(count >= 8 for count in rois.values())
        rounded = all(value is not None and value >= 0.5 for value in measured_radii)
        return tile_count > 0 and degenerate_count == 0 and covered and rounded

    radius_ok = all(value is not None and value >= 0.5 for value in radii)
    coverage_ok = all(count >= 8 for count in roi_counts.values())
    missing_roi = dict(roi_counts)
    missing_roi["monolith_foot"] = 0
    sharp_edge = list(radii)
    sharp_edge[0] = 0.25
    controls_rejected = (
        not surface_pass(tiles, bad_area, missing_roi, radii)
        and not surface_pass(tiles, bad_area, roi_counts, sharp_edge)
        and not surface_pass(0, 0, roi_counts, radii)
    )
    passed = surface_pass(tiles, bad_area, roi_counts, radii) and controls_rejected
    for name, count in roi_counts.items():
        print(f"[MotifSurfaceROI] view={args.view} name={name} triangles={count} "
              f"whole_surface=LOD0 status={'PASS' if count >= 8 else 'FAIL'}")
    if all(value is not None for value in radii):
        print(f"[MotifFilletMeasure] view={args.view} foot_radius_m={foot_radius:.3f} "
              f"cap_radius_m={cap_radius:.3f} platform_radius_m={platform_radius:.3f} "
              f"fit_points={fit_counts} target_min_m=0.500")
    else:
        print(f"[MotifFilletMeasure] view={args.view} foot_radius_m={foot_radius} "
              f"cap_radius_m={cap_radius} platform_radius_m={platform_radius} "
              f"fit_points={fit_counts} target_min_m=0.500")
    print(f"[MotifSurfaceAudit] view={args.view} meshes={len(mesh_rows)} lod0_tiles={tiles} triangles={len(triangles)} "
        f"degenerate={bad_area} feature_rois={len(roi_counts)} covered={coverage_ok} "
        f"fillets_ge_0.5m={radius_ok} missing_roi_control={'REJECTED' if not surface_pass(tiles, bad_area, missing_roi, radii) else 'FAILED'} "
        f"sharp_fillet_control={'REJECTED' if not surface_pass(tiles, bad_area, roi_counts, sharp_edge) else 'FAILED'} "
        f"empty_mesh_control={'REJECTED' if not surface_pass(0, 0, roi_counts, radii) else 'FAILED'} "
        f"negative_controls={'REJECTED' if controls_rejected else 'FAILED'} "
          f"status={'PASS' if passed else 'FAIL'}")
    return 0 if passed else 1


if __name__ == "__main__":
    try:
        raise SystemExit(main())
    except Exception as error:
        print(f"[MotifSurfaceAudit] status=FAIL error={error}", file=sys.stderr)
        raise SystemExit(1)