# lena_leafbikini lane, stage 07: BARBIE-FILL the crotch — melt the briefs leaves into a
# smooth featureless surface. The bust holes stay cut open (stage 05 behaviour).
#
#   blender --background --factory-startup --python 07_fill_crotch.py -- \
#           <pristine.glb> <leaf_mask.npz> <out.glb> [--z-split 0.615] [--free-rings 2]
#           [--puff-mm 0.0] [--blend <out.blend>]
#
# TARGET. "Barbie doll anatomy": a completely smooth, undifferentiated pelvic surface — no
# cleft, no features, a taut convex continuation of belly into inner thighs. That is exactly
# what a bi-harmonic membrane produces: solve L²x = 0 with the surrounding skin held fixed.
# The rim supplies POSITION (the "distance between the two sides"), the collar behind it
# supplies SLOPE through the second Laplacian application (the "angle"), and a bi-harmonic
# surface cannot invent detail — no crease, no cleft, by construction.
#
# WHY MELT, NOT FILL. This stage was first written as hole-filling on the stage-05 cut, and it
# failed twice, instructively:
#
#   attempt 1 — triangle_fill each rim + densify + membrane. The briefs rim is ONE ~2,500-vert
#     loop snaking front -> between the legs -> back; beauty triangulation of a loop that long
#     and that non-convex connects the WRONG BANKS at every bend. The membrane then faithfully
#     smooths garbage into rippled sheets. Bonus failure: the rim itself still carried leaf-root
#     remnants, and a membrane anchors position AND slope to its rim, so it reproduced the
#     crumple (the mesh-repair playbook's "collar on CLEAN skin" lesson, re-learned).
#   attempt 2 — rim erosion (3 rings) + interleaved Delaunay flips + double solve. Better rims,
#     same disease: flips are local, the mis-bridging is global. Long sliver strands shot off
#     the hips where chords bridged front rim to back rim, and triangle_fill did not even close
#     the pinched loop (2,177 faces where ~2,471 were needed; 819 boundary edges left).
#
#   The fix is to stop inventing topology. THE LEAF SHELL IS the disk that spans the hole — the
#   scan's own manifold surface, connected to the true rim at every point, no bank ever bridged
#   wrongly. So in the crotch band the leaves are not deleted at all: their vertices are FREED,
#   a --free-rings collar of surrounding skin is freed with them (this erases the under-leaf rim
#   crease, the same move as the nude lane's fair_rim_band), and the membrane collapses the
#   whole shell onto the smooth spanning surface. Folded flaps famously resist melting by
#   ITERATIVE flow (the 06*-era lesson) — but L²x = 0 is linear with a unique solution, and PCG
#   run to convergence lands on it regardless of where the folds start.
#
# THE REMNANT SWEEP. Near the melt zone the colour key is re-run without the stage-04 mask's
# blind spots (min-comp speckle filter, value gate): any fixed vertex within 6 rings of the
# primary free set that is greenish (hue 46..200, no gates) or blown-white (val >= 0.78,
# sat <= 0.25 — nothing on her actual skin is that colour) joins the melt, grown 2 rings.
# Freeing a few honest skin verts by accident is harmless — the membrane returns them almost
# in place; a pinned leaf fragment is not.
#
# THE BALLOON EXCISION. The melt's first run still left half a dozen smooth raised nubs, and a
# debug bake proved they were FREE verts — melted, converged, and still bulging. That is not a
# solver bug, it is what L²x = 0 does to a PENDANT BALLOON: a fully-masked leaf is a closed
# shell attached along its root line, its excess surface area has nowhere to go, and the
# bi-harmonic solution — smooth in GRAPH terms, with no maximum principle — parks it as a
# rounded mound. Two remedies failed before this one worked:
#   * proudness detection (60-sweep, then 400-sweep Taubin reference): a smoothing-built
#     reference partially FOLLOWS any bump wider than its radius, so visibly 4 mm nubs measured
#     1.7 mm and thresholds caught only their tips;
#   * Laplacian deflation of what it did catch: flattened tips, kept the wave.
# What a wad cannot hide is its AREA: several layers of surface over one spot put several times
# the vertices of honest membrane into the same cell of a 4 mm grid. So: detect wads by vertex
# density, EXCISE every face touching one, and refill the scars with triangle_fill + densify +
# a small membrane solve — which is exactly the right tool at this scale (round holes a
# centimetre or two across; its failure mode was only ever the giant winding channel).
# CAVEAT the run exposed: the whole melted shell lies ~4 layers deep (cell median 37), so
# density above median finds the FLAT piles — worth excising, they would z-fight — but the
# visible INFLATED caps sit at ~2 layers, BELOW median. Density cannot see them either.
#
# THE FLOW FINISH — the step that actually guarantees a smooth result, with no detector at
# all: damped pure-Laplacian flow over the entire changed region, fixed skin held, weight
# ramping 0 -> 1 over the first 8 rings so the membrane's C1 rim blend survives. The maximum
# principle does what every detector could not promise: a raised cap has strictly nowhere to
# go but down, while the broad pubic web barely moves (flow erases features at ~1/size², and
# the web is 5-10x wider than any cap). A short Taubin polish follows to undo the slight
# overall shrink. This ordering — bi-harmonic for shape, flow for guarantees — is the recipe.
#
# TEXTURE. Melted faces keep leaf texels, so every face whose vertices are all masked gets ONE
# donor texel — an old thigh-band vertex whose albedo is closest to the band's median skin tone
# and whose normal-map texel is nearest neutral. Interpolating UVs instead is meaningless here:
# the atlas is Tripo chart soup (the sibling Lena mesh had 5,870 charts). Flat is correct —
# Barbie plastic has no albedo detail either. The baked leaf contact shadows still darken the
# surviving skin just outside the melt; that is an albedo problem for a later stage.
#
# TOPOLOGY GROUND RULES (as stages 04/05): the glTF importer splits every UV seam, so the mesh
# is welded (exact duplicates, 1e-5) before boundaries or adjacency mean anything — UVs live
# per face corner and survive the weld. Custom split normals do not survive the bmesh round
# trip; on the WELDED mesh "clear + shade smooth" is seamless, which is the second reason the
# weld comes first.
import bpy, bmesh, sys, os, time, argparse
from collections import deque
import numpy as np


