# Stage 50 (v02): body-only unwrap; the head KEEPS its original Tripo UVs.
#
# Chain of failures that led here, all measured on v02/48_headsafe.blend:
#   - MINIMUM_STRETCH over everything: the undecimated head carries Tripo's eyelash and brow
#     shells — thousands of needle slivers — and SLIM COLLAPSED those charts to points
#     (half of all loops at one UV coordinate, coverage 0.0%). v01 never saw this because
#     decimation had already crushed the lashes before any unwrap ran.
#   - ANGLE_BASED over everything: survived the slivers but packed at 17.3% coverage,
#     1.15 px/mm — half of v01's texel density, on the version whose whole point is quality.
#
# The escape: the head was deliberately never edited, so its ORIGINAL Tripo UV charts are
# still valid — and they are artist-grade (the face was the best-mapped region of the source
# atlas, ~2.6 px/mm). So: SLIM on the body only (proven at this density in v01),
# average-island-scale to equalise everything, boost the head charts by HEAD_BOOST for the
# quality this version exists for, then pack the lot together.
# Stage 49 (v02): 24_seams.py with ONE fix — a density-aware hip landmark.
#
# The v01 hip rule walked 0.01-thick midline slices and declared the crotch at the first slice
# holding <3 verts. On the head-protected v02 body (27.5k body verts vs 32k) that fired at
# u=0.610 — mid-belly. A hip ring above the real crotch turns the leg region into two joined
# tubes with too few cuts, and the minimum-stretch solver blows one chart up and packs the rest
# to nothing (coverage 0.0%). Slices are 2x thicker here and the threshold scales with the
# body's actual vertex density. Copied rather than edited: 24_seams.py is the recipe that
# regenerates the SHIPPED v01, and its behaviour must not drift.
# Stage 24: a proper human UV atlas — ~10 anatomical charts instead of Tripo's 5,870 blobs.
#
#   blender --background --python 24_seams.py -- <mesh.glb|blend> <out_dir> [decimate_ratio]
#
# Every seam is placed off a MEASURED landmark, not a guessed threshold, so the same rules port
# between the hires sculpt and the decimated game body:
#   neck / wrist / ankle  = local minimum of cross-section radius (the narrow part)
#   shoulder              = smallest |s| whose slice is short in u (arm, not torso)
#   hip                   = lowest slice that still has vertices on the midline (the crotch)
# Lengthwise cuts open each tube flat, hidden where nobody looks:
#   torso + head -> back midline;  arms -> back of the arm;  legs -> inner side.
# Hands / feet / head need no lengthwise cut: a tube with one closed end is already a disc.
import bpy, bmesh, sys, os, time
import numpy as np

argv = sys.argv[sys.argv.index("--") + 1:]
SRC = argv[0]
OUTDIR = os.path.abspath(argv[1])
RATIO = float(argv[2]) if len(argv) > 2 else 0.0
os.makedirs(OUTDIR, exist_ok=True)
t0 = time.time()


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


if SRC.lower().endswith(".glb"):
    bpy.ops.wm.read_homefile(use_empty=True)
    bpy.ops.import_scene.gltf(filepath=SRC)
else:
    bpy.ops.wm.open_mainfile(filepath=SRC)
ob = max([o for o in bpy.data.objects if o.type == 'MESH'], key=lambda o: len(o.data.vertices))
bpy.context.view_layer.objects.active = ob
for o in bpy.data.objects:
    o.select_set(o is ob)
bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
log(f"body '{ob.name}' {len(ob.data.vertices)}v {len(ob.data.polygons)}f")

if 0 < RATIO < 1:
    m = ob.modifiers.new("dec", 'DECIMATE')
    m.ratio = RATIO
    bpy.ops.object.modifier_apply(modifier=m.name)
    log(f"decimated -> {len(ob.data.vertices)}v {len(ob.data.polygons)}f")

