"""Mesh Little ghost: one object replacing the 15-prim ghost (ghost.json). Run headless:

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

Design frame = the prim ghost's root frame: origin at the old head sphere's centre, +x is where it looks, z up.
A gumdrop sheet (head dome flowing into a bell skirt) with soft folds that end in a wavy hem, the underside open a
little like a real sheet, two arm nubs. Eyes and mouth are thin decals that follow the head surface, 4 mm proud, as a
second material. Signed distance fields -> surface nets -> decimate (mirrored in y) -> snapped back to the field,
normals from the field. Exported Z-up (`export_yup=False`) and centred on its bounding box, like walker_mesh.py.
Faces: 0 sheet (#f4f6fb), 1 face (#1b1b24).
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, 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 ""

SHEET = "#f4f6fb"
FACE = "#1b1b24"


# ---------------------------------------------------------------- distance fields (P: (N, 3) array)

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


LOBES = 8            # hem lobes round the skirt; one hangs at the front centre, under the face
HEM_Z = -0.335       # cusp height (where the hem rises between lobes)
HEM_DROP = 0.045     # how much lower each lobe hangs
SHELL = 0.022        # sheet thickness at the open hem


def hem_z(theta):
    return HEM_Z - HEM_DROP * np.abs(np.cos(LOBES / 2 * theta)) ** 0.8


def blob(P):
    """Head + bell, no hem: the sheet's outer surface, extended below the hem."""
    x, y, z = P[:, 0], P[:, 1], P[:, 2]
    th = np.arctan2(y, x)
    rho = np.hypot(x, y)
    # folds: the skirt swells where the lobes hang, more toward the hem
    amp = 0.075 * sstep(-0.08, -0.36, z)
    s = 1 + amp * np.cos(LOBES * th)
    Q = np.stack([x / s, y / s, z], 1)
    head = ell(P, (0, 0, 0.0), (0.158, 0.168, 0.155))
    bell = rcone(Q, (0, 0, -0.03), (0, 0, -0.29), 0.155, 0.222)
    d = smin(head, bell * np.minimum(s, 1), 0.06)
    return d


def arms(P):
    return mirror(lambda Q: rcone(Q, (0.0, 0.12, -0.095), (0.035, 0.218, -0.105), 0.05, 0.04))(P)


def sheet(P):
    d = smin(blob(P), arms(P), 0.045)
    th = np.arctan2(P[:, 1], P[:, 0])
    # open underneath: hollow out the bell from below, leaving a SHELL-thick hem
    inner = smax(blob(P) + SHELL, P[:, 2] + 0.2, 0.05)
    d = smax(d, -inner, 0.008)
    return smax(d, hem_z(th) - P[:, 2], 0.012)


# face decals: an outline in the y-z plane, as seen from the front, extruded along x and clipped to a thin layer
# over the head surface (4 mm proud, 6 mm buried)
EYE_C, EYE_R = (0.058, -0.016), (0.024, 0.034)
MOUTH_C, MOUTH_R = (0.0, -0.07), (0.015, 0.018)


def ellipse2(u, v, c, r):
    """Approximate signed distance to an ellipse in the plane (good enough near the outline)."""
    q = np.stack([(u - c[0]) / r[0], (v - c[1]) / r[1]], 1)
    k0 = np.linalg.norm(q, axis=1)
    k1 = np.linalg.norm(q / np.array(r), axis=1)
    return k0 * (k0 - 1.0) / np.maximum(k1, 1e-9)


def layer(P, outline):
    h = blob(P)
    d = smax(outline, h - 0.004, 0.002)
    d = smax(d, -(h + 0.006), 0.002)
    return smax(d, 0.05 - P[:, 0], 0.002)        # front half only


def face(P):
    eyes = layer(P, mirror(lambda Q: ellipse2(Q[:, 1], Q[:, 2], EYE_C, EYE_R))(P))
    mouth = layer(P, ellipse2(P[:, 1], P[:, 2], MOUTH_C, MOUTH_R))
    return np.minimum(eyes, mouth)


# ---------------------------------------------------------------- meshing (as walker_mesh.py)

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


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):
    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)))
    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)
    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]
    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]


def srgb(h):
    c = [int(h[i:i + 2], 16) / 255 for i in (1, 3, 5)]
    return tuple(v / 12.92 if v <= 0.04045 else ((v + 0.055) / 1.055) ** 2.4 for v in c)


def material(name, hexcol, rough=0.55, emit=0.0):
    m = bpy.data.materials.get(name) or bpy.data.materials.new(name)
    m.use_nodes = True
    b = m.node_tree.nodes["Principled BSDF"]
    b.inputs["Base Color"].default_value = (*srgb(hexcol), 1)
    b.inputs["Roughness"].default_value = rough
    if emit:
        b.inputs["Emission Color"].default_value = (*srgb(hexcol), 1)
        b.inputs["Emission Strength"].default_value = emit
    return m


def assemble(name, parts):
    Vs, Fs, Ns, Ms, off = [], [], [], [], 0
    for i, (V, F, f, mat) in enumerate(parts):
        g = grad(f, V)
        Ns.append(g / np.linalg.norm(g, axis=1, keepdims=True))
        Vs.append(V)
        Fs.append(F + off)
        Ms += [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])
    print(f"{name}: {(np.einsum('ij,ij->i', fn, gc) < 0).sum()} faces against the field, "
          f"{(np.linalg.norm(fn, axis=1) < 1e-8).sum()} degenerate")
    centre = (V.min(0) + V.max(0)) / 2
    V = V - centre
    me = bpy.data.meshes.new(name)
    me.from_pydata(V.tolist(), [], F.tolist())
    for _, _, _, mat in parts:
        me.materials.append(mat)
    me.polygons.foreach_set("material_index", Ms)
    me.polygons.foreach_set("use_smooth", [True] * len(F))
    uv = me.uv_layers.new(name="UVMap")
    fn = np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]])
    dom = np.abs(fn).argmax(1)
    lv = np.empty(len(me.loops), int)
    me.loops.foreach_get("vertex_index", lv)
    keep = np.array([[1, 2], [0, 2], [0, 1]])[np.repeat(dom, 3)]
    uvs = np.take_along_axis(V[lv], keep, 1) + 0.5
    uv.data.foreach_set("uv", uvs.ravel())
    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)
    return ob, centre, V


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)

