"""Mesh Burger: body, head and tail for the walking shiba (shibawalk.lsl), replacing the sphere/torus prims. Run:

  "C:/Program Files/Blender Foundation/Blender 5.2/blender.exe" -b -P gen_mesh.py -- [--out DIR] [--render]
      [--tris SCALE --suffix _lite] [--probe]

Design frame = the body (root) frame of shibawalk.lsl: origin at the body centre, +x forward, +z up. Each part is a
signed distance field -> surface nets -> decimated -> snapped back to the field, normals from the field (the
walker/ghost pipeline). Then a near-seamless UV layout (smooth charts folded at y = 0; see "UV layout (v2)"), and a
baseColor texture painted texel by texel from the 3D point under it: orange/cream urajiro regions, fur strands (3D noise smeared along a fur-flow direction, so seams don't show),
soft baked occlusion. Exported Z-up, x forward, each centred on its bounding box (SL puts the prim origin there),
so in the linkset every part sits at its bbox centre with rotation 0. Positions are printed and in README.md.
Claude Opus 5.5 (subagent for Tessel), 2026-10-11."""

import math
import os
import sys

import bmesh
import bpy
import numpy as np
from mathutils import Euler, Quaternion, Vector

argv = sys.argv[sys.argv.index("--") + 1:] if "--" in sys.argv else []
OUT = argv[argv.index("--out") + 1] if "--out" in argv else os.path.dirname(os.path.abspath(__file__))
RENDER = "--render" in argv
PROBE = "--probe" in argv
TRIS = float(argv[argv.index("--tris") + 1]) if "--tris" in argv else 1.0
SUFFIX = argv[argv.index("--suffix") + 1] if "--suffix" in argv else ""

HEAD_PIVOT = np.array([0.19, 0.0, 0.09])
TAIL_PIVOT = np.array([-0.14, 0.0, 0.06])
HIPS = [(0.13, 0.065, -0.05), (0.13, -0.065, -0.05), (-0.14, 0.065, -0.05), (-0.14, -0.065, -0.05)]


# ---------------------------------------------------------------- distance-field kit (P: (N, 3))

def ell(P, c, r):
    q = (P - c) / r
    k0 = np.linalg.norm(q, axis=1)
    k1 = np.linalg.norm(q / r, axis=1)
    return k0 * (k0 - 1.0) / np.maximum(k1, 1e-9)


def rcone(P, a, b, r1, r2):
    """Round cone from sphere (a, r1) to sphere (b, r2)."""
    a = np.asarray(a, float)
    ba = np.asarray(b, float) - a
    l2 = ba @ ba
    rr = r1 - r2
    a2 = l2 - rr * rr
    il2 = 1.0 / l2
    pa = P - a
    y = pa @ ba
    z = y - l2
    w = pa * l2 - np.outer(y, ba)
    x2 = np.einsum("ij,ij->i", w, w)
    y2 = y * y * l2
    z2 = z * z * l2
    k = np.sign(rr) * rr * rr * x2
    d3 = (np.sqrt(np.maximum(x2 * a2 * il2, 0)) + y * rr) * il2 - r1
    d1 = np.sqrt(x2 + z2) * il2 - r2
    d2 = np.sqrt(x2 + y2) * il2 - r1
    return np.where(np.sign(z) * a2 * z2 > k, d1, np.where(np.sign(y) * a2 * y2 < k, d2, d3))


def smin(a, b, k):
    h = np.clip(0.5 + 0.5 * (b - a) / k, 0, 1)
    return b * (1 - h) + a * h - k * h * (1 - h)


def smax(a, b, k):
    return -smin(-a, -b, k)


def sstep(e0, e1, x):
    t = np.clip((x - e0) / (e1 - e0), 0, 1)
    return t * t * (3 - 2 * t)


def mirror(f):
    def g(P):
        Q = P.copy()
        Q[:, 1] = np.abs(Q[:, 1])
        return f(Q)
    return g


def unit(v):
    v = np.asarray(v, float)
    return v / np.linalg.norm(v, axis=-1, keepdims=True)


def grad(f, P, e=1e-4):
    g = np.empty_like(P)
    for i in range(3):
        o = np.zeros(3)
        o[i] = e
        g[:, i] = (f(P + o) - f(P - o)) / (2 * e)
    return g


# ---------------------------------------------------------------- value noise

def hash3(ix, iy, iz, seed):
    h = (ix * 73856093) ^ (iy * 19349663) ^ (iz * 83492791) ^ (seed * 2654435761)
    h &= 0xFFFFFFFF
    h = ((h ^ (h >> 13)) * 1274126177) & 0xFFFFFFFF
    h ^= h >> 16
    return (h & 0xFFFF) / 32767.5 - 1.0


def vnoise(P, h, seed=0):
    Q = P / h
    I = np.floor(Q).astype(np.int64)
    f = Q - I
    s = f * f * (3 - 2 * f)
    out = np.zeros(len(P))
    for dx in (0, 1):
        wx = s[:, 0] if dx else 1 - s[:, 0]
        for dy in (0, 1):
            wy = s[:, 1] if dy else 1 - s[:, 1]
            for dz in (0, 1):
                wz = s[:, 2] if dz else 1 - s[:, 2]
                out += wx * wy * wz * hash3(I[:, 0] + dx, I[:, 1] + dy, I[:, 2] + dz, seed)
    return out


def strands(P, D, width, length, seed):
    """Fine noise smeared along D (line integral convolution in 3D): short streaks, about `width` across."""
    acc = np.zeros(len(P))
    wsum = 0.0
    for s in np.linspace(-0.5, 0.5, 9):
        w = 1.0 - abs(s) * 1.4
        acc += w * vnoise(P + D * (s * length), width, seed)
        wsum += w
    acc /= wsum
    return acc / (acc.std() + 1e-9)


# ---------------------------------------------------------------- shapes (body frame)

def body(P):
    d = ell(P, (-0.01, 0, 0.0), (0.195, 0.096, 0.086))                    # torso
    d = smin(d, ell(P, (0.125, 0, -0.008), (0.088, 0.094, 0.094)), 0.05)   # deep chest
    d = smin(d, ell(P, (-0.125, 0, 0.0), (0.088, 0.094, 0.086)), 0.04)    # rump
    sym = lambda Q: smin(ell(Q, (-0.135, 0.052, -0.035), (0.072, 0.05, 0.068)),      # haunches
                         ell(Q, (0.125, 0.05, -0.04), (0.06, 0.046, 0.064)), 0.02)   # shoulders
    d = smin(d, mirror(sym)(P), 0.035)
    d = smin(d, rcone(P, (0.135, 0, 0.02), (0.175, 0, 0.056), 0.06, 0.045), 0.04)  # neck root (inside the ruff)
    # a slight tuck behind the ribs: the underside rises toward the waist
    x = P[:, 0]
    floor_z = -0.1 + 0.022 * np.exp(-((x + 0.04) / 0.075) ** 2)
    d = smax(d, floor_z - P[:, 2], 0.03)
    return d