# The source mesh is non-manifold (1,649 edges with >2 faces) and decimation turns that into a
# scatter of degenerate/duplicate triangles. Each one becomes its own UV island — that is where
# the confetti comes from, not the seam rules. Clean it here, before any of it reaches the unwrap.
def clean_mesh(tag):
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.mesh.remove_doubles(threshold=1e-5)
    bpy.ops.mesh.dissolve_degenerate(threshold=1e-6)
    bpy.ops.mesh.delete_loose(use_verts=True, use_edges=True, use_faces=False)
    # (Splitting the 1,649 non-manifold edges was tried here and rejected: it cost +3,512 verts
    # and made the collapsed-face count worse, because the collapse is a solver problem, not a
    # topology one — see the unwrap method below.)
    bpy.ops.mesh.normals_make_consistent(inside=False)
    bpy.ops.object.mode_set(mode='OBJECT')
    bm_ = bmesh.new(); bm_.from_mesh(ob.data)
    nm = sum(1 for e in bm_.edges if len(e.link_faces) > 2)
    dg = sum(1 for f in bm_.faces if f.calc_area() < 1e-12)
    bm_.free()
    log(f"clean[{tag}]: {len(ob.data.vertices)}v {len(ob.data.polygons)}f  "
        f"nonmanifold_edges={nm} degenerate_faces={dg}")


# SKIP_CLEAN=1: the transplanted head's pristine look depends on its authored custom normals,
# and clean_mesh's normals_make_consistent re-winds nonmanifold shells — exactly the lash/eye
# geometry — which is how the master's eyes got broken in the first place. The cost of skipping
# is only cosmetic confetti islands from the body's 1,649 nonmanifold edges.
if os.environ.get("SKIP_CLEAN") != "1":
    clean_mesh("after decimate")
else:
    log("clean_mesh SKIPPED (SKIP_CLEAN=1) — protecting authored head normals")

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)
lo, hi = co.min(axis=0), co.max(axis=0)
span = hi - lo
UP = int(np.argmax(span))
LR = int(np.argmax(np.where(np.arange(3) == UP, -1, span)))
FB = 3 - UP - LR
u = (co[:, UP] - lo[UP]) / span[UP]
s = co[:, LR] - 0.5 * (lo[LR] + hi[LR])
d = co[:, FB] - 0.5 * (lo[FB] + hi[FB])
HALF = 0.5 * span[LR]
sn = s / HALF                                    # left-right, normalised to +-1
print(f"axes up={'xyz'[UP]} lr={'xyz'[LR]} fb={'xyz'[FB]}   height {span[UP]:.4f} units")

# ---------------------------------------------------------------- landmarks
def radius_profile(mask, coord, lo_c, hi_c, nb):
    """mean in-slice radius vs coord, over `nb` bins — the narrow parts are the joints."""
    ed = np.linspace(lo_c, hi_c, nb + 1)
    mid = 0.5 * (ed[:-1] + ed[1:])
    out = np.full(nb, np.nan)
    for i in range(nb):
        m = mask & (coord >= ed[i]) & (coord < ed[i + 1])
        if m.sum() < 30:
            continue
        A = np.stack([s[m], d[m], u[m]], axis=1)
        A = np.delete(A, 0 if coord is sn else 2, axis=1) if False else A
        # radius measured in the two axes perpendicular to `coord`
        if coord is u:
            P = np.stack([s[m], d[m]], axis=1)
        else:
            P = np.stack([d[m], (u[m] - u[m].mean()) * span[UP]], axis=1)
        out[i] = np.linalg.norm(P - P.mean(axis=0), axis=1).mean()
    return mid, out


def local_min(mid, prof, lo_c, hi_c):
    m = (mid >= lo_c) & (mid <= hi_c) & np.isfinite(prof)
    if not m.any():
        return 0.5 * (lo_c + hi_c)
    return float(mid[m][np.argmin(prof[m])])


def taper_end(mid, prof, lo_c, hi_c, from_high, tol=0.08):
    """Where the limb stops tapering, taken from the EXTREMITY side.

    The global minimum of the radius profile is not the joint: on the leg it sits up the shin,
    well above where the foot ends. The joint is the end of the taper nearest the extremity —
    the first bin (walking in from that side) that reaches within `tol` of the minimum.
    """
    m = (mid >= lo_c) & (mid <= hi_c) & np.isfinite(prof)
    if not m.any():
        return 0.5 * (lo_c + hi_c)
    mm, pp = mid[m], prof[m]
    close = np.nonzero(pp <= pp.min() * (1.0 + tol))[0]
    return float(mm[close[-1] if from_high else close[0]])


