animation/characters/work/lena/22_unwrap_test.py

# Stage 22 (feasibility test, writes only into its out dir): can we replace Tripo's 5,870-chart
# soup with a proper human atlas — a dozen anatomical charts, seams hidden where a character
# artist would put them (back midline, inner arm, inner leg, wrist/ankle/neck rings)?
#
#   blender --background --python 22_unwrap_test.py -- <body.glb|blend> <out_dir> [decimate_ratio]
#
# Reports the same metrics as 20_atlas_probe so old and new atlas are directly comparable:
# island count, atlas coverage, texel-density spread, and per-chart stretch.
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"[uw {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)
log(f"body '{ob.name}' {len(ob.data.vertices)}v {len(ob.data.polygons)}f")

# glTF import leaves the object rotated/parented; work in world space
bpy.ops.object.transform_apply(location=True, rotation=True, scale=True)
me = ob.data
co = np.empty(len(me.vertices) * 3); me.vertices.foreach_get("co", co); co = co.reshape(-1, 3)
lo, hi = co.min(axis=0), co.max(axis=0)
print(f"BBOX x {lo[0]:.3f}..{hi[0]:.3f}  y {lo[1]:.3f}..{hi[1]:.3f}  z {lo[2]:.3f}..{hi[2]:.3f}")
span = hi - lo
UP = int(np.argmax(span))
LR = int(np.argmax(np.where(np.arange(3) == UP, -1, span)))
FB = 3 - UP - LR
print(f"axes: up={'xyz'[UP]} span {span[UP]:.3f} | left-right={'xyz'[LR]} span {span[LR]:.3f} "
      f"| front-back={'xyz'[FB]} span {span[FB]:.3f}")

if RATIO > 0 and 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")
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)

# normalise to a body frame: u = up in 0..1, s = signed left-right, d = signed front-back
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]

# where the limbs are, as fractions of the left-right span
print(f"|s| percentiles: " + " ".join(f"p{p}={np.percentile(np.abs(s),p)/HALF:.2f}"
                                      for p in (50, 80, 90, 95, 99)))
print(f"u percentiles:  " + " ".join(f"p{p}={np.percentile(u,p):.2f}" for p in (5, 25, 50, 75, 95)))

# ---- seam rules, all in body-frame fractions so they port between densities ----
ARM_IN = 0.34 * HALF        # shoulder: |s| beyond this is arm
WRIST = 0.80 * HALF
NECK_U = 0.855              # above this is head
ANKLE_U = 0.055
CROTCH_U = 0.46

bm = bmesh.new()
bm.from_mesh(me)
bm.verts.ensure_lookup_table()
V = np.array([v.co for v in bm.verts])
uu = (V[:, UP] - lo[UP]) / span[UP]
ss = V[:, LR] - 0.5 * (lo[LR] + hi[LR])
dd = V[:, FB] - 0.5 * (lo[FB] + hi[FB])
is_arm = np.abs(ss) > ARM_IN
is_head = uu > NECK_U
is_leg = uu < CROTCH_U

for e in bm.edges:
    e.seam = False
n_seam = 0
for e in bm.edges:
    a, b = e.verts[0].index, e.verts[1].index
    # 1. rings: neck, wrists, ankles, crotch band -> separate head / hands / feet charts
    if (uu[a] > NECK_U) != (uu[b] > NECK_U):
        e.seam = True; n_seam += 1; continue
    if (np.abs(ss[a]) > WRIST) != (np.abs(ss[b]) > WRIST):
        e.seam = True; n_seam += 1; continue
    if (uu[a] > ANKLE_U) != (uu[b] > ANKLE_U):
        e.seam = True; n_seam += 1; continue
    # 2. shoulder ring: torso | arm
    if (np.abs(ss[a]) > ARM_IN) != (np.abs(ss[b]) > ARM_IN):
        e.seam = True; n_seam += 1; continue
    # 3. lengthwise cut so each tube can open flat:
    #    torso/head -> back midline; arms -> underside; legs -> inner side
    if is_arm[a] and is_arm[b]:
        if (dd[a] > 0) != (dd[b] > 0):          # back of the arm
            e.seam = True; n_seam += 1; continue
    elif is_leg[a] and is_leg[b]:
        sgn = 1.0 if ss[a] + ss[b] >= 0 else -1.0
        if ((ss[a] * sgn) < (ss[b] * sgn)) and dd[a] > 0 and dd[b] > 0:
            pass                                 # handled by the midline test below
        if (dd[a] > 0) != (dd[b] > 0) and np.abs(ss[a]) < 0.30 * HALF:
            e.seam = True; n_seam += 1; continue
    else:
        if dd[a] > 0 and dd[b] > 0 and (ss[a] > 0) != (ss[b] > 0):
            e.seam = True; n_seam += 1; continue
log(f"marked {n_seam} seam edges")
bm.to_mesh(me)
bm.free()

bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='SELECT')
bpy.ops.uv.unwrap(method='ANGLE_BASED', margin=0.002)
log("unwrapped")
try:
    bpy.ops.uv.pack_islands(rotate=True, margin=0.004, scale=True)