def interior_edges(faces):
    """Edges whose every adjacent face is a patch face — the only ones safe to subdivide."""
    return list({e for f in faces for e in f.edges
                 if all(lf in faces for lf in e.link_faces)})


def refresh(faces, *rets):
    """Re-collect live patch faces after a bmesh op invalidated / created some."""
    out = {f for f in faces if f.is_valid}
    for ret in rets:
        for key in ("geom", "geom_inner", "faces"):
            for g in ret.get(key, ()):
                if isinstance(g, bmesh.types.BMFace) and g.is_valid:
                    out.add(g)
    return out

argv = sys.argv[sys.argv.index("--") + 1:] if "--" in sys.argv else []
ap = argparse.ArgumentParser()
ap.add_argument("glb")
ap.add_argument("mask")
ap.add_argument("out")
ap.add_argument("--z-split", type=float, default=0.615,
                help="fraction of body height separating briefs (melted) from bust (cut open); "
                     "briefs mask tops out at 0.594, bust starts at 0.631")
ap.add_argument("--free-rings", type=int, default=2,
                help="rings of surrounding skin freed with the leaves, to erase the rim crease")
ap.add_argument("--puff-mm", type=float, default=0.0,
                help="optional outward dome on top of the membrane, peak amplitude in mm")
ap.add_argument("--blend", default="")
A = ap.parse_args(argv)
GLB, MASK, OUT = os.path.abspath(A.glb), os.path.abspath(A.mask), os.path.abspath(A.out)
os.makedirs(os.path.dirname(OUT), exist_ok=True)
t0 = time.time()


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


def smoothstep(x):
    x = np.clip(x, 0.0, 1.0)
    return x * x * (3.0 - 2.0 * x)


# =============================================================================================
# load original + mask, split the mask at the waist
# =============================================================================================
bpy.ops.wm.read_factory_settings(use_empty=True)
bpy.ops.import_scene.gltf(filepath=GLB)
body = max([o for o in bpy.data.objects if o.type == 'MESH'], key=lambda o: len(o.data.vertices))
me = body.data
bpy.context.view_layer.objects.active = body
body.select_set(True)
n0 = len(me.vertices)
log(f"in : '{body.name}' {n0}v {len(me.polygons)}f")

z = np.load(MASK)
inv, mk = z["inv"].astype(np.int64), z["mask"]
if len(inv) != n0:
    raise SystemExit(f"[melt] FATAL: mask was built for {len(inv)} verts, this GLB has {n0}")
vm = mk[inv]                                   # per raw vertex: is it leaf?

co = np.empty(n0 * 3); me.vertices.foreach_get("co", co); P0 = co.reshape(-1, 3)
Z0, H = float(P0[:, 2].min()), float(P0[:, 2].max() - P0[:, 2].min())
MM = 1000.0 * 1.777 / H
zf0 = (P0[:, 2] - Z0) / H
crotch = vm & (zf0 <= A.z_split)
bust = vm & (zf0 > A.z_split)
log(f"mask: {vm.sum()} leaf verts -> {crotch.sum()} briefs (melt), {bust.sum()} bust (cut)")

# per-vertex albedo HSV for the remnant sweep — sampled NOW, on the raw import, because UV
# indexing goes stale the moment bmesh touches the topology
base_img = None
for mat in [m_ for m_ in me.materials if m_]:
    bsdf = next((x for x in mat.node_tree.nodes if x.type == 'BSDF_PRINCIPLED'), None)
    lnk = bsdf and bsdf.inputs["Base Color"].links
    if lnk:
        nd = lnk[0].from_node
        while nd.type != 'TEX_IMAGE':
            up = [i for i in nd.inputs if i.links]
            if not up:
                break
            nd = up[0].links[0].from_node
        if nd.type == 'TEX_IMAGE':
            base_img = nd.image
if base_img is None:
    raise SystemExit("[melt] FATAL: no base-colour image")
nl0 = len(me.loops)
lv0 = np.empty(nl0, dtype=np.int32); me.loops.foreach_get("vertex_index", lv0)
uv0 = np.empty(nl0 * 2); me.uv_layers.active.data.foreach_get("uv", uv0); uv0 = uv0.reshape(-1, 2)
vuv0 = np.zeros((n0, 2)); vuv0[lv0[::-1]] = uv0[::-1]
w_, h_ = base_img.size
buf = np.empty(w_ * h_ * 4, dtype=np.float32); base_img.pixels.foreach_get(buf)
px = buf.reshape(h_, w_, 4)[:, :, :3]; del buf
xi = np.clip((vuv0[:, 0] * (w_ - 1)).astype(np.int64), 0, w_ - 1)
yi = np.clip((vuv0[:, 1] * (h_ - 1)).astype(np.int64), 0, h_ - 1)
C = px[yi, xi].astype(np.float64); del px
S = np.clip(np.where(C <= 0.0031308, C * 12.92, 1.055 * np.maximum(C, 0) ** (1 / 2.4) - 0.055), 0, 1)
R, G, B = S[:, 0], S[:, 1], S[:, 2]
mx = S.max(1); mn = S.min(1); dd = mx - mn
hue = np.zeros(n0)
nz = dd > 1e-6
im = np.argmax(S, axis=1)
sel = nz & (im == 0); hue[sel] = 60 * (((G[sel] - B[sel]) / dd[sel]) % 6)
sel = nz & (im == 1); hue[sel] = 60 * ((B[sel] - R[sel]) / dd[sel] + 2)
sel = nz & (im == 2); hue[sel] = 60 * ((R[sel] - G[sel]) / dd[sel] + 4)
sat = np.where(mx > 1e-6, dd / np.maximum(mx, 1e-6), 0.0)

