"""Blender import + objective topology/UV metrics + neutral 4-view renders for each benchmark model.glb. No visual scoring.

What it measures (per GLB, after `import_scene.gltf`):
  import: seconds, objects/meshes/materials/images imported, importer warnings caught
  topology: triangles/quads/ngons, non-manifold edges, boundary edges, loose vertices, degenerate (zero-area) faces,
            duplicate vertices (find_doubles, 1e-6), connected mesh parts (shells)
  uv: layer count, faces without UVs, UV verts outside [0,1], summed UV area (≈ 0–1 square coverage ignoring overlap),
      overlapping faces (Blender's uv.select_overlap, C implementation) as % of faces, texel-density consistency
      (coefficient of variation of sqrt(uvArea/faceArea)), island count (skipped above 300k faces)
  topology is reported twice: as delivered (glTF splits vertices along UV/normal seams, so seam edges look non-manifold)
      and after welding coincident vertices (remove_doubles 1e-6) — the welded numbers are the real holes/loose parts.
  bounds: size, centre, lowest point (pivot on floor?)
  render: 4 fixed world-axis views (camera at −Y, +Y, −X, +X; files named front/back/left/right in that order — engines
      differ in which axis they face, so read the labels as camera positions, not the model's front) with EEVEE, the file's
      own PBR materials, a neutral 3-point-ish world light — for human comparison only.

Usage (Blender 4.2+ / 5.x, background):
  blender -b -P mesh_quality_blender.py -- <glb-dir> <out-dir> [--no-render] [--render-only] [--only stem1,stem2]
Writes <out-dir>/mesh-quality.json (merged, resumable) and <out-dir>/renders/<stem>__{front,back,left,right}.png
"""
import bpy, bmesh, json, os, sys, time, math
from mathutils import Vector

argv = sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else []
GLB_DIR, OUT_DIR = argv[0], argv[1]
NO_RENDER = '--no-render' in argv
RENDER_ONLY = '--render-only' in argv
ONLY = set(argv[argv.index('--only') + 1].split(',')) if '--only' in argv else None
RENDER_DIR = os.path.join(OUT_DIR, 'renders'); os.makedirs(RENDER_DIR, exist_ok=True)
OUT = os.path.join(OUT_DIR, 'mesh-quality.json')

def reset_scene():
    bpy.ops.wm.read_factory_settings(use_empty=True)

def poly_area_2d(pts):
    a = 0.0
    for i in range(len(pts)):
        x1, y1 = pts[i]; x2, y2 = pts[(i + 1) % len(pts)]
        a += x1 * y2 - x2 * y1
    return abs(a) / 2

def uv_overlap_faces(meshes):
    """Faces whose UVs overlap other faces, via Blender's uv.select_overlap (edit mode, UV sync select)."""
    total = 0
    for o in meshes:
        if not o.data.uv_layers: continue
        bpy.ops.object.select_all(action='DESELECT'); o.select_set(True); bpy.context.view_layer.objects.active = o
        bpy.ops.object.mode_set(mode='EDIT'); bpy.ops.mesh.select_all(action='SELECT')
        bpy.context.scene.tool_settings.use_uv_select_sync = True
        bpy.ops.uv.select_overlap()
        bm = bmesh.from_edit_mesh(o.data); total += sum(1 for f in bm.faces if f.select)
        bpy.ops.object.mode_set(mode='OBJECT')
    return total

def welded_topology(meshes):
    """Non-manifold/boundary/shells after merging coincident vertices — seam splits no longer count."""
    nonmanifold = boundary = shells = 0
    for o in meshes:
        bm = bmesh.new(); bm.from_mesh(o.data)
        bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=1e-6)
        for e in bm.edges:
            if not e.is_manifold: nonmanifold += 1
            if e.is_boundary: boundary += 1
        seen = set()
        for v in bm.verts:
            if v.index in seen: continue
            shells += 1; stack = [v]
            while stack:
                x = stack.pop()
                if x.index in seen: continue
                seen.add(x.index)
                for e in x.link_edges:
                    for y in e.verts:
                        if y.index not in seen: stack.append(y)
        bm.free()
    return {'nonManifoldEdges': nonmanifold, 'boundaryEdges': boundary, 'shells': shells}

