# ============================================================================================
# REJECTED 2026-08-06 — DO NOT RUN. Its premise is false and running it DAMAGES the mesh.
#
# It assumes the scan-panel seams are disconnected vertex runs that need welding. They are not.
# 16c_topology.py measured the opposite: the whole mesh has only 830 boundary edges in 689 loops
# of 3-5 edges, and for kink verts the nearest vertex outside the 3-ring sits at 1.58x the local
# edge length (median) versus 1.73x for control skin — no crack network exists.
#
# The default eps of 0.0012 units is about the 1.85 mm MEAN EDGE LENGTH, so "twins" were mostly
# ordinary neighbours. Welding them collapsed real triangles and tore the mesh:
#     boundary edges   830 -> 24,114
#     non-manifold     1,649 -> 39,854
# The output was discarded. The lines are GEOMETRY, not topology — see 26_finish.py part A, which
# removes them with a ring-median filter.
#
# Kept only as the record of the dead end. Everything below argues confidently for a fix that
# does not work.
# ============================================================================================
# Stage 17: TRUE panel weld + hole closure, self-detected on the CURRENT mesh.
#
#   blender --background --python 17_panel_weld.py -- <in.blend> <out.blend> [eps_mm_units]
#
# WHY THIS EXISTS (stage 15 already tried to weld and the lines survived)
# Stage 15 chose its merge set by mapping the RAW glb's boundary verts onto the working mesh at
# 2.5 mm. By then the sculpt + five membrane passes had moved those verts further than the
# tolerance, so most of the panel network was never selected. Measured after stage 15
# (16_diagnose): 58.7% of shading-kink verts STILL have a non-adjacent twin within 0.8 mm, i.e.
# the two sides of each seam are separate vertex runs that shade independently. A crack survives
# any amount of vertex MOVEMENT, which is why 11-14's membranes could not remove it.
#
# So this stage never consults the raw mesh. It finds the defect where it actually is:
#   candidates = shading-kink verts ∪ boundary verts, grown 2 rings
#   a PAIR is two candidates within eps that are NOT edge-adjacent and whose normals agree
#     (dot > 0.5). The normal test is what makes this safe: two sides of one seam face the same
#     way, whereas two surfaces that merely come close (inner thighs, armpit) face opposite ways
#     and are never paired.
#   pairs -> union-find -> bmesh.ops.weld_verts with an explicit targetmap.
# weld_verts (not remove_doubles) because the targetmap is exact: only vertices this script has
# validated get merged, so no collateral merge is possible inside the eps ball.
#
# Loops keep their own UVs through a weld, so the baked atlas is unaffected — a welded vertex
# simply carries two UV corners, which is what every UV seam already is.
#
# Then holes: the body should be watertight below the chin. The largest boundary cluster is the
# bra-bow excision at the sternum (z 0.648-0.729, 329 edges) — that hole is the black gash that
# reads as a broken cleavage. Filled with bmesh triangle_fill and smoothed by the next stage.
# Head openings (z > HEAD_Z: mouth, eyes, nostrils) are left alone: they are supposed to be open.
import bpy, bmesh, sys, time, math
import numpy as np
from mathutils import Vector
from mathutils.kdtree import KDTree

argv = sys.argv[sys.argv.index("--") + 1:]
BLEND, OUT = argv[0], argv[1]
EPS = float(argv[2]) if len(argv) > 2 else 0.0012      # mesh units (~2.2 mm real: 1 unit=1.815 m)
t0 = time.time()

KINK_DEG = 4.0          # generous: anything that could read as a line
DEV_MIN = 0.00008       # or a small-scale bump/groove this deep (mesh units)
NORM_DOT = 0.5          # same-facing test that makes the weld safe
GROW = 2
Z_LO, Z_HI = 0.04, 0.90   # below the chin, above the soles
X_MAX = 0.36              # excludes hands/wrists so fingers can never weld together
HEAD_Z = 0.90             # holes above this are real openings (mouth/eyes/nostrils)


def log(m):
    print(f"[weld {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)
log(f"in: {n_v}v {len(me.polygons)}f custom_normals={me.has_custom_normals}")

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)
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):
    acc = np.add.reduceat(X[n_s], ptr[:-1], axis=0)
    acc[np.diff(ptr) == 0] = X[np.diff(ptr) == 0]
    return acc / cnt[:, None]