# ride the mask + HSV through the weld as attributes (per-vertex custom data survives
# remove_doubles on the surviving vertex of each duplicate cluster)
for name, arr in (("melt_m", crotch), ("bust_m", bust)):
    at = me.attributes.new(name=name, type='INT', domain='POINT')
    at.data.foreach_set("value", arr.astype(np.int32))
for name, arr in (("hsv_h", hue), ("hsv_s", sat), ("hsv_v", mx)):
    at = me.attributes.new(name=name, type='FLOAT', domain='POINT')
    at.data.foreach_set("value", arr.astype(np.float32))

# =============================================================================================
# weld, cut the bust open, free the briefs
# =============================================================================================
bm = bmesh.new()
bm.from_mesh(me)
bmesh.ops.remove_doubles(bm, verts=list(bm.verts), dist=1e-5)
bm.verts.ensure_lookup_table()
lm = bm.verts.layers.int["melt_m"]
lb = bm.verts.layers.int["bust_m"]
log(f"welded: {len(bm.verts)}v")

# bust: delete fully-masked faces (stage 05 cut rule), then despike the new rim
kill = [f for f in bm.faces if all(v[lb] for v in f.verts)]
bmesh.ops.delete(bm, geom=kill, context='FACES')
log(f"bust cut: -{len(kill)} faces")
for it in range(4):
    spikes = [f for f in bm.faces
              if sum(1 for e in f.edges if len(e.link_faces) == 1) >= 2
              and (sum(v.co.z for v in f.verts) / len(f.verts) - Z0) / H > A.z_split]
    if not spikes:
        break
    bmesh.ops.delete(bm, geom=spikes, context='FACES')
    log(f"bust despike pass {it+1}: -{len(spikes)} dangling faces")
loose = [v for v in bm.verts if not v.link_faces]
if loose:
    bmesh.ops.delete(bm, geom=loose, context='VERTS')

# briefs: free = leaf verts + a skin collar, grown over true (welded) adjacency
free_set = {v for v in bm.verts if v[lm]}
for _ in range(A.free_rings):
    free_set |= {o for v in free_set for e in v.link_edges for o in e.verts}
log(f"melt set: {len(free_set)} free verts (leaves + {A.free_rings}-ring skin collar)")

# the remnant sweep (see header): re-key the fixed verts near the melt without the mask's
# speckle filter or value gate, so missed leaf fragments melt too instead of pinning welts
lh = bm.verts.layers.float["hsv_h"]
lsat = bm.verts.layers.float["hsv_s"]
lval = bm.verts.layers.float["hsv_v"]
near = set(free_set)
for _ in range(6):
    near |= {o for v in near for e in v.link_edges for o in e.verts}
adds = {v for v in near - free_set
        if (v.co.z - Z0) / H <= A.z_split + 0.01
        and ((46.0 <= v[lh] <= 200.0) or (v[lval] >= 0.78 and v[lsat] <= 0.25))}
grown = set(adds)
for _ in range(2):
    grown |= {o for v in grown for e in v.link_edges for o in e.verts}
free_set |= grown
log(f"remnant sweep: +{len(adds)} keyed (+{len(grown - adds)} ring growth) "
    f"-> {len(free_set)} free verts")

# mark free verts and melted-face texels via flags that survive to_mesh
for v in bm.verts:
    v.select_set(False)
for v in free_set:
    v.select_set(True)
texel = {v for v in bm.verts if v[lm]} | grown   # leaf faces + swept remnants, not the collar
for f in bm.faces:
    f.select_set(all(v in texel for v in f.verts))
bm.to_mesh(me)
bm.free()
me.update()
n = len(me.vertices)
log(f"topology done: {n}v {len(me.polygons)}f")

# =============================================================================================
# the melt: matrix-free PCG on L²x = 0, everything but the briefs held fixed
# =============================================================================================
co = np.empty(n * 3); me.vertices.foreach_get("co", co); P = co.reshape(-1, 3).copy()
P_orig = P.copy()                               # for the fixed-verts gate at the end
vsel = np.empty(n, dtype=bool); me.vertices.foreach_get("select", vsel)
free = vsel.copy()
log(f"free verts: {free.sum()}")

ev = np.empty(len(me.edges) * 2, dtype=np.int32); me.edges.foreach_get("vertices", ev)
ea, eb = ev[0::2].astype(np.int64), ev[1::2].astype(np.int64)
zfw = (P[:, 2] - Z0) / H                        # fixed verts never move, so this stays valid


def grow_np(mask, rings):
    out = mask.copy()
    for _ in range(rings):
        hit = np.zeros(n, dtype=bool)
        m = out[ea] | out[eb]
        hit[ea[m]] = True
        hit[eb[m]] = True
        out |= hit
    return out


# full-mesh adjacency in CSR form, built once — the solver restricts it per pass
fsrc = np.concatenate([ea, eb]); fdst = np.concatenate([eb, ea])
fo = np.argsort(fsrc, kind='stable'); fsrc, fdst = fsrc[fo], fdst[fo]
fdeg = np.bincount(fsrc, minlength=n).astype(np.float64)
fptr = np.concatenate([[0], np.cumsum(fdeg)]).astype(np.int64)


