# Stage 16c (diagnosis only): settle whether the visible lines are TOPOLOGY or GEOMETRY, and
# inventory the real holes. Stage 15 asserted "the two sides of every seam are disconnected
# vertex runs"; stage 17's weld attempt on that premise tore the mesh (830 -> 24k boundary
# edges), which is itself evidence the premise is wrong.
#
#   blender --background --python 16c_topology.py -- <blend>
#
# Tests
# 1. LOCAL EDGE LENGTH — the scale everything else must be judged against. A "twin at 0.8 mm"
#    means nothing until you know the mesh spacing is ~1.2 mm; at that scale a non-adjacent
#    vertex 0.8 mm away is just a 2-ring neighbour, not a crack. This is the control the earlier
#    split-panel test lacked.
# 2. TRUE TWIN TEST — nearest vertex that is outside the 3-ring topological neighbourhood,
#    normalised by local edge length. A real crack gives a spike at ratio << 1; ordinary mesh
#    gives a distribution centred near/above 1.
# 3. HOLE INVENTORY — boundary edges grouped into loops, with size and location, so filling can
#    be per-loop instead of one net (which is what produced 39k non-manifold edges).
# 4. RELIEF SCALE-SPACE — for the kink verts, small-scale vs large-scale normal offset. A scan
#    seam is thin (small-scale relief, no large-scale relief); real anatomy (underbust fold,
#    gluteal fold, clavicle) has both. This is the discriminator a heal mask must use so it
#    erases seams without erasing her.
import bpy, bmesh, sys, time
import numpy as np
from mathutils import Vector
from mathutils.kdtree import KDTree
from collections import deque

argv = sys.argv[sys.argv.index("--") + 1:]
BLEND = argv[0]
t0 = time.time()
UNIT_MM = 1815.0     # 1 mesh unit = 1.815 m (body is 0.979 units for 1.777 m)


def log(m):
    print(f"[topo {time.time()-t0:6.1f}s] {m}", flush=True)


bpy.ops.wm.open_mainfile(filepath=BLEND)
ob = max([o for o in bpy.data.objects if o.type == 'MESH'],
         key=lambda o: len(o.data.vertices))
me = ob.data
n_v = len(me.vertices)
co = np.empty(n_v * 3)
me.vertices.foreach_get("co", co)
co = co.reshape(-1, 3)
nrm = np.empty(n_v * 3)
me.vertices.foreach_get("normal", nrm)
nrm = nrm.reshape(-1, 3)
ev = np.empty(len(me.edges) * 2, dtype=np.int32)
me.edges.foreach_get("vertices", ev)
ev = ev.reshape(-1, 2)
log(f"{n_v}v {len(me.polygons)}f {len(me.edges)}e")

# ---- 1. local edge length ----
elen = np.linalg.norm(co[ev[:, 0]] - co[ev[:, 1]], axis=1)
acc = np.zeros(n_v)
cntv = np.zeros(n_v)
np.add.at(acc, ev[:, 0], elen)
np.add.at(acc, ev[:, 1], elen)
np.add.at(cntv, ev[:, 0], 1.0)
np.add.at(cntv, ev[:, 1], 1.0)
L = acc / np.maximum(cntv, 1)
print(f"EDGE LENGTH mesh units: mean {elen.mean():.6f} ({elen.mean()*UNIT_MM:.2f} real mm)  "
      f"p05 {np.percentile(elen,5):.6f}  p95 {np.percentile(elen,95):.6f}")

order = np.concatenate([ev[:, 0], ev[:, 1]])
nbr = np.concatenate([ev[:, 1], ev[:, 0]])
srt = np.argsort(order, kind="stable")
o_s, n_s = order[srt], nbr[srt]
ptr = np.searchsorted(o_s, np.arange(n_v + 1))
cnt = np.maximum(np.diff(ptr), 1)


def nbr_mean(X):
    a = np.add.reduceat(X[n_s], ptr[:-1], axis=0)
    a[np.diff(ptr) == 0] = X[np.diff(ptr) == 0]
    return a / cnt[:, None]


N = nrm.copy()
for _ in range(5):
    N = nbr_mean(N)
    N /= np.maximum(np.linalg.norm(N, axis=1, keepdims=True), 1e-12)
ang = np.degrees(np.arccos(np.clip((nrm * N).sum(axis=1), -1, 1)))
torso = (co[:, 2] > 0.28) & (co[:, 2] < 0.90)
kink = torso & (ang > 12.0)
log(f"kink>12deg in torso: {kink.sum()}")