M_SHEET, M_FACE = material("sheet", SHEET, 0.6), material("face", FACE, 0.35)

sV, sF = decimated(sheet, (-0.3, -0.3, -0.42), (0.3, 0.3, 0.2), 0.004, 1000 * TRIS, "sheet")
fV, fF = decimated(face, (0.08, -0.1, -0.1), (0.18, 0.1, 0.05), 0.0015, 220 * TRIS, "face")
ghost, centre, Vf = assemble("Little ghost", [(sV, sF, sheet, M_SHEET), (fV, fF, face, M_FACE)])
lo, hi = Vf.min(0), Vf.max(0)
print(f"ghost: {len(sF)} + {len(fF)} = {len(ghost.data.polygons)} tris, {len(ghost.data.vertices)} verts, "
      f"size {tuple(round(float(v), 3) for v in hi - lo)}, bbox centre {tuple(round(float(v), 4) for v in centre)} "
      f"from the old root (head centre)")
f0 = lambda p: blob(np.atleast_2d(np.asarray(p, float)))[0]
xs = np.linspace(0.25, 0, 2501)
for nm, (y, z) in (("eye", EYE_C), ("mouth", MOUTH_C)):
    sx = xs[np.argmax(blob(np.stack([xs, np.full_like(xs, y), np.full_like(xs, z)], 1)) < 0)]
    print(f"  {nm} centre on the head surface at x {sx:+.3f}")

if PROBE:
    sys.exit(0)

os.makedirs(OUT, exist_ok=True)
export(ghost, os.path.join(OUT, f"ghost{SUFFIX}.glb"))

if not RENDER:
    sys.exit(0)

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

scene = bpy.context.scene
scene.render.engine = "BLENDER_EEVEE"
scene.eevee.taa_render_samples = 24
scene.render.resolution_x, scene.render.resolution_y = 720, 540
scene.view_settings.view_transform = "Standard"
scene.world = bpy.data.worlds.new("w")
scene.world.use_nodes = True
bg = scene.world.node_tree.nodes["Background"]
bg.inputs["Color"].default_value = (*srgb("#7E8FA6"), 1)
bg.inputs["Strength"].default_value = 0.7

M_GROUND = material("ground", "#4F5E3A", 0.95)
bpy.ops.mesh.primitive_plane_add(size=30)
ground = bpy.context.active_object
ground.data.materials.append(M_GROUND)

HOVER = 0.6          # old root (head centre) height above the ground in the renders
ghost.location = Vector(centre) + Vector((0, 0, HOVER))

# the old prim ghost from ghost.json, beside it for scale
import json
spec = json.load(open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "ghost.json")))
OLD_AT = Vector((0, -0.8, HOVER))
old = []
M_OSHEET, M_OFACE = material("old sheet", SHEET, 0.6), material("old face", FACE, 0.35)
for p in spec["prims"]:
    loc = OLD_AT + Vector(p.get("pos", (0, 0, 0)))
    mat = M_OFACE if p["color"] == FACE else M_OSHEET
    if p["shape"] == "sphere":
        bpy.ops.mesh.primitive_uv_sphere_add(segments=24, ring_count=12, radius=0.5)
        o = bpy.context.active_object
        o.scale = p["size"]
    else:  # the tapered hollow sheet: an open cone frustum
        bpy.ops.mesh.primitive_cone_add(vertices=32, radius1=0.5 * p["size"][0], radius2=0.5 * p["size"][0] *
                                        p["top_size"][0], depth=p["size"][2], end_fill_type="NOTHING")
        o = bpy.context.active_object
        o.modifiers.new("s", "SOLIDIFY").thickness = 0.01
    o.rotation_euler = Euler([math.radians(a) for a in p.get("rot", (0, 0, 0))])
    o.location = loc
    o.data.materials.append(mat)
    for f in o.data.polygons:
        f.use_smooth = True
    old.append(o)


def camera(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
    return cam


sun = bpy.data.objects.new("sun", bpy.data.lights.new("sun", "SUN"))
sun.data.energy = 2.6
sun.data.angle = math.radians(10)
sun.rotation_euler = Euler((math.radians(55), 0, math.radians(-50)))
scene.collection.objects.link(sun)
fill = bpy.data.objects.new("fill", bpy.data.lights.new("fill", "SUN"))
fill.data.energy = 0.7
fill.rotation_euler = Euler((math.radians(60), 0, math.radians(140)))
scene.collection.objects.link(fill)


def shot(fname, target, az, el, dist, lens=50):
    cam = camera(target, az, el, dist, lens)
    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, 0, HOVER - 0.1)
shot("ghost_front.png", T, 15, 12, 1.6)
shot("ghost_34.png", T, 40, 22, 1.8)
shot("ghost_side.png", T, 90, 5, 1.7)
shot("ghost_back.png", T, 200, 20, 1.8)
ground.hide_render = True
shot("ghost_below.png", T, 30, -35, 1.5)
ground.hide_render = False
shot("ghost_far.png", T, 25, 25, 4.0)            # a few metres away, standing eye height-ish
shot("compare_old_new.png", (0, -0.4, HOVER - 0.1), 20, 18, 2.6)