def solve_membrane(free_mask, tag, warm_harmonic):
    """Bi-harmonic solve for the free verts; reads and writes me's positions in place.
    Returns (ridx, Fl, deg, ptr, dst) so the puff step can reuse the last pass's graph."""
    co_ = np.empty(n * 3); me.vertices.foreach_get("co", co_); Pv = co_.reshape(-1, 3).copy()
    region = grow_np(free_mask, 3)              # ring1 enters L, ring2 enters L², ring3 margin
    ridx = np.nonzero(region)[0]
    loc = np.full(n, -1, dtype=np.int64)
    loc[ridx] = np.arange(len(ridx))
    m = region[ea] & region[eb]
    ra, rb = loc[ea[m]], loc[eb[m]]
    src = np.concatenate([ra, rb]); dst = np.concatenate([rb, ra])
    o = np.argsort(src, kind='stable'); src, dst = src[o], dst[o]
    deg = np.bincount(src, minlength=len(ridx)).astype(np.float64)
    ptr = np.concatenate([[0], np.cumsum(deg)]).astype(np.int64)
    Fl = free_mask[ridx]
    log(f"solve[{tag}]: {len(ridx)} region verts ({Fl.sum()} free), {len(src)//2} edges")

    def Lap(Xv):
        s = np.add.reduceat(Xv[dst], ptr[:-1], axis=0)
        s[deg == 0] = 0.0
        return s - deg[:, None] * Xv

    X = Pv[ridx].copy()

    def T_free(XF):
        Y = np.zeros_like(X)
        Y[Fl] = XF
        return Lap(Lap(Y))[Fl]

    if warm_harmonic:
        # Lx = 0 converges in ~diameter sweeps and lands within a crease of the bi-harmonic
        # answer, cutting the expensive solve's iterations roughly in half. Only worth it when
        # starting from the raw leaf shell — a re-solve already sits near the answer.
        Y = X.copy()
        for it in range(3000):
            d = Lap(Y)
            Y[Fl] += 0.9 / np.maximum(deg[Fl], 1.0)[:, None] * d[Fl]
            if it % 500 == 499 and float(np.abs(d[Fl]).max()) * MM < 1e-4:
                break
        X[Fl] = Y[Fl]
        log(f"  harmonic warm start: {it+1} sweeps")

    r = -Lap(Lap(X))[Fl]
    Mjac = (deg[Fl] ** 2 + deg[Fl])[:, None]    # diag(L²) = deg² + deg
    zv = r / Mjac
    p = zv.copy()
    rz = float((r * zv).sum())
    b0 = float(np.linalg.norm(-Lap(Lap(np.where(Fl[:, None], 0.0, X)))[Fl])) + 1e-30
    xF = X[Fl].copy()
    rn = float(np.linalg.norm(r))
    for it in range(20000):
        Ap = T_free(p)
        alpha = rz / (float((p * Ap).sum()) + 1e-300)
        xF += alpha * p
        r -= alpha * Ap
        rn = float(np.linalg.norm(r))
        if rn / b0 < 3e-7:
            break
        zv = r / Mjac
        rz2 = float((r * zv).sum())
        p = zv + (rz2 / rz) * p
        rz = rz2
        if it % 2000 == 1999:
            log(f"  PCG iter {it+1}: residual {rn/b0:.2e}")
    X[Fl] = xF
    log(f"solve[{tag}]: {it+1} iters (residual {rn/b0:.2e}), "
        f"max move {np.linalg.norm(X[Fl]-Pv[ridx][Fl],axis=1).max()*MM:.1f} mm")
    Pv[ridx] = X
    me.vertices.foreach_set("co", Pv.ravel())
    me.update()
    return ridx, Fl, deg, ptr, dst


ridx, Fl, deg, ptr, dst = solve_membrane(free, "melt", warm_harmonic=True)


# THE SOLID-WELT PASS. Some leaf roots are not shells at all — they are SOLID ridges sculpted
# into the body surface and painted in skin tones, which is why every shell-hunting detector
# (density, occlusion, fin normals) returned almost nothing while four caps sat in plain view.
# A solid bump is honest single surface, so the ORIGINAL remedy is the right one: free it and
# let the membrane pull it down — no excess area, no pendant balloon. Detector exactly as the
# debug probe validated it: 400-sweep Taubin reference over the WHOLE band (everything moves,
# so there is no anchored-strip chord and no frozen-zone blindness), proud along the normal
# > 0.8 mm, fixed verts only. At that reference the smooth melt web reads ~0.2 mm (p95) and
# the caps read 1-2.4 mm.
def full_nbmean(Xv):
    s = np.add.reduceat(Xv[fdst], fptr[:-1], axis=0)
    s[fdeg == 0] = Xv[fdeg == 0]
    return s / np.maximum(fdeg, 1.0)[:, None]


