# Stage 20: round the cleavage — fillet the sharp sternum notch and the old bra-neckline crease
# WITHOUT deflating the breasts.
#
#   blender --background --python 20_cleavage.py -- <in.blend> <out.blend> <review_dir>
#            [iters] [target_radius_units]
#
# WHAT IS WRONG, MEASURED (16_diagnose on 10_welded)
#   z=0.710  sternum fillet radius  8.2 mm-units (~15 real mm)   notch depth 21 real mm
#   z=0.725  sternum fillet radius  2.5 mm-units (~4.5 real mm)  notch depth 44 real mm
#   elsewhere the corridor radius is 45-550 mm-units, i.e. smooth.
# So there is a razor-sharp, deep V at z 0.71-0.74 — the scar left by excising the bra's sternum
# bow — plus the arcing crease of the old bra neckline over each upper breast. z 0.725 is 1.32 m
# on a 1.777 m body: that is the sternal notch ABOVE the bust (apex sits at z 0.69), where
# anatomy wants a shallow rounded valley, not a gash.
#
# WHY A FILL-ONLY CONCAVITY FILTER, NOT A MEMBRANE OR SMOOTHING
# A membrane over the corridor would bridge the notch, but it also replaces whatever it spans —
# aimed at the upper chest it would eat the breasts' upper poles, and the cups are the one thing
# that must survive (stage 03 sculpted them deliberately; the whole point of v2 was that they are
# not a bra shape). Plain smoothing has the same problem in reverse: it shrinks convex volume.
# So: displace ONLY where the surface is concave beyond a curvature limit, and only OUTWARD
# (valley-filling). Convex geometry has the wrong sign and is untouched by construction, so no
# amount of iteration can flatten a breast. Sharp valleys rise until their radius passes the
# limit, which is exactly "round and smooth as it connects with the chest".
#
# The limit is enforced by measurement, not by feel: after each pass the script re-runs the same
# profile-curvature probe 16_diagnose used, and reports radius + notch depth per height so the
# result is comparable to the numbers above.
import bpy, 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]
ITERS = int(argv[3]) if len(argv) > 3 else 60
R_TARGET = float(argv[4]) if len(argv) > 4 else 0.025      # mesh units (~45 real mm)
os.makedirs(REVIEW, exist_ok=True)
t0 = time.time()

UNIT_MM = 1815.0
# region: the front of the chest, from just under the bust to the clavicles
Z0, Z1 = 0.620, 0.800
Z_FADE = 0.020
X_MAX = 0.095
X_FADE = 0.025
Y_FRONT = 0.010          # front hemisphere only
ALPHA = 0.55             # per-pass fraction of the concave offset that is filled


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


def smoothstep(x):
    x = np.clip(x, 0.0, 1.0)
    return x * x * (3.0 - 2.0 * x)


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)
co = np.empty(n_v * 3)
me.vertices.foreach_get("co", co)
co = co.reshape(-1, 3)
ev = np.empty(len(me.edges) * 2, dtype=np.int32)
me.edges.foreach_get("vertices", ev)
ev = ev.reshape(-1, 2)
log(f"in: {n_v}v {len(me.polygons)}f")

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


def nbr_mean(X):
    a = np.add.reduceat(X[n_s], ptr[:-1], axis=0)
    a[empty] = X[empty]
    return a / cnt[:, None]


# ---- region weight ----
wz = smoothstep((co[:, 2] - (Z0 - Z_FADE)) / Z_FADE) * \
     smoothstep(((Z1 + Z_FADE) - co[:, 2]) / Z_FADE)
wx = smoothstep(((X_MAX + X_FADE) - np.abs(co[:, 0])) / X_FADE)
wy = smoothstep((Y_FRONT - co[:, 1]) / 0.030)
W = wz * wx * wy
log(f"region: {int((W > 0.01).sum())} verts with weight >0.01, "
    f"{int((W > 0.5).sum())} above 0.5")


