VoxelForge/Standalone/Research/WorldAssemblerM3/Linux37/RunMouthFragmentAuditL37.py

#!/usr/bin/env python3
"""Compare the FI fade-ring crease against L35 using the same surface samples."""
import argparse
import importlib.util
import math
import sys
from pathlib import Path

ROOT = Path(__file__).resolve().parents[4]
RING_PATH = ROOT / "Standalone/Research/WorldAssemblerM3/Linux32/AuditRingTerraceL32.py"
SPEC = importlib.util.spec_from_file_location("ring_profile_l32", RING_PATH)
RING = importlib.util.module_from_spec(SPEC)
SPEC.loader.exec_module(RING)

CENTER = (6.298233, 1.137233)
INNER, OUTER, STEP = 5.5, 9.5, 0.125
HEIGHT_MIN, HEIGHT_MAX, PREFERRED = -80.0, -40.0, -56.0
# Chosen before the L37 mouth edit: a 20-degree normal jump over one 0.125 m
# sample interval is visible at the standard camera distance.
NORMAL_LIMIT_DEG = 20.0
COVERAGE_MIN = 0.90
FRAGMENT_MAX_TRIANGLES = 64
FRAGMENT_MIN_AREA_M2 = 0.0001
FRAGMENT_MAX_AREA_M2 = 0.05
FRAGMENT_MAX_DISTANCE_M = 12.0


def small_components(triangles, center=None):
    parent = list(range(len(triangles)))
    sizes = [1] * len(triangles)
    owners = {}

    def find(index):
        while parent[index] != index:
            parent[index] = parent[parent[index]]
            index = parent[index]
        return index

    def union(a, b):
        a, b = find(a), find(b)
        if a == b:
            return
        if sizes[a] < sizes[b]:
            a, b = b, a
        parent[b] = a
        sizes[a] += sizes[b]

    for index, (_, _, vertices) in enumerate(triangles):
        for point, _ in vertices:
            key = tuple(round(value, 5) for value in point)
            other = owners.setdefault(key, index)
            union(index, other)

    components = {}
    for index, (_, _, vertices) in enumerate(triangles):
        root = find(index)
        item = components.setdefault(root, {"count": 0, "area": 0.0, "sum": [0.0] * 3})
        item["count"] += 1
        points = [vertex[0] for vertex in vertices]
        for point in points:
            for axis in range(3):
                item["sum"][axis] += point[axis]
        ab = tuple(points[1][axis] - points[0][axis] for axis in range(3))
        ac = tuple(points[2][axis] - points[0][axis] for axis in range(3))
        cross = (ab[1] * ac[2] - ab[2] * ac[1],
                 ab[2] * ac[0] - ab[0] * ac[2],
                 ab[0] * ac[1] - ab[1] * ac[0])
        item["area"] += 0.5 * math.sqrt(sum(value * value for value in cross))

    cx, cy = center or CENTER
    small = []
    below_floor = []
    for item in components.values():
        centroid = tuple(value / (3 * item["count"]) for value in item["sum"])
        distance = math.hypot(centroid[0] - cx, centroid[1] - cy)
        if (item["count"] <= FRAGMENT_MAX_TRIANGLES
                and item["area"] <= FRAGMENT_MAX_AREA_M2
                and distance <= FRAGMENT_MAX_DISTANCE_M):
            record = (item["count"], item["area"], distance, centroid)
            (below_floor if item["area"] < FRAGMENT_MIN_AREA_M2 else small).append(record)
    return len(components), small, below_floor