mvW = (zfw >= 0.40) & (zfw <= A.z_split + 0.02)
for wpass in range(2):                          # re-reference and re-detect: the first fix
    cow = np.empty(n * 3)                       # exposes whatever its 3-ring growth missed
    me.vertices.foreach_get("co", cow)
    Pw = cow.reshape(-1, 3).copy()
    Qw = Pw.copy()
    for _ in range(400):
        Qw[mvW] += 0.50 * (full_nbmean(Qw) - Qw)[mvW]
        Qw[mvW] += -0.53 * (full_nbmean(Qw) - Qw)[mvW]
    nrw = np.empty(n * 3); me.vertices.foreach_get("normal", nrw)
    proudW = np.einsum('ij,ij->i', Pw - Qw, nrw.reshape(-1, 3)) * MM
    # 0.5 mm, not 0.8: the caps are ~1.8 mm PLATEAUS with sharp edges (the hard shadows in
    # clay renders oversell their height), and 0.8 clipped 220 crown verts while the body of
    # each plateau survived. Honest skin reads p95 +0.13 mm against this reference — 0.5 mm
    # is still 4x above the noise floor.
    weltS = ~free & mvW & (proudW > 0.5)
    log(f"solid-welt pass {wpass+1}: {weltS.sum()} proud fixed verts "
        f"(band fixed p95 {np.percentile(proudW[~free & mvW], 95):+.2f} mm, "
        f"max {proudW[~free & mvW].max():+.2f} mm)")
    if not weltS.sum():
        break
    weltG = grow_np(weltS, 3) & ~free
    free |= weltG
    solve_membrane(weltG, f"weltfix{wpass+1}", warm_harmonic=False)
    # their paint is leaf-root shadow, not skin: hand their faces to the donor texel too
    fselW = np.empty(len(me.polygons), dtype=bool); me.polygons.foreach_get("select", fselW)
    lvW = np.empty(len(me.loops), dtype=np.int32); me.loops.foreach_get("vertex_index", lvW)
    lsW = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_start", lsW)
    ltW = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_total", ltW)
    allin = np.add.reduceat(weltG[lvW].astype(np.int32), lsW.astype(np.int64)) == ltW
    me.polygons.foreach_set("select", fselW | allin)
    log(f"solid-welt pass {wpass+1}: freed {weltG.sum()} verts, "
        f"{int(allin.sum())} faces to donor texel")

# THE BALLOON DEFLATION (see header). Wide reference: 400 Taubin sweeps over the melt zone —
# diffusion radius ~sqrt(400) = 20 rings (~26 mm here), wide enough that a centimetre nub reads
# as fully proud instead of being absorbed into its own reference.
def full_nbmean(Xv):
    s = np.add.reduceat(Xv[fdst], fptr[:-1], axis=0)
    s[fdeg == 0] = Xv[fdeg == 0]
    return s / np.maximum(fdeg, 1.0)[:, None]


co2 = np.empty(n * 3); me.vertices.foreach_get("co", co2); Pm = co2.reshape(-1, 3).copy()
# the reference zone must extend PAST anything the detectors are asked to judge: outside `mv`
# the smoothed copy equals the input and proudness is identically zero by construction — a
# 4-ring halo silently blinded the welt detector to caps sitting 5+ rings out
mv = grow_np(free, 30) & (zfw <= A.z_split + 0.03)
# Detection is by DENSITY, not proudness: a smoothing-built reference partially follows any
# bump wider than its radius (a 400-sweep probe read the visibly 4 mm nubs at 1.7 mm), but a
# wad cannot hide its area — multiple layers over one spot of surface put several times the
# verts of honest membrane into the same cell of a 3D grid.
CELL = 4.0 / MM
key3 = np.floor(Pm[free] / CELL).astype(np.int64)
_, cinv, ccnt = np.unique(key3, axis=0, return_inverse=True, return_counts=True)
per_vert_cnt = ccnt[cinv]
med_cnt = float(np.median(per_vert_cnt))
hot_thr = max(3.0 * med_cnt, 18.0)
balloon = np.zeros(n, dtype=bool)
balloon[np.nonzero(free)[0][per_vert_cnt > hot_thr]] = True
balloon = grow_np(balloon, 1)
log(f"balloon pass: cell median {med_cnt:.0f} verts, threshold {hot_thr:.0f} -> "
    f"{balloon.sum()} wad verts to excise")

# THE FOLD DETECTOR — for the rim-attached flaps density cannot see (~2 layers, BELOW the
# piled median) and the flow finish cannot reach (they live in the rim-damped zone). A folded
# flap betrays itself by its NORMALS: its flanks and underside disagree with the smoothed
# reference field by 80-180 degrees, which honest skin never does — even the walls of a deep
# concave crease stay within ~70 degrees of a 150-sweep reference. Restricted to FREE verts,
# so fixed anatomy can never be excised no matter how it folds.
ltF = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_total", ltF)
if (ltF == 3).all():
    lvF = np.empty(len(me.loops), dtype=np.int32); me.loops.foreach_get("vertex_index", lvF)
    lsF = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_start", lsF)

    def vnormals(Pts):
        va, vb, vc = lvF[lsF], lvF[lsF + 1], lvF[lsF + 2]
        fn = np.cross(Pts[vb] - Pts[va], Pts[vc] - Pts[va])
        acc = np.zeros_like(Pts)
        for idx in (va, vb, vc):
            np.add.at(acc, idx, fn)
        return acc / np.maximum(np.linalg.norm(acc, axis=1, keepdims=True), 1e-12)

    Qf = Pm.copy()
    for _ in range(400):
        Qf[mv] += 0.50 * (full_nbmean(Qf) - Qf)[mv]
        Qf[mv] += -0.53 * (full_nbmean(Qf) - Qf)[mv]
    ncur = vnormals(Pm)
    dotn = np.einsum('ij,ij->i', ncur, vnormals(Qf))
    fold = free & (dotn < 0.15)
    log(f"fold pass: {fold.sum()} inverted-normal verts "
        f"(free dot p05 {np.percentile(dotn[free], 5):+.2f})")
    balloon |= grow_np(fold, 1)

    # THE FIXED-WELT PASS — by RAY-CAST OCCLUSION, the playbook's own move, after every
    # reference-surface detector failed for a structural reason worth recording:
    #   * proudness vs a Taubin reference: Taubin is shape-PRESERVING — a 1.5 cm cap sits
    #     inside its passband, so the reference reproduces the cap and P-Q reads ~0 forever;
    #   * proudness at all: a FIN's wall normals are perpendicular to its height, the dot
    #     is ~0 no matter how far it sticks out;
    #   * diffusion references: run wide they chord across convex anatomy (+1.9 mm on honest
    #     hips), run narrow they cannot see the cap tops standing 15+ rings out.
    # A flap needs no reference: it stands OVER surface, so a short ray cast INWARD from it
    # hits geometry within millimetres — its own opposite wall (fins are 1-2 mm thick) or the
    # web below — while honest skin's inward ray travels centimetres of flesh before exiting.
    # The gluteal crease is safe by construction: its walls' inward rays point into the flesh,
    # AWAY from each other; only a +n ray could cross the crease gap, and none is cast.
    import mathutils
    deps = bpy.context.evaluated_depsgraph_get()
    bvh = mathutils.bvhtree.BVHTree.FromObject(body, deps)
    nearF = grow_np(free, 20) & ~free & (zfw >= 0.40) & (zfw <= A.z_split + 0.01)
    EPS, DMAX = 0.4 / MM, 4.5 / MM
    widx = []
    for vi in np.nonzero(nearF)[0]:
        p, nv = Pm[vi], ncur[vi]
        o = mathutils.Vector((p[0] - nv[0] * EPS, p[1] - nv[1] * EPS, p[2] - nv[2] * EPS))
        if bvh.ray_cast(o, mathutils.Vector((-nv[0], -nv[1], -nv[2])), DMAX)[0] is not None:
            widx.append(int(vi))
    welt = np.zeros(n, dtype=bool)
    welt[widx] = True
    log(f"fixed-welt pass: {welt.sum()} occluded (flap) verts of {nearF.sum()} candidates")
    balloon |= grow_np(welt, 1)