EYE_YZ = (0.043, 0.166)
NOSE_C, NOSE_R = (0.364, 0.0, 0.121), (0.0135, 0.0175, 0.012)
EAR_BASE, EAR_TIP = np.array([0.224, 0.05, 0.196]), np.array([0.238, 0.073, 0.284])


def head_core(P):
    d = ell(P, (0.245, 0, 0.148), (0.084, 0.088, 0.078))                          # cranium
    d = smin(d, ell(P, (0.314, 0, 0.111), (0.054, 0.035, 0.03)), 0.025)            # muzzle
    d = smin(d, ell(P, (0.302, 0, 0.088), (0.03, 0.024, 0.016)), 0.015)            # chin
    d = smin(d, mirror(lambda Q: ell(Q, (0.272, 0.052, 0.112), (0.048, 0.042, 0.04)))(P), 0.03)   # cheek puffs
    d = smin(d, np.linalg.norm(P - HEAD_PIVOT, axis=1) - 0.068, 0.035)             # neck ruff, a ball about the pivot
    return d


def ear_frame():
    a = unit(EAR_TIP - EAR_BASE)
    t = np.array([1.0, 0.35, 0.0])
    t = unit(t - a * (a @ t))           # thickness axis: the ear's front faces forward and a little out
    w = np.cross(a, t)
    return a, t, w, float(np.linalg.norm(EAR_TIP - EAR_BASE))


def ear_local(P):
    a, t, w, L = ear_frame()
    D = P - EAR_BASE
    return D @ a, D @ w, D @ t, L


def ear(P):
    u, v, q, L = ear_local(P)
    s = np.clip(u / L, 0, 1)
    hw = 0.031 * (1 - s) + 0.0035
    ht = 0.013 * (1 - s) + 0.003
    side = (np.sqrt((v / hw) ** 2 + (q / ht) ** 2) - 1) * ht
    d = smax(side, smax(-0.02 - u, u - L, 0.002), 0.003)    # runs on 2 cm below the base, buried in the skull
    # the open front of the ear: a shallow scoop
    hi = hw * 0.62
    ti = ht * 0.75
    inner = (np.sqrt((v / hi) ** 2 + ((q - ht * 0.95) / ti) ** 2) - 1) * ti
    inner = smax(inner, smax(0.12 * L - u, u - 0.82 * L, 0.002), 0.002)
    return smax(d, -inner, 0.0025)


def eye_frame():
    y, z = EYE_YZ
    xs = np.linspace(0.36, 0.2, 1601)
    pts = np.stack([xs, np.full_like(xs, y), np.full_like(xs, z)], 1)
    S = pts[np.argmax(head_core(pts) < 0)]
    n = unit(grad(head_core, S[None])[0])
    th = unit(np.cross([0, 0, 1.0], n))
    tv = np.cross(n, th)
    return S, n, th, tv


def eye(P):
    S, n, th, tv = EYE
    c = S - n * 0.0035
    D = P - c
    Q = np.stack([D @ th, D @ tv, D @ n], 1)
    return ell(Q, (0, 0, 0), (0.0135, 0.013, 0.0075))


def head(P):
    return smin(head_core(P), mirror(ear)(P), 0.014)


def eyes(P):
    return mirror(eye)(P)


def nose(P):
    return ell(P, NOSE_C, NOSE_R)


# tail: a fat curl rising off the rump, rolling forward over the back; the tip comes to rest a little to the left
TAIL_PTS = np.array([(-0.150, 0.000, 0.055), (-0.188, 0.000, 0.105), (-0.192, 0.002, 0.160), (-0.158, 0.006, 0.198),
                     (-0.104, 0.012, 0.196), (-0.070, 0.018, 0.160), (-0.072, 0.024, 0.122), (-0.100, 0.030, 0.108),
                     (-0.124, 0.034, 0.124)])
TAIL_R = np.array([0.024, 0.029, 0.033, 0.034, 0.033, 0.030, 0.025, 0.019, 0.012])
CURL_C = np.array([-0.135, 0.01, 0.148])


def catmull(pts, n):
    P = np.vstack([2 * pts[0] - pts[1], pts, 2 * pts[-1] - pts[-2]])
    out = []
    for i in range(1, len(P) - 2):
        p0, p1, p2, p3 = P[i - 1], P[i], P[i + 1], P[i + 2]
        for t in np.linspace(0, 1, n, endpoint=False):
            out.append(0.5 * (2 * p1 + (-p0 + p2) * t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t * t
                              + (-p0 + 3 * p1 - 3 * p2 + p3) * t ** 3))
    out.append(pts[-1])
    return np.array(out)


TAIL_C = catmull(TAIL_PTS, 4)
TAIL_RC = np.interp(np.linspace(0, len(TAIL_R) - 1, len(TAIL_C)), np.arange(len(TAIL_R)), TAIL_R)


def tail(P):
    d = np.full(len(P), 1e9)
    for i in range(len(TAIL_C) - 1):
        d = smin(d, rcone(P, TAIL_C[i], TAIL_C[i + 1], TAIL_RC[i], TAIL_RC[i + 1]), 0.004)
    # tufts: a soft lumpy outline so it reads as fluff, not a hose
    return d - 0.0035 * vnoise(P, 0.012, 5)


def scene_sdf(P):
    return np.minimum.reduce([body(P), head(P), eyes(P), nose(P), tail(P)])


# ---------------------------------------------------------------- paint

def srgb_lin(c):
    c = np.asarray(c, float)
    return np.where(c <= 0.04045, c / 12.92, ((c + 0.055) / 1.055) ** 2.4)


def hexc(h):
    return np.array([int(h[i:i + 2], 16) / 255 for i in (1, 3, 5)])


ORANGE, SADDLE, LIGHT_OR = hexc("#d9772b"), hexc("#b65c1d"), hexc("#e9a35e")
CREAM, CREAM_SHADE = hexc("#f1dbb5"), hexc("#dcc097")
INNER_EAR = hexc("#f7e2cc")
NOSE_BLACK, EYE_DARK, RIM, MOUTH = hexc("#1b1513"), hexc("#140d0b"), hexc("#2b1a12"), hexc("#4a2c1c")


def mix(a, b, t):
    t = np.asarray(t)[:, None]
    return a * (1 - t) + b * t


def tangent(D0, N):
    D = D0 - N * np.einsum("ij,ij->i", D0, N)[:, None]
    n = np.linalg.norm(D, axis=1, keepdims=True)
    return np.where(n > 1e-6, D / np.maximum(n, 1e-9), np.array([1.0, 0, 0]))