def audit(path, scene_data=None, center=None, height_min=None, height_max=None,
          preferred=None, inner=None, outer=None, normal_limit=None,
          coverage_limit=None, fragment_gate=True):
    cx, cy = center or CENTER
    inner = INNER if inner is None else inner
    outer = OUTER if outer is None else outer
    normal_limit = NORMAL_LIMIT_DEG if normal_limit is None else normal_limit
    coverage_limit = COVERAGE_MIN if coverage_limit is None else coverage_limit
    bounds = (cx - outer, cx + outer, cy - outer, cy + outer)
    if scene_data is None:
        scene_data = RING.read_scene(str(path), 0,
            HEIGHT_MIN if height_min is None else height_min,
            HEIGHT_MAX if height_max is None else height_max, bounds)
    voxel_um, triangles, tiles = scene_data
    component_count, fragments, below_floor = small_components(triangles, (cx, cy))
    print(f"[MouthFragmentComponents] mesh={path.name} components={component_count} "
          f"small_near_mouth={len(fragments)} subthreshold_near_mouth={len(below_floor)} "
          f"limits=(triangles<={FRAGMENT_MAX_TRIANGLES},"
          f"area={FRAGMENT_MIN_AREA_M2:.4f}..{FRAGMENT_MAX_AREA_M2:.3f}m2,"
          f"distance<={FRAGMENT_MAX_DISTANCE_M:.1f}m) "
          f"first={fragments[0] if fragments else None} "
          f"subthreshold={below_floor[0] if below_floor else None} "
          f"gate={'on' if fragment_gate else 'diagnostic'} "
          f"status={'PASS' if not fragments or not fragment_gate else 'FAIL'}")
    preferred = PREFERRED if preferred is None else preferred
    cell = max(2.0, 16.0 * voxel_um / 1_000_000.0)
    buckets = RING.build_buckets(triangles, True, bounds, cell)
    lines = []
    count = round((outer - inner) / STEP)
    for degrees in (135.0, 180.0, 225.0, 270.0, 315.0):
        angle = math.radians(degrees)
        lines.append((f"radial-{degrees:.0f}", [
            (cx + radius * math.cos(angle), cy + radius * math.sin(angle))
            for radius in (inner + i * STEP for i in range(count + 1))]))
    tangent_count = round(8.0 / STEP)
    for degrees in (225.0, 315.0):
        angle = math.radians(degrees)
        normal = (math.cos(angle), math.sin(angle))
        tangent = (-normal[1], normal[0])
        lines.append((f"tangent-{degrees:.0f}", [
            (cx + 6.5 * normal[0] + (-4.0 + i * STEP) * tangent[0],
             cy + 6.5 * normal[1] + (-4.0 + i * STEP) * tangent[1])
            for i in range(tangent_count + 1)]))

    max_angle = 0.0
    max_pair = None
    min_coverage = 1.0
    for name, points in lines:
        samples = [RING.sample_surface(buckets, x, y, True, cell, preferred)
                   for x, y in points]
        coverage = sum(value is not None for value in samples) / len(samples)
        min_coverage = min(min_coverage, coverage)
        line_max = 0.0
        for index in range(1, len(samples)):
            prior, current = samples[index - 1], samples[index]
            if prior is None or current is None:
                continue
            dot = max(-1.0, min(1.0, sum(a * b for a, b in zip(prior[1], current[1]))))
            delta = math.degrees(math.acos(dot))
            if delta > line_max:
                line_max = delta
            if delta > max_angle:
                max_angle = delta
                max_pair = (name, points[index - 1], prior[0], points[index], current[0])
        print(f"[MouthFadeLine] mesh={path.name} line={name} coverage={coverage:.3f} "
              f"max_normal_jump_deg={line_max:.6f}")
    normal_pass = max_angle <= normal_limit and min_coverage >= coverage_limit
    passed = normal_pass and (not fragment_gate or not fragments)
    print(f"[MouthFadeAudit] mesh={path.name} tiles={tiles} lines={len(lines)} "
          f"coverage_min={min_coverage:.3f} coverage_limit={coverage_limit:.2f} "
          f"max_normal_jump_deg={max_angle:.6f} normal_limit_deg={normal_limit:.1f} "
          f"fragments={len(fragments)} "
          f"worst_pair={max_pair} status={'PASS' if passed else 'FAIL'}")
    return passed, max_angle


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--control", type=Path,
                        default=ROOT / "Saved/Linux/l35/owner-fi-ground-mouth.vfmesh")
    parser.add_argument("--prior", type=Path,
                        default=ROOT / "Saved/Linux/l36/owner-fi-ground-mouth.vfmesh")
    parser.add_argument("--current", type=Path,
                        default=ROOT / "Saved/Linux/l37/owner-fi-ground-mouth.vfmesh")
    args = parser.parse_args()
    control = args.control if args.control.is_absolute() else ROOT / args.control
    prior = args.prior if args.prior.is_absolute() else ROOT / args.prior
    current = args.current if args.current.is_absolute() else ROOT / args.current
    control_ok, control_max = audit(control)
    prior_ok, prior_max = audit(prior)
    current_ok, current_max = audit(current)
    controls = (not control_ok and not prior_ok
                and control_max > NORMAL_LIMIT_DEG and prior_max > NORMAL_LIMIT_DEG)
    expected = controls and current_ok and current_max <= NORMAL_LIMIT_DEG
    print(f"[MouthFadeNegativeControl] L35_max_deg={control_max:.6f} "
          f"L35_detected={not control_ok} L36_max_deg={prior_max:.6f} "
          f"L36_detected={not prior_ok} L37_max_deg={current_max:.6f} "
          f"normal_limit_deg={NORMAL_LIMIT_DEG:.1f} current_pass={current_ok} "
          f"status={'PASS' if expected else 'FAIL'}")
    raise SystemExit(0 if expected else 1)


if __name__ == "__main__":
    main()