# build_optics.py — models, textures, rigs and animates a pair of acetate
# spectacles (plus an eye chart and a display puck) in Blender, saves the
# .blend, exports a GLB for the three.js viewer on demo.grn.dk/optics/ and
# optionally renders a 1080p hero still with Cycles.
#
#   blender -b -noaudio -P build_optics.py -- --out /path/to/outdir [--still] [--samples N]
#
# The film frames are rendered separately from the saved .blend
# (blender -b -noaudio optics.blend -a) so the long job can run in the background.
#
# Timeline (24 fps, 192 frames = 8 s):
#   1-60    assembly: lenses, rims, bridge, pads and temples fly together
#   60-100  the temples unfold
#   100-193 one full turn of the display pivot
# Lens objects carry two shape keys, 'plus' and 'minus', exported as glTF
# morph targets so the viewer can dial a prescription.
import bpy, bmesh, math, sys, os
from mathutils import Vector

argv = sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else []
OUT = os.path.abspath(argv[argv.index('--out') + 1]) if '--out' in argv else os.getcwd()
STILL = '--still' in argv
SAMPLES = int(argv[argv.index('--samples') + 1]) if '--samples' in argv else 48
os.makedirs(OUT, exist_ok=True)

FPS, FRAMES = 24, 192
TAU = math.tau
F_ASSEMBLED, F_OPEN = 60, 100

bpy.ops.wm.read_factory_settings(use_empty=True)
scene = bpy.context.scene
scene.render.fps = FPS
scene.frame_start, scene.frame_end = 1, FRAMES
prefs = bpy.context.preferences.edit
prefs.keyframe_new_interpolation_type = 'BEZIER'

MODEL = bpy.data.collections.new('Optics')      # exported
STAGE = bpy.data.collections.new('Stage')       # render only
scene.collection.children.link(MODEL)
scene.collection.children.link(STAGE)

# ---------------------------------------------------------------- materials
def principled(name, color, metallic=0.0, rough=0.5, coat=0.0, transmission=0.0, ior=1.45,
               emission=None, strength=0.0, alpha=1.0):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    b = m.node_tree.nodes['Principled BSDF']
    b.inputs['Base Color'].default_value = (*color, 1)
    b.inputs['Metallic'].default_value = metallic
    b.inputs['Roughness'].default_value = rough
    b.inputs['Coat Weight'].default_value = coat
    b.inputs['Transmission Weight'].default_value = transmission
    b.inputs['IOR'].default_value = ior
    b.inputs['Alpha'].default_value = alpha
    if emission:
        b.inputs['Emission Color'].default_value = (*emission, 1)
        b.inputs['Emission Strength'].default_value = strength
    return m

def tortoise(name):
    """Procedural tortoiseshell acetate for the Cycles render (the viewer paints
    its own canvas texture; procedural nodes do not survive glTF)."""
    m = principled(name, (0.30, 0.13, 0.03), 0.0, 0.12, coat=1.0)
    nt = m.node_tree; n = nt.nodes; b = n['Principled BSDF']
    tex = n.new('ShaderNodeTexNoise'); tex.inputs['Scale'].default_value = 26
    tex.inputs['Detail'].default_value = 2.0; tex.inputs['Roughness'].default_value = 0.5
    ramp = n.new('ShaderNodeValToRGB')
    ramp.color_ramp.elements[0].position = 0.40; ramp.color_ramp.elements[0].color = (0.035, 0.014, 0.004, 1)
    ramp.color_ramp.elements[1].position = 0.60; ramp.color_ramp.elements[1].color = (0.70, 0.32, 0.06, 1)
    e = ramp.color_ramp.elements.new(0.50); e.color = (0.18, 0.06, 0.015, 1)
    nt.links.new(tex.outputs['Fac'], ramp.inputs['Fac'])
    nt.links.new(ramp.outputs['Color'], b.inputs['Base Color'])
    return m

