# Stage 25: the black dashes/speckles. Diagnose FIRST, then fix only what the numbers show.
#
#   blender --background --python 25_speckle.py -- <in.blend> <out.blend> <review_dir> [--apply]
#
# WHY NOT "FILL THE HOLES" (that has now failed twice)
# The mesh has 830 boundary edges in 689 loops of 3-5 edges, plus 1649 non-manifold edges. Both
# bmesh.ops.holes_fill (per loop) and the edit-mode mesh.fill_holes operator left the count at
# exactly 830 — they refused every one. A 3-edge run of boundary that cannot be filled is not a
# closed triangular hole; it is an OPEN CHAIN, i.e. these are dangling flaps and slivers hanging
# off the surface, not perforations. Roughly 2 non-manifold edges per boundary edge fits that
# reading. So filling is the wrong verb; the candidates are flipped winding (a backfacing triangle
# renders black in EEVEE, which is exactly what a "black dash" looks like) and tiny stray shards.
#
# Tests, in order, all reported before anything is modified:
#   1. FLIPPED FACES  — face normal against the locally smoothed vertex normal. >90 deg apart
#      means the triangle faces inward and will render black.
#   2. STRAY SHARDS   — connected components of the face graph. The body is one component; small
#      components are debris and can be deleted outright.
#   3. SLIVERS        — near-zero-area and extreme-aspect triangles, which shade unpredictably.
# --apply then: recalculate consistent winding, delete shards under SHARD_MAX faces, dissolve
# degenerate slivers. Nothing here moves a vertex, so the sculpt and the heals are untouched.
import bpy, bmesh, sys, os, math, time
import numpy as np
from mathutils import Vector

argv = sys.argv[sys.argv.index("--") + 1:]
BLEND, OUT, REVIEW = argv[0], argv[1], argv[2]
APPLY = "--apply" in argv
os.makedirs(REVIEW, exist_ok=True)
t0 = time.time()
SHARD_MAX = 200          # faces; the body itself is ~1.7M so this is unambiguous debris
UNIT_MM = 1815.0


def log(m):
    print(f"[spk {time.time()-t0:6.1f}s] {m}", flush=True)


bpy.ops.wm.open_mainfile(filepath=BLEND)
ob = max([o for o in bpy.data.objects if o.type == 'MESH'],
         key=lambda o: len(o.data.vertices))
me = ob.data
n_v = len(me.vertices)
n_f = len(me.polygons)
log(f"in: {n_v}v {n_f}f custom_normals={me.has_custom_normals}")

co = np.empty(n_v * 3)
me.vertices.foreach_get("co", co)
co = co.reshape(-1, 3)
vn = np.empty(n_v * 3)
me.vertices.foreach_get("normal", vn)
vn = vn.reshape(-1, 3)
fn = np.empty(n_f * 3)
me.polygons.foreach_get("normal", fn)
fn = fn.reshape(-1, 3)
fc = np.empty(n_f * 3)
me.polygons.foreach_get("center", fc)
fc = fc.reshape(-1, 3)
ev = np.empty(len(me.edges) * 2, dtype=np.int32)
me.edges.foreach_get("vertices", ev)
ev = ev.reshape(-1, 2)

# smoothed vertex normal field, as the reference for "which way is out"
order = np.concatenate([ev[:, 0], ev[:, 1]])
nbr = np.concatenate([ev[:, 1], ev[:, 0]])
srt = np.argsort(order, kind="stable")
o_s, n_s = order[srt], nbr[srt]
ptr = np.searchsorted(o_s, np.arange(n_v + 1))
cnt = np.maximum(np.diff(ptr), 1)
empty = np.diff(ptr) == 0
N = vn.copy()
for _ in range(8):
    a = np.add.reduceat(N[n_s], ptr[:-1], axis=0)
    a[empty] = N[empty]
    N = a / cnt[:, None]
    N /= np.maximum(np.linalg.norm(N, axis=1, keepdims=True), 1e-12)