def throat(P, N):
    """Cream bib down the front of the neck and chest; >0 inside. Shared by body and head so the seam matches."""
    x, y, z = P.T
    front = x - (0.165 + 0.4 * np.clip(z, 0, None))              # the bib's back edge leans forward as it rises
    width = 0.042 - np.abs(y) + 0.15 * np.clip(-z, 0, None)       # narrow at the throat, wider down the chest
    return np.minimum.reduce([front, width, 0.115 - z])


def body_flow(P, N):
    x, y, z = P.T
    side = np.clip(np.abs(y) / 0.08, 0, 1)
    D0 = np.stack([-np.ones_like(x), 0.35 * np.sign(y) * side, -0.75 * side], 1)
    front = sstep(0.12, 0.2, x) * np.clip(N[:, 0], 0, 1)
    D0 = D0 * (1 - front[:, None]) + np.array([0.2, 0, -1.0]) * front[:, None]
    return tangent(D0, N)


def body_paint(P, N, s, clump):
    x, y, z = P.T
    fz = 0.0045 * s + 0.004 * clump                               # furry edges between the colours
    belly = -0.052 - z                                            # under the lower flanks
    rump = np.minimum(-0.168 - x, -0.0 - z)                       # cream under the tail and behind the thighs
    cream = np.maximum.reduce([belly, throat(P, N), rump]) + fz
    m = sstep(-0.004, 0.004, cream)
    top = sstep(0.02, 0.09, z) * np.clip(N[:, 2], 0, 1)           # a slightly darker saddle along the spine
    low = sstep(0.0, -0.05, z)                                     # orange warms and lightens toward the cream
    orange = mix(mix(ORANGE, SADDLE, 0.55 * top), LIGHT_OR, 0.45 * low)
    cream_c = mix(CREAM, CREAM_SHADE, 0.35 * sstep(-0.06, -0.1, z))
    return mix(orange, cream_c, m)


def head_flow(P, N):
    x, y, z = P.T
    D0 = np.stack([-np.ones_like(x), np.zeros_like(x), -0.15 * np.ones_like(x)], 1)
    cheek = sstep(0.03, 0.07, np.abs(y)) * sstep(0.15, 0.11, z)
    D0 += np.stack([0.4 * np.ones_like(x), 1.3 * np.sign(y), -0.4 * np.ones_like(x)], 1) * cheek[:, None]
    neck = sstep(0.2, 0.16, x) * sstep(0.12, 0.08, z)
    D0 += np.array([0.3, 0, -1.2]) * neck[:, None]
    u, v, q, L = ear_local(np.stack([x, np.abs(y), z], 1))
    in_ear = sstep(-0.004, 0.008, u) * sstep(0.02, 0.01, np.sqrt(v * v + q * q) - 0.6 * np.clip(L - u, 0, L) * 0.4)
    a = ear_frame()[0]
    ear_dir = np.stack([np.full_like(x, a[0]), np.sign(y) * a[1], np.full_like(x, a[2])], 1)
    D0 = D0 * (1 - in_ear[:, None]) + ear_dir * 2 * in_ear[:, None]
    return tangent(D0, N)


def head_paint(P, N, s, clump):
    x, y, z = P.T
    ay = np.abs(y)
    fz = 0.003 * s + 0.003 * clump
    # urajiro: everything below a line that runs just under the eyes and sweeps back and down to the throat, except
    # an orange stripe down the bridge of the muzzle that narrows to the nose
    under = 0.152 - 0.6 * np.clip(0.3 - x, 0, None) - z
    bridge_w = np.clip(0.018 - 0.16 * (x - 0.29), 0.006, 0.018)
    bridge = np.minimum.reduce([bridge_w - ay, z - 0.122 + 0.3 * np.clip(x - 0.335, 0, None), x - 0.27])
    cream = np.maximum(np.minimum(under, -bridge), throat(P, N)) + fz
    # eyebrow dots
    for c in BROWS:
        cream = np.maximum(cream, 0.0072 - np.linalg.norm(P - c, axis=1) + 0.5 * fz)
    m = sstep(-0.003, 0.003, cream)
    top = sstep(0.17, 0.22, z) * np.clip(N[:, 2], 0, 1)
    nearc = sstep(0.0, -0.012, -cream)                     # orange lightens next to the cream, and along the bridge
    orange = mix(mix(ORANGE, SADDLE, 0.4 * top), LIGHT_OR, 0.35 * nearc + 0.4 * sstep(0.3, 0.35, x))
    col = mix(orange, CREAM, m)
    # ears: cream scoop inside, darker tips
    u, v, q, L = ear_local(np.stack([x, ay, z], 1))
    a, t, w, _ = ear_frame()
    s_ = np.clip(u / L, 0, 1)
    hw = 0.031 * (1 - s_) + 0.0035
    t_side = np.stack([np.full_like(x, t[0]), np.sign(y) * t[1], np.full_like(x, t[2])], 1)
    front = np.einsum("ij,ij->i", N, t_side) > 0.2
    on_ear = (mirror(ear)(P) < 0.002) & (head_core(P) > 0.004)
    inner = (u > 0.08 * L) & (u < 0.86 * L) & (q > 0.0) & front & on_ear
    ear_in = inner * sstep(-0.004, 0.002, hw * 0.7 - np.abs(v) + 0.0015 * s)
    col = mix(col, INNER_EAR, ear_in)
    tipdark = sstep(0.7, 1.0, s_) * (u > 0) * (np.sqrt(v * v + q * q) < 0.02)
    col = mix(col, SADDLE * 0.8, 0.5 * tipdark)
    # dark rims round the eyes, and a little mouth under the nose
    rim = sstep(0.0035, 0.0012, eyes(P))
    col = mix(col, RIM, rim)
    col = mix(col, MOUTH, mouth(P, N))
    return col


def mouth(P, N):
    """A small 'w': a philtrum line down from the nose, splitting into two curves. Painted on the muzzle front."""
    x, y, z = P.T
    nb = NOSE_C[2] - NOSE_R[2] + 0.001
    ys = np.linspace(0, 0.019, 30)
    zs = nb - 0.009 - 0.004 * np.sin(np.pi * ys / 0.019) + 0.003 * (ys / 0.019) ** 2
    pts = np.concatenate([np.stack([np.zeros(8), np.linspace(nb, nb - 0.009, 8)], 1),
                          np.stack([ys, zs], 1)])
    yz = np.stack([np.abs(y), z], 1)
    d = np.min(np.linalg.norm(yz[:, None, :] - pts[None], axis=2), axis=1)
    on = (x > 0.33) & (N[:, 0] > 0.2)
    return sstep(0.0011, 0.0005, d) * on


