"""Grouper, round 3: the topology and lines check (plain Python 3 + numpy; reads build/checks/skin_struct.json, which the build
writes, and build/build_report.json). The goals the Administrator set after the user's "fix the topology and lines needs to be
smooth":
  A  every design line a fair curve: per vertex of its loop on the mesh, the turning angle between consecutive edges in the
     surface's tangent plane <= 4 deg (<= 6 deg inside designed round corners), changing by <= 2 deg from one vertex to the next;
  B  every free edge loop turns <= 8 deg per vertex in the surface; neighbouring quads' areas within 2.5x; the spacing between
     stations changes <= 1.3x from one interval to the next; no whole-ring station pairs closer than 2 cm;
  C  quads only, no folds, twist <= 15 deg (<= 10 deg for 99 %), no slivers (aspect <= 8, walls and chamfers aside); the
     silhouettes in the orthographic views turn <= 5 deg per vertex.
Writes build/checks/topology.json (B, C) and build/checks/lines.json (A, with every line's per-vertex turning angles).
    python build/gr_topo_check.py
"""
import json, math, os, collections
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
CHK = os.path.join(HERE, "checks")
D = json.load(open(os.path.join(CHK, "skin_struct.json")))
P = np.array(D["pos"], float)
PB = np.array(D["base"], float)
N = np.array(D["nrm"], float)
F = np.array(D["faces"], int)
FT = D["ftag"]
meta = D["meta"]
used = np.zeros(len(P), bool)
used[F.ravel()] = True


def tangent_turn(Pts, Nrm):
    """signed turning angle (deg) at each interior vertex of a polyline, in the tangent plane of the given normals"""
    out = []
    for i in range(1, len(Pts) - 1):
        n = Nrm[i] / max(1e-12, np.linalg.norm(Nrm[i]))
        a = Pts[i] - Pts[i - 1]
        b = Pts[i + 1] - Pts[i]
        a = a - n * a.dot(n)
        b = b - n * b.dot(n)
        la, lb = np.linalg.norm(a), np.linalg.norm(b)
        if la < 1e-9 or lb < 1e-9:
            out.append(0.0)
            continue
        c = max(-1.0, min(1.0, a.dot(b) / (la * lb)))
        s = np.sign(np.cross(a, b).dot(n))
        out.append(math.degrees(math.acos(c)) * (s if s != 0 else 1.0))
    return np.array(out)


def chart_normals(Pts):
    """normals for a free 3D outline (a trim's): the nearest skin vertex's normal"""
    out = []
    for p in Pts:
        d = np.einsum("ij,ij->i", PB - p, PB - p)
        d[~used] = 1e9
        out.append(N[int(np.argmin(d))])
    return np.array(out)


# ================================================================== A. the design lines
LINES = {}
CORNER = {}
for name, ids in D["lines"].items():
    if len(ids) < 3:
        continue
    # (the stern ring crosses the grooves and the band: its offset positions step by design; its curve is the base surface's.
    # FRAME_BOT's lines start on the band's front wall, a designed step: they are judged from the next station on)
    if name in ("FB0", "FB1") and len(ids) > 4:
        ids = ids[1:-1] if ids[0] != ids[-1] else ids
    pts = (PB if name == "stern_ring" else P)[ids]
    nrm = N[ids]
    closed = ids[0] == ids[-1]
    if closed:
        pts = np.vstack([pts[-2:-1], pts, pts[1:2]])[:-1]
        nrm = np.vstack([nrm[-2:-1], nrm, nrm[1:2]])[:-1]
    LINES[name] = (pts, nrm)
    CORNER_IDS = ids
    # designed round corners: the U-turn of the upper groove and of FRAME_BOT round the nose (their run through the cap)
    cap_ids = set(v for row in D["cap"] for v in row)
    CORNER[name] = np.array([v in cap_ids for v in CORNER_IDS])[1:-1] if not closed else np.zeros(len(CORNER_IDS) - 1, bool)
for name, pts in D.get("trims", {}).items():
    pts = np.array(pts, float)
    pts = np.vstack([pts[-1:], pts, pts[:1]])
    tn = D.get("trim_normals", {}).get(name)
    if tn is not None:                                    # the outline's own surface normals (the build's)
        tn = np.array(tn, float)
        nrm = np.vstack([tn[-1:], tn, tn[:1]])
    else:
        nrm = chart_normals(pts)
    LINES["trim_" + name] = (pts, nrm)
    # corners of an outline: where it turns faster than a 0.25 m radius
    seg = np.linalg.norm(np.diff(pts, axis=0), axis=1)
    t = tangent_turn(pts, nrm)
    k = np.radians(np.abs(t)) / np.maximum(1e-6, 0.5 * (seg[:-1] + seg[1:]))
    kk = np.convolve(k, np.ones(5) / 5, mode="same")
    CORNER["trim_" + name] = kk > 1 / 0.25