# ---- candidates: shading kinks + small-scale relief + existing boundary ----
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)))
sm = co.copy()
for _ in range(12):
    sm = nbr_mean(sm)
dev = np.abs(((co - sm) * N).sum(axis=1))

bm = bmesh.new()
bm.from_mesh(me)
bm.verts.ensure_lookup_table()
bnd_v = np.zeros(n_v, dtype=bool)
for e in bm.edges:
    if len(e.link_faces) == 1:
        bnd_v[e.verts[0].index] = True
        bnd_v[e.verts[1].index] = True
n_bnd0 = len([e for e in bm.edges if len(e.link_faces) == 1])
n_nm0 = len([e for e in bm.edges if len(e.link_faces) > 2])
log(f"before: boundary_edges={n_bnd0} nonmanifold_edges={n_nm0} boundary_verts={bnd_v.sum()}")

zone = (co[:, 2] > Z_LO) & (co[:, 2] < Z_HI) & (np.abs(co[:, 0]) < X_MAX)
cand = zone & ((ang > KINK_DEG) | (dev > DEV_MIN) | bnd_v)
for _ in range(GROW):
    hit = cand[ev[:, 0]] | cand[ev[:, 1]]
    c2 = cand.copy()
    c2[ev[:, 0]] |= hit
    c2[ev[:, 1]] |= hit
    cand = c2 & zone
cidx = np.nonzero(cand)[0]
log(f"candidates: {len(cidx)} verts (kink>{KINK_DEG}deg {(zone&(ang>KINK_DEG)).sum()}, "
    f"dev {(zone&(dev>DEV_MIN)).sum()}, boundary {(zone&bnd_v).sum()}, +{GROW} rings)")

# adjacency lookup restricted to candidates
adj = {}
sel_mask = cand
mE = sel_mask[ev[:, 0]] & sel_mask[ev[:, 1]]
for a, b in ev[mE]:
    adj.setdefault(int(a), set()).add(int(b))
    adj.setdefault(int(b), set()).add(int(a))

kd = KDTree(len(cidx))
for j, i in enumerate(cidx):
    kd.insert(Vector(co[i]), j)
kd.balance()
log("KD built over candidates")

# ---- pair up twins ----
pairs = []
dists = []
for j, i in enumerate(cidx):
    ai = adj.get(int(i), ())
    for (_, k, d) in kd.find_range(Vector(co[i]), EPS):
        o = int(cidx[k])
        if o <= int(i) or o in ai:
            continue
        if float(nrm[i] @ nrm[o]) < NORM_DOT:
            continue
        pairs.append((int(i), o))
        dists.append(d)
log(f"pairs: {len(pairs)}")
if dists:
    dh = np.array(dists)
    print("PAIR-DISTANCE histogram (mesh units, 1 unit = 1.815 m):")
    hist, edges = np.histogram(dh, bins=np.linspace(0, EPS, 9))
    for c, lo, hi in zip(hist, edges[:-1], edges[1:]):
        print(f"   {lo*1000:5.3f}-{hi*1000:5.3f} mm-units: {c:7d}  "
              f"({lo*1815:5.2f}-{hi*1815:5.2f} real mm)")

# ---- union-find ----
parent = {}


def find(x):
    parent.setdefault(x, x)
    while parent[x] != x:
        parent[x] = parent[parent[x]]
        x = parent[x]
    return x


def union(a, b):
    ra, rb = find(a), find(b)
    if ra != rb:
        parent[min(ra, rb)] = min(ra, rb)
        parent[max(ra, rb)] = min(ra, rb)


for a, b in pairs:
    union(a, b)
clusters = {}
for v in list(parent):
    clusters.setdefault(find(v), []).append(v)
sizes = np.array([len(c) for c in clusters.values()])
log(f"clusters: {len(clusters)} covering {int(sizes.sum())} verts "
    f"(max {sizes.max() if len(sizes) else 0}, mean {sizes.mean() if len(sizes) else 0:.2f})")

# guard: a huge cluster would mean the eps ball is chaining across a whole region
if len(sizes) and sizes.max() > 40:
    log(f"WARNING: largest cluster {sizes.max()} verts — chaining suspected; "
        f"clusters >40 verts are SKIPPED")