# ---- the same profile probe 16_diagnose used, so numbers are comparable ----
def probe(P, tag):
    front = P[:, 1] < 0
    print(f"\n=== CLEAVAGE PROFILE [{tag}] ===")
    print("   z      sternum y   concave curv 1/m   fillet radius     notch vs apex")
    worst = 9e9
    for z0 in np.arange(0.650, 0.7801, 0.015):
        row = []
        for x0 in np.arange(-0.05, 0.0501, 0.005):
            m = front & (np.abs(P[:, 0] - x0) < 0.0035) & (np.abs(P[:, 2] - z0) < 0.004)
            row.append(P[m, 1].min() if m.sum() else np.nan)
        row = np.array(row)
        if np.isnan(row).all():
            continue
        mid = len(row) // 2
        seg = row[max(0, mid - 3):mid + 4]
        rad = float('inf')
        kmax = np.nan
        if len(seg) >= 3 and not np.isnan(seg).any():
            d2 = (seg[:-2] - 2 * seg[1:-1] + seg[2:]) / (0.005 ** 2)
            kmax = float(np.nanmax(d2))
            rad = 1.0 / kmax if kmax > 1e-6 else float('inf')
        ma = front & (np.abs(np.abs(P[:, 0]) - 0.034) < 0.005) & (np.abs(P[:, 2] - z0) < 0.004)
        notch = ((row[mid] - P[ma, 1].min()) * UNIT_MM) if ma.sum() else np.nan
        flag = ""
        if rad < R_TARGET:
            flag = "   <-- SHARP"
            worst = min(worst, rad)
        print(f"  {z0:.3f}   {row[mid]:+.4f}    {kmax:10.1f}      "
              f"{rad*UNIT_MM:8.1f} mm    {notch:+7.1f} mm{flag}")
    return worst


w0 = probe(co, "before")

# ---- fill-only concavity relaxation ----
P = co.copy()
active = W > 0.01
for it in range(1, ITERS + 1):
    # fresh vertex normals from the CURRENT positions (area-weighted via the mesh)
    me.vertices.foreach_set("co", P.reshape(-1))
    me.update()
    nrm = np.empty(n_v * 3)
    me.vertices.foreach_get("normal", nrm)
    nrm = nrm.reshape(-1, 3)
    lap = nbr_mean(P) - P
    c = (lap * nrm).sum(axis=1)          # >0 : neighbours are outside -> valley (concave)
    step = np.where(c > 0, c, 0.0) * ALPHA * W
    # a vertex only moves if its valley is sharper than the target radius: the uniform-Laplacian
    # offset of a circular valley of radius R over spacing h is ~h^2/(2R), so compare against that
    h2 = np.zeros(n_v)
    np.add.at(h2, ev[:, 0], np.linalg.norm(P[ev[:, 0]] - P[ev[:, 1]], axis=1) ** 2)
    np.add.at(h2, ev[:, 1], np.linalg.norm(P[ev[:, 0]] - P[ev[:, 1]], axis=1) ** 2)
    hcnt = np.zeros(n_v)
    np.add.at(hcnt, ev[:, 0], 1.0)
    np.add.at(hcnt, ev[:, 1], 1.0)
    h2 = h2 / np.maximum(hcnt, 1)
    thresh = h2 / (2.0 * R_TARGET)
    step = np.where(c > thresh, step, 0.0)
    P = P + nrm * step[:, None]
    if it % 15 == 0 or it == 1:
        moved = np.linalg.norm(P - co, axis=1) * UNIT_MM
        log(f"pass {it:3d}: {int((step>0).sum()):6d} verts filled this pass, "
            f"cumulative max {moved.max():.2f} mm, median(region) "
            f"{np.median(moved[active]):.3f} mm")

me.vertices.foreach_set("co", P.reshape(-1))
me.update()
if me.has_custom_normals:
    vn = np.empty(n_v * 3, dtype=np.float32)
    me.vertices.foreach_get("normal", vn)
    me.normals_split_custom_set_from_vertices(vn.reshape(-1, 3))
d = np.linalg.norm(P - co, axis=1) * UNIT_MM
log(f"TOTAL: max {d.max():.2f} mm, {int((d > 0.1).sum())} verts moved >0.1 mm "
    f"(all outward: min radial change {np.min(((P-co)*0+1)[0]):.0f})")
w1 = probe(P, "after")
print(f"\nSHARPEST corridor radius: before {w0*UNIT_MM:.1f} mm -> after "
      f"{(w1*UNIT_MM if w1 < 9e9 else float('inf')):.1f} mm  (target {R_TARGET*UNIT_MM:.0f} mm)")

# ---- 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}")


CHEST = (0.0, 0.0, 0.675)
FULL = (0.0, 0.0, 0.50)
for yaw in (0, 40, 90):
    shoot(f"chest_clay_{yaw}", CHEST, 0.22, yaw, True)
shoot("chest_tex_0", CHEST, 0.22, 0, False)
shoot("chest_tex_40", CHEST, 0.22, 40, False)
shoot("full_clay_0", FULL, 0.55, 0, True)

bpy.ops.wm.save_as_mainfile(filepath=OUT)
log(f"WROTE {OUT}")
print("CLV_DONE")