torso_side = np.abs(sn) < 0.30
mid_u, prof_u = radius_profile(torso_side, u, 0.0, 1.0, 60)
NECK_U = local_min(mid_u, prof_u, 0.80, 0.93)
# The ankle must be measured on ONE leg. Over both, the "radius" is really the gap between
# them, which falls monotonically from the feet up and has no minimum to find.
one_leg = (sn > 0.05) & (u < 0.35)
mid_l, prof_l = radius_profile(one_leg, u, 0.0, 0.35, 35)
print("  single-leg radius profile (u -> radius), the ankle is the narrow point:")
print("    " + "  ".join(f"{m:.2f}:{r*1000:.0f}" for m, r in zip(mid_l, prof_l)
                         if np.isfinite(r)))
ANKLE_U = taper_end(mid_l, prof_l, 0.02, 0.14, from_high=False)

# shoulder: the smallest |s| whose slice is SHORT in u (an arm), scanning outward
absn = np.abs(sn)
ARM_IN = 0.35
for cut in np.arange(0.15, 0.60, 0.01):
    m = (absn >= cut) & (absn < cut + 0.03)
    if m.sum() < 30:
        continue
    if (u[m].max() - u[m].min()) < 0.13:
        ARM_IN = float(cut)
        break
mid_a, prof_a = radius_profile(absn > ARM_IN, absn, ARM_IN, 1.0, 40)
print("  arm radius profile (|s| -> radius), the wrist is the narrow point before the hand:")
print("    " + "  ".join(f"{m:.2f}:{r*1000:.0f}" for m, r in zip(mid_a, prof_a)
                         if np.isfinite(r)))
# The profile reads: radius falls to the wrist, bulges again over the palm, then tapers down
# the fingers. Search below the palm bulge or the "wrist" lands among the fingers.
WRIST = taper_end(mid_a, prof_a, 0.60, 0.80, from_high=True)

# hip: the crotch is the HIGHEST strictly-empty midline slice. Measured on this mesh: the
# |sn|<0.025 strip is zero below u=0.42 and populated above — except one freak slice at
# 0.50-0.52 holding a single vertex, which is what tripped both threshold-based rules (v01's
# "<3 verts" and the density-scaled version) into calling mid-belly the crotch. No stray
# vertex can fake a STRICTLY empty slice, and there is no surface between the legs to put
# one there.
HIP_U = 0.45
for lv in np.arange(0.60, 0.20, -0.005):
    m = (u >= lv) & (u < lv + 0.02) & (absn < 0.025)
    if m.sum() == 0:
        HIP_U = float(lv + 0.02)
        break
print(f"LANDMARKS  neck_u {NECK_U:.3f}  hip_u {HIP_U:.3f}  ankle_u {ANKLE_U:.3f}  "
      f"arm_in {ARM_IN:.3f}  wrist {WRIST:.3f}   (fractions of height / half-span)")

is_arm = absn > ARM_IN
is_hand = absn > WRIST
is_head = u > NECK_U
is_leg = (u < HIP_U) & ~is_arm
is_foot = u < ANKLE_U
# centre lines for the lengthwise cuts
# A constant centre works only for a perfectly axis-aligned limb. Hers droop, so a level set of
# u wanders off the arm and cuts it twice. Use a measured centre LINE: median u per |s| bin for
# the arms, median |s| per u bin for each leg, linearly interpolated.
def centre_line(mask, along, of, nb, lo_a, hi_a):
    ed = np.linspace(lo_a, hi_a, nb + 1)
    mid = 0.5 * (ed[:-1] + ed[1:])
    val = np.full(nb, np.nan)
    for i in range(nb):
        m = mask & (along >= ed[i]) & (along < ed[i + 1])
        if m.sum() >= 20:
            val[i] = np.median(of[m])
    ok = np.isfinite(val)
    if ok.sum() < 2:
        return lambda q: np.full_like(q, np.nanmedian(of[mask]) if mask.any() else 0.0)
    mid_o, val_o = mid[ok], val[ok]
    return lambda q: np.interp(np.asarray(q), mid_o, val_o)


arm_u_of = centre_line(is_arm & ~is_hand, absn, u, 24, ARM_IN, WRIST)
leg_s_of = {}
for sg in (-1, 1):
    m = is_leg & ~is_foot & (np.sign(sn) == sg)
    leg_s_of[sg] = centre_line(m, u, absn, 20, ANKLE_U, HIP_U)