def eye_paint(P, N, s, clump):
    col = mix(np.tile(EYE_DARK, (len(P), 1)), hexc("#3a2418"), 0.5 * sstep(0.004, -0.008, (P[:, 2] - EYE_YZ[1])))
    # catchlights: a big one up and to the front-outside, a tiny one below; a painted highlight reads at any angle
    L1 = unit(np.array([0.75, 0.0, 0.66]))
    L2 = unit(np.array([0.9, 0.0, -0.4]))
    sideways = np.stack([np.zeros(len(P)), np.sign(P[:, 1]) * 0.35, np.zeros(len(P))], 1)
    c1 = np.einsum("ij,ij->i", N, unit(L1 + sideways))
    c2 = np.einsum("ij,ij->i", N, unit(L2 + sideways * 0.3))
    col = mix(col, np.array([1.0, 0.98, 0.95]), sstep(0.955, 0.975, c1))
    col = mix(col, np.array([0.8, 0.75, 0.72]), sstep(0.988, 0.994, c2))
    return col


def nose_paint(P, N, s, clump):
    col = mix(np.tile(NOSE_BLACK, (len(P), 1)), hexc("#3a2c28"), 0.6 * sstep(0.3, 0.9, N[:, 2]))
    # nostrils: two small dark commas on the front
    for sy in (1, -1):
        d = np.linalg.norm((P - (NOSE_C[0] + 0.012, sy * 0.0075, NOSE_C[2] - 0.001)) * (1, 1, 0.8), axis=1)
        col = mix(col, np.array([0.03, 0.02, 0.02]), sstep(0.0042, 0.0028, d))
    return col


def tail_near(P):
    d = np.linalg.norm(P[:, None, :] - TAIL_C[None], axis=2)
    i = d.argmin(1)
    i0 = np.clip(i - 1, 0, len(TAIL_C) - 1)
    i1 = np.clip(i + 1, 0, len(TAIL_C) - 1)
    T = unit(TAIL_C[i1] - TAIL_C[i0])
    R = P - TAIL_C[i]
    R = unit(R - T * np.einsum("ij,ij->i", R, T)[:, None])
    return i, T, R


def tail_flow(P, N):
    i, T, R = tail_near(P)
    return tangent(T + 0.7 * R, N)


def tail_paint(P, N, s, clump):
    i, T, R = tail_near(P)
    inward = unit(CURL_C - TAIL_C[i])
    inward = unit(inward - T * np.einsum("ij,ij->i", inward, T)[:, None])
    k = np.einsum("ij,ij->i", R, inward)                   # +1 facing the middle of the curl
    under = -N[:, 2] * sstep(3, 0, i)                       # the underside of the rising base
    cream = np.maximum(k + 0.25, under - 0.3) + 0.06 * s + 0.04 * clump
    m = sstep(-0.12, 0.12, cream)
    tip = sstep(len(TAIL_C) - 6, len(TAIL_C) - 1, i)        # the tip goes pale
    orange = mix(mix(ORANGE, SADDLE, 0.35 * np.clip(-k, 0, 1)), LIGHT_OR, 0.5 * tip)
    return mix(orange, CREAM, np.maximum(m, 0.6 * tip))


# ---------------------------------------------------------------- meshing

def surface_nets(f, lo, hi, h):
    xs, ys, zs = (np.arange(lo[i] - 2 * h, hi[i] + 2 * h, h) for i in range(3))
    G = np.stack(np.meshgrid(xs, ys, zs, indexing="ij"), -1)
    F = f(G.reshape(-1, 3)).reshape(G.shape[:3])
    for ax in range(3):
        assert F.take([0, -1], axis=ax).min() > 0, "shape reaches the grid edge"
    n = np.array(F.shape)
    cdims = n - 1
    cid = lambda i, j, k: (i * cdims[1] + j) * cdims[2] + k
    acc = np.zeros((cdims.prod(), 3))
    cnt = np.zeros(cdims.prod())
    quads = []
    for ax in range(3):
        s0 = [slice(None)] * 3
        s1 = [slice(None)] * 3
        s0[ax] = slice(0, -1)
        s1[ax] = slice(1, None)
        f0, f1 = F[tuple(s0)], F[tuple(s1)]
        cross = (f0 < 0) != (f1 < 0)
        idx = np.argwhere(cross)
        a, b = f0[cross], f1[cross]
        t = a / (a - b)
        pts = np.stack([xs[0] + idx[:, 0] * h, ys[0] + idx[:, 1] * h, zs[0] + idx[:, 2] * h], 1)
        pts[:, ax] += t * h
        u, v = [x for x in range(3) if x != ax]
        cells = []
        for du, dv in ((-1, -1), (0, -1), (0, 0), (-1, 0)):
            c = idx.copy()
            c[:, u] += du
            c[:, v] += dv
            cells.append(c)
        for c in cells:
            ok = np.all((c >= 0) & (c < cdims), axis=1)
            ci = cid(c[ok, 0], c[ok, 1], c[ok, 2])
            np.add.at(acc, ci, pts[ok])
            np.add.at(cnt, ci, 1)
        q = np.stack([cid(c[:, 0], c[:, 1], c[:, 2]) for c in cells], 1)
        flip = (a < 0) if ax != 1 else (a >= 0)
        q[~flip] = q[~flip][:, ::-1]
        quads.append(q)
    quads = np.concatenate(quads)
    used = np.nonzero(cnt)[0]
    remap = -np.ones(cdims.prod(), int)
    remap[used] = np.arange(len(used))
    V = acc[used] / cnt[used, None]
    return V, remap[quads]


def project(f, V, iters=4):
    for _ in range(iters):
        g = grad(f, V)
        V = V - (f(V) / np.maximum(np.einsum("ij,ij->i", g, g), 1e-9))[:, None] * g
    return V