lines_out = {}
A_ok = True
for name, (pts, nrm) in LINES.items():
    t = tangent_turn(pts, nrm)
    corner = CORNER.get(name, np.zeros(len(t), bool))
    if len(corner) != len(t):
        corner = np.zeros(len(t), bool)
    lim = np.where(corner, 6.0, 4.0)
    over = np.abs(t) > lim + 1e-6
    dt = np.abs(np.diff(t))
    rec = {"vertices": int(len(t)), "max_turn_deg": round(float(np.abs(t).max()), 2) if len(t) else 0.0,
           "max_turn_outside_corners": round(float(np.abs(t[~corner]).max()), 2) if (~corner).any() else 0.0,
           "max_turn_in_corners": round(float(np.abs(t[corner]).max()), 2) if corner.any() else None,
           "over_limit": int(over.sum()), "max_change_deg": round(float(dt.max()), 2) if len(dt) else 0.0,
           "change_over_2": int((dt > 2.0 + 1e-6).sum()),
           "turn_deg": [round(float(v), 2) for v in t]}
    rec["pass"] = bool(rec["over_limit"] == 0 and rec["change_over_2"] == 0)
    A_ok &= rec["pass"]
    lines_out[name] = rec

# ================================================================== B / C. the grid
nv = len(P)
edges = collections.defaultdict(list)                 # edge -> faces
for fi, f in enumerate(F):
    for q in range(4):
        a, b = int(f[q]), int(f[(q + 1) % 4])
        edges[(min(a, b), max(a, b))].append(fi)
vert_edges = collections.defaultdict(set)
for (a, b) in edges:
    vert_edges[a].add((a, b)); vert_edges[b].add((a, b))
valence = {v: len(es) for v, es in vert_edges.items()}
boundary_v = set()
for e, fs in edges.items():
    if len(fs) == 1:
        boundary_v.update(e)
faces_of_edge = edges


def opposite_edge(v, e):
    """the edge across the vertex v from edge e (a regular interior vertex: the one sharing no face with e)"""
    if valence.get(v, 0) != 4 or v in boundary_v:
        return None
    fe = set(faces_of_edge[e])
    for e2 in vert_edges[v]:
        if e2 == e:
            continue
        if not (set(faces_of_edge[e2]) & fe):
            return e2
    return None


line_edges = set()
for name, ids in D["lines"].items():
    for a, b in zip(ids, ids[1:]):
        line_edges.add((min(a, b), max(a, b)))
wall_faces = set(i for i, t in enumerate(FT) if t[1] == "wall")
wall_edges = set()
for fi in wall_faces:
    f = F[fi]
    for q in range(4):
        a, b = int(f[q]), int(f[(q + 1) % 4])
        wall_edges.add((min(a, b), max(a, b)))
seen = set()
loops = []
for e0 in edges:
    if e0 in seen:
        continue
    seq = [e0]
    seen.add(e0)
    for direction in (0, 1):
        e = e0
        v = e0[direction]
        while True:
            e2 = opposite_edge(v, e)
            if e2 is None or e2 in seen:
                break
            seen.add(e2)
            if direction == 0:
                seq.insert(0, e2)
            else:
                seq.append(e2)
            v = e2[0] if e2[1] == v else e2[1]
            e = e2
    loops.append(seq)


def loop_vertices(seq):
    if len(seq) == 1:
        return list(seq[0])
    vs = []
    a, b = seq[0]
    first = a if a not in seq[1] else b
    vs.append(first)
    cur = first
    for e in seq:
        nxt = e[1] if e[0] == cur else e[0]
        vs.append(nxt)
        cur = nxt
    return vs