ACETATE = tortoise('Acetate tortoise')
BLACK   = principled('Acetate black', (0.01, 0.01, 0.012), 0.0, 0.12, coat=1.0)
STEEL   = principled('Steel hinge', (0.62, 0.62, 0.64), 1.0, 0.28)
GOLD    = principled('Gold trim', (0.95, 0.72, 0.32), 1.0, 0.22)
LENS    = principled('Lens glass', (1.0, 1.0, 1.0), 0.0, 0.02, transmission=1.0, ior=1.52)
ACRYLIC = principled('Acrylic stand', (0.95, 0.97, 1.0), 0.0, 0.05, transmission=0.95, ior=1.49)
PAD     = principled('Silicone pad', (0.88, 0.86, 0.82), 0.0, 0.35, transmission=0.4, ior=1.41)
PUCK    = principled('Display puck', (0.02, 0.02, 0.024), 0.0, 0.18, coat=1.0)
BOARD   = principled('Chart board', (0.78, 0.76, 0.70), 0.0, 0.7, emission=(1, 1, 1), strength=0.06)
INK     = principled('Chart ink', (0.02, 0.02, 0.025), 0.0, 0.6)

# ---------------------------------------------------------------- helpers
def link(obj, coll=MODEL):
    for c in obj.users_collection:
        c.objects.unlink(obj)
    coll.objects.link(obj)
    return obj

def smooth(obj, on=True):
    obj.data.polygons.foreach_set('use_smooth', [on] * len(obj.data.polygons))
    obj.data.update()

def bevel(obj, width=0.0006, segments=3, angle=math.radians(30)):
    m = obj.modifiers.new('Bevel', 'BEVEL')
    m.width, m.segments, m.limit_method, m.angle_limit = width, segments, 'ANGLE', angle
    m.harden_normals = True
    return m

def from_bmesh(name, bm, mat, loc=(0, 0, 0), coll=MODEL, do_smooth=True):
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces)
    me = bpy.data.meshes.new(name); bm.to_mesh(me); bm.free()
    o = bpy.data.objects.new(name, me); o.location = loc
    me.materials.append(mat)
    if do_smooth: smooth(o)
    coll.objects.link(o)
    return o

def cylinder(name, r, h, loc, mat, verts=32, bev=0.0004, axis='Z', coll=MODEL):
    bpy.ops.mesh.primitive_cylinder_add(vertices=verts, radius=r, depth=h, location=loc)
    o = bpy.context.object; o.name = name; o.data.materials.append(mat)
    if axis == 'X': o.rotation_euler = (0, math.pi / 2, 0)
    if axis == 'Y': o.rotation_euler = (math.pi / 2, 0, 0)
    smooth(o)
    if bev: bevel(o, min(bev, h * 0.3, r * 0.3))
    return link(o, coll)