lv = np.empty(len(me.loops), dtype=np.int32)
me.loops.foreach_get("vertex_index", lv)
ls = np.empty(n_f, dtype=np.int32)
me.polygons.foreach_get("loop_start", ls)
lt = np.empty(n_f, dtype=np.int32)
me.polygons.foreach_get("loop_total", lt)
# reference normal per face = mean of its verts' smoothed normals
acc = np.zeros((n_f, 3))
for k in range(int(lt.max())):
    sel = lt > k
    acc[sel] += N[lv[ls[sel] + k]]
acc /= np.maximum(np.linalg.norm(acc, axis=1, keepdims=True), 1e-12)
fdot = (fn * acc).sum(axis=1)
flipped = fdot < 0.0
print(f"FLIPPED FACES: {int(flipped.sum())} of {n_f} "
      f"({100.0*flipped.mean():.4f}%)  [dot<0 vs smoothed field]")
print(f"  nearly-perpendicular (|dot|<0.2): {int((np.abs(fdot)<0.2).sum())}")
if flipped.sum():
    zs = fc[flipped, 2]
    hist, edges = np.histogram(zs, bins=12)
    print("  flipped-face z histogram:")
    for c, lo, hi in zip(hist, edges[:-1], edges[1:]):
        if c:
            print(f"     z {lo:.3f}-{hi:.3f}: {c}")

# ---- stray shards ----
bm = bmesh.new()
bm.from_mesh(me)
bm.faces.ensure_lookup_table()
seen = np.zeros(len(bm.faces), dtype=bool)
comps = []
from collections import deque
for f0 in bm.faces:
    if seen[f0.index]:
        continue
    q = deque([f0])
    seen[f0.index] = True
    size = 0
    members = []
    while q:
        f = q.popleft()
        size += 1
        members.append(f)
        for e in f.edges:
            for g in e.link_faces:
                if not seen[g.index]:
                    seen[g.index] = True
                    q.append(g)
    comps.append((size, members))
comps.sort(key=lambda t: -t[0])
print(f"FACE COMPONENTS: {len(comps)} (largest {[c[0] for c in comps[:5]]})")
shards = [c for c in comps if c[0] <= SHARD_MAX]
print(f"  shards <= {SHARD_MAX} faces: {len(shards)} components, "
      f"{sum(c[0] for c in shards)} faces total")
for size, mem in shards[:10]:
    ctr = Vector((0, 0, 0))
    for f in mem:
        ctr += f.calc_center_median()
    ctr /= len(mem)
    print(f"     shard {size:4d} faces at ({ctr.x:+.3f},{ctr.y:+.3f},{ctr.z:+.3f})")

areas = np.array([f.calc_area() for f in bm.faces])
tiny = areas < (1e-5 ** 2)
print(f"SLIVERS: {int(tiny.sum())} faces with area < (0.01 mm-unit)^2; "
      f"min area {areas.min():.3e}, p01 {np.percentile(areas,1):.3e}")