except TypeError:
    bpy.ops.uv.pack_islands(margin=0.004)
log("packed")
bpy.ops.object.mode_set(mode='OBJECT')

# ---- metrics (same definitions as 20_atlas_probe) ----
me = ob.data
n_l = len(me.loops)
loops_v = np.empty(n_l, dtype=np.int32); me.loops.foreach_get("vertex_index", loops_v)
uv = np.empty(n_l * 2); me.uv_layers.active.data.foreach_get("uv", uv); uv = uv.reshape(-1, 2)
n_f = len(me.polygons)
l_start = np.empty(n_f, dtype=np.int32); me.polygons.foreach_get("loop_start", l_start)
l_tot = np.empty(n_f, dtype=np.int32); me.polygons.foreach_get("loop_total", l_tot)
tri = l_tot == 3
li = l_start[tri]

Q = 1 << 20
key = (loops_v.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(key, return_inverse=True)
n_uvv = uvv.max() + 1
T = np.stack([uvv[li], uvv[li + 1], uvv[li + 2]], axis=1)
parent = np.arange(n_uvv, dtype=np.int64)


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


for a, b, c in T:
    ra, rb, rc = find(a), find(b), find(c)
    if ra != rb:
        parent[rb] = ra
    if ra != rc:
        parent[rc] = ra
_, isl = np.unique(np.array([find(i) for i in range(n_uvv)]), return_inverse=True)
n_isl = isl.max() + 1

Puv = np.clip(uv, 0, 1)
tuv = np.stack([Puv[li], Puv[li + 1], Puv[li + 2]], axis=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]))
co = np.empty(n_v * 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)
fisl = isl[T[:, 0]]
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)

# real-world scale: body height in metres / up-span in mesh units
HEIGHT_M = 1.777
UNITM = HEIGHT_M / span[UP]
W = 4096
print(f"\n=== NEW ATLAS ===")
print(f"islands: {n_isl}   (Tripo atlas: 5870)")
print(f"atlas UV area used: {isl_auv.sum()*100:.1f}%   (Tripo: 62.0%)")
print(f"uv-vertices {n_uvv} for {n_v} verts -> {100.0*(n_uvv-n_v)/n_v:+.1f}% duplication")
o = np.argsort(-isl_auv)
cum = np.cumsum(isl_auv[o]) / max(isl_auv.sum(), 1e-12)
print(f"  {int(np.searchsorted(cum,0.99))+1} islands cover 99% of the used area  (Tripo: 84)")
print(f"  islands with <20 faces: {int((isl_nf<20).sum())}  (Tripo: 5763)")
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)
big = o[:int(np.searchsorted(cum, 0.99)) + 1]
db = dens[big]; db = db[np.isfinite(db)]
print(f"texel density over 99%-area islands: min {db.min():.2f} median {np.median(db):.2f} "
      f"max {db.max():.2f} px/mm -> {db.max()/max(db.min(),1e-9):.1f}x spread  (Tripo: 1.8x)")
print("\ntop 16 islands:  #   faces   uv_area%   3D cm2    px/mm   uv centre")
for i in o[:16]:
    m = fisl == i
    cu = tuv[m].reshape(-1, 2).mean(axis=0)
    print(f"  {i:5d} {isl_nf[i]:7d}  {isl_auv[i]*100:8.3f} {isl_a3[i]*UNITM*UNITM*1e4:9.1f}  "
          f"{dens[i]:6.2f}   ({cu[0]:.3f},{cu[1]:.3f})")

# per-triangle stretch: how anisotropic is the mapping (1.0 = conformal)
ok = (a3 > 1e-12) & (auv > 1e-14)
sc = np.sqrt(auv[ok] / a3[ok])
sc /= np.median(sc)
print(f"\narea-scale ratio vs median, per triangle: p05 {np.percentile(sc,5):.2f} "
      f"p50 {np.percentile(sc,50):.2f} p95 {np.percentile(sc,95):.2f} "
      f"p99 {np.percentile(sc,99):.2f}")

# ---- picture of the new layout ----
R = 1024
rng = np.random.RandomState(3)
pal = rng.rand(n_isl, 3) * 0.75 + 0.2
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
    dd_ = ((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(dd_) < 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])) / dd_
    bb = ((P[2, 1] - P[0, 1]) * (gx - P[2, 0]) + (P[0, 0] - P[2, 0]) * (gy - P[2, 1])) / dd_
    cc = 1.0 - aa - bb
    ins = (aa >= 0) & (bb >= 0) & (cc >= 0)
    if ins.any():
        IS[gy[ins], gx[ins]] = pal[fisl[fi]]
img = bpy.data.images.new("newislands", 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, "new_islands.png")
img.file_format = 'PNG'
img.filepath_raw = p
img.save(filepath=p)
log(f"wrote {p}")
bpy.context.preferences.filepaths.save_version = 0  # no .blend1 autosave
bpy.ops.wm.save_as_mainfile(filepath=os.path.join(OUTDIR, "unwrapped.blend"))
print("UNWRAP_TEST_DONE")