# Stage 6e: melt open-slit rims in the crotch box.
# The remaining flap and the dotted briefs-edge lines both live on boundary edges (the mesh
# is not watertight). Roughness misses them — a folded flap is locally smooth. Select verts
# on boundary edges inside the box, grow 3 rings, melt hard; also stitch: each boundary vert
# pairs with its nearest non-neighbour boundary vert within 2.5 mm and both move to the
# midpoint, closing the slit gap positionally (topology untouched, masks.npz stays valid).
#   blender --background --python 06e_slit_melt.py -- <in.blend> <out.blend>
import bpy, sys, time
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]
t0 = time.time()

BOX = lambda co: (np.abs(co[:, 0]) < 0.075) & (co[:, 2] > 0.340) & (co[:, 2] < 0.480)
MELT_ITERS = 100
STITCH_R = 0.0025


def log(m):
    print(f"[slit {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)

n_e = len(me.edges)
ev = np.empty(n_e * 2, dtype=np.int32)
me.edges.foreach_get("vertices", ev)
ev = ev.reshape(-1, 2)

# boundary edges: adjacent to exactly one face
l_tot = np.empty(len(me.polygons), dtype=np.int32)
me.polygons.foreach_get("loop_total", l_tot)
l_start = np.empty(len(me.polygons), dtype=np.int32)
me.polygons.foreach_get("loop_start", l_start)
l_v = np.empty(len(me.loops), dtype=np.int32)
me.loops.foreach_get("vertex_index", l_v)
ecount = {}
for fs, ft in zip(l_start, l_tot):
    idxs = l_v[fs:fs + ft]
    for k in range(ft):
        a, b = idxs[k], idxs[(k + 1) % ft]
        key = (a, b) if a < b else (b, a)
        ecount[key] = ecount.get(key, 0) + 1
bnd_v = np.zeros(n_v, dtype=bool)
for (a, b), c in ecount.items():
    if c == 1:
        bnd_v[a] = bnd_v[b] = True
log(f"boundary verts total: {bnd_v.sum()}")

box = BOX(co)
hot = box & bnd_v
log(f"slit verts in box: {hot.sum()}")

# stitch pass: pair each hot vert with nearest hot vert that is not a mesh neighbour
order = np.concatenate([ev[:, 0], ev[:, 1]])
nbr = np.concatenate([ev[:, 1], ev[:, 0]])
srt = np.argsort(order, kind="stable")
o_s = order[srt]
n_s = nbr[srt]
ptr = np.searchsorted(o_s, np.arange(n_v + 1))
cnt = np.maximum(ptr[1:] - ptr[:-1], 1)

Q = co.copy()
hidx = np.nonzero(hot)[0]
if len(hidx):
    tree = KDTree(len(hidx))
    for j, i in enumerate(hidx):
        tree.insert(Vector(Q[i]), j)
    tree.balance()
    neigh_sets = {int(i): set(int(x) for x in n_s[ptr[i]:ptr[i + 1]]) for i in hidx}
    stitched = 0
    done = set()
    for j, i in enumerate(hidx):
        if j in done:
            continue
        best = None
        for (_, k, dist) in tree.find_range(Vector(Q[i]), STITCH_R):
            if k == j or k in done:
                continue
            ik = int(hidx[k])
            if ik in neigh_sets[int(i)]:
                continue
            if best is None or dist < best[1]:
                best = (k, dist)
        if best is not None:
            ik = int(hidx[best[0]])
            mid = 0.5 * (Q[i] + Q[ik])
            Q[i] = mid
            Q[ik] = mid
            done.add(j)
            done.add(best[0])
            stitched += 1
    log(f"stitched {stitched} slit pairs")

    m = hot.copy()
    for _ in range(3):
        h = m[ev[:, 0]] | m[ev[:, 1]]
        m2 = m.copy()
        m2[ev[:, 0]] |= h
        m2[ev[:, 1]] |= h
        m = m2
    sidx = np.nonzero(m)[0]
    for _ in range(MELT_ITERS):
        acc = np.zeros_like(Q)
        np.add.at(acc, o_s, Q[n_s])
        mean = acc / cnt[:, None]
        Q[sidx] = 0.5 * Q[sidx] + 0.5 * mean[sidx]
    d = np.linalg.norm(Q - co, axis=1)
    log(f"melted {len(sidx)} rim verts, max move {d.max():.4f}")

me.vertices.foreach_set("co", Q.reshape(-1))
me.update()
if me.has_custom_normals:
    vn = np.empty(n_v * 3, dtype=np.float32)
    me.vertices.foreach_get("normal", vn)
    me.normals_split_custom_set_from_vertices(vn.reshape(-1, 3))
bpy.context.preferences.filepaths.save_version = 0  # no .blend1 autosave
bpy.ops.wm.save_as_mainfile(filepath=OUT)
log(f"WROTE {OUT}")
print("SLIT_DONE")