if APPLY:
    # delete shards
    if shards:
        geom = [f for _, mem in shards for f in mem]
        bmesh.ops.delete(bm, geom=geom, context='FACES')
        log(f"deleted {len(geom)} shard faces")
    # drop degenerate slivers
    bm.faces.ensure_lookup_table()
    dgn = [f for f in bm.faces if f.calc_area() < 1e-12]
    if dgn:
        bmesh.ops.delete(bm, geom=dgn, context='FACES')
        log(f"deleted {len(dgn)} degenerate faces")
    bmesh.ops.delete(bm, geom=[v for v in bm.verts if not v.link_faces], context='VERTS')
    bm.to_mesh(me)
    bm.free()
    me.update()
    # consistent winding — this is what actually removes backfacing black dashes
    bpy.context.view_layer.objects.active = ob
    ob.select_set(True)
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.mesh.normals_make_consistent(inside=False)
    bpy.ops.object.mode_set(mode='OBJECT')
    me.update()
    log(f"recalculated winding; now {len(me.vertices)}v {len(me.polygons)}f")
    # re-measure
    n_f2 = len(me.polygons)
    fn2 = np.empty(n_f2 * 3)
    me.polygons.foreach_get("normal", fn2)
    fn2 = fn2.reshape(-1, 3)
    lv2 = np.empty(len(me.loops), dtype=np.int32)
    me.loops.foreach_get("vertex_index", lv2)
    ls2 = np.empty(n_f2, dtype=np.int32)
    me.polygons.foreach_get("loop_start", ls2)
    lt2 = np.empty(n_f2, dtype=np.int32)
    me.polygons.foreach_get("loop_total", lt2)
    n_v2 = len(me.vertices)
    vn2 = np.empty(n_v2 * 3)
    me.vertices.foreach_get("normal", vn2)
    vn2 = vn2.reshape(-1, 3)
    acc2 = np.zeros((n_f2, 3))
    for k in range(int(lt2.max())):
        sel = lt2 > k
        acc2[sel] += vn2[lv2[ls2[sel] + k]]
    acc2 /= np.maximum(np.linalg.norm(acc2, axis=1, keepdims=True), 1e-12)
    fl2 = ((fn2 * acc2).sum(axis=1) < 0)
    print(f"FLIPPED FACES after make_consistent: {int(fl2.sum())} of {n_f2}")

    # make_consistent propagates orientation across shared edges, so the 1649 non-manifold edges
    # block it — it only got 1511 down to 1198. The remainder are individually wrong relative to
    # their own neighbourhood, so reverse exactly those: it moves no vertex and each reversal
    # brings that triangle into agreement with the faces around it.
    for rnd in range(3):
        n_f3 = len(me.polygons)
        fn3 = np.empty(n_f3 * 3)
        me.polygons.foreach_get("normal", fn3)
        fn3 = fn3.reshape(-1, 3)
        lv3 = np.empty(len(me.loops), dtype=np.int32)
        me.loops.foreach_get("vertex_index", lv3)
        ls3 = np.empty(n_f3, dtype=np.int32)
        me.polygons.foreach_get("loop_start", ls3)
        lt3 = np.empty(n_f3, dtype=np.int32)
        me.polygons.foreach_get("loop_total", lt3)
        nv3 = len(me.vertices)
        vv3 = np.empty(nv3 * 3)
        me.vertices.foreach_get("normal", vv3)
        vv3 = vv3.reshape(-1, 3)
        ac3 = np.zeros((n_f3, 3))
        for k in range(int(lt3.max())):
            sel = lt3 > k
            ac3[sel] += vv3[lv3[ls3[sel] + k]]
        ac3 /= np.maximum(np.linalg.norm(ac3, axis=1, keepdims=True), 1e-12)
        bad = np.nonzero((fn3 * ac3).sum(axis=1) < 0)[0]
        if len(bad) == 0:
            log(f"reversal round {rnd}: none left")
            break
        bm2 = bmesh.new()
        bm2.from_mesh(me)
        bm2.faces.ensure_lookup_table()
        bmesh.ops.reverse_faces(bm2, faces=[bm2.faces[int(i)] for i in bad])
        bm2.to_mesh(me)
        bm2.free()
        me.update()
        log(f"reversal round {rnd}: reversed {len(bad)} faces")
    n_f4 = len(me.polygons)
    fn4 = np.empty(n_f4 * 3)
    me.polygons.foreach_get("normal", fn4)
    fn4 = fn4.reshape(-1, 3)
    lv4 = np.empty(len(me.loops), dtype=np.int32)
    me.loops.foreach_get("vertex_index", lv4)
    ls4 = np.empty(n_f4, dtype=np.int32)
    me.polygons.foreach_get("loop_start", ls4)
    lt4 = np.empty(n_f4, dtype=np.int32)
    me.polygons.foreach_get("loop_total", lt4)
    vv4 = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("normal", vv4)
    vv4 = vv4.reshape(-1, 3)
    ac4 = np.zeros((n_f4, 3))
    for k in range(int(lt4.max())):
        sel = lt4 > k
        ac4[sel] += vv4[lv4[ls4[sel] + k]]
    ac4 /= np.maximum(np.linalg.norm(ac4, axis=1, keepdims=True), 1e-12)
    print(f"FLIPPED FACES final: {int(((fn4*ac4).sum(axis=1) < 0).sum())} of {n_f4}")
    n_v2 = len(me.vertices)
    vnn = np.empty(n_v2 * 3, dtype=np.float32)
    me.vertices.foreach_get("normal", vnn)
    me.normals_split_custom_set_from_vertices(vnn.reshape(-1, 3))
    # SAVE BEFORE RENDERING. The render helper swaps clay into the material slots, and 18/20/25
    # originally saved afterwards — which is how 23_cleavage.blend lost its texture wiring and
    # made 21_tone.py abort. Saving first means the file on disk always keeps the real material.
    bpy.ops.wm.save_as_mainfile(filepath=OUT)
    log(f"WROTE {OUT}")