# A foot cut at the ankle is an L (ankle-heel-toes) and a hand is a flat paddle; neither
# flattens as one chart. Split each along its silhouette, exactly where an artist would:
# the foot into upper/sole at its mid-height, the hand into back/palm at its mid-thickness.
U_SOLE = float(np.median(u[is_foot])) if is_foot.sum() else ANKLE_U * 0.5
D_PALM = float(np.median(d[is_hand])) if is_hand.sum() else 0.0
print(f"  foot split at u {U_SOLE:.3f} (upper|sole)   hand split at d {D_PALM:+.4f} (back|palm)")

qa = np.linspace(ARM_IN, WRIST, 5)
print(f"  arm centre line u at |s|={np.round(qa,2).tolist()}: "
      f"{np.round(arm_u_of(qa), 3).tolist()}")
ql = np.linspace(ANKLE_U, HIP_U, 5)
print(f"  leg(+) centre line |s| at u={np.round(ql,2).tolist()}: "
      f"{np.round(leg_s_of[1](ql), 3).tolist()}")

# The torso is a tube with FOUR holes (neck, two armholes, hip). The back midline joins neck to
# hip; two more cuts are needed or the chart stays multiply-connected and the unwrap stretches it
# badly. Run a relief cut across the BACK at armpit height, from each armhole in to the midline.
# upper body only — at |s| = ARM_IN the feet also splay out that far, and they win a p5
m_pit = (np.abs(absn - ARM_IN) < 0.03) & (u > 0.5)
U_PIT = float(np.percentile(u[m_pit], 5)) if m_pit.sum() > 30 else 0.62
print(f"  armpit u {U_PIT:.3f}  -> relief cut across the back at that height")

# ---------------------------------------------------------------- mark seams
bm = bmesh.new()
bm.from_mesh(me)
tally = {}


def hit(k):
    tally[k] = tally.get(k, 0) + 1
    return True


for e in bm.edges:
    e.seam = False
for e in bm.edges:
    a, b = e.verts[0].index, e.verts[1].index
    # rings, in order of priority
    if is_head[a] != is_head[b]:
        e.seam = hit("neck ring"); continue
    if is_hand[a] != is_hand[b]:
        e.seam = hit("wrist rings"); continue
    if is_foot[a] != is_foot[b]:
        e.seam = hit("ankle rings"); continue
    if is_arm[a] != is_arm[b]:
        e.seam = hit("shoulder rings"); continue
    if (u[a] < HIP_U) != (u[b] < HIP_U):
        e.seam = hit("hip ring"); continue
    # lengthwise cuts
    if is_foot[a] and is_foot[b]:
        if (u[a] > U_SOLE) != (u[b] > U_SOLE):
            e.seam = hit("foot sole line"); continue
    elif is_hand[a] and is_hand[b]:
        if (d[a] > D_PALM) != (d[b] > D_PALM):
            e.seam = hit("hand palm line"); continue
    elif is_arm[a] and is_arm[b] and not (is_hand[a] or is_hand[b]):
        ca, cb = arm_u_of(absn[a]), arm_u_of(absn[b])
        if d[a] > 0 and d[b] > 0 and (u[a] > ca) != (u[b] > cb):
            e.seam = hit("arm back line"); continue
    elif is_leg[a] and is_leg[b] and not (is_foot[a] or is_foot[b]):
        sg = int(1 if sn[a] + sn[b] >= 0 else -1)
        if absn[a] < leg_s_of[sg](u[a]) and absn[b] < leg_s_of[sg](u[b]) \
                and (d[a] > 0) != (d[b] > 0):
            e.seam = hit("leg inner line"); continue
    elif not (is_arm[a] or is_arm[b] or is_foot[a] or is_foot[b]):
        # torso and head share one continuous back midline
        if d[a] > 0 and d[b] > 0 and (sn[a] > 0) != (sn[b] > 0):
            e.seam = hit("back midline"); continue
        # armpit relief: back only, from each armhole inward to the midline
        if d[a] > 0 and d[b] > 0 and not is_head[a] and not is_head[b] \
                and (u[a] > U_PIT) != (u[b] > U_PIT):
            e.seam = hit("armpit relief"); continue
for k in sorted(tally):
    print(f"  seam '{k}': {tally[k]} edges")
bm.to_mesh(me)
bm.free()
log(f"marked {sum(tally.values())} seam edges")

# ---------------------------------------------------------------- unwrap
UNWRAP_METHOD = os.environ.get("UNWRAP_METHOD", "MINIMUM_STRETCH")