def decimated(f, lo, hi, h, tris, name, sym=True):
    """Surface nets -> collapse decimate to ~tris -> (mirrored in y if sym) -> snapped onto the field."""
    V, Q = surface_nets(f, lo, hi, h)
    me = bpy.data.meshes.new(name + "_raw")
    me.from_pydata(V.tolist(), [], Q.tolist())
    ob = bpy.data.objects.new(name + "_raw", me)
    bpy.context.scene.collection.objects.link(ob)
    m = ob.modifiers.new("dec", "DECIMATE")
    m.ratio = min(1.0, tris / (2 * len(Q)))     # quads -> tris; with sym, half is kept then mirrored back
    m.use_collapse_triangulate = True
    bm = bmesh.new()
    bm.from_object(ob, bpy.context.evaluated_depsgraph_get())
    bpy.data.objects.remove(ob)
    bpy.data.meshes.remove(me)
    if sym:
        bmesh.ops.bisect_plane(bm, geom=bm.verts[:] + bm.edges[:] + bm.faces[:], plane_co=(0, 0, 0),
                               plane_no=(0, 1, 0), clear_inner=True)
    bmesh.ops.triangulate(bm, faces=bm.faces[:])
    V = np.array([v.co[:] for v in bm.verts])
    bm.verts.index_update()
    F = np.array([[v.index for v in fc.verts] for fc in bm.faces])
    bm.free()
    area2 = np.linalg.norm(np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]]), axis=1)
    F = F[area2 > 1e-10]
    if sym:
        rim = np.nonzero(np.abs(V[:, 1]) < 1e-6)[0]
        V[rim, 1] = 0
        n = len(V)
        other = np.arange(n) + n
        other[rim] = rim
        V = np.concatenate([V, V * (1, -1, 1)])
        F = np.concatenate([F, other[F[:, ::-1]]])
    V = project(f, V)
    area2 = np.linalg.norm(np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]]), axis=1)
    F = F[area2 > 1e-10]
    used = np.unique(F)
    remap = -np.ones(len(V), int)
    remap[used] = np.arange(len(used))
    return V[used], remap[F]


# ---------------------------------------------------------------- texture

def raster(UV, T3, W):
    """UV: (T, 3, 2) per-triangle uvs; T3: (T, 3, 3) the triangles' 3D points. -> per-texel 3D point, triangle id."""
    P = np.zeros((W, W, 3))
    tri = -np.ones((W, W), int)
    for t in range(len(UV)):
        a, b, c = UV[t] * W - 0.5                  # texel (i, j) centre is at pixel coord (i, j)
        lo = np.floor(np.minimum.reduce([a, b, c])).astype(int)
        hi = np.ceil(np.maximum.reduce([a, b, c])).astype(int)
        lo = np.clip(lo, 0, W - 1)
        hi = np.clip(hi, 0, W - 1)
        if (hi < lo).any():
            continue
        I, J = np.meshgrid(np.arange(lo[0], hi[0] + 1), np.arange(lo[1], hi[1] + 1), indexing="ij")
        q = np.stack([I.ravel(), J.ravel()], 1).astype(float)
        v0, v1 = b - a, c - a
        den = v0[0] * v1[1] - v1[0] * v0[1]
        if abs(den) < 1e-12:
            continue
        r = q - a
        l1 = (r[:, 0] * v1[1] - v1[0] * r[:, 1]) / den
        l2 = (v0[0] * r[:, 1] - r[:, 0] * v0[1]) / den
        l0 = 1 - l1 - l2
        ok = (l0 >= -1e-6) & (l1 >= -1e-6) & (l2 >= -1e-6)
        ii, jj = I.ravel()[ok], J.ravel()[ok]
        P[jj, ii] = l0[ok, None] * T3[t, 0] + l1[ok, None] * T3[t, 1] + l2[ok, None] * T3[t, 2]
        tri[jj, ii] = t
    return P, tri


def dilate(img, filled):
    img = img.copy()
    filled = filled.copy()
    for _ in range(len(img)):
        if filled.all():
            break
        acc = np.zeros_like(img)
        cnt = np.zeros(filled.shape)
        for dy, dx in ((0, 1), (0, -1), (1, 0), (-1, 0), (1, 1), (1, -1), (-1, 1), (-1, -1)):
            f = np.roll(np.roll(filled, dy, 0), dx, 1)
            acc += np.roll(np.roll(img, dy, 0), dx, 1) * f[..., None]
            cnt += f
        grow = (~filled) & (cnt > 0)
        img[grow] = acc[grow] / cnt[grow, None]
        filled |= grow
    return img


def occlusion(P, N, steps=5, delta=0.007):
    ao = np.zeros(len(P))
    for k in range(1, steps + 1):
        ao += (k * delta - scene_sdf(P + N * (k * delta))) / 2 ** k
    return np.clip(ao / delta, 0, 1)


# UV layout (v2). SL's uploader makes the lower LODs with meshoptimizer, and when it can't reach the triangle target
# while keeping UV seams it welds across them (or goes "sloppy"), which drags UVs from one island into another: the
# likely cause of v1's calico patches in SL (Smart UV Project, dozens of islands; the files themselves round-trip
# fine, check_glb.py). So v2 has almost no seams: each chart is a smooth map of
# the surface, folded at the y = 0 plane (left and right halves share texels, and a fold is not a seam). Only the
# ears are cut out (cylinder maps that wrap round the full texture width, so REPEAT sampling hides the wrap).
# Returns per-corner uvs, which triangles paint the texture (one side of each fold), and which ones wrap in u.

M = 0.012          # margin kept free round charts that touch the texture edge


def fit(x, lo, hi):
    return lo + (hi - lo) * x


def uv_body(C):
    Q = C - np.array([0.0, 0.0, -0.02])                     # axis low enough that both poles land in cream
    yy, zz = np.abs(Q[..., 1]) / 0.105, Q[..., 2] / 0.095
    u = np.arctan2(yy, zz) / np.pi                          # 0 along the spine, 1 along the belly midline
    v = np.arctan2(Q[..., 0] / 0.22, np.hypot(yy, zz)) / np.pi + 0.5
    return np.stack([fit(u, M, 1 - M), fit(v, M, 1 - M)], -1)


HEAD_C = np.array([0.245, 0.0, 0.13])
U_CUT = 0.012      # ear triangles start this far up the ear axis from EAR_BASE
HEAD_V, EAR_V, BOX_V = (M, 0.70), (0.72, 0.86), 0.875


def uv_head_fur(C):
    Q = (C - HEAD_C) / np.array([0.09, 0.09, 0.085])
    u = np.arctan2(np.abs(Q[..., 1]), Q[..., 0]) / np.pi   # 0 at the nose, 1 at the nape
    v = np.arctan2(Q[..., 2], np.hypot(Q[..., 0], Q[..., 1])) / np.pi + 0.5
    return np.stack([fit(u, M, 1 - M), fit(v, *HEAD_V)], -1)


def ear_coords(C):
    F = C.copy()
    F[..., 1] = np.abs(F[..., 1])
    u, v, q, L = ear_local(F.reshape(-1, 3))
    s_ = np.clip(u / L, 0, 1)
    hw = 0.031 * (1 - s_) + 0.0035
    ht = 0.013 * (1 - s_) + 0.003
    sh = C.shape[:-1]
    return u.reshape(sh), (v / hw).reshape(sh), (q / ht).reshape(sh), L


def uv_ear(C):
    u, vn, qn, L = ear_coords(C)
    phi = np.arctan2(qn, vn)                               # round the ear; +-pi on the rim nearest the head's midline
    U = (phi + np.pi) / (2 * np.pi)
    V = fit(np.clip((u - U_CUT) / (L - U_CUT), 0, 1), *EAR_V)
    return np.stack([U, V], -1)