else:
    bm.free()

# ---- renders ----
scn = bpy.context.scene
wd = bpy.data.worlds.new("W")
wd.color = (0.22, 0.22, 0.24)
scn.world = wd
key = bpy.data.objects.new("Key", bpy.data.lights.new("Key", 'SUN'))
key.data.energy = 3.0
key.data.use_shadow = False
bpy.context.collection.objects.link(key)
fl_ = bpy.data.objects.new("Fill", bpy.data.lights.new("Fill", 'SUN'))
fl_.data.energy = 1.0
fl_.data.use_shadow = False
bpy.context.collection.objects.link(fl_)
cam = bpy.data.objects.new("Cam", bpy.data.cameras.new("Cam"))
cam.data.lens = 85
bpy.context.collection.objects.link(cam)
scn.camera = cam
scn.render.engine = 'BLENDER_EEVEE' if bpy.app.version >= (4, 2) else 'BLENDER_EEVEE_NEXT'
scn.render.resolution_x = scn.render.resolution_y = 1000
clay = bpy.data.materials.new("Clay")
clay.use_nodes = True
clay.node_tree.nodes["Principled BSDF"].inputs["Base Color"].default_value = (0.62, 0.60, 0.58, 1)
clay.node_tree.nodes["Principled BSDF"].inputs["Roughness"].default_value = 0.45
orig = [ms.material for ms in ob.material_slots]


def shoot(tag, ctr, span, yaw_deg, use_clay=True):
    for i, ms in enumerate(ob.material_slots):
        ms.material = clay if use_clay else orig[i]
    yaw = math.radians(yaw_deg)
    dist = span * 3.0
    cam.location = Vector(ctr) + Vector((math.sin(yaw) * dist, -math.cos(yaw) * dist, 0.02))
    cam.rotation_euler = (Vector(ctr) - cam.location).to_track_quat('-Z', 'Y').to_euler()
    key.rotation_euler = (math.radians(62), 0, math.radians(35 + yaw_deg))
    fl_.rotation_euler = (math.radians(75), 0, math.radians(yaw_deg - 110))
    scn.render.filepath = os.path.abspath(os.path.join(REVIEW, f"{tag}.png"))
    bpy.ops.render.render(write_still=True)
    log(f"render {tag}")


shoot("chest_clay_0", (0.0, 0.0, 0.675), 0.22, 0)
shoot("chest_clay_40", (0.0, 0.0, 0.675), 0.22, 40)
shoot("chest_tex_0", (0.0, 0.0, 0.675), 0.22, 0, False)
shoot("hip_clay_0", (0.0, 0.0, 0.53), 0.22, 0)
shoot("full_clay_0", (0.0, 0.0, 0.50), 0.55, 0)
shoot("full_tex_0", (0.0, 0.0, 0.50), 0.55, 0, False)
print("SPK_DONE")