else:
    log("fold pass: SKIPPED (non-triangle faces present)")

if balloon.sum():
    # excise every face touching a wad vert (kills whole balloons, leaves no orphan shells),
    # tidy the scar, and refill: at this scale — round holes a centimetre or two across —
    # triangle_fill is exactly the right tool; its failure mode was the giant winding channel
    bm = bmesh.new(); bm.from_mesh(me)
    bm.verts.ensure_lookup_table()
    bidx = set(np.nonzero(balloon)[0].tolist())
    kill = [f for f in bm.faces if any(v.index in bidx for v in f.verts)]
    bmesh.ops.delete(bm, geom=kill, context='FACES')
    for it in range(4):
        spikes = [f for f in bm.faces
                  if sum(1 for e in f.edges if len(e.link_faces) == 1) >= 2
                  and (sum(v.co.z for v in f.verts) / len(f.verts) - Z0) / H <= A.z_split]
        if not spikes:
            break
        bmesh.ops.delete(bm, geom=spikes, context='FACES')
    loose = [v for v in bm.verts if not v.link_faces]
    if loose:
        bmesh.ops.delete(bm, geom=loose, context='VERTS')
    log(f"excision: -{len(kill)} wad faces (+{len(loose)} loose verts)")

    # zip the pre-existing Tripo slits where they meet the scars (nude lane, close_rim_slits):
    # a scar boundary that runs into a slit is a RIBBON, not a closed loop, and both fill
    # operators refuse it — these were the 10 holes that survived three sweeps untouched
    sl = [v for v in bm.verts
          if Z0 + 0.40 * H <= v.co.z <= Z0 + (A.z_split + 0.01) * H
          and any(len(e.link_faces) == 1 for e in v.link_edges)]
    v0 = len(bm.verts)
    bmesh.ops.remove_doubles(bm, verts=sl, dist=0.8 / MM)
    bm.verts.ensure_lookup_table()
    log(f"slit weld: {len(sl)} band boundary verts, {v0 - len(bm.verts)} merged at 0.8 mm")
    # the weld leaves zero-area faces and zero-length edges; feeding those to triangle_fill
    # took Blender down with an access violation, not an exception — clean them first
    dg = [e for e in bm.edges
          if Z0 + 0.38 * H <= (e.verts[0].co.z + e.verts[1].co.z) / 2 <= Z0 + 0.64 * H]
    bmesh.ops.dissolve_degenerate(bm, dist=1e-5, edges=dg)
    bm.verts.ensure_lookup_table()
    log(f"degenerate dissolve: {len(bm.verts)}v {len(bm.faces)}f")

    pre2 = {(round(v.co.x, 6), round(v.co.y, 6), round(v.co.z, 6)) for v in bm.verts}
    patch_faces = set()
    for sweep in range(3):                      # re-detect after filling: one triangle_fill
        bedges2 = [e for e in bm.edges          # failing silently must not leave a pinhole
                   if len(e.link_faces) == 1
                   and Z0 + 0.42 * H <= (e.verts[0].co.z + e.verts[1].co.z) / 2
                   <= Z0 + A.z_split * H]
        v2b2 = {}
        for e in bedges2:
            for v in e.verts:
                v2b2.setdefault(v, []).append(e)
        seen2, scars = set(), []
        for e0 in bedges2:
            if e0 in seen2:
                continue
            comp, q = [], deque([e0])
            seen2.add(e0)
            while q:
                e = q.popleft()
                comp.append(e)
                for v in e.verts:
                    for e2 in v2b2[v]:
                        if e2 not in seen2:
                            seen2.add(e2)
                            q.append(e2)
            if len(comp) >= 3:
                scars.append(comp)
        if not scars:
            break
        log(f"scar sweep {sweep+1}: {len(scars)} holes "
            f"(sizes {sorted(len(c) for c in scars)[::-1][:10]})")
        sweep_faces = set()                     # densify THIS sweep's fills only — letting a
        for comp in scars:                      # later sweep's sliver scars re-densify earlier
            if len(comp) > 900:                 # patches to their microscopic target once
                log(f"  REFUSING {len(comp)}-edge boundary tangle (fill would crash/garble)")
                continue
            live = [e for e in comp if e.is_valid]
            try:                                # (3.3M-face / 17 h lesson)
                ret = bmesh.ops.triangle_fill(bm, use_beauty=True, use_dissolve=False,
                                              edges=live)
                newf = [g for g in ret["geom"] if isinstance(g, bmesh.types.BMFace)]
            except Exception:
                newf = []
            if not newf:
                try:
                    ret = bmesh.ops.holes_fill(bm, edges=live, sides=0)
                    newf = list(ret["faces"])
                except Exception:
                    newf = []
            sweep_faces.update(newf)
        tgt = max(1.5 * float(np.median([e.calc_length() for c in scars for e in c])),
                  2.4 / MM)                     # floored at 2.4 mm: sliver rims must not set it
        for it in range(5):
            sweep_faces = {f for f in sweep_faces if f.is_valid}
            if len(sweep_faces) > 120000:
                log(f"  densify CAPPED at {len(sweep_faces)} faces")
                break
            longe = [e for e in interior_edges(sweep_faces) if e.calc_length() > 1.45 * tgt]
            if not longe:
                break
            r1 = bmesh.ops.subdivide_edges(bm, edges=longe, cuts=1, use_grid_fill=True)
            sweep_faces = refresh(sweep_faces, r1)
            r2 = bmesh.ops.triangulate(bm, faces=list(sweep_faces))
            sweep_faces = refresh(set(), r2)
            r3 = bmesh.ops.beautify_fill(bm, faces=list(sweep_faces),
                                         edges=interior_edges(sweep_faces))
            sweep_faces = refresh(sweep_faces, r3)
        patch_faces = {f for f in patch_faces if f.is_valid} | sweep_faces
    log(f"scar fill: {len(patch_faces)} patch faces")
    for f in patch_faces:                       # patches join the donor-texel set
        if f.is_valid:
            f.select_set(True)
    bm.to_mesh(me)
    bm.free()
    me.update()

    # topology changed: rebuild the globals the solver reads, then relax the patches
    n = len(me.vertices)
    ev = np.empty(len(me.edges) * 2, dtype=np.int32); me.edges.foreach_get("vertices", ev)
    ea, eb = ev[0::2].astype(np.int64), ev[1::2].astype(np.int64)
    co2 = np.empty(n * 3); me.vertices.foreach_get("co", co2)
    Pn = co2.reshape(-1, 3)
    patch_free = np.array([tuple(k) not in pre2 for k in np.round(Pn, 6)], dtype=bool)
    log(f"patch verts: {patch_free.sum()}")
    if patch_free.sum():
        solve_membrane(patch_free, "patch", warm_harmonic=False)

co3 = np.empty(n * 3); me.vertices.foreach_get("co", co3); Pout = co3.reshape(-1, 3).copy()
# every vertex is now either at a pristine post-cut position or it is part of the melt/patch;
# `changed` is the union of moved and newly created — the gate + donor selection key off it
orig_keys = {tuple(k) for k in np.round(P_orig, 6)}
changed = np.array([tuple(k) not in orig_keys for k in np.round(Pout, 6)], dtype=bool)
free = changed

# THE FLOW FINISH. Whatever pendant caps survived every detector above die here, and nothing
# has to find them first: damped pure-Laplacian flow over the ENTIRE changed region, fixed skin
# held. Laplacian flow obeys the maximum principle — no point can move outside the hull of its
# neighbours — so a raised cap has strictly nowhere to go but down, while the broad pubic web
# barely moves (flow erases features at a rate ~1/size², and the web is 5-10x wider than any
# cap). The weight ramps from 0 at the fixed rim to 1 by ring 8, so the bi-harmonic C1 blend
# earned by the membrane is untouched where it matters.
fsrc = np.concatenate([ea, eb]); fdst = np.concatenate([eb, ea])
fo = np.argsort(fsrc, kind='stable'); fsrc, fdst = fsrc[fo], fdst[fo]
fdeg = np.bincount(fsrc, minlength=n).astype(np.float64)
fptr = np.concatenate([[0], np.cumsum(fdeg)]).astype(np.int64)


def nbmean2(Xv):
    s = np.add.reduceat(Xv[fdst], fptr[:-1], axis=0)
    s[fdeg == 0] = Xv[fdeg == 0]
    return s / np.maximum(fdeg, 1.0)[:, None]


depth = np.zeros(n)
reach = ~free
d = 0
while not reach.all() and d < 200:
    d += 1
    nxt = reach.copy()
    hit = np.zeros(n, dtype=bool)
    m2 = reach[fsrc]
    hit[fdst[m2]] = True
    nxt |= hit
    ring = nxt & ~reach
    if not ring.any():
        break
    depth[ring] = d
    reach = nxt
# ramp over 3 rings, NOT 8. A pixel-ray probe finally identified the last "caps" as the melt
# web itself bridging taut over the inguinal hollow — and the hip-vine channel is only
# ~10-20 rings wide, so an 8-ring ramp kept essentially the whole strip in the damped zone
# and the flow never engaged exactly where the bridge needed pulling down. Three rings still
# protects the immediate C1 blend; everything past it flows.
w = smoothstep(depth / 3.0)[:, None]
log(f"flow finish: max depth {int(depth.max())} rings")
for _ in range(300):
    Pout[free] += (0.55 * w[free]) * (nbmean2(Pout) - Pout)[free]
for _ in range(8):                              # Taubin polish: undo the slight flow shrink
    Pout[free] += (0.55 * w[free]) * (nbmean2(Pout) - Pout)[free]
    Pout[free] += (-0.58 * w[free]) * (nbmean2(Pout) - Pout)[free]
moved_fin = np.linalg.norm(Pout[free] - co3.reshape(-1, 3)[free], axis=1)
log(f"flow finish: moved p50 {np.percentile(moved_fin,50)*MM:.2f} mm, "
    f"max {moved_fin.max()*MM:.1f} mm")
me.vertices.foreach_set("co", Pout.ravel())
me.update()

# optional Barbie dome: outward along the membrane normal, smoothstep of rim distance,
# zero value AND zero slope at the rim so the C1 blend survives
if A.puff_mm > 0:
    dist = np.zeros(len(ridx))
    unv = set(np.nonzero(Fl)[0].tolist())
    cur = set(np.nonzero(~Fl)[0].tolist())
    d = 0
    while unv and cur:
        d += 1
        nxt = set()
        for c in cur:
            for j in range(int(ptr[c]), int(ptr[c + 1])):
                nb = int(dst[j])
                if nb in unv:
                    unv.discard(nb)
                    dist[nb] = d
                    nxt.add(nb)
        cur = nxt
    t = dist / max(dist.max(), 1.0)
    nrm = np.empty(n * 3); me.vertices.foreach_get("normal", nrm); nrm = nrm.reshape(-1, 3)
    Pout[ridx] += nrm[ridx] * (smoothstep(t) * (A.puff_mm / MM))[:, None] * Fl[:, None]
    me.vertices.foreach_set("co", Pout.ravel())
    me.update()
    log(f"puff: +{A.puff_mm} mm dome over {int(dist.max())} rings")

# =============================================================================================
# texture the melt: one clean donor texel on the leaf-texel faces
# =============================================================================================
imgs = {}
for mat in [m_ for m_ in me.materials if m_]:
    for nd in mat.node_tree.nodes:
        if nd.type == 'TEX_IMAGE' and nd.image:
            for out_ in nd.outputs:
                for lnk in out_.links:
                    if lnk.to_socket.name == 'Base Color':
                        imgs['base'] = nd.image
            if 'normal' in nd.image.name.lower():
                imgs['normal'] = nd.image


def sample(img, uvs):
    w, h = img.size
    buf = np.empty(w * h * 4, dtype=np.float32)
    img.pixels.foreach_get(buf)
    px = buf.reshape(h, w, 4)[:, :, :3]
    xi = np.clip((uvs[:, 0] * (w - 1)).astype(np.int64), 0, w - 1)
    yi = np.clip((uvs[:, 1] * (h - 1)).astype(np.int64), 0, h - 1)
    out = px[yi, xi].copy()
    del buf, px
    return out


nl = len(me.loops)
lv = np.empty(nl, dtype=np.int32); me.loops.foreach_get("vertex_index", lv)
uvb = np.empty(nl * 2); me.uv_layers.active.data.foreach_get("uv", uvb); uvb = uvb.reshape(-1, 2)
vuv = np.zeros((n, 2)); vuv[lv[::-1]] = uvb[::-1]

zf = (Pout[:, 2] - Z0) / H
cand = np.nonzero(~free & (zf > 0.30) & (zf < 0.42))[0][::37]      # thigh band, thinned
cb = sample(imgs['base'], vuv[cand])
med = np.median(cb, axis=0)
score = np.linalg.norm(cb - med, axis=1)
if 'normal' in imgs:
    cn = sample(imgs['normal'], vuv[cand])
    score += 2.0 * np.linalg.norm(cn - np.array([0.5, 0.5, 1.0]), axis=1)
best = int(np.argmin(score))
donor = cand[best]
log(f"donor texel: vert {donor} zf={zf[donor]:.3f} albedo={np.round(cb[best],3)} "
    f"(band median {np.round(med,3)})")

fsel = np.empty(len(me.polygons), dtype=bool); me.polygons.foreach_get("select", fsel)
ls = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_start", ls)
lt = np.empty(len(me.polygons), dtype=np.int32); me.polygons.foreach_get("loop_total", lt)
duv = vuv[donor]
touched = 0
for fi in np.nonzero(fsel)[0]:
    for li in range(ls[fi], ls[fi] + lt[fi]):
        uvb[li] = duv
        touched += 1
me.uv_layers.active.data.foreach_set("uv", uvb.ravel())
log(f"UVs: {touched} loops on {int(fsel.sum())} melted faces -> donor texel")

# =============================================================================================
# normals, gates, export
# =============================================================================================
if me.has_custom_normals:
    bpy.ops.mesh.customdata_custom_splitnormals_clear()
me.polygons.foreach_set("use_smooth", np.ones(len(me.polygons), dtype=bool))
me.update()

# the excision/refill renumbers vertices, so "fixed didn't move" is asserted by position:
# every vertex is either bit-identical to a pristine post-cut position, or it is changed —
# and everything changed must live inside the crotch band
bm = bmesh.new(); bm.from_mesh(me)
band_open = sum(1 for e in bm.edges if len(e.link_faces) == 1
                and (0.5 * (e.verts[0].co.z + e.verts[1].co.z) - Z0) / H <= A.z_split)
bm.free()
mz = zf[changed]
# +0.04, not +0.01: the solid-welt pass detects up to z_split+0.02 and grows 3 rings, so its
# legitimate reach is a little above the split — the gate must allow what the recipe declares
in_band_ok = bool((mz.min() >= 0.42) and (mz.max() <= A.z_split + 0.04))
log(f"gate changed geometry confined to band: z {mz.min():.3f}..{mz.max():.3f} "
    f"({'PASS' if in_band_ok else 'FAIL'})")
log(f"gate crotch band boundary edges: {band_open} (pre-existing Tripo slits only)")
if not in_band_ok:
    raise SystemExit("[melt] FATAL: geometry changed outside the crotch band")

if A.blend:
    bpy.ops.wm.save_as_mainfile(filepath=os.path.abspath(A.blend))
    log(f"WROTE {A.blend}")
for o in bpy.data.objects:
    o.select_set(True)
bpy.ops.export_scene.gltf(filepath=OUT, export_format='GLB', use_selection=True,
                          export_yup=True, export_skins=False, export_animations=False,
                          export_apply=False, export_image_format='AUTO',
                          export_tangents=False, export_normals=True)
log(f"WROTE {OUT} ({os.path.getsize(OUT)/1e6:.2f} MB)")