free_max = []
free_over = 0
worst_free = []
for seq in loops:
    if len(seq) < 3:
        continue
    if sum(1 for e in seq if e in line_edges) > 0.5 * len(seq):
        continue
    if sum(1 for e in seq if e in wall_edges) > 0.5 * len(seq):
        continue
    vs = loop_vertices(seq)
    t = tangent_turn(P[vs], N[vs])
    # vertices on a design line or a wall belong to their features; the free loop is judged at its own vertices
    on_line = np.array([any(((min(vs[i], vs[i + d]), max(vs[i], vs[i + d])) in (line_edges | wall_edges)) for d in (-1, 1))
                        for i in range(1, len(vs) - 1)])
    tf = np.abs(t[~on_line]) if (~on_line).any() else np.zeros(1)
    m = float(tf.max())
    free_max.append(m)
    free_over += int((tf > 8.0).sum())
    if m > 8.0:
        i = int(np.argmax(np.where(~on_line, np.abs(t), 0)))
        worst_free.append((round(m, 1), [round(c, 3) for c in P[vs[i + 1]]]))
worst_free.sort(reverse=True)

# quads: twist, folds, aspect, areas
def quad_metrics(f):
    a, b, c, d = P[f]
    n1, n2 = np.cross(b - a, c - a), np.cross(c - a, d - a)
    n3, n4 = np.cross(b - a, d - a), np.cross(c - b, d - b)
    def ang(u, v):
        lu, lv = np.linalg.norm(u), np.linalg.norm(v)
        if lu < 1e-14 or lv < 1e-14:
            return 0.0
        return math.degrees(math.acos(max(-1.0, min(1.0, u.dot(v) / (lu * lv)))))
    tw = max(ang(n1, n2), ang(n3, n4))
    area = 0.5 * (np.linalg.norm(n1) + np.linalg.norm(n2))
    sides = [np.linalg.norm(b - a), np.linalg.norm(c - b), np.linalg.norm(d - c), np.linalg.norm(a - d)]
    asp = max(sides) / max(1e-9, min(sides))
    fn = n1 + n2
    sn = N[f].sum(0)
    fold = fn.dot(sn) < 0 or tw > 90.0
    return tw, area, asp, fold


tw, area, asp, folds = [], [], [], 0
for fi, f in enumerate(F):
    t_, a_, s_, fo = quad_metrics(f)
    tw.append(t_); area.append(a_); asp.append(s_)
    folds += int(fo)
tw = np.array(tw); area = np.array(area); asp = np.array(asp)
flat = np.array([i not in wall_faces for i in range(len(F))])
ring_faces = np.array([t[0] in ("ring", "mirror") for t in FT])
ratios = []
for e, fs in edges.items():
    if len(fs) == 2 and flat[fs[0]] and flat[fs[1]]:
        a1, a2 = area[fs[0]], area[fs[1]]
        ratios.append(max(a1, a2) / max(1e-12, min(a1, a2)))
ratios = np.array(ratios)

# stations: the planar rings' spacing along the hull (image x of the planar stations); the transition rings' median gap
STX = meta["STX"]
dx = np.diff(STX)
sr = np.maximum(dx[1:] / dx[:-1], dx[:-1] / dx[1:])
rings = D["rings"]
gaps = [float(np.median(np.linalg.norm(P[rings[i + 1]] - P[rings[i]], axis=1))) for i in range(len(rings) - 1)]
gr = [max(a / b, b / a) for a, b in zip(gaps, gaps[1:])]
pairs_lt_2cm = int(sum(1 for g in gaps if g < 0.02))

# silhouettes of the skin in the orthographic views. SIDE and TOP: the outline is the extreme vertex of each ring of the hull
# (the crest and keel / the widest point), forward of the rings the cap's own (its centre line, its widest vertex per row);
# FRONT and BACK: the envelope of every ring's polygon, cast from the section's centre in 0.02 deg steps (each step's turning
# is the turning at the outline's vertex it passes)
def _turns(poly):
    out = []
    for i in range(1, len(poly) - 1):
        u, v = poly[i] - poly[i - 1], poly[i + 1] - poly[i]
        lu, lv = np.hypot(*u), np.hypot(*v)
        if lu < 1e-7 or lv < 1e-7:
            continue
        out.append(math.degrees(math.acos(max(-1.0, min(1.0, u.dot(v) / (lu * lv))))))
    return np.array(out) if out else np.zeros(1)


def sil_chain(vals):
    t = _turns(np.array(vals))
    return {"vertices": len(vals), "max_turn_deg": round(float(t.max()), 2), "over_5": int((t > 5.0).sum())}


