#!/usr/bin/env python3
"""Measure axial normal ripple in a camera-matched tunnel mesh window."""
import importlib.util
import math
import re
import sys
from pathlib import Path


ROOT = Path(__file__).resolve().parents[4]
RESEARCH = Path(__file__).resolve().parent
RING_AUDIT = ROOT / "Standalone/Research/WorldAssemblerM3/Linux32/AuditRingTerraceL32.py"
SPEC = importlib.util.spec_from_file_location("ring_profile_l32", RING_AUDIT)
PROFILE = importlib.util.module_from_spec(SPEC)
SPEC.loader.exec_module(PROFILE)
CAMERA_LOG = ROOT / "Saved/Linux/l34/run-capture-standard-set-final-radius3.log"
CAMERA_RE = re.compile(
    r"^\[L34CaptureCamera\] view=(\S+).*?camera_m=\(([^)]*)\) "
    r"target_m=\(([^)]*)\)", re.M)
MAXIMUM_AXIAL_NORMAL_AMPLITUDE = 0.0255


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


def fit(samples, spacing=0.25):
    if len(samples) < 32:
        return 0.0, 0.0
    best_period, best_amplitude = 0.0, 0.0
    for index in range(8, 129):
        period = index * 0.0625
        if period > 8.0:
            break
        omega = 2.0 * math.pi / period
        cc = ss = cs = yc = ys = 0.0
        for i, value in enumerate(samples):
            phase = omega * (i * spacing)
            c, s = math.cos(phase), math.sin(phase)
            cc += c * c
            ss += s * s
            cs += c * s
            yc += value * c
            ys += value * s
        determinant = cc * ss - cs * cs
        if determinant <= 1.0e-9:
            continue
        cosine = (yc * ss - ys * cs) / determinant
        sine = (ys * cc - yc * cs) / determinant
        amplitude = math.hypot(cosine, sine)
        if amplitude > best_amplitude:
            best_period, best_amplitude = period, amplitude
    return best_period, best_amplitude


def main():
    if len(sys.argv) != 3:
        raise SystemExit("usage: AuditTunnelMeshNormalsL35.py MESH L34_VIEW_NAME")
    mesh = Path(sys.argv[1])
    view_name = sys.argv[2]
    cameras = {name: (tuple(map(float, camera.split(","))),
                      tuple(map(float, target.split(","))))
               for name, camera, target in CAMERA_RE.findall(
                   CAMERA_LOG.read_text(encoding="utf-8", errors="replace"))}
    if view_name not in cameras:
        raise SystemExit(f"missing L34 camera: {view_name}")
    camera, target = cameras[view_name]
    axis = normalize(tuple(target[i] - camera[i] for i in range(3)))
    down = (0.0, 0.0, -1.0)
    axis_down = sum(axis[i] * down[i] for i in range(3))
    floor_dir = normalize(tuple(down[i] - axis_down * axis[i] for i in range(3)))
    x_min, x_max = min(camera[0], target[0]) - 8.0, max(camera[0], target[0]) + 8.0
    y_min, y_max = min(camera[1], target[1]) - 8.0, max(camera[1], target[1]) + 8.0
    z_min, z_max = min(camera[2], target[2]) - 8.0, max(camera[2], target[2]) + 8.0
    voxel_um, triangles, tiles = PROFILE.read_scene(
        str(mesh), None, z_min, z_max, (x_min, x_max, y_min, y_max))
    spacing = 0.25
    bin_count = 48
    sums = [0.0] * bin_count
    counts = [0] * bin_count
    for _, _, verts in triangles:
        centroid = tuple(sum(v[0][axis_index] for v in verts) / 3.0
                         for axis_index in range(3))
        delta = tuple(centroid[i] - camera[i] for i in range(3))
        along = sum(delta[i] * axis[i] for i in range(3))
        if not 0.0 <= along < bin_count * spacing:
            continue
        radial = tuple(delta[i] - along * axis[i] for i in range(3))
        radius = math.sqrt(sum(value * value for value in radial))
        if not 0.5 <= radius <= 4.0:
            continue
        radial_unit = tuple(value / radius for value in radial)
        if sum(radial_unit[i] * floor_dir[i] for i in range(3)) < 0.25:
            continue
        normal = normalize(tuple(sum(v[1][axis_index] for v in verts) / 3.0
                                 for axis_index in range(3)))
        axial = abs(sum(normal[i] * axis[i] for i in range(3)))
        bin_index = int(along / spacing)
        sums[bin_index] += axial
        counts[bin_index] += 1
    valid = [(i, sums[i] / counts[i]) for i in range(bin_count) if counts[i]]
    if len(valid) < 32:
        print(f"[TunnelMeshNormalAudit] mesh={mesh} view={view_name} tiles={tiles} "
              f"bins={len(valid)}/{bin_count} status=INSUFFICIENT")
        raise SystemExit(1)
    values = [sums[i] / counts[i] if counts[i] else 0.0 for i in range(bin_count)]
    mean = sum(value for _, value in valid) / len(valid)
    values = [value - mean for value in values]
    period, amplitude = fit(values, spacing)
    passed = amplitude < MAXIMUM_AXIAL_NORMAL_AMPLITUDE
    print(f"[TunnelMeshNormalAudit] mesh={mesh} view={view_name} voxel_um={voxel_um} "
          f"tiles={tiles} bins={len(valid)}/{bin_count} floor_sector=0.25 "
          f"period_m={period:.3f} axial_normal_amplitude={amplitude:.5f} "
          f"axial_normal_amplitude_deg={math.degrees(math.asin(min(1.0, amplitude))):.3f} "
          f"threshold={MAXIMUM_AXIAL_NORMAL_AMPLITUDE:.4f} "
          f"status={'PASS' if passed else 'FAIL'}")
    if not passed:
        raise SystemExit(1)


if __name__ == "__main__":
    main()