def analyze(path):
    reset_scene()
    t0 = time.time()
    bpy.ops.import_scene.gltf(filepath=path)
    import_s = round(time.time() - t0, 2)
    meshes = [o for o in bpy.data.objects if o.type == 'MESH']
    rec = {'importSeconds': import_s, 'objects': len(bpy.data.objects), 'meshObjects': len(meshes),
           'materials': len(bpy.data.materials), 'images': len([i for i in bpy.data.images if i.name != 'Render Result']),
           'imageSizes': sorted({(i.size[0], i.size[1]) for i in bpy.data.images if i.name != 'Render Result' and i.size[0]}, reverse=True)}
    tris = quads = ngons = 0; nonmanifold = boundary = loose = degenerate = 0; doubles = 0; shells = 0
    faces_total = 0; faces_no_uv = 0; uv_layers = 0
    uv_out = 0; uv_verts = 0; ratios = []; uv_area_sum = 0.0
    bbox_min = Vector((1e9,) * 3); bbox_max = Vector((-1e9,) * 3)
    for o in meshes:
        me = o.data; mw = o.matrix_world
        for v in me.vertices:
            p = mw @ v.co
            bbox_min = Vector(map(min, bbox_min, p)); bbox_max = Vector(map(max, bbox_max, p))
        bm = bmesh.new(); bm.from_mesh(me); bm.transform(mw)
        uv_layers = max(uv_layers, len(me.uv_layers))
        uv = bm.loops.layers.uv.active
        for f in bm.faces:
            n = len(f.verts); faces_total += 1
            if n == 3: tris += 1
            elif n == 4: quads += 1
            else: ngons += 1
            area = f.calc_area()
            if area < 1e-12: degenerate += 1
            if uv is None: faces_no_uv += 1; continue
            pts = [tuple(l[uv].uv) for l in f.loops]
            for (u, v) in pts:
                uv_verts += 1
                if u < -1e-4 or u > 1 + 1e-4 or v < -1e-4 or v > 1 + 1e-4: uv_out += 1
            ua = poly_area_2d(pts); uv_area_sum += ua
            if area > 1e-12 and ua > 1e-12: ratios.append(math.sqrt(ua / area))
        for e in bm.edges:
            if not e.is_manifold: nonmanifold += 1
            if e.is_boundary: boundary += 1
        loose += sum(1 for v in bm.verts if not v.link_edges)
        res = bmesh.ops.find_doubles(bm, verts=bm.verts, dist=1e-6); doubles += len(res['targetmap'])
        # connected parts (shells)
        seen = set()
        for v in bm.verts:
            if v.index in seen: continue
            shells += 1; stack = [v]
            while stack:
                x = stack.pop()
                if x.index in seen: continue
                seen.add(x.index)
                for e in x.link_edges:
                    for y in e.verts:
                        if y.index not in seen: stack.append(y)
        bm.free()
    size = bbox_max - bbox_min
    rec.update({
        'topology': {'faces': faces_total, 'triangles': tris, 'quads': quads, 'ngons': ngons, 'nonManifoldEdges': nonmanifold,
                     'boundaryEdges': boundary, 'looseVertices': loose, 'degenerateFaces': degenerate, 'duplicateVertices': doubles, 'shells': shells,
                     'welded': welded_topology(meshes)},
        'bounds': {'size': [round(c, 4) for c in size], 'centre': [round(c, 4) for c in (bbox_min + bbox_max) / 2],
                   'minY': round(bbox_min.y, 4), 'floorContact': abs(bbox_min.y) <= 0.01 * max(size) if max(size) else None},
    })
    if faces_total and faces_no_uv < faces_total:
        overlap_faces = uv_overlap_faces(meshes)
        mean = sum(ratios) / len(ratios) if ratios else 0
        cv = (math.sqrt(sum((r - mean) ** 2 for r in ratios) / len(ratios)) / mean) if ratios and mean else None
        rec['uv'] = {'layers': uv_layers, 'facesWithoutUv': faces_no_uv, 'uvVerts': uv_verts, 'uvVertsOutside01': uv_out,
                     'uvAreaSum': round(uv_area_sum, 4), 'overlappingFaces': overlap_faces,
                     'overlappingFacesPct': round(100 * overlap_faces / faces_total, 2),
                     'texelDensityCV': round(cv, 3) if cv is not None else None}
        if faces_total <= 300_000:
            # UV islands: faces connected through shared UV coordinates
            rec['uv']['islands'] = uv_islands(meshes)
        else:
            rec['uv']['islands'] = None; rec['uv']['islandsNote'] = 'skipped (>300k faces)'
    else:
        rec['uv'] = {'layers': uv_layers, 'facesWithoutUv': faces_no_uv}
    return rec

def uv_islands(meshes):
    total = 0
    for o in meshes:
        bm = bmesh.new(); bm.from_mesh(o.data); uv = bm.loops.layers.uv.active
        if uv is None: bm.free(); continue
        parent = {}
        def find(x):
            while parent.get(x, x) != x:
                parent[x] = parent.get(parent[x], parent[x]); x = parent[x]
            return x
        def union(a, b):
            ra, rb = find(a), find(b)
            if ra != rb: parent[ra] = rb
        key = {}
        for f in bm.faces:
            for l in f.loops:
                k = (l.vert.index, round(l[uv].uv.x, 5), round(l[uv].uv.y, 5))
                if k in key: union(key[k], f.index)
                else: key[k] = f.index; parent.setdefault(f.index, f.index)
            parent.setdefault(f.index, f.index)
        total += len({find(f.index) for f in bm.faces})
        bm.free()
    return total