# The silhouettes, judged on the base surface (the grooves' and the frames' steps are drawn features, not facets). SIDE: the
# crest and the keel of every planar ring (the extreme vertex up and down) and the cap's centre column round the nose; TOP: the
# widest vertex of every planar ring; FRONT and BACK: their outlines are made of the rings' own polygons, so each planar ring's
# turning per vertex in its plane (the cross-section's facets) is reported, with where it is worst.
NPL = len(meta["STX"])
rings_ = D["rings"]
crest = [PB[r][np.argmax(PB[r][:, 2])][[1, 2]] for r in rings_[:NPL]]
keel = [PB[r][np.argmin(PB[r][:, 2])][[1, 2]] for r in rings_[:NPL]]
cap_c = D["cap"][len(D["cap"]) // 2]
nose_profile = [PB[v][[1, 2]] for v in cap_c[::-1]]
widest = [PB[r][np.argmax(PB[r][:, 0])][[1, 0]] for r in rings_[:NPL]]


def ring_facets():
    worst, over, allt = [], 0, []
    for i, r in enumerate(rings_[:NPL]):
        Q = PB[r][:, [0, 2]]
        Q = np.vstack([Q[-1:], Q, Q[:1]])
        t = _turns(Q)
        allt.append(t)
        over += int((t > 5.0).sum())
        k = int(np.argmax(t))
        worst.append((round(float(t[k]), 1), round(float(meta["STX"][i]), 1), k))
    allt = np.concatenate(allt)
    worst.sort(reverse=True)
    return {"rings": NPL, "vertices": int(len(allt)), "max_turn_deg": round(float(allt.max()), 2), "over_5": over,
            "share_le_5": round(float((allt <= 5.0).mean()), 4), "p95": round(float(np.percentile(allt, 95)), 2),
            "worst_rings_x_col": worst[:8]}


sil = {"side_crest": sil_chain(crest), "side_keel": sil_chain(keel), "side_nose_profile": sil_chain(nose_profile),
       "top_widest": sil_chain(widest), "front_back_ring_facets": ring_facets()}

rep = json.load(open(os.path.join(HERE, "build_report.json")))
hull_audit = [a for a in rep.get("audit", []) if a.get("name") == "SM_Grouper_Hull"]
topo = {
    "A_lines_pass": bool(A_ok),
    "A_summary": {k: {kk: v[kk] for kk in ("max_turn_deg", "max_turn_outside_corners", "max_turn_in_corners", "over_limit",
                                           "max_change_deg", "change_over_2", "pass")} for k, v in lines_out.items()},
    "B_free_loops": {"loops": len(free_max), "max_turn_deg": round(max(free_max), 2) if free_max else 0.0,
                     "vertices_over_8": free_over, "worst": worst_free[:12]},
    "B_area_ratio": {"pairs": int(len(ratios)), "max": round(float(ratios.max()), 2), "over_2.5": int((ratios > 2.5).sum()),
                     "p99": round(float(np.percentile(ratios, 99)), 2)},
    "B_station_spacing": {"planar_max_ratio": round(float(sr.max()), 3), "planar_over_1.3": int((sr > 1.3).sum()),
                          "ring_gap_min_m": round(min(gaps), 4), "ring_pairs_under_2cm": pairs_lt_2cm,
                          "ring_gap_max_ratio": round(max(gr), 3), "ring_gap_ratio_over_1.3": int(sum(1 for g in gr if g > 1.3))},
    "C_quads": {"skin_quads": int(len(F)), "skin_triangles_split_in_finish": rep.get("twisted_quads_split", {}).get("skin"),
                "hull_mesh": hull_audit[0] if hull_audit else None, "folds": int(folds),
                "twist_max_deg": round(float(tw.max()), 2), "twist_over_15": int((tw > 15).sum()),
                "twist_le_10_share": round(float((tw <= 10).mean()), 4),
                "aspect_max_flat": round(float(asp[flat].max()), 2), "aspect_over_8_flat": int((asp[flat] > 8).sum()),
                "aspect_over_8_flat_rings": int((asp[flat & ring_faces] > 8).sum())},
    "C_silhouettes": sil,
    "valence": dict(collections.Counter(v for k, v in valence.items() if k not in boundary_v)),
}
json.dump(topo, open(os.path.join(CHK, "topology.json"), "w"), indent=1)
json.dump(lines_out, open(os.path.join(CHK, "lines.json"), "w"), indent=1)
print("TOPO A lines pass:", A_ok)
for k, v in topo["A_summary"].items():
    print("TOPO A", k, v)
for k in ("B_free_loops", "B_area_ratio", "B_station_spacing", "C_quads", "C_silhouettes", "valence"):
    print("TOPO", k, json.dumps(topo[k]))