"""The encoding's own error, layer by layer: the descriptor's rings at a layer centre against the exact mesh slice there
(two-way: ring vertices to the slice boundary, slice boundary to the rings), % of the diagonal. No points, no query.

    env -u PYTHONPATH .venv/bin/python scripts/ring_error.py <cached .fcp> <object> [layer_stride]
"""
import sys
from pathlib import Path

import numpy as np
from scipy.spatial import cKDTree

sys.path.insert(0, str(Path(__file__).parent))
from common import CACHE, fcp, load_mesh  # noqa: E402
from fcpq import FCPQuery  # noqa: E402

LPAD = 1.1  # fcp._slice.L


def resample(poly, step):
    seg = np.diff(poly, axis=0)
    ln = np.linalg.norm(seg, axis=1)
    out = [poly[:1]]
    for a, s, l in zip(poly[:-1], seg, ln):
        k = max(1, int(np.ceil(l / step)))
        out.append(a + s * (np.arange(1, k + 1)[:, None] / k))
    return np.vstack(out)


def main():
    f, name = sys.argv[1], sys.argv[2]
    stride = int(sys.argv[3]) if len(sys.argv) > 3 else 8
    d = fcp.FCP.from_bytes((CACHE / f).read_bytes())
    m = load_mesh(name)
    diag = np.linalg.norm(m.extents)
    V = (np.asarray(m.vertices) - d.pose.center) @ d.pose.rotation  # aligned, world units
    import trimesh
    ma = trimesh.Trimesh(V, m.faces, process=False)
    q = FCPQuery(d, nt=1024)
    sc = d.pose.scale
    L = d.layers
    rows = []
    for layer in range(0, L, stride):
        zc = -LPAD + (layer + 0.5) * 2 * LPAD / L
        ps = [p for p in range(q.P) if q.h0[p] <= layer <= q.h1[p]]
        sec = ma.section(plane_origin=[0, 0, zc * sc[2]], plane_normal=[0, 0, 1])
        if not ps and sec is None:
            continue
        ex = np.vstack([resample(np.asarray(c)[:, :2], diag * 1e-5) for c in sec.discrete]) if sec is not None else np.zeros((0, 2))
        fr = np.vstack([q.rings[q.first[p] + layer - q.h0[p]] * sc[:2] for p in ps]) if ps else np.zeros((0, 2))
        if len(ex) == 0 or len(fr) == 0:
            rows.append((layer, zc, len(ps), np.inf, np.inf))
            continue
        a = cKDTree(ex).query(fr)[0].max() / diag * 100  # descriptor ring off the true boundary
        b = cKDTree(resample(np.vstack([np.vstack([q.rings[q.first[p] + layer - q.h0[p]], q.rings[q.first[p] + layer - q.h0[p]][:1]]) * sc[:2] for p in ps[:1]]), 1) if False else fr).query(ex)[0].max() / diag * 100  # true boundary not covered
        rows.append((layer, zc, len(ps), a, b))
    r = np.array(rows, float)
    fin = np.isfinite(r[:, 3])
    print(f"{f}: {len(r)} layers checked, parts {q.P}")
    print(f"  ring -> slice max %diag: median {np.median(r[fin, 3]):.4f}  p95 {np.percentile(r[fin, 3], 95):.4f}  max {r[fin, 3].max():.4f}")
    print(f"  slice -> ring max %diag: median {np.median(r[fin, 4]):.4f}  p95 {np.percentile(r[fin, 4], 95):.4f}  max {r[fin, 4].max():.4f}")
    print(f"  layers with a part but no slice, or a slice but no part: {int((~fin).sum())}")
    worst = r[fin][np.argsort(-np.maximum(r[fin, 3], r[fin, 4]))[:12]]
    for w in worst:
        print(f"    layer {int(w[0]):5d} z {w[1]:+.4f} parts {int(w[2])}: ring->slice {w[3]:.4f}  slice->ring {w[4]:.4f}")
    np.save(CACHE / f"ringerr_{Path(f).stem}.npy", r)


if __name__ == "__main__":
    main()