# Stage 17: heal the cut/scan lines — pinhole fill + thin-relief membrane. NO topology weld.
#
#   blender --background --python 17_line_heal.py -- <in.blend> <out.blend> [rounds]
#
# WHAT THE MEASUREMENTS ACTUALLY SAY (16c_topology.py, and they contradict stages 15 + 17-weld)
#   * The lines are NOT cracks. For kink verts the nearest vertex outside the 3-ring sits at
#     1.58x the local edge length (median); only 3.2% are closer than 0.35x, against 0.1% for
#     control skin. There is no disconnected-panel network to weld. Welding on that false premise
#     (eps 1.2 mm ~= the 1.85 mm mean edge length) merged ordinary neighbours and tore the mesh:
#     830 -> 24k boundary edges. Do not try it again.
#   * The mesh is nearly closed: 830 boundary edges in 689 loops of 3-5 edges each — scattered
#     PINHOLES, which is what the dotted black dashes in clay renders are. There is no sternum
#     hole; that gash is a sharp crease (see 18_cleavage.py), not an opening.
#   * The lines are thin GEOMETRIC relief: along kink verts |offset from the locally smooth
#     surface| is 0.22 mm median / 1.55 mm p95, versus 0.05 mm on quiet skin.
#
# So: fill the pinholes per-loop (never as one edge net — one net is what produced 39k
# non-manifold edges), then remove the thin relief with a collar-fixed bi-harmonic membrane.
#
# WHY A MEMBRANE AND NOT SMOOTHING. Taubin (stage 11) sheds high frequencies but PRESERVES low
# ones, and a 1.5 mm ridge riding on a curved torso is not purely high-frequency — that is why
# 30 Taubin pairs left the lines legible. A bi-harmonic membrane with two fixed collar rings
# matches position AND slope at the band edge, so the broad shape (underbust fold, hip curve,
# clavicle) is reproduced exactly while the thin ridge inside is replaced by the interpolant.
#
# WHY THE MASK IS SELF-DETECTED. Stages 11/14/15 all chose their masks by mapping the raw glb's
# boundary verts onto the working mesh; by then the sculpt had moved those verts past the
# tolerance, so most of the line network was never in the mask — the membranes were solving over
# the wrong verts, which is the real reason "the lines survived every pass". Detection here runs
# on the CURRENT geometry: very-small-scale relief (4-ring) AND a shading kink. Broad anatomy has
# low relief at that scale by construction, so it is excluded automatically rather than by hand.
import bpy, bmesh, sys, time
import numpy as np

argv = sys.argv[sys.argv.index("--") + 1:]
BLEND, OUT = argv[0], argv[1]
ROUNDS = int(argv[2]) if len(argv) > 2 else 3
t0 = time.time()

UNIT_MM = 1815.0            # 1 mesh unit = 1.815 m
DEV_MM = 0.20               # thin-relief threshold, real mm
KINK_DEG = 8.0
GROW = 2                    # band half-width, rings
COLLAR = 2                  # fixed rings outside the band (position + slope BC)
Z_LO, Z_HI = 0.04, 0.90     # below the chin: the face keeps its own detail
X_MAX = 0.36                # exclude hands/wrists
PIN_MAX_EDGES = 8           # a "pinhole" — bigger openings are left for a human to look at
CLAMP_MM = 2.5              # cap per-vertex displacement, real mm

# WIDTH IS SETTLED — DO NOT WIDEN THE BAND. 19_wide_heal.py measured the seam cross-section: the
# disturbance is one vertex wide (|offset| 0.23 mm median at the centre, already at the 0.074 mm
# background level by ring 2). Trials at GROW=4 and GROW=8 were no better than this narrow band
# (kink>12deg 6674 and 6751 vs 6994) while moving material up to 38 mm. Narrow is correct.
#
# CLAMP_MM exists because the solver is free to do something dramatic wherever a band happens to
# span a real feature rather than a seam — the unclamped run made a 17 mm excursion. Seam relief
# tops out around 1.5 mm, so 2.5 mm is generous for the defect and still forbids reshaping her.


def log(m):
    print(f"[line {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
log(f"in: {len(me.vertices)}v {len(me.polygons)}f")

# =============================================================================
# 1. pinhole fill — per loop
# =============================================================================
bm = bmesh.new()
bm.from_mesh(me)
open_e = [e for e in bm.edges if len(e.link_faces) == 1]
eset = {e.index for e in open_e}
emap = {e.index: e for e in open_e}
from collections import deque
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:
        cur = emap[q.popleft()]
        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)
sizes = sorted((len(lp) for lp in loops), reverse=True)
log(f"boundary loops: {len(loops)} (edges: {len(open_e)}, largest {sizes[:6]})")
bm.free()

# bmesh.ops.holes_fill silently refused every one of these loops (830 -> 839 boundary edges, no
# new faces), so use the edit-mode operator, which handles the small non-manifold fans these
# pinholes actually are.
bpy.context.view_layer.objects.active = ob
ob.select_set(True)
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='DESELECT')
bpy.ops.mesh.select_mode(type='EDGE')
bpy.ops.mesh.select_non_manifold(extend=False, use_boundary=True, use_wire=True,
                                 use_multi_face=False, use_non_contiguous=False, use_verts=False)
bpy.ops.mesh.fill_holes(sides=PIN_MAX_EDGES)
bpy.ops.object.mode_set(mode='OBJECT')
me.update()
bm = bmesh.new()
bm.from_mesh(me)
n_open_after = len([e for e in bm.edges if len(e.link_faces) == 1])
n_nm_after = 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')
    bm.to_mesh(me)
    me.update()
    log(f"deleted {len(dgn)} degenerate faces")