def sphere(name, r, loc, mat, seg=32, scale=(1, 1, 1), coll=MODEL):
    bpy.ops.mesh.primitive_uv_sphere_add(segments=seg, ring_count=seg // 2, radius=r, location=loc)
    o = bpy.context.object; o.name = name; o.scale = scale; o.data.materials.append(mat); smooth(o)
    return link(o, coll)

def box(name, size, loc, mat, bev=0.0006, coll=MODEL):
    bpy.ops.mesh.primitive_cube_add(size=1, location=loc)
    o = bpy.context.object; o.name = name; o.scale = size
    bpy.ops.object.transform_apply(scale=True)
    o.data.materials.append(mat); smooth(o); bevel(o, bev, 3)
    return link(o, coll)

def empty(name, loc, coll=MODEL):
    o = bpy.data.objects.new(name, None); o.location = loc; o.empty_display_size = 0.01
    coll.objects.link(o); return o

def parent(child, par):
    bpy.context.view_layer.update()      # fresh matrix_world, else the inverse is stale (identity)
    child.parent = par
    child.matrix_parent_inverse = par.matrix_world.inverted()

def key_loc(obj, frame, loc):
    obj.location = loc; obj.keyframe_insert('location', frame=frame)

def key_rot(obj, frame, rot):
    obj.rotation_mode = 'XYZ'; obj.rotation_euler = rot; obj.keyframe_insert('rotation_euler', frame=frame)

# ---------------------------------------------------------------- lens outline
def outline(cx, n=96, a=0.0265, b=0.0195, mirror=False):
    """Wayfarer-ish lens shape in the XZ plane around x=cx: keystone (wider at
    the top), flatter brow, rounder bottom. Returns [(x, z)] counter-clockwise."""
    pts = []
    for i in range(n):
        t = i * TAU / n
        x, z = a * math.cos(t), b * math.sin(t)
        x *= 1 + 0.14 * (z / b)               # keystone
        if z > 0: z *= 0.86                    # flatter brow
        x += 0.004 * (z / b) ** 2 * (1 if mirror else -1) * 0.5   # slight outer-corner lift
        x = -x if mirror else x
        pts.append((cx + x, z))
    if mirror: pts.reverse()
    return pts

def offset(pts, d):
    """Offset a closed 2-D polygon outwards by d."""
    n = len(pts); out = []
    for i in range(n):
        x0, z0 = pts[i - 1]; x1, z1 = pts[i]; x2, z2 = pts[(i + 1) % n]
        tx, tz = x2 - x0, z2 - z0
        l = math.hypot(tx, tz); nx, nz = tz / l, -tx / l
        out.append((x1 + nx * d, z1 + nz * d))
    # make sure we went outwards (ccw polygon → right-hand normal is outwards)
    def area(p): return sum(p[i][0] * p[(i + 1) % n][1] - p[(i + 1) % n][0] * p[i][1] for i in range(n)) / 2
    if abs(area(out)) < abs(area(pts)):
        out = [(x1 - (ox - x1), z1 - (oz - z1)) for (x1, z1), (ox, oz) in zip(pts, out)]
    return out

def rim(name, pts, width=0.0036, depth=0.0046, mat=ACETATE):
    inner, outer = pts, offset(pts, width)
    bm = bmesh.new(); n = len(pts); h = depth / 2
    of = [bm.verts.new((x, -h, z)) for x, z in outer]
    inf = [bm.verts.new((x, -h, z)) for x, z in inner]
    ob = [bm.verts.new((x, h, z)) for x, z in outer]
    ib = [bm.verts.new((x, h, z)) for x, z in inner]
    for i in range(n):
        j = (i + 1) % n
        bm.faces.new((of[i], of[j], inf[j], inf[i]))
        bm.faces.new((ib[i], ib[j], ob[j], ob[i]))
        bm.faces.new((of[j], of[i], ob[i], ob[j]))
        bm.faces.new((inf[i], inf[j], ib[j], ib[i]))
    o = from_bmesh(name, bm, mat)
    bevel(o, 0.0009, 3)
    return o

def lens(name, pts, K=10, t=0.0018, curve=0.0022, sag=0.0032):
    """Meniscus lens filling the outline, front face towards -Y. Shape keys
    'plus' (thick centre) and 'minus' (thick edge)."""
    cx = sum(x for x, _ in pts) / len(pts); cz = sum(z for _, z in pts) / len(pts)
    n = len(pts)
    bm = bmesh.new()
    rings_f, rings_b, rho_of = [], [], {}
    def surf(y0, sign):
        rings = []
        for k in range(K, 0, -1):
            s = k / K
            ring = []
            for x, z in pts:
                px, pz = cx + (x - cx) * s, cz + (z - cz) * s
                v = bm.verts.new((px, y0 + sign * 0, pz)); rho_of[v] = s; ring.append(v)
            rings.append(ring)
        c = bm.verts.new((cx, y0, cz)); rho_of[c] = 0.0
        return rings, c
    rf, cf = surf(-t / 2, -1)
    rb, cb = surf(t / 2, 1)
    for v, s in rho_of.items():
        v.co.y -= curve * (1 - s * s)        # meniscus: both surfaces bow forward
    for rings, c, flip in ((rf, cf, False), (rb, cb, True)):
        for k in range(K - 1):
            A, B = rings[k], rings[k + 1]
            for i in range(n):
                j = (i + 1) % n
                f = (A[i], A[j], B[j], B[i])
                bm.faces.new(f[::-1] if flip else f)
        last = rings[-1]
        for i in range(n):
            j = (i + 1) % n
            f = (last[i], last[j], c)
            bm.faces.new(f[::-1] if flip else f)
    for i in range(n):
        j = (i + 1) % n
        bm.faces.new((rf[0][j], rf[0][i], rb[0][i], rb[0][j]))
    verts = list(bm.verts)
    rho = [rho_of[v] for v in verts]
    front = [v in set(v2 for r in rf for v2 in r) or v is cf for v in verts]
    o = from_bmesh(name, bm, LENS)
    me = o.data
    o.shape_key_add(name='Basis', from_mix=False)
    plus = o.shape_key_add(name='plus', from_mix=False)
    minus = o.shape_key_add(name='minus', from_mix=False)
    for i, (s, isf) in enumerate(zip(rho, front)):
        if isf: plus.data[i].co.y -= sag * (1 - s * s)    # front bulges more
        else:   minus.data[i].co.y += sag * s * s          # back thickens at the edge
    return o

def sweep(name, path, side, hw, hh, mat, m=14, cap=True):
    """Loft a rounded-rectangle section along a polyline. hw/hh are functions
    of u in [0,1] giving half-width (along `side`) and half-height."""
    bm = bmesh.new(); rings = []
    N = len(path); S = Vector(side)
    for i, p in enumerate(path):
        t = (path[min(i + 1, N - 1)] - path[max(i - 1, 0)]).normalized()
        e1 = (S - t * S.dot(t)).normalized(); e2 = t.cross(e1).normalized()
        u = i / (N - 1); w, h = hw(u), hh(u)
        ring = []
        for k in range(m):
            a = k * TAU / m
            cx, cz = math.cos(a), math.sin(a)
            px = math.copysign(abs(cx) ** 0.55, cx) * w; pz = math.copysign(abs(cz) ** 0.55, cz) * h
            ring.append(bm.verts.new(p + e1 * px + e2 * pz))
        rings.append(ring)
    for A, B in zip(rings, rings[1:]):
        for k in range(m):
            bm.faces.new((A[k], A[(k + 1) % m], B[(k + 1) % m], B[k]))
    if cap:
        bm.faces.new(rings[0][::-1]); bm.faces.new(rings[-1])
    return from_bmesh(name, bm, mat)

def arc_path(p0, p1, bulge, n=16):
    """Points on a parabola from p0 to p1 bowing by `bulge` (a Vector) at the middle."""
    return [Vector(p0).lerp(Vector(p1), u) + Vector(bulge) * (4 * u * (1 - u)) for u in (i / (n - 1) for i in range(n))]

# ---------------------------------------------------------------- the frame
display = empty('Display pivot', (0, 0, 0))
LX, RX = -0.0335, 0.0335
outR = outline(RX); outL = outline(LX, mirror=True)
parts_explode = []   # (object, exploded local offset)

for side, pts, sx in (('R', outR, 1), ('L', outL, -1)):
    r = rim(f'Rim {side}', pts); parent(r, display); parts_explode.append((r, Vector((sx * 0.035, 0, 0))))
    l = lens(f'Lens {side}', pts); parent(l, display); parts_explode.append((l, Vector((0, -0.07, 0))))
    # end piece + hinge on the outer edge
    ex = sx * 0.0625
    ep = box(f'End piece {side}', (0.0085, 0.0046, 0.0095), (ex - sx * 0.0025, 0.0015, 0.0055), ACETATE, bev=0.0012)
    parent(ep, r)
    hinge = cylinder(f'Hinge {side}', 0.0015, 0.0095, (ex + sx * 0.0005, 0.0050, 0.0055), STEEL, verts=24, bev=0.0003)
    parent(hinge, r)
    for k, z in enumerate((0.0015, 0.0095)):
        s = sphere(f'Screw {side}{k}', 0.0011, (ex + sx * 0.0005, 0.0050, z), GOLD, seg=16)
        parent(s, r)
    # nose pad arm + pad
    px = sx * 0.0125
    arm = sweep(f'Pad arm {side}', arc_path((sx * 0.0105, 0.001, 0.000), (px, 0.0105, -0.0045), (sx * 0.002, 0.002, 0)), (0, 0, 1),
                lambda u: 0.0007, lambda u: 0.0007, STEEL, m=10)
    pad = sphere(f'Nose pad {side}', 0.0032, (px + sx * 0.001, 0.0125, -0.0055), PAD, seg=24, scale=(0.55, 0.45, 1.25))
    pad.rotation_euler = (math.radians(-10), math.radians(sx * 22), 0)
    for o in (arm, pad):
        parent(o, display); parts_explode.append((o, Vector((0, 0, -0.035))))
    # temple on its hinge pivot
    H = Vector((ex + sx * 0.0005, 0.0050, 0.0055))
    piv = empty(f'Temple pivot {side}', H); parent(piv, display)
    parts_explode.append((piv, Vector((sx * 0.045, 0.045, 0))))
    path = [Vector((0, y, 0)) for y in (i * 0.095 / 22 for i in range(23))]
    R = 0.022; c = Vector((0, 0.095, -R))
    for i in range(1, 15):
        a = math.radians(95) * i / 14
        path.append(c + Vector((0, R * math.sin(a), R * math.cos(a))))
    end = path[-1]; tdir = (path[-1] - path[-2]).normalized()
    for i in range(1, 6): path.append(end + tdir * 0.005 * i)
    # slight inward toe of the temple so the pair does not look parallel
    for p in path: p.x = -sx * 0.06 * p.y
    temple = sweep(f'Temple {side}', path, (1, 0, 0),
                   lambda u: 0.0021 - 0.0006 * u, lambda u: 0.0058 - 0.0028 * u, ACETATE)
    temple.location = H; parent(temple, piv)
    tip = sweep(f'Temple tip {side}', [Vector((0, y, 0)) for y in (i * 0.012 / 4 for i in range(5))], (1, 0, 0),
                lambda u: 0.0024, lambda u: 0.0050, STEEL, m=10)
    tip.location = H + Vector((-sx * 0.06 * 0.002, 0.002, 0)); parent(tip, piv)
    tip.hide_render = True; tip.hide_viewport = True   # hinge cover, kept simple: hidden
    # fold animation
    fold = (0, 0, sx * math.radians(88))
    key_rot(piv, 1, fold); key_rot(piv, F_ASSEMBLED, fold); key_rot(piv, F_OPEN, (0, 0, 0))

bridge = sweep('Bridge', arc_path((-0.0085, -0.0005, 0.0075), (0.0085, -0.0005, 0.0075), (0, -0.0015, 0.0045)), (0, 1, 0),
               lambda u: 0.0022, lambda u: 0.0030, ACETATE)
parent(bridge, display); parts_explode.append((bridge, Vector((0, 0, 0.04))))
brow = sweep('Brow bar', arc_path((-0.0085, 0.0020, 0.0110), (0.0085, 0.0020, 0.0110), (0, 0.0005, 0.0030)), (0, 1, 0),
             lambda u: 0.0007, lambda u: 0.0007, GOLD, m=10)
parent(brow, display); parts_explode.append((brow, Vector((0, 0, 0.04))))

# stand: acrylic post from the puck up to the bridge, turning with the glasses
stand = cylinder('Acrylic stand', 0.0045, 0.052, (0, 0.006, -0.026), ACRYLIC, verts=32, bev=0.0006)
parent(stand, display)
puck = cylinder('Display puck', 0.085, 0.008, (0, 0, -0.056), PUCK, verts=96, bev=0.0015)
ringg = bpy.ops.mesh.primitive_torus_add(major_radius=0.085, minor_radius=0.0012, major_segments=96, minor_segments=12, location=(0, 0, -0.052))
ringo = bpy.context.object; ringo.name = 'Puck ring'; ringo.data.materials.append(GOLD); smooth(ringo); link(ringo)

# ---------------------------------------------------------------- assembly keyframes
for o, off in parts_explode:
    home = o.location.copy()
    key_loc(o, 1, home + off); key_loc(o, F_ASSEMBLED, home)
    if o.type == 'MESH':
        o.rotation_mode = 'XYZ'
        key_rot(o, 1, (math.radians(35) * (1 if off.x >= 0 else -1), 0, math.radians(20)))
        key_rot(o, F_ASSEMBLED, (0, 0, 0))

prefs.keyframe_new_interpolation_type = 'LINEAR'
key_rot(display, F_OPEN, (0, 0, 0)); key_rot(display, FRAMES + 1, (0, 0, TAU))
prefs.keyframe_new_interpolation_type = 'BEZIER'

# ---------------------------------------------------------------- eye chart
board = box('Chart board', (0.44, 0.004, 0.44), (0, 0.30, 0.02), BOARD, bev=0.002)
ROWS = [('E', 0.060), ('F P', 0.045), ('T O Z', 0.034), ('L P E D', 0.026), ('P E C F D', 0.020), ('E D F C Z P', 0.015), ('F E L O P Z D', 0.011)]
z = 0.185
for i, (txt, size) in enumerate(ROWS):
    cu = bpy.data.curves.new(f'Chart row {i}', 'FONT')
    cu.body = txt; cu.size = size; cu.align_x = 'CENTER'; cu.extrude = 0.0004; cu.space_character = 1.35
    o = bpy.data.objects.new(f'Chart row {i}', cu); o.location = (0, 0.2975, z - size * 0.9)
    o.rotation_euler = (math.pi / 2, 0, 0); cu.materials.append(INK)
    MODEL.objects.link(o)
    z -= size * 1.55 + 0.006

# ---------------------------------------------------------------- stage (render only)
bpy.ops.mesh.primitive_plane_add(size=6, location=(0, 0, -0.06))
floor = bpy.context.object; floor.name = 'Studio floor'
floor.data.materials.append(principled('Studio floor', (0.03, 0.032, 0.036), 0.0, 0.22, coat=0.9))
link(floor, STAGE)

def area(name, loc, energy, color, size, target=(0, 0, 0), coll=STAGE, shape='SQUARE', size_y=None):
    ld = bpy.data.lights.new(name, 'AREA'); ld.energy, ld.color, ld.size = energy, color, size
    if size_y: ld.shape = 'RECTANGLE'; ld.size_y = size_y
    lo = bpy.data.objects.new(name, ld); lo.location = loc; coll.objects.link(lo)
    d = Vector(target) - Vector(loc)
    lo.rotation_euler = d.to_track_quat('-Z', 'Y').to_euler()
    return lo
area('Key',     ( 0.45, -0.55, 0.55), 55, (1.0, 0.95, 0.88), 0.9, size_y=0.35)
area('Rim',     (-0.60,  0.25, 0.40), 30, (0.80, 0.88, 1.0), 0.5)
area('Fill',    (-0.40, -0.60, 0.10), 12, (1.0, 1.0, 1.0), 0.8)
area('Softbox', ( 0.00, -0.30, 0.90), 25, (1.0, 1.0, 1.0), 1.4, size_y=0.25)
world = bpy.data.worlds.new('World'); scene.world = world; world.use_nodes = True
bg = world.node_tree.nodes['Background']; bg.inputs['Color'].default_value = (0.012, 0.013, 0.016, 1); bg.inputs['Strength'].default_value = 1.0

cam_data = bpy.data.cameras.new('Camera'); cam_data.lens = 70
cam_data.dof.use_dof = True; cam_data.dof.aperture_fstop = 2.8; cam_data.dof.focus_distance = 0.40
cam = bpy.data.objects.new('Camera', cam_data); STAGE.objects.link(cam)
cam.location = (0.05, -0.40, 0.045)
cam.rotation_euler = (Vector((0, 0, -0.002)) - Vector(cam.location)).to_track_quat('-Z', 'Y').to_euler()
scene.camera = cam

# ---------------------------------------------------------------- render settings
scene.render.engine = 'CYCLES'
scene.cycles.device = 'CPU'
scene.cycles.samples = SAMPLES
scene.cycles.use_denoising = True
scene.cycles.use_adaptive_sampling = True
scene.cycles.max_bounces = 8; scene.cycles.transmission_bounces = 8; scene.cycles.glossy_bounces = 4
scene.cycles.caustics_reflective = False; scene.cycles.caustics_refractive = False
scene.render.resolution_x, scene.render.resolution_y = 1280, 720
scene.render.image_settings.file_format = 'PNG'
scene.render.filepath = os.path.join(OUT, 'frames', 'optics_')
scene.view_settings.view_transform = 'AgX'
scene.view_settings.look = 'AgX - Punchy'

# ---------------------------------------------------------------- save, export, still
blend = os.path.join(OUT, 'optics.blend')
bpy.ops.wm.save_as_mainfile(filepath=blend)

# Freeze modifiers and text into plain meshes for export (the exporter drops
# shape keys on objects it has to apply modifiers to; the lenses have none, but
# freezing everything keeps the GLB independent of exporter behaviour).
dg = bpy.context.evaluated_depsgraph_get()
for o in list(MODEL.all_objects):
    if o.type == 'FONT':
        me = bpy.data.meshes.new_from_object(o.evaluated_get(dg)); me.materials.append(INK)
        n = bpy.data.objects.new(o.name + ' mesh', me); n.matrix_world = o.matrix_world; n.parent = o.parent
        MODEL.objects.link(n); bpy.data.objects.remove(o)
    elif o.type == 'MESH' and o.modifiers and not o.data.shape_keys:
        me = bpy.data.meshes.new_from_object(o.evaluated_get(dg))
        for m in o.data.materials: me.materials.append(m)
        o.modifiers.clear(); o.data = me; smooth(o)

bpy.ops.object.select_all(action='DESELECT')
for o in MODEL.all_objects:
    if not o.hide_render: o.select_set(True)
bpy.ops.export_scene.gltf(
    filepath=os.path.join(OUT, 'optics.glb'), export_format='GLB',
    use_selection=True, export_apply=False, export_yup=True,
    export_animations=True, export_animation_mode='SCENE', export_force_sampling=True,
    export_morph=True, export_morph_normal=True,
    export_cameras=False, export_lights=False, export_extras=False,
)
tris = sum(len(o.data.polygons) for o in MODEL.all_objects if o.type == 'MESH')
print(f'OPTICS: {len(list(MODEL.all_objects))} objects, {tris} polygons exported')

if STILL:
    bpy.ops.wm.open_mainfile(filepath=blend)
    scene = bpy.context.scene
    scene.frame_set(150)
    scene.render.resolution_x, scene.render.resolution_y = 1920, 1080
    scene.cycles.samples = max(SAMPLES, 128)
    scene.render.filepath = os.path.join(OUT, 'optics-hero.png')
    bpy.ops.render.render(write_still=True)
    print('OPTICS: hero still written')