targetmap = {}
merged_verts = 0
for root, members in clusters.items():
    if len(members) < 2 or len(members) > 40:
        continue
    members = sorted(members)
    keep = members[0]
    ctr = co[members].mean(axis=0)
    bm.verts[keep].co = Vector(ctr)
    for m in members[1:]:
        targetmap[bm.verts[m]] = bm.verts[keep]
        merged_verts += 1
log(f"targetmap: merging {merged_verts} verts into {len(set(targetmap.values()))} survivors")

if targetmap:
    bmesh.ops.weld_verts(bm, targetmap=targetmap)
    log("weld_verts done")

# ---- close holes below the chin ----
bm.edges.ensure_lookup_table()
open_e = [e for e in bm.edges if len(e.link_faces) == 1]
body_e = [e for e in open_e
          if max(v.co.z for v in e.verts) < HEAD_Z and min(v.co.z for v in e.verts) > Z_LO * 0.5]
log(f"open edges after weld: {len(open_e)} total, {len(body_e)} below the chin -> filling")
if body_e:
    res = bmesh.ops.triangle_fill(bm, use_beauty=True, use_dissolve=False, edges=body_e)
    log(f"triangle_fill created {len(res.get('geom', []))} elements")
    # anything still open: try holes_fill as a second pass
    bm.edges.ensure_lookup_table()
    still = [e for e in bm.edges if len(e.link_faces) == 1
             and max(v.co.z for v in e.verts) < HEAD_Z]
    if still:
        bmesh.ops.holes_fill(bm, edges=still, sides=0)
        log(f"holes_fill on {len(still)} remaining open edges")

bm.edges.ensure_lookup_table()
n_bnd1 = len([e for e in bm.edges if len(e.link_faces) == 1])
n_nm1 = len([e for e in bm.edges if len(e.link_faces) > 2])
dgn = [f for f in bm.faces if f.calc_area() < 1e-12]
if dgn:
    bmesh.ops.delete(bm, geom=dgn, context='FACES')
    log(f"deleted {len(dgn)} degenerate faces")
bm.to_mesh(me)
bm.free()
me.update()
n_after = len(me.vertices)
log(f"after: {n_after}v (-{n_v - n_after}) boundary_edges={n_bnd1} nonmanifold_edges={n_nm1}")

# ---- shared vertex normals across the now-shared seams ----
vn = np.empty(n_after * 3, dtype=np.float32)
me.vertices.foreach_get("normal", vn)
me.normals_split_custom_set_from_vertices(vn.reshape(-1, 3))
log("custom split normals reset from vertex normals")

# ---- re-measure the lines ----
co2 = np.empty(n_after * 3)
me.vertices.foreach_get("co", co2)
co2 = co2.reshape(-1, 3)
nr2 = np.empty(n_after * 3)
me.vertices.foreach_get("normal", nr2)
nr2 = nr2.reshape(-1, 3)
ev2 = np.empty(len(me.edges) * 2, dtype=np.int32)
me.edges.foreach_get("vertices", ev2)
ev2 = ev2.reshape(-1, 2)
o2 = np.concatenate([ev2[:, 0], ev2[:, 1]])
n2 = np.concatenate([ev2[:, 1], ev2[:, 0]])
s2 = np.argsort(o2, kind="stable")
o2s, n2s = o2[s2], n2[s2]
p2 = np.searchsorted(o2s, np.arange(n_after + 1))
c2 = np.maximum(np.diff(p2), 1)
N2 = nr2.copy()
for _ in range(5):
    acc = np.add.reduceat(N2[n2s], p2[:-1], axis=0)
    acc[np.diff(p2) == 0] = N2[np.diff(p2) == 0]
    N2 = acc / c2[:, None]
    N2 /= np.maximum(np.linalg.norm(N2, axis=1, keepdims=True), 1e-12)
ang2 = np.degrees(np.arccos(np.clip((nr2 * N2).sum(axis=1), -1, 1)))
torso2 = (co2[:, 2] > 0.28) & (co2[:, 2] < 0.90)
for thr in (8.0, 12.0, 20.0):
    print(f"KINK >{thr:4.1f}deg after weld: {int((torso2 & (ang2 > thr)).sum()):6d} verts")

bpy.context.preferences.filepaths.save_version = 0  # no .blend1 autosave
bpy.ops.wm.save_as_mainfile(filepath=OUT)
log(f"WROTE {OUT}")
print("WELD_DONE")