def render_views(stem):
    scene = bpy.context.scene
    scene.render.engine = 'BLENDER_EEVEE'
    scene.render.resolution_x = scene.render.resolution_y = 512; scene.render.film_transparent = False
    scene.eevee.taa_render_samples = 16
    scene.view_settings.view_transform = 'Standard'
    world = bpy.data.worlds.new('W'); scene.world = world; world.use_nodes = True
    bg = world.node_tree.nodes.get('Background')
    if bg: bg.inputs[0].default_value = (0.82, 0.82, 0.84, 1.0); bg.inputs[1].default_value = 1.0
    meshes = [o for o in bpy.data.objects if o.type == 'MESH']
    mn = Vector((1e9,) * 3); mx = Vector((-1e9,) * 3)
    for o in meshes:
        for c in o.bound_box:
            p = o.matrix_world @ Vector(c); mn = Vector(map(min, mn, p)); mx = Vector(map(max, mx, p))
    centre = (mn + mx) / 2; radius = max((mx - mn).length / 2, 1e-3)
    # key + fill + rim sun lights so metallic/dark materials still read
    for i, (d, e) in enumerate(((Vector((-1, -1, 1.5)), 3.0), (Vector((1, -0.5, 0.8)), 1.5), (Vector((0.3, 1, 1)), 2.0))):
        ld = bpy.data.lights.new(f'L{i}', 'SUN'); ld.energy = e; ld.angle = 0.6
        lo = bpy.data.objects.new(f'L{i}', ld); scene.collection.objects.link(lo)
        lo.rotation_euler = (-d).to_track_quat('-Z', 'Y').to_euler()
    cam_data = bpy.data.cameras.new('C'); cam_data.type = 'ORTHO'; cam_data.ortho_scale = radius * 2.3
    cam = bpy.data.objects.new('C', cam_data); scene.collection.objects.link(cam); scene.camera = cam
    views = {'front': Vector((0, -1, 0)), 'back': Vector((0, 1, 0)), 'left': Vector((-1, 0, 0)), 'right': Vector((1, 0, 0))}
    out = {}
    for name, d in views.items():
        cam.location = centre + d * radius * 4
        cam.rotation_euler = (d * -1).to_track_quat('-Z', 'Y').to_euler()
        scene.render.filepath = os.path.join(RENDER_DIR, f'{stem}__{name}.png')
        bpy.ops.render.render(write_still=True)
        out[name] = os.path.relpath(scene.render.filepath, OUT_DIR)
    return out

if __name__ == '__main__':
    data = json.load(open(OUT)) if os.path.exists(OUT) else {'blender': bpy.app.version_string, 'runs': {}}
    stems = sorted(f[:-len('.model.glb')] for f in os.listdir(GLB_DIR) if f.endswith('.model.glb'))
    for stem in stems:
        if ONLY and stem not in ONLY: continue
        path = os.path.join(GLB_DIR, stem + '.model.glb')
        if RENDER_ONLY:
            # re-render with the current renderer settings, keep the existing metrics
            reset_scene(); bpy.ops.import_scene.gltf(filepath=path)
            rec = data['runs'].get(stem, {}); rec['renders'] = render_views(stem); rec['renderer'] = 'EEVEE'
            data['runs'][stem] = rec; json.dump(data, open(OUT, 'w'), indent=1); print('rendered', stem, flush=True); continue
        if stem in data['runs'] and (NO_RENDER or data['runs'][stem].get('renders')): continue
        t0 = time.time()
        try:
            rec = analyze(path)
            if not NO_RENDER: rec['renders'] = render_views(stem); rec['renderer'] = 'EEVEE'
            rec['analysisSeconds'] = round(time.time() - t0, 1)
        except Exception as e:
            rec = {'error': str(e)[:300]}
        data['runs'][stem] = rec
        json.dump(data, open(OUT, 'w'), indent=1)
        t = rec.get('topology', {}); u = rec.get('uv', {})
        w = t.get('welded', {})
        print(f"{stem:40s} imp {rec.get('importSeconds')}s faces={t.get('faces')} nonmanifold={t.get('nonManifoldEdges')}/welded {w.get('nonManifoldEdges')} boundary(w)={w.get('boundaryEdges')} shells={t.get('shells')}/w {w.get('shells')} uvOverlap={u.get('overlappingFacesPct')}% uvArea={u.get('uvAreaSum')} cv={u.get('texelDensityCV')} islands={u.get('islands')} {rec.get('analysisSeconds')}s err={rec.get('error')}", flush=True)
    print('wrote', OUT)