def uv_eye(C):
    S, n, th, tv = EYE
    F = C.copy()
    F[..., 1] = np.abs(F[..., 1])
    Q = F - (S - n * 0.0035)
    a, b = Q @ th, Q @ tv                                   # front view of the eye (planar; the back is buried)
    return np.stack([0.07 + 0.06 * a / 0.0145, BOX_V + 0.06 + 0.06 * b / 0.0145], -1), Q @ n


def uv_nose(C):
    Q = (C - np.array(NOSE_C)) / np.array(NOSE_R)
    u = np.arctan2(np.abs(Q[..., 1]), Q[..., 0]) / np.pi
    v = np.arctan2(Q[..., 2], np.hypot(Q[..., 0], Q[..., 1])) / np.pi + 0.5
    return np.stack([fit(u, 0.15, 0.27), fit(v, BOX_V, BOX_V + 0.12)], -1)


SEG_A, SEG_B = TAIL_C[:-1], TAIL_C[1:]
SEG_L = np.linalg.norm(SEG_B - SEG_A, axis=1)
SEG_T = (SEG_B - SEG_A) / SEG_L[:, None]
SEG_S0 = np.concatenate([[0], np.cumsum(SEG_L)])[:-1]
TAIL_LEN = SEG_L.sum()


def tail_param(P):
    """Arc length of the nearest point on the tail's centre line (running on past the ends), the side direction
    (unit, from the centre line), the outward direction (away from the middle of the curl) and the tangent."""
    rel = P[:, None, :] - SEG_A[None]
    s = np.einsum("nsk,sk->ns", rel, SEG_T)
    sc = np.clip(s, 0, SEG_L[None])
    d = np.linalg.norm(rel - sc[..., None] * SEG_T[None], axis=2)
    i = d.argmin(1)
    si = s[np.arange(len(P)), i]
    last = len(SEG_L) - 1
    si = np.where(i == 0, np.minimum(si, SEG_L[0]), np.where(i == last, np.maximum(si, 0), np.clip(si, 0, SEG_L[i])))
    t = SEG_S0[i] + si
    foot = SEG_A[i] + np.clip(si, 0, SEG_L[i])[:, None] * SEG_T[i]
    T = SEG_T[i]
    R = P - foot
    R = unit(R - T * np.einsum("ij,ij->i", R, T)[:, None] + 1e-12)
    o = foot - CURL_C
    o = unit(o - T * np.einsum("ij,ij->i", o, T)[:, None])
    return t, R, o, T


def uv_tail(C):
    t, R, o, T = tail_param(C.reshape(-1, 3))
    U = np.arccos(np.clip(np.einsum("ij,ij->i", R, o), -1, 1)) / np.pi     # 0 outside of the curl, 1 inside
    side = np.einsum("ij,ij->i", R, np.cross(T, o))
    V = np.clip(t / TAIL_LEN, -0.02, 1.02)
    uv = np.stack([fit(U, M, 1 - M), fit(V, 0.03, 0.97)], -1)
    return uv.reshape(C.shape[:-1] + (2,)), side.reshape(C.shape[:-1])


def layout(name, T3, pid):
    """-> uv (T, 3, 2), paint (T,) bool, wrap (T,) bool, chart (T,) int for the report."""
    nT = len(T3)
    cen = T3.mean(1)
    plus = cen[:, 1] >= 0
    wrap = np.zeros(nT, bool)
    chart = np.zeros(nT, int)
    if name == "body":
        return uv_body(T3), plus, wrap, chart
    if name == "tail":
        uv, side = uv_tail(T3)
        return uv, side.mean(1) >= 0, wrap, chart
    uv = np.zeros((nT, 3, 2))
    paint = plus.copy()
    fur = pid == 0
    u, vn, qn, L = ear_coords(cen)
    is_ear = fur & (u > U_CUT) & (np.hypot(vn, qn) < 2.0)
    hf = fur & ~is_ear
    uv[hf] = uv_head_fur(T3[hf])
    e = uv_ear(T3[is_ear])
    span = e[..., 0].max(1) - e[..., 0].min(1)
    fix = span > 0.5                                        # straddles the wrap: unwrap past 1, REPEAT does the rest
    e[fix, :, 0] = np.where(e[fix, :, 0] < 0.5, e[fix, :, 0] + 1, e[fix, :, 0])
    uv[is_ear] = e
    wrap[np.nonzero(is_ear)[0][fix]] = True
    chart[is_ear] = 1
    ey = pid == 1
    uv[ey], depth = uv_eye(T3[ey])
    paint[ey] &= depth.mean(1) > 0                          # paint from the visible front; the buried back follows
    chart[ey] = 2
    no = pid == 2
    uv[no] = uv_nose(T3[no])
    chart[no] = 3
    return uv, paint, wrap, chart


def paint_texture(name, ob, parts, W, strand_w):
    """Unwrap ob and paint a W x W baseColor PNG from each triangle's part (sdf, flow, paint)."""
    me = ob.data
    nl = len(me.loops)
    lv = np.empty(nl, int)
    me.loops.foreach_get("vertex_index", lv)
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    co = co.reshape(-1, 3) + OB_CENTRE[name]               # back to the body frame
    T3 = co[lv].reshape(-1, 3, 3)
    tri_part = PART_OF_TRI[name]
    uv, paint_tri, wrap, chart = layout(name, T3, tri_part)
    me.uv_layers.active.data.foreach_set("uv", uv.ravel())
    me.update()
    # each painted side of a fold should keep one orientation; count triangles that flip (folded over in UV)
    area = np.cross(uv[:, 1] - uv[:, 0], uv[:, 2] - uv[:, 0])
    for c in np.unique(chart):
        sel = paint_tri & (chart == c)
        if sel.any():
            pos = (area[sel] > 0).sum()
            print(f"  {name} chart {c}: {sel.sum()} painted tris, {min(pos, sel.sum() - pos)} flipped, "
                  f"{wrap[chart == c].sum()} wrapping")
    sel = np.nonzero(paint_tri)[0]
    extra = sel[wrap[sel]]
    UVp = np.concatenate([uv[sel], uv[extra] - np.array([1.0, 0])])
    T3p = np.concatenate([T3[sel], T3[extra]])
    owner = np.concatenate([sel, extra])
    P, tri = raster(UVp, T3p, W)
    filled = tri >= 0
    tri = np.where(filled, owner[np.maximum(tri, 0)], -1)
    img = np.zeros((W, W, 3))
    Pf, tf = P[filled], tri[filled]
    pid = tri_part[tf]
    out = np.zeros((len(Pf), 3))
    for k, (sdf, flow, paint) in enumerate(parts):
        sel = pid == k
        if not sel.any():
            continue
        X = Pf[sel]
        X = project(sdf, X, 2)
        Nn = unit(grad(sdf, X))
        s = strands(X, flow(X, Nn), strand_w, strand_w * 7, 11 + k) if flow else np.zeros(len(X))
        clump = vnoise(X, 0.007, 3 + k) + 0.5 * vnoise(X, 0.0035, 4 + k)
        col = paint(X, Nn, s, clump)
        fur = 1 + 0.07 * np.clip(s, -2.5, 2.5) + 0.045 * clump
        if paint in (eye_paint, nose_paint):
            fur = np.ones(len(X))
        ao = occlusion(X, Nn)
        if paint is head_paint:                # the ear scoops are shallow; full occlusion there reads as grey
            ao *= np.where(mirror(ear)(X) < 0.002, 0.35, 1.0)
        out[sel] = np.clip(col * fur[:, None] * (1 - 0.32 * ao)[:, None], 0, 1)
    img[filled] = out
    img = dilate(img, filled)
    im = bpy.data.images.new(name + "_tex", W, W, alpha=False)
    rgba = np.concatenate([img, np.ones((W, W, 1))], 2)
    im.pixels.foreach_set(rgba.ravel().astype(np.float32))
    path = os.path.join(OUT, f"{name}{SUFFIX}.png")
    im.filepath_raw = path
    im.file_format = "PNG"
    im.save()
    bpy.data.images.remove(im)
    print("wrote", path, f"({filled.mean() * 100:.0f}% of texels covered)")
    return bpy.data.images.load(path)