# ---- 2. true twin test (outside the 3-ring) ----
kidx = np.nonzero(kink)[0]
sample = kidx[::max(1, len(kidx) // 3000)]
ctrl = np.nonzero(torso & (ang < 3.0))[0]
ctrl = ctrl[::max(1, len(ctrl) // 3000)]      # control: quiet skin, same test
kd = KDTree(n_v)
for i in range(n_v):
    kd.insert(Vector(co[i]), i)
kd.balance()
adj = {}
for a, b in ev:
    adj.setdefault(int(a), set()).add(int(b))
    adj.setdefault(int(b), set()).add(int(a))


def ring3(i):
    seen = {i}
    frontier = {i}
    for _ in range(3):
        nxt = set()
        for v in frontier:
            nxt |= adj.get(v, set())
        nxt -= seen
        seen |= nxt
        frontier = nxt
    return seen


def twin_ratio(idxs):
    out = []
    for i in idxs:
        excl = ring3(int(i))
        best = None
        for (_, j, d) in kd.find_range(Vector(co[i]), L[i] * 2.0):
            if int(j) in excl:
                continue
            if float(nrm[i] @ nrm[j]) < 0.0:
                continue      # opposite-facing surface (thigh vs thigh) is not a seam twin
            if best is None or d < best:
                best = d
        out.append(best / L[i] if best is not None else np.nan)
    return np.array(out)


tr_k = twin_ratio(sample)
tr_c = twin_ratio(ctrl)
log("twin test done")
for nm, tr in (("KINK verts", tr_k), ("CONTROL quiet skin", tr_c)):
    v = tr[~np.isnan(tr)]
    print(f"TWIN-RATIO {nm}: n={len(v)}/{len(tr)}  "
          f"p05={np.percentile(v,5):.2f} p25={np.percentile(v,25):.2f} "
          f"median={np.median(v):.2f} p75={np.percentile(v,75):.2f}"
          if len(v) else f"TWIN-RATIO {nm}: no hits")
    if len(v):
        print(f"    fraction with a non-3-ring vertex closer than 0.35x edge length: "
              f"{100.0*(v<0.35).mean():.1f}%   (<0.6x: {100.0*(v<0.6).mean():.1f}%)")

# ---- 3. hole inventory ----
bm = bmesh.new()
bm.from_mesh(me)
open_e = [e for e in bm.edges if len(e.link_faces) == 1]
print(f"\nHOLES: {len(open_e)} boundary edges, "
      f"{len([e for e in bm.edges if len(e.link_faces) > 2])} non-manifold edges")
eset = set(e.index for e in open_e)
emap = {e.index: e for e in open_e}
seen = set()
loops = []
for e in open_e:
    if e.index in seen:
        continue
    q = deque([e.index])
    seen.add(e.index)
    comp = []
    while q:
        ei = q.popleft()
        cur = emap[ei]
        comp.append(cur)
        for v in cur.verts:
            for e2 in v.link_edges:
                if e2.index in eset and e2.index not in seen:
                    seen.add(e2.index)
                    q.append(e2.index)
    loops.append(comp)
loops.sort(key=len, reverse=True)
print(f"HOLES: {len(loops)} separate boundary loops")
for i, lp in enumerate(loops[:14]):
    zs = [v.co.z for e in lp for v in e.verts]
    xs = [v.co.x for e in lp for v in e.verts]
    ys = [v.co.y for e in lp for v in e.verts]
    print(f"   loop{i:2d}: {len(lp):5d} edges  z {min(zs):.3f}-{max(zs):.3f}  "
          f"x {min(xs):+.3f}..{max(xs):+.3f}  y {min(ys):+.3f}..{max(ys):+.3f}")
small = [lp for lp in loops if len(lp) <= 60]
print(f"HOLES: {len(small)} loops <=60 edges (safe per-loop fills), "
      f"{len(loops)-len(small)} larger")
bm.free()

# ---- 4. relief scale-space ----
def smooth_n(X, k):
    Y = X.copy()
    for _ in range(k):
        Y = nbr_mean(Y)
    return Y


sm_s = smooth_n(co, 8)
sm_l = smooth_n(co, 60)
dev_s = ((co - sm_s) * N).sum(axis=1)
dev_l = ((co - sm_l) * N).sum(axis=1)
print("\nRELIEF SCALE-SPACE (real mm along the smooth normal)")
for nm, m in (("kink>12deg", kink),
              ("quiet skin", torso & (ang < 3.0))):
    if not m.sum():
        continue
    print(f"  {nm}: |small-scale| median {np.median(np.abs(dev_s[m]))*UNIT_MM:.3f} mm, "
          f"p95 {np.percentile(np.abs(dev_s[m]),95)*UNIT_MM:.3f} mm | "
          f"|large-scale| median {np.median(np.abs(dev_l[m]))*UNIT_MM:.3f} mm, "
          f"p95 {np.percentile(np.abs(dev_l[m]),95)*UNIT_MM:.3f} mm")
band = np.abs(dev_s) * UNIT_MM
bandpass = torso & (band > 0.35) & (np.abs(dev_l) * UNIT_MM < 1.6)
print(f"  BAND-PASS candidate mask (thin relief >0.35 mm, broad relief <1.6 mm): "
      f"{int(bandpass.sum())} verts")
hist, edges = np.histogram(co[bandpass, 2], bins=14)
for c, lo, hi in zip(hist, edges[:-1], edges[1:]):
    if c > 50:
        print(f"     z {lo:.3f}-{hi:.3f}: {c:6d}")
print("TOPO_DONE")