def do_unwrap(pack):
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    # Angle-based (ABF) is conformal: it preserves angles and is free to crush area. At her folds
    # and flaps that means whole patches land under one texel and come out untextured. The
    # minimum-stretch (SLIM) solver optimises area distortion instead, which is what a texture
    # transfer actually needs.
    try:
        bpy.ops.uv.unwrap(method=UNWRAP_METHOD, margin=0.0)
    except TypeError:
        log(f"unwrap method {UNWRAP_METHOD} unavailable — falling back to ANGLE_BASED")
        bpy.ops.uv.unwrap(method='ANGLE_BASED', margin=0.0)
    bpy.ops.uv.average_islands_scale()
    if pack:
        try:
            bpy.ops.uv.pack_islands(rotate=True, margin=0.003, scale=True)
        except TypeError:
            bpy.ops.uv.pack_islands(margin=0.003)
    bpy.ops.object.mode_set(mode='OBJECT')


def get_islands():
    m_ = ob.data
    n_l_, n_f_ = len(m_.loops), len(m_.polygons)
    lv = np.empty(n_l_, dtype=np.int32); m_.loops.foreach_get("vertex_index", lv)
    uv_ = np.empty(n_l_ * 2); m_.uv_layers.active.data.foreach_get("uv", uv_)
    uv_ = uv_.reshape(-1, 2)
    ls = np.empty(n_f_, dtype=np.int32); m_.polygons.foreach_get("loop_start", ls)
    lt = np.empty(n_f_, dtype=np.int32); m_.polygons.foreach_get("loop_total", lt)
    li_ = ls[lt == 3]
    Q = 1 << 20
    k = (lv.astype(np.int64) * Q * Q
         + np.round(np.clip(uv_[:, 0], 0, 1) * (Q - 1)).astype(np.int64) * Q
         + np.round(np.clip(uv_[:, 1], 0, 1) * (Q - 1)).astype(np.int64))
    _, uvv_ = np.unique(k, return_inverse=True)
    n_uvv_ = uvv_.max() + 1
    T_ = np.stack([uvv_[li_], uvv_[li_ + 1], uvv_[li_ + 2]], axis=1)
    par = np.arange(n_uvv_, dtype=np.int64)

    def find(x):
        r = x
        while par[r] != r:
            r = par[r]
        while par[x] != r:
            par[x], x = r, par[x]
        return r

    for a_, b_, c_ in T_:
        ra, rb, rc = find(a_), find(b_), find(c_)
        if ra != rb:
            par[rb] = ra
        if ra != rc:
            par[rc] = ra
    _, isl_ = np.unique(np.array([find(i) for i in range(n_uvv_)]), return_inverse=True)
    return isl_, isl_[T_[:, 0]], li_, lv, uv_, isl_.max() + 1


# Keep the source UVs on a second layer. Stage 31 needs them to resample the existing maps into
# the new layout — unwrapping over them in place would throw the textures away.
src_layer = ob.data.uv_layers.active
src_layer.name = "UVMap_tripo"
new_layer = ob.data.uv_layers.new(name="UVMap_atlas", do_init=True)
ob.data.uv_layers.active = new_layer
for i, l in enumerate(ob.data.uv_layers):
    if l is new_layer:
        ob.data.uv_layers.active_index = i
log(f"uv layers: {[l.name for l in ob.data.uv_layers]}  active="
    f"{ob.data.uv_layers.active.name}")

HEAD_BOOST = float(os.environ.get("HEAD_BOOST", "1.5"))
# EXACTLY the seam threshold: with any offset, faces between the neck seam and the selection
# line end up in islands that are half re-unwrapped and half Tripo — incoherent charts.
NECK_FACE = NECK_U

# UVMap_atlas was created with do_init=True, i.e. a copy of the Tripo layer — so faces we do
# NOT unwrap keep their original mapping. Select body faces only, by every-corner test.
import bmesh as _bm
bpy.ops.object.mode_set(mode='EDIT')
# FACE select mode, or per-face .select is flushed away by vertex-mode sync and the unwrap
# silently runs on a leftover fraction of the body (measured: 9,757 of ~56k faces).
bpy.ops.mesh.select_mode(type='FACE')
bpy.ops.mesh.select_all(action='DESELECT')
bmm = _bm.from_edit_mesh(ob.data)
bmm.verts.ensure_lookup_table()
zs = [v.co.z for v in bmm.verts]
zlo_, zhi_ = min(zs), max(zs)
vhead = [((v.co.z - zlo_) / (zhi_ - zlo_)) > NECK_FACE for v in bmm.verts]
nbody = 0
for f in bmm.faces:
    sel = not any(vhead[v.index] for v in f.verts)
    f.select = sel
    nbody += int(sel)
_bm.update_edit_mesh(ob.data)
log(f"body-only unwrap: {nbody} faces selected, head keeps Tripo UVs")
try:
    bpy.ops.uv.unwrap(method='MINIMUM_STRETCH', margin=0.0)
except TypeError:
    bpy.ops.uv.unwrap(method='ANGLE_BASED', margin=0.0)
bpy.ops.mesh.select_all(action='SELECT')
bpy.ops.uv.average_islands_scale()
bpy.ops.object.mode_set(mode='OBJECT')

# boost the head charts before packing: scale each majority-head island about its centroid
isl, fisl, li, loops_v, uv, n_isl = get_islands()
me2 = ob.data
co2 = np.empty(len(me2.vertices) * 3); me2.vertices.foreach_get("co", co2); co2 = co2.reshape(-1, 3)
u2 = (co2[:, 2] - co2[:, 2].min()) / (co2[:, 2].max() - co2[:, 2].min())
tri_head = u2[loops_v[li]] > NECK_FACE
head_isl = set()
for i in range(n_isl):
    m = fisl == i
    if m.sum() and tri_head[m].mean() > 0.5:
        head_isl.add(i)
uvl2 = me2.uv_layers["UVMap_atlas"]
buf = np.empty(len(me2.loops) * 2); uvl2.data.foreach_get("uv", buf)
buf = buf.reshape(-1, 2)
# island id per loop: quantized (vertex,uv) identity -> island, same construction as get_islands
loop_isl = np.full(len(me2.loops), -1, dtype=np.int64)
loop_isl[li] = fisl
loop_isl[li + 1] = fisl
loop_isl[li + 2] = fisl
boosted = 0
for i in head_isl:
    m = loop_isl == i
    c = buf[m].mean(axis=0)
    buf[m] = c + (buf[m] - c) * HEAD_BOOST
    boosted += 1
uvl2.data.foreach_set("uv", buf.reshape(-1))
log(f"boosted {boosted} head islands x{HEAD_BOOST}")

bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='SELECT')
try:
    bpy.ops.uv.pack_islands(rotate=True, margin=0.003, scale=True)
except TypeError:
    bpy.ops.uv.pack_islands(margin=0.003)
bpy.ops.object.mode_set(mode='OBJECT')
isl, fisl, li, loops_v, uv, n_isl = get_islands()
log("unwrapped + packed (body SLIM, head Tripo)")

tuv = np.stack([np.clip(uv, 0, 1)[li], np.clip(uv, 0, 1)[li + 1], np.clip(uv, 0, 1)[li + 2]], 1)
auv = 0.5 * np.abs((tuv[:, 1, 0] - tuv[:, 0, 0]) * (tuv[:, 2, 1] - tuv[:, 0, 1])
                   - (tuv[:, 2, 0] - tuv[:, 0, 0]) * (tuv[:, 1, 1] - tuv[:, 0, 1]))
me = ob.data
co = np.empty(len(me.vertices) * 3); me.vertices.foreach_get("co", co); co = co.reshape(-1, 3)
P3 = co[np.stack([loops_v[li], loops_v[li + 1], loops_v[li + 2]], axis=1)]
a3 = 0.5 * np.linalg.norm(np.cross(P3[:, 1] - P3[:, 0], P3[:, 2] - P3[:, 0]), axis=1)
isl_auv = np.bincount(fisl, weights=auv, minlength=n_isl)
isl_a3 = np.bincount(fisl, weights=a3, minlength=n_isl)
isl_nf = np.bincount(fisl, minlength=n_isl)
UNITM = 1.777 / span[UP]
W = 4096
dens = np.where(isl_a3 > 0, np.sqrt(np.maximum(isl_auv, 0) / np.maximum(isl_a3, 1e-12))
                * W / (UNITM * 1000.0), np.nan)
o = np.argsort(-isl_auv)
cum = np.cumsum(isl_auv[o]) / max(isl_auv.sum(), 1e-12)
n99 = int(np.searchsorted(cum, 0.99)) + 1
print(f"\n=== NEW ATLAS  (Tripo baseline in brackets) ===")
print(f"islands {n_isl}  [5870]      99%-of-area islands {n99}  [84]")
print(f"coverage {isl_auv.sum()*100:.1f}%  [62.0]      confetti <20 faces {int((isl_nf<20).sum())}  [5763]")
db = dens[o[:n99]]; db = db[np.isfinite(db)]
print(f"texel density {db.min():.2f}..{db.max():.2f} px/mm -> {db.max()/max(db.min(),1e-9):.2f}x "
      f"spread  [1.8x]")
# distortion, measured only over the real charts — the confetti islands are degenerate
# triangles whose ratios are meaningless and would own the tail
real = np.isin(fisl, o[:n99])
ok = real & (a3 > 1e-12) & (auv > 1e-14)
sc = np.sqrt(auv[ok] / a3[ok]); sc /= np.median(sc)
print(f"per-triangle area-scale vs median, real charts only: p05 {np.percentile(sc,5):.2f}  "
      f"p95 {np.percentile(sc,95):.2f}  p99 {np.percentile(sc,99):.2f}   (1.0 = undistorted)")
print("\nthe charts:      #    faces   uv_area%    3D cm2  stretch_p95   where")
for i in o[:16]:
    if isl_auv[i] * 100 < 0.05:
        break
    m = fisl == i
    mo = m & (a3 > 1e-12) & (auv > 1e-14)
    st = np.sqrt(auv[mo] / a3[mo])
    st = st / np.median(st) if len(st) else np.array([1.0])
    vs = np.unique(np.stack([loops_v[li], loops_v[li + 1], loops_v[li + 2]], 1)[m].ravel())
    lab = []
    for nm, msk in (("head", is_head), ("hand", is_hand), ("foot", is_foot),
                    ("arm", is_arm & ~is_hand), ("leg", is_leg & ~is_foot)):
        if msk[vs].mean() > 0.6:
            lab.append(nm)
    side = "L" if sn[vs].mean() > 0.15 else ("R" if sn[vs].mean() < -0.15 else "mid")
    print(f"  {i:6d} {isl_nf[i]:8d} {isl_auv[i]*100:9.3f} {isl_a3[i]*UNITM*UNITM*1e4:9.1f} "
          f"{np.percentile(st,95):11.2f}   {'+'.join(lab) or 'torso'} {side}")

# picture
R = 1024
rng = np.random.RandomState(3)
pal = rng.rand(n_isl, 3) * 0.7 + 0.25
IS = np.zeros((R, R, 3))
tp = tuv * (R - 1)
for fi in range(len(tp)):
    P = tp[fi]
    x0, x1 = int(P[:, 0].min()), int(np.ceil(P[:, 0].max()))
    y0, y1 = int(P[:, 1].min()), int(np.ceil(P[:, 1].max()))
    if x1 < x0 or y1 < y0 or x1 - x0 > 64 or y1 - y0 > 64:
        continue
    dt = ((P[1, 1] - P[2, 1]) * (P[0, 0] - P[2, 0]) + (P[2, 0] - P[1, 0]) * (P[0, 1] - P[2, 1]))
    if abs(dt) < 1e-12:
        continue
    gx, gy = np.meshgrid(np.arange(x0, min(x1, R - 1) + 1), np.arange(y0, min(y1, R - 1) + 1))
    aa = ((P[1, 1] - P[2, 1]) * (gx - P[2, 0]) + (P[2, 0] - P[1, 0]) * (gy - P[2, 1])) / dt
    bb = ((P[2, 1] - P[0, 1]) * (gx - P[2, 0]) + (P[0, 0] - P[2, 0]) * (gy - P[2, 1])) / dt
    ins = (aa >= 0) & (bb >= 0) & (1 - aa - bb >= 0)
    if ins.any():
        IS[gy[ins], gx[ins]] = pal[fisl[fi]]
img = bpy.data.images.new("isl", R, R, alpha=False)
A = np.ones((R, R, 4), dtype=np.float32); A[:, :, :3] = IS
img.pixels.foreach_set(A.reshape(-1))
p = os.path.join(OUTDIR, "charts.png")
img.file_format = 'PNG'; img.filepath_raw = p; img.save(filepath=p)
bpy.ops.wm.save_as_mainfile(filepath=os.path.join(OUTDIR, "seamed.blend"))
log(f"wrote {p} + seamed.blend")
print("SEAMS_DONE")