# ---------------------------------------------------------------- objects

def material(name, image, rough=0.8):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    nt = m.node_tree
    b = nt.nodes["Principled BSDF"]
    b.inputs["Roughness"].default_value = rough
    tex = nt.nodes.new("ShaderNodeTexImage")
    tex.image = image
    nt.links.new(tex.outputs["Color"], b.inputs["Base Color"])
    return m


OB_CENTRE, PART_OF_TRI = {}, {}


def build(name, parts):
    """parts: [(V, F, sdf)] in the body frame -> one bbox-centred object (centre recorded), normals from the fields."""
    Vs, Fs, Ns, pid, off = [], [], [], [], 0
    for i, (V, F, f) in enumerate(parts):
        Ns.append(unit(grad(f, V)))
        Vs.append(V)
        Fs.append(F + off)
        pid += [i] * len(F)
        off += len(V)
    V, F, N = np.concatenate(Vs), np.concatenate(Fs), np.concatenate(Ns)
    fn = np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]])
    gc = np.concatenate([grad(f, V_[F_].mean(1)) for (V_, F_, f) in parts])
    bad = (np.einsum("ij,ij->i", fn, gc) < 0).sum()
    centre = (V.min(0) + V.max(0)) / 2
    me = bpy.data.meshes.new(name)
    me.from_pydata((V - centre).tolist(), [], F.tolist())
    me.polygons.foreach_set("use_smooth", [True] * len(F))
    me.uv_layers.new(name="UVMap")
    me.normals_split_custom_set_from_vertices(N.tolist())
    me.update()
    ob = bpy.data.objects.new(name, me)
    bpy.context.scene.collection.objects.link(ob)
    OB_CENTRE[name] = centre
    PART_OF_TRI[name] = np.array(pid)
    lo, hi = V.min(0), V.max(0)
    print(f"{name}: {len(F)} tris, {len(V)} verts, {bad} faces against the field; size "
          f"{tuple(round(float(v), 3) for v in hi - lo)}, bbox centre {tuple(round(float(v), 4) for v in centre)}")
    return ob


def finish(ob, mats, face_of_part):
    for m in mats:
        ob.data.materials.append(m)
    ob.data.polygons.foreach_set("material_index", [face_of_part[p] for p in PART_OF_TRI[ob.name]])
    ob.data.update()


def export(ob, path):
    for o in bpy.context.scene.objects:
        o.select_set(o == ob)
    bpy.context.view_layer.objects.active = ob
    bpy.ops.export_scene.gltf(filepath=path, export_format="GLB", use_selection=True, export_yup=False,
                              export_apply=True)
    print("wrote", path)


# ---------------------------------------------------------------- build

for o in list(bpy.data.objects):
    bpy.data.objects.remove(o)

EYE = eye_frame()
BROWS = []
for sy in (1, -1):
    y, z = 0.036 * sy, 0.192
    xs = np.linspace(0.36, 0.2, 1601)
    pts = np.stack([xs, np.full_like(xs, y), np.full_like(xs, z)], 1)
    BROWS.append(pts[np.argmax(head_core(pts) < 0)])

# fit checks: hips (thigh tops) should be well inside the body; the head ball about the pivot
fb = lambda p: float(body(np.atleast_2d(np.asarray(p, float)))[0])
for hp in HIPS[:3:2]:
    print(f"  hip {hp}: body sdf {fb(hp):+.3f} (thigh radius 0.0275)")
print(f"  tail pivot {tuple(TAIL_PIVOT)}: body sdf {fb(TAIL_PIVOT):+.3f}, tail sdf {float(tail(TAIL_PIVOT[None])[0]):+.3f}")
print(f"  eye surface point {np.round(EYE[0], 3)}, normal {np.round(EYE[1], 2)}")

bV, bF = decimated(body, (-0.24, -0.13, -0.13), (0.24, 0.13, 0.14), 0.003, 1600 * TRIS, "body")
hV, hF = decimated(head, (0.1, -0.13, 0.0), (0.39, 0.13, 0.31), 0.0018, 1300 * TRIS, "head")
eV, eF = decimated(eyes, (0.28, -0.07, 0.14), (0.34, 0.07, 0.19), 0.0008, 120 * TRIS, "eyes")
nV, nF = decimated(nose, (0.34, -0.03, 0.1), (0.38, 0.03, 0.14), 0.001, 90 * TRIS, "nose")
tV, tF = decimated(tail, (-0.24, -0.06, 0.02), (-0.02, 0.09, 0.25), 0.0022, 1300 * TRIS, "tail", sym=False)

body_ob = build("body", [(bV, bF, body)])
head_ob = build("head", [(hV, hF, head), (eV, eF, eyes), (nV, nF, nose)])
tail_ob = build("tail", [(tV, tF, tail)])

if PROBE:
    sys.exit(0)

os.makedirs(OUT, exist_ok=True)
img_b = paint_texture("body", body_ob, [(body, body_flow, body_paint)], 512, 0.0016)
img_h = paint_texture("head", head_ob, [(head, head_flow, head_paint), (eyes, None, eye_paint),
                                        (nose, None, nose_paint)], 512, 0.0011)