bm.free()
log(f"after fill: {len(me.vertices)}v {len(me.polygons)}f "
    f"boundary_edges={n_open_after} nonmanifold={n_nm_after}")

# =============================================================================
# 2. thin-relief membrane, re-detected each round
# =============================================================================
n_v = len(me.vertices)
co = np.empty(n_v * 3)
me.vertices.foreach_get("co", co)
co = co.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)
empty = np.diff(ptr) == 0


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


def smooth_n(X, k):
    Y = X.copy()
    for _ in range(k):
        Y = nbr_mean(Y)
    return Y


def grow(mask, rings):
    m = mask.copy()
    for _ in range(rings):
        hit = m[ev[:, 0]] | m[ev[:, 1]]
        m2 = m.copy()
        m2[ev[:, 0]] |= hit
        m2[ev[:, 1]] |= hit
        m = m2
    return m


zone = (co[:, 2] > Z_LO) & (co[:, 2] < Z_HI) & (np.abs(co[:, 0]) < X_MAX)
navel = (np.abs(co[:, 0]) < 0.022) & (co[:, 2] > 0.495) & (co[:, 2] < 0.555) & (co[:, 1] < 0)
protect = navel
log(f"zone {int(zone.sum())} verts, protected {int(protect.sum())}")


def detect(P):
    """Thin relief + shading kink, both measured on P."""
    nrm = np.empty(n_v * 3)
    me.vertices.foreach_get("normal", nrm)
    nrm = nrm.reshape(-1, 3)
    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)))
    dev = ((P - smooth_n(P, 4)) * N).sum(axis=1) * UNIT_MM       # real mm, very local
    m = zone & ~protect & (np.abs(dev) > DEV_MM) & (ang > KINK_DEG)
    return m, ang, dev


def bilaplacian_solve(P, free_m):
    """Collar-fixed bi-harmonic membrane over free_m, solved with CG on the S-subgraph."""
    collar = grow(free_m, COLLAR) & ~free_m
    S = np.nonzero(free_m | collar)[0]
    in_S = np.zeros(n_v, dtype=bool)
    in_S[S] = True
    glb = np.full(n_v, -1, dtype=np.int64)
    glb[S] = np.arange(len(S))
    se = ev[in_S[ev].all(axis=1)]
    a_ = glb[se[:, 0]]
    b_ = glb[se[:, 1]]
    deg = np.zeros(len(S))
    np.add.at(deg, a_, 1.0)
    np.add.at(deg, b_, 1.0)
    free = free_m[S]

    def Ls(X):
        out = deg[:, None] * X
        np.add.at(out, a_, -X[b_])
        np.add.at(out, b_, -X[a_])
        return out

    def A_op(U):
        X = np.zeros((len(S), 3))
        X[free] = U
        return Ls(Ls(X))[free]

    Xc = np.zeros((len(S), 3))
    Xc[~free] = P[S[~free]]
    rhs = -Ls(Ls(Xc))[free]
    U = P[S[free]].copy()
    r = rhs - A_op(U)
    p = r.copy()
    rs = (r * r).sum()
    rs0 = max(rs, 1e-30)
    it = 0
    for it in range(200000):
        Ap = A_op(p)
        den = (p * Ap).sum()
        if abs(den) < 1e-30:
            break
        al = rs / den
        U += al * p
        r -= al * Ap
        rs2 = (r * r).sum()
        if rs2 < 1e-20 or rs2 < rs0 * 1e-14:
            rs = rs2
            break
        p = r + (rs2 / rs) * p
        rs = rs2
    Q = P.copy()
    Q[S[free]] = U
    return Q, len(S), int(free.sum()), it, rs / rs0


P = co.copy()
for rnd in range(1, ROUNDS + 1):
    mask, ang, dev = detect(P)
    if not mask.sum():
        log(f"round {rnd}: nothing left to heal")
        break
    band = grow(mask, GROW) & zone & ~protect
    Q, ns, nf, it, rel = bilaplacian_solve(P, band)
    # clamp: a seam is <=1.5 mm of relief, so anything larger is the solver reshaping anatomy
    disp = Q - P
    dmag = np.linalg.norm(disp, axis=1) * UNIT_MM
    over = dmag > CLAMP_MM
    if over.any():
        disp[over] *= (CLAMP_MM / dmag[over])[:, None]
        Q = P + disp
        log(f"round {rnd}: clamped {int(over.sum())} verts to {CLAMP_MM} mm "
            f"(largest pre-clamp {dmag.max():.1f} mm)")
    d = np.linalg.norm(Q - P, axis=1) * UNIT_MM
    log(f"round {rnd}: seed {int(mask.sum())} -> band {nf} (S={ns})  CG it={it} rel={rel:.1e}  "
        f"moved max {d.max():.3f} mm, median(band) {np.median(d[band]):.3f} mm")
    P = Q
    me.vertices.foreach_set("co", P.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))
    m2, ang2, dev2 = detect(P)
    torso = (co[:, 2] > 0.28) & (co[:, 2] < 0.90)
    log(f"   after: seed mask {int(m2.sum())}, kink>12deg {int((torso & (ang2 > 12)).sum())}, "
        f"kink>20deg {int((torso & (ang2 > 20)).sum())}")

# final normals + save
me.vertices.foreach_set("co", P.reshape(-1))
me.update()
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))
d_tot = np.linalg.norm(P - co, axis=1) * UNIT_MM
log(f"TOTAL displacement: max {d_tot.max():.3f} mm, "
    f"{int((d_tot > 0.05).sum())} verts moved >0.05 mm")
bpy.context.preferences.filepaths.save_version = 0  # no .blend1 autosave
bpy.ops.wm.save_as_mainfile(filepath=OUT)
log(f"WROTE {OUT}")
print("LINE_DONE")