img_t = paint_texture("tail", tail_ob, [(tail, tail_flow, tail_paint)], 256, 0.0018)
finish(body_ob, [material("fur", img_b)], [0])
finish(head_ob, [material("head fur", img_h), material("eyes nose", img_h, 0.15)], [0, 1, 1])
finish(tail_ob, [material("tail fur", img_t)], [0])
for ob in (body_ob, head_ob, tail_ob):
    export(ob, os.path.join(OUT, f"{ob.name}{SUFFIX}.glb"))

print("PLACEMENT: bbox centres in the design frame, then local positions relative to the mesh body (the root):")
for ob in (body_ob, head_ob, tail_ob):
    c, rel = OB_CENTRE[ob.name], OB_CENTRE[ob.name] - OB_CENTRE["body"]
    print(f"  {ob.name}: centre <{', '.join(f'{v:.4f}' for v in c)}>, local <{', '.join(f'{v:.3f}' for v in rel)}>")

if not RENDER:
    sys.exit(0)

# ---------------------------------------------------------------- renders

scene = bpy.context.scene
scene.render.engine = "BLENDER_EEVEE"
scene.eevee.taa_render_samples = 16
scene.render.resolution_x, scene.render.resolution_y = 640, 480
scene.view_settings.view_transform = "Standard"
scene.world = bpy.data.worlds.new("w")
scene.world.use_nodes = True
scene.world.node_tree.nodes["Background"].inputs["Color"].default_value = (*srgb_lin(hexc("#9fb0c4")), 1)
scene.world.node_tree.nodes["Background"].inputs["Strength"].default_value = 0.8

H = 0.305


def flat(name, hexcol, rough=0.85):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    b = m.node_tree.nodes["Principled BSDF"]
    b.inputs["Base Color"].default_value = (*srgb_lin(hexc(hexcol)), 1)
    b.inputs["Roughness"].default_value = rough
    return m


def fur_mat(name, png):
    p = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", png)
    return material(name, bpy.data.images.load(p))


bpy.ops.mesh.primitive_plane_add(size=6, location=(0, 0, 0))
bpy.context.active_object.data.materials.append(flat("ground", "#7d8a6a", 0.95))

M_LEG_OR, M_LEG_CR = fur_mat("leg or", "fur_orange.png"), fur_mat("leg cr", "fur_cream.png")
root = bpy.data.objects.new("root", None)
scene.collection.objects.link(root)
root.location = (0, 0, H)
for ob in (body_ob, head_ob, tail_ob):
    ob.parent = root
    ob.location = Vector(OB_CENTRE[ob.name])

L1 = L2 = 0.14
HOMES = [(0.15, 0.075), (0.15, -0.075), (-0.12, 0.075), (-0.12, -0.075)]


def prim(kind, size, loc, rot, mat):
    if kind == "sphere":
        bpy.ops.mesh.primitive_uv_sphere_add(segments=24, ring_count=12, radius=0.5)
    else:
        bpy.ops.mesh.primitive_cylinder_add(vertices=20, radius=0.5, depth=1)
    o = bpy.context.active_object
    o.scale = size
    o.rotation_mode = "QUATERNION"
    o.rotation_quaternion = rot
    o.location = loc
    o.data.materials.append(mat)
    for p in o.data.polygons:
        p.use_smooth = True
    o.parent = root
    return o


for hip, (hx, hy) in zip(HIPS, HOMES):
    hip = Vector(hip)
    foot = Vector((hx, hy, -H))
    d = foot - hip
    dist = min(d.length, L1 + L2 - 0.004)
    dr = d.normalized()
    c = max(-1, min(1, (L1 * L1 + dist * dist - L2 * L2) / (2 * L1 * dist)))
    a = math.acos(c)
    pole = Vector((-1, 0, 0))
    p = (pole - dr * pole.dot(dr)).normalized()
    knee = hip + dr * (L1 * math.cos(a)) + p * (L1 * math.sin(a))
    end = hip + dr * dist
    up = Vector((0, 0, 1))
    prim("cyl", (0.055, 0.055, 0.14), (hip + knee) / 2, up.rotation_difference((knee - hip).normalized()), M_LEG_OR)
    prim("cyl", (0.04, 0.04, 0.14), (knee + end) / 2, up.rotation_difference((end - knee).normalized()), M_LEG_CR)
    prim("sphere", (0.055, 0.045, 0.032), end + Vector((0.012, 0, 0.012)), Quaternion(), M_LEG_CR)

sun = bpy.data.objects.new("sun", bpy.data.lights.new("sun", "SUN"))
sun.data.energy = 3.0
sun.data.angle = math.radians(8)
sun.rotation_euler = Euler((math.radians(50), 0, math.radians(-35)))
scene.collection.objects.link(sun)
fill = bpy.data.objects.new("fill", bpy.data.lights.new("fill", "SUN"))
fill.data.energy = 0.8
fill.rotation_euler = Euler((math.radians(65), 0, math.radians(150)))
scene.collection.objects.link(fill)


def pose(head_yaw=0.0, head_tilt=0.0, wag=0.0):
    w = Euler((head_tilt, 0, head_yaw)).to_quaternion()
    for ob, piv, q in ((head_ob, HEAD_PIVOT, w), (tail_ob, TAIL_PIVOT, Euler((0, 0, wag)).to_quaternion())):
        ob.rotation_mode = "QUATERNION"
        ob.rotation_quaternion = q
        ob.location = Vector(piv) + q @ (Vector(OB_CENTRE[ob.name]) - Vector(piv))


def shot(fname, target, az, el, dist, lens=50):
    cam = bpy.data.objects.new("cam", bpy.data.cameras.new("cam"))
    scene.collection.objects.link(cam)
    cam.data.lens = lens
    t = Vector(target)
    d = Vector((math.cos(math.radians(el)) * math.cos(math.radians(az)),
                math.cos(math.radians(el)) * math.sin(math.radians(az)), math.sin(math.radians(el))))
    cam.location = t + d * dist
    cam.rotation_euler = (-d).to_track_quat("-Z", "Y").to_euler()
    scene.camera = cam
    scene.render.filepath = os.path.join(OUT, fname.replace(".png", SUFFIX + ".png"))
    bpy.ops.render.render(write_still=True)
    bpy.data.objects.remove(cam)
    print("rendered", fname)


T = (0.02, 0, H - 0.02)
shot("preview.png", T, 35, 18, 1.35)
shot("preview_side.png", T, 90, 8, 1.4)
shot("preview_back.png", T, 215, 25, 1.35)
shot("preview_face.png", (0.27, 0, H + 0.14), 12, 8, 0.75)
shot("preview_head.png", (0.25, 0, H + 0.15), 55, 20, 0.7)
pose(0.75, 0.35, 0.5)
shot("preview_look.png", T, 40, 22, 1.35)