Fr0zkaandClaude Opus 5.5
9458217d34
Core.Caves.PrimitiveField v1, opt-in per strate (legacy TunnelNetwork assets stay legacy): solid host rock and one smooth union of two elliptical spine tubes and an ellipsoid room, each intersected with its rising floor by cubic smooth-max (W 2.4 m, measured fillet 0.70 m, no trailing floor cap). Bound in the L40 field table and L42 capabilities, BVH feature candidates, isTileShapeFree aware, exact interval boxes (0 misses). Capsule flood fill arrival -> departure passes (headroom >= 0.3 m, grade <= 15 deg); the generalized whole-surface audit passes with its controls failing. Cold whole view 0.727 s, near ring 0.19 s. Graph ordering/serialization changes are scoped to native caves (the flagged L47 demo is unchanged); verify tiers write captures under --out only. fast/medium/full PASS, Clang and GCC. Research R23 (deepened): per-tile LOD0 kernel speed plan. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016dhxM43dnWvUu2B5Ph1iG1
729 lines
35 KiB
Python
729 lines
35 KiB
Python
#!/usr/bin/env python3
|
|
"""Audit every LOD0 mesh triangle across a resolved multi-feature cave graph."""
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import array
|
|
import json
|
|
import math
|
|
from pathlib import Path
|
|
import struct
|
|
import sys
|
|
|
|
|
|
ROOT = Path(__file__).resolve().parents[4]
|
|
CELL_M = 1.0
|
|
CORRIDOR_MARGIN_M = 3.0
|
|
OPEN_GROUND_EXTENSION_M = 5.0
|
|
WELD_M = 0.0001
|
|
|
|
|
|
class Triangle:
|
|
__slots__ = ("verts", "center", "area", "normal", "smooth", "quality",
|
|
"route_distance", "route_lateral", "route_s", "floor_z",
|
|
"route_radius", "route_height", "route_edge", "in_region", "reference")
|
|
|
|
def __init__(self, verts, area, normal, smooth, quality):
|
|
self.verts = verts
|
|
self.center = tuple(sum(v[i] for v in verts) / 3.0 for i in range(3))
|
|
self.area = area
|
|
self.normal = normal
|
|
self.smooth = smooth
|
|
self.quality = quality
|
|
self.route_distance = 0.0
|
|
self.route_lateral = 0.0
|
|
self.route_s = 0.0
|
|
self.floor_z = 0.0
|
|
self.route_radius = 0.0
|
|
self.route_height = 0.0
|
|
self.route_edge = -1
|
|
self.in_region = False
|
|
self.reference = False
|
|
|
|
|
|
def normalize(value):
|
|
length = math.sqrt(sum(item * item for item in value))
|
|
return tuple(item / length for item in value) if length > 1.0e-12 else (0.0, 0.0, 0.0)
|
|
|
|
|
|
def dot(a, b):
|
|
return sum(a[i] * b[i] for i in range(3))
|
|
|
|
|
|
def angle_degrees(a, b):
|
|
return math.degrees(math.acos(max(-1.0, min(1.0, dot(a, b)))))
|
|
|
|
|
|
def quantile(values, fraction):
|
|
if not values:
|
|
return 0.0
|
|
ordered = sorted(values)
|
|
return ordered[min(len(ordered) - 1, max(0, int(math.ceil(fraction * len(ordered))) - 1))]
|
|
|
|
|
|
def read_mesh(path: Path):
|
|
triangles = []
|
|
selected_tiles = 0
|
|
lod_counts = {}
|
|
with path.open("rb") as stream:
|
|
prefix = stream.read(60)
|
|
if len(prefix) != 60 or prefix[:8] != b"VFMSH001":
|
|
raise ValueError(f"{path}: invalid VFMSH001 mesh header")
|
|
version, voxel_um, tile_count = struct.unpack_from("<3I", prefix, 8)
|
|
if version != 1 or not voxel_um or tile_count > 1_000_000:
|
|
raise ValueError(f"{path}: unsupported mesh version or tile count")
|
|
for _ in range(tile_count):
|
|
header = stream.read(36)
|
|
if len(header) != 36:
|
|
raise ValueError(f"{path}: truncated tile header")
|
|
key = struct.unpack_from("<4i", header, 0)
|
|
surface_count, _surface_aux, vertex_count, index_count, mask_count = \
|
|
struct.unpack_from("<5I", header, 16)
|
|
lod_counts[key[3]] = lod_counts.get(key[3], 0) + 1
|
|
position_bytes = stream.read(vertex_count * 12)
|
|
normal_bytes = stream.read(vertex_count * 12)
|
|
index_bytes = stream.read(index_count * 4)
|
|
masks = stream.read(mask_count)
|
|
if (len(position_bytes) != vertex_count * 12
|
|
or len(normal_bytes) != vertex_count * 12
|
|
or len(index_bytes) != index_count * 4
|
|
or len(masks) != mask_count):
|
|
raise ValueError(f"{path}: truncated tile payload")
|
|
if key[3] != 0:
|
|
continue
|
|
selected_tiles += 1
|
|
positions = array.array("f")
|
|
positions.frombytes(position_bytes)
|
|
normals = array.array("f")
|
|
normals.frombytes(normal_bytes)
|
|
indices = array.array("I")
|
|
indices.frombytes(index_bytes)
|
|
if sys.byteorder != "little":
|
|
positions.byteswap()
|
|
normals.byteswap()
|
|
indices.byteswap()
|
|
for offset in range(0, index_count - 2, 3):
|
|
if masks and masks[offset // 3]:
|
|
continue
|
|
ids = indices[offset:offset + 3]
|
|
if max(ids) >= surface_count:
|
|
continue
|
|
verts = tuple(tuple(float(positions[index * 3 + axis]) for axis in range(3))
|
|
for index in ids)
|
|
edge_a = tuple(verts[1][axis] - verts[0][axis] for axis in range(3))
|
|
edge_b = tuple(verts[2][axis] - verts[0][axis] for axis in range(3))
|
|
cross = (edge_a[1] * edge_b[2] - edge_a[2] * edge_b[1],
|
|
edge_a[2] * edge_b[0] - edge_a[0] * edge_b[2],
|
|
edge_a[0] * edge_b[1] - edge_a[1] * edge_b[0])
|
|
cross_length = math.sqrt(dot(cross, cross))
|
|
if cross_length <= 1.0e-12:
|
|
continue
|
|
face_normal = tuple(value / cross_length for value in cross)
|
|
area = cross_length * 0.5
|
|
smooth_normal = normalize(tuple(
|
|
sum(float(normals[index * 3 + axis]) for index in ids) for axis in range(3)))
|
|
if dot(face_normal, smooth_normal) < 0.0:
|
|
face_normal = tuple(-value for value in face_normal)
|
|
edge_squares = 0.0
|
|
for a, b in ((verts[0], verts[1]), (verts[1], verts[2]),
|
|
(verts[2], verts[0])):
|
|
edge_squares += sum((a[i] - b[i]) ** 2 for i in range(3))
|
|
quality = 4.0 * math.sqrt(3.0) * area / max(edge_squares, 1.0e-15)
|
|
triangles.append(Triangle(verts, area, face_normal, smooth_normal, quality))
|
|
if stream.read(1):
|
|
raise ValueError(f"{path}: trailing bytes after mesh tiles")
|
|
if not triangles:
|
|
raise ValueError(f"{path}: no LOD0 surface triangles")
|
|
return triangles, voxel_um, selected_tiles, lod_counts
|
|
|
|
|
|
def bezier(points, t):
|
|
one = 1.0 - t
|
|
weights = (one ** 3, 3.0 * one * one * t, 3.0 * one * t * t, t ** 3)
|
|
return tuple(sum(weights[i] * points[i][axis] for i in range(4)) for axis in range(3))
|
|
|
|
|
|
def make_route(layout, phase_m=0.0):
|
|
edges = layout.get("cavern_edges", [])
|
|
rooms = layout.get("cavern_rooms", [])
|
|
if not edges or not rooms:
|
|
raise ValueError("cavern audit needs at least one resolved room and edge")
|
|
segments = []
|
|
routes = []
|
|
for edge_index, edge in enumerate(sorted(edges, key=lambda item: item["stable_id"])):
|
|
points = [[value / 1_000_000.0 for value in row]
|
|
for row in edge["control_points_um"]]
|
|
bore_radius = edge["bore_radius_um"] / 1_000_000.0
|
|
bore_height = edge["bore_height_um"] / 1_000_000.0
|
|
samples = [bezier(points, index / 160.0) for index in range(161)]
|
|
horizontal_s = [0.0]
|
|
for a, b in zip(samples, samples[1:]):
|
|
horizontal_s.append(horizontal_s[-1] + math.hypot(b[0] - a[0], b[1] - a[1]))
|
|
first_dx = samples[1][0] - samples[0][0]
|
|
first_dy = samples[1][1] - samples[0][1]
|
|
first_len = max(1.0e-9, math.hypot(first_dx, first_dy))
|
|
extension = (samples[0][0] - first_dx / first_len * OPEN_GROUND_EXTENSION_M,
|
|
samples[0][1] - first_dy / first_len * OPEN_GROUND_EXTENSION_M,
|
|
samples[0][2] - (samples[1][2] - samples[0][2]) / first_len
|
|
* OPEN_GROUND_EXTENSION_M)
|
|
route_points = [extension, *samples]
|
|
route_s = [-OPEN_GROUND_EXTENSION_M, *horizontal_s]
|
|
for i, (a, b) in enumerate(zip(route_points, route_points[1:])):
|
|
segments.append((a, b, route_s[i], route_s[i + 1],
|
|
bore_radius, bore_height, edge_index))
|
|
routes.append({"edge": edge, "points": samples, "radius": bore_radius,
|
|
"height": bore_height, "length": horizontal_s[-1], "index": edge_index})
|
|
resolved_rooms = [{"room": room,
|
|
"center": tuple(value / 1_000_000.0 for value in room["center_um"]),
|
|
"half": tuple(value / 1_000_000.0 for value in room["half_extents_um"]),
|
|
"floor": room["floor_um"] / 1_000_000.0}
|
|
for room in sorted(rooms, key=lambda item: item["stable_id"])]
|
|
return {
|
|
"segments": segments,
|
|
"edges": routes,
|
|
"rooms": resolved_rooms,
|
|
"phase_m": phase_m,
|
|
}
|
|
|
|
|
|
def route_nearest(point, route):
|
|
best = None
|
|
for a, b, s0, s1, radius, height, edge_index in route["segments"]:
|
|
dx, dy = b[0] - a[0], b[1] - a[1]
|
|
denominator = dx * dx + dy * dy
|
|
t = 0.0 if denominator <= 1.0e-12 else max(0.0, min(1.0,
|
|
((point[0] - a[0]) * dx + (point[1] - a[1]) * dy) / denominator))
|
|
x, y = a[0] + dx * t, a[1] + dy * t
|
|
distance = math.hypot(point[0] - x, point[1] - y)
|
|
if best is None or distance < best[0]:
|
|
segment_length = math.sqrt(max(denominator, 1.0e-12))
|
|
lateral = (dx * (point[1] - y) - dy * (point[0] - x)) / segment_length
|
|
best = (distance, s0 + (s1 - s0) * t,
|
|
a[2] + (b[2] - a[2]) * t + route["phase_m"],
|
|
lateral, radius, height, edge_index,
|
|
route["edges"][edge_index]["length"])
|
|
return best
|
|
|
|
|
|
def mark_region(triangles, route):
|
|
for tri in triangles:
|
|
x, y, z = tri.center
|
|
nearest = route_nearest(tri.center, route)
|
|
distance, distance_s, floor_z, lateral, radius, height, edge_index, edge_length = nearest
|
|
tri.route_distance = distance
|
|
tri.route_lateral = lateral
|
|
tri.route_s = distance_s
|
|
tri.floor_z = floor_z
|
|
tri.route_radius = radius
|
|
tri.route_height = height
|
|
tri.route_edge = edge_index
|
|
corridor = (-OPEN_GROUND_EXTENSION_M <= distance_s <= edge_length
|
|
and distance <= radius + CORRIDOR_MARGIN_M
|
|
and floor_z - CORRIDOR_MARGIN_M <= z
|
|
<= floor_z + height + CORRIDOR_MARGIN_M)
|
|
room_box = any(abs(x - room["center"][0]) <= room["half"][0] + CORRIDOR_MARGIN_M
|
|
and abs(y - room["center"][1]) <= room["half"][1] + CORRIDOR_MARGIN_M
|
|
and room["floor"] - CORRIDOR_MARGIN_M <= z
|
|
<= room["floor"] + 2.0 * room["half"][2] + CORRIDOR_MARGIN_M
|
|
for room in route["rooms"])
|
|
tri.in_region = corridor or room_box
|
|
inside_room = any(abs(x - room["center"][0]) <= room["half"][0] + 0.45
|
|
and abs(y - room["center"][1]) <= room["half"][1] + 0.45
|
|
for room in route["rooms"])
|
|
outside_bore = distance > radius + 0.45 and not inside_room
|
|
# Keep the reference on native, ground-facing terrain. The outer mesh
|
|
# also contains steep cavern walls and room edges, which are features.
|
|
tri.reference = tri.in_region and outside_bore and tri.smooth[2] >= 0.45
|
|
|
|
|
|
def vertex_key(point):
|
|
return tuple(int(round(value / WELD_M)) for value in point)
|
|
|
|
|
|
def make_edges(triangles, indices):
|
|
edge_map = {}
|
|
for index in indices:
|
|
tri = triangles[index]
|
|
points = tri.verts
|
|
for a, b in ((0, 1), (1, 2), (2, 0)):
|
|
ka, kb = vertex_key(points[a]), vertex_key(points[b])
|
|
if ka == kb:
|
|
continue
|
|
key = (ka, kb) if ka < kb else (kb, ka)
|
|
delta = tuple(points[a][axis] - points[b][axis] for axis in range(3))
|
|
length = math.sqrt(dot(delta, delta))
|
|
midpoint = tuple((points[a][axis] + points[b][axis]) * 0.5 for axis in range(3))
|
|
edge_map.setdefault(key, []).append((index, length, midpoint))
|
|
pairs = []
|
|
for incident in edge_map.values():
|
|
unique = {}
|
|
for index, length, midpoint in incident:
|
|
unique[index] = (length, midpoint)
|
|
incident = list(unique.items())
|
|
if len(incident) < 2:
|
|
continue
|
|
if len(incident) > 4:
|
|
incident = incident[:4]
|
|
for left in range(len(incident)):
|
|
for right in range(left + 1, len(incident)):
|
|
i, (length, midpoint) = incident[left]
|
|
j, _ = incident[right]
|
|
angle = max(angle_degrees(triangles[i].normal, triangles[j].normal),
|
|
0.75 * angle_degrees(triangles[i].smooth, triangles[j].smooth))
|
|
pairs.append((i, j, angle, length, midpoint))
|
|
return pairs
|
|
|
|
|
|
def cell_of(point):
|
|
return tuple(math.floor(value / CELL_M) for value in point)
|
|
|
|
|
|
def score_surface(triangles, considered, pairs, edge_cutoff, rough_cutoff,
|
|
ref_cell_density):
|
|
considered_set = set(considered)
|
|
variation = {index: 0.0 for index in considered}
|
|
max_jump = {index: 0.0 for index in considered}
|
|
severe_count = {index: 0.0 for index in considered}
|
|
for i, j, angle, length, _midpoint in pairs:
|
|
if i not in considered_set or j not in considered_set:
|
|
continue
|
|
weighted = angle * length * 0.5
|
|
variation[i] += weighted
|
|
variation[j] += weighted
|
|
max_jump[i] = max(max_jump[i], angle)
|
|
max_jump[j] = max(max_jump[j], angle)
|
|
if angle >= edge_cutoff:
|
|
severe_count[i] += 0.5
|
|
severe_count[j] += 0.5
|
|
cells = {}
|
|
for index in considered:
|
|
tri = triangles[index]
|
|
key = cell_of(tri.center)
|
|
cell = cells.setdefault(key, {"area": 0.0, "variation": 0.0, "severe": 0.0,
|
|
"sliver_area": 0.0, "slivers": 0, "triangles": 0, "max_jump": 0.0,
|
|
"reference_area": 0.0, "reference_variation": 0.0,
|
|
"reference_severe": 0.0, "reference_sliver_area": 0.0,
|
|
"reference_slivers": 0})
|
|
cell["area"] += tri.area
|
|
cell["variation"] += variation[index]
|
|
cell["severe"] += severe_count[index]
|
|
cell["triangles"] += 1
|
|
cell["max_jump"] = max(cell["max_jump"], max_jump[index])
|
|
sliver = tri.quality < 0.07 and tri.area < 0.12
|
|
if sliver:
|
|
cell["sliver_area"] += tri.area
|
|
cell["slivers"] += 1
|
|
if tri.reference:
|
|
cell["reference_area"] += tri.area
|
|
cell["reference_variation"] += variation[index]
|
|
cell["reference_severe"] += severe_count[index]
|
|
if sliver:
|
|
cell["reference_sliver_area"] += tri.area
|
|
cell["reference_slivers"] += 1
|
|
ref_rough = []
|
|
ref_density = []
|
|
ref_sliver_fraction = []
|
|
for cell in cells.values():
|
|
if cell["reference_area"] >= max(0.15, cell["area"] * 0.65):
|
|
ref_rough.append(cell["reference_variation"] / max(cell["reference_area"], 1.0e-9))
|
|
ref_density.append(cell["reference_severe"] / max(cell["reference_area"], 1.0e-9))
|
|
ref_sliver_fraction.append(cell["reference_sliver_area"]
|
|
/ max(cell["reference_area"], 1.0e-9))
|
|
reference_median = quantile(ref_rough, 0.5)
|
|
reference_mad = quantile([abs(value - reference_median) for value in ref_rough], 0.5)
|
|
cell_rough_limit = max(10.0, quantile(ref_rough, 0.995) * 1.5,
|
|
reference_median + 8.0 * reference_mad)
|
|
edge_density_limit = max(1.5, quantile(ref_density, 0.995) * 1.7,
|
|
ref_cell_density * 1.7)
|
|
sliver_limit = max(0.22, quantile(ref_sliver_fraction, 0.995) + 0.18)
|
|
hot = set()
|
|
for key, cell in cells.items():
|
|
area = max(cell["area"], 1.0e-9)
|
|
roughness = cell["variation"] / area
|
|
edge_density = cell["severe"] / area
|
|
sliver_fraction = cell["sliver_area"] / area
|
|
rough_cluster = roughness > cell_rough_limit and cell["severe"] >= 2.0
|
|
edge_cluster = cell["severe"] >= 4.0 and edge_density > edge_density_limit
|
|
sliver_cluster = cell["slivers"] >= 5 and sliver_fraction > sliver_limit
|
|
if rough_cluster or edge_cluster or sliver_cluster:
|
|
hot.add(key)
|
|
cell["roughness_deg_per_m"] = roughness
|
|
cell["severe_edge_density_m2"] = edge_density
|
|
cell["sliver_area_fraction"] = sliver_fraction
|
|
cell["reasons"] = [name for name, active in (("normal-variation", rough_cluster),
|
|
("edge-density", edge_cluster), ("slivers", sliver_cluster)) if active]
|
|
clusters = []
|
|
flagged_cells = set()
|
|
remaining = set(hot)
|
|
while remaining:
|
|
stack = [remaining.pop()]
|
|
members = []
|
|
while stack:
|
|
current = stack.pop()
|
|
members.append(current)
|
|
x, y, z = current
|
|
neighbors = [(x + dx, y + dy, z + dz)
|
|
for dx in (-1, 0, 1) for dy in (-1, 0, 1) for dz in (-1, 0, 1)
|
|
if dx or dy or dz]
|
|
for neighbor in neighbors:
|
|
if neighbor in remaining:
|
|
remaining.remove(neighbor)
|
|
stack.append(neighbor)
|
|
cells_in_cluster = [cells[key] for key in members]
|
|
cluster_area = sum(cell["area"] for cell in cells_in_cluster)
|
|
cluster_severe = sum(cell["severe"] for cell in cells_in_cluster)
|
|
cluster_slivers = sum(cell["slivers"] for cell in cells_in_cluster)
|
|
accepted = len(members) >= 2 and cluster_area >= 0.25 and \
|
|
(cluster_severe >= 8.0 or cluster_slivers >= 12)
|
|
if accepted:
|
|
flagged_cells.update(members)
|
|
clusters.append({"cells": len(members), "area_m2": cluster_area,
|
|
"severe_edges": cluster_severe, "sliver_triangles": cluster_slivers,
|
|
"bounds_m": [[min(key[i] for key in members) * CELL_M,
|
|
(max(key[i] for key in members) + 1) * CELL_M]
|
|
for i in range(3)],
|
|
"reasons": sorted({reason for cell in cells_in_cluster
|
|
for reason in cell["reasons"]})})
|
|
clusters.sort(key=lambda cluster: (cluster["severe_edges"], cluster["area_m2"]), reverse=True)
|
|
flagged_triangles = [index for index in considered
|
|
if cell_of(triangles[index].center) in flagged_cells]
|
|
return {"cells": cells, "hot_cells": hot, "flagged_cells": flagged_cells,
|
|
"clusters": clusters, "hot_triangles": flagged_triangles,
|
|
"thresholds": {"edge_angle_deg": edge_cutoff,
|
|
"roughness_deg_per_m": cell_rough_limit,
|
|
"severe_edge_density_m2": edge_density_limit,
|
|
"sliver_area_fraction": sliver_limit,
|
|
"reference_cell_roughness_p995_deg_per_m": quantile(ref_rough, 0.995),
|
|
"reference_cell_roughness_median_deg_per_m": reference_median}}
|
|
|
|
|
|
def floor_profile_audit(triangles, considered, region_indices, route):
|
|
"""Measure the whole floor patch as area-weighted bins, not selected rails."""
|
|
if route is None:
|
|
return {"present": False, "detected": False, "triangles": 0, "bins": 0,
|
|
"lateral_bands": 0, "ripple_amplitude_m": 0.0,
|
|
"max_adjacent_step_m": 0.0}
|
|
along_step = 0.25
|
|
lateral_step = 0.5
|
|
bins = {}
|
|
floor_triangles = 0
|
|
wall_triangles = 0
|
|
wall_bins = set()
|
|
for index in region_indices:
|
|
tri = triangles[index]
|
|
height = tri.center[2] - tri.floor_z
|
|
if (0.0 <= tri.route_s <= min(10.0,
|
|
route["edges"][tri.route_edge]["length"])
|
|
and tri.route_radius - 0.9 <= tri.route_distance
|
|
<= tri.route_radius + 1.0
|
|
and 0.45 <= height <= tri.route_height + 0.5
|
|
and abs(tri.smooth[2]) <= 0.45):
|
|
wall_triangles += 1
|
|
wall_bins.add(int(tri.route_s))
|
|
for index in considered:
|
|
tri = triangles[index]
|
|
if (tri.route_s < 0.0 or tri.route_s > min(10.0,
|
|
route["edges"][tri.route_edge]["length"])
|
|
or tri.route_distance > min(1.9, tri.route_radius - 0.1)
|
|
or abs(tri.center[2] - tri.floor_z) > 0.60
|
|
or tri.smooth[2] < 0.55):
|
|
continue
|
|
floor_triangles += 1
|
|
along_bin = int(math.floor(tri.route_s / along_step))
|
|
lateral_bin = int(math.floor((tri.route_lateral + 2.0) / lateral_step))
|
|
key = (lateral_bin, along_bin)
|
|
row = bins.setdefault(key, [0.0, 0.0, 0])
|
|
residual = tri.center[2] - tri.floor_z
|
|
row[0] += residual * tri.area
|
|
row[1] += tri.area
|
|
row[2] += 1
|
|
by_lateral = {}
|
|
for (lateral_bin, along_bin), (weighted_residual, area, count) in bins.items():
|
|
if area < 0.002:
|
|
continue
|
|
by_lateral.setdefault(lateral_bin, {})[along_bin] = weighted_residual / area
|
|
amplitudes = []
|
|
steps = []
|
|
valid_bins = 0
|
|
for values in by_lateral.values():
|
|
ordered = sorted(values.items())
|
|
run = []
|
|
runs = []
|
|
for item in ordered:
|
|
if run and item[0] != run[-1][0] + 1:
|
|
runs.append(run)
|
|
run = []
|
|
run.append(item)
|
|
if run:
|
|
runs.append(run)
|
|
for segment in runs:
|
|
if len(segment) < 8:
|
|
continue
|
|
valid_bins += len(segment)
|
|
raw = [value for _, value in segment]
|
|
highpass = []
|
|
for position, value in enumerate(raw):
|
|
nearby = raw[max(0, position - 4):min(len(raw), position + 5)]
|
|
highpass.append(value - sum(nearby) / len(nearby))
|
|
amplitudes.append(max(highpass) - min(highpass))
|
|
steps.extend(abs(raw[i] - raw[i - 1]) for i in range(1, len(raw)))
|
|
max_ripple = max(amplitudes, default=0.0) * 0.5
|
|
max_step = max(steps, default=0.0)
|
|
present = (floor_triangles >= 24 and valid_bins >= 8
|
|
and wall_triangles >= 24 and len(wall_bins) >= 3)
|
|
detected = present and (max_ripple > 0.015 or max_step > 0.050)
|
|
return {"present": present, "detected": detected, "triangles": floor_triangles,
|
|
"bins": valid_bins, "lateral_bands": len(by_lateral),
|
|
"wall_triangles": wall_triangles, "wall_bins": len(wall_bins),
|
|
"ripple_amplitude_m": max_ripple, "max_adjacent_step_m": max_step,
|
|
"thresholds": {"ripple_amplitude_m": 0.015, "adjacent_step_m": 0.050}}
|
|
|
|
|
|
def wall_band_measure(pairs, route, edge_cutoff):
|
|
"""Locate high dihedrals around the bore-foot and the L47 0.25 m band edge."""
|
|
if route is None:
|
|
return {}
|
|
bins = {}
|
|
for _i, _j, angle, _length, midpoint in pairs:
|
|
distance, route_s, floor_z, _lateral, radius, _height, _edge, _edge_length = \
|
|
route_nearest(midpoint, route)
|
|
height = midpoint[2] - floor_z
|
|
if (0.0 <= route_s <= 10.0 and abs(distance - radius) <= 0.75
|
|
and -0.25 <= height <= 0.875):
|
|
index = int(math.floor((height + 0.25) / 0.125))
|
|
entry = bins.setdefault(index, [])
|
|
entry.append(angle)
|
|
result = {}
|
|
for index, angles in sorted(bins.items()):
|
|
center_height = (index + 0.5) * 0.125 - 0.25
|
|
result[f"{center_height:.3f}"] = {
|
|
"edges": len(angles),
|
|
"over_reference_cutoff": sum(value >= edge_cutoff for value in angles),
|
|
"p95_dihedral_deg": quantile(angles, 0.95),
|
|
"max_dihedral_deg": max(angles, default=0.0),
|
|
}
|
|
return result
|
|
|
|
|
|
def radial_wall_measure(pairs, route, edge_cutoff):
|
|
"""Measure dihedrals across the hard bore-radius gate at the entrance foot."""
|
|
if route is None:
|
|
return {}
|
|
bins = {}
|
|
for _i, _j, angle, _length, midpoint in pairs:
|
|
distance, route_s, floor_z, _lateral, radius, _height, _edge, _edge_length = \
|
|
route_nearest(midpoint, route)
|
|
height = midpoint[2] - floor_z
|
|
if (-5.0 <= route_s <= 3.0 and -1.0 <= height <= 0.75
|
|
and abs(distance - radius) <= 1.0):
|
|
index = int(math.floor((distance - (radius - 1.0)) / 0.125))
|
|
bins.setdefault(index, []).append(angle)
|
|
return {f"{radius - 1.0 + (index + 0.5) * 0.125:.3f}": {
|
|
"edges": len(angles),
|
|
"over_reference_cutoff": sum(value >= edge_cutoff for value in angles),
|
|
"p95_dihedral_deg": quantile(angles, 0.95),
|
|
"max_dihedral_deg": max(angles, default=0.0),
|
|
} for index, angles in sorted(bins.items())}
|
|
|
|
|
|
def svg_for(path, triangles, route, hot_triangles, label):
|
|
selected = [triangles[index] for index in hot_triangles]
|
|
all_region = [tri for tri in triangles if tri.in_region]
|
|
points = [point for tri in all_region for point in tri.verts]
|
|
if not points:
|
|
return
|
|
xmin, xmax = min(p[0] for p in points), max(p[0] for p in points)
|
|
ymin, ymax = min(p[1] for p in points), max(p[1] for p in points)
|
|
width, height, margin = 1200.0, 900.0, 35.0
|
|
sx = (width - 2 * margin) / max(1.0e-6, xmax - xmin)
|
|
sy = (height - 2 * margin) / max(1.0e-6, ymax - ymin)
|
|
scale = min(sx, sy)
|
|
draw_w, draw_h = (xmax - xmin) * scale, (ymax - ymin) * scale
|
|
ox, oy = (width - draw_w) * 0.5, (height - draw_h) * 0.5
|
|
def project(point):
|
|
return ox + (point[0] - xmin) * scale, height - oy - (point[1] - ymin) * scale
|
|
lines = [f'<svg xmlns="http://www.w3.org/2000/svg" width="{int(width)}" height="{int(height)}" viewBox="0 0 {int(width)} {int(height)}">',
|
|
'<rect width="100%" height="100%" fill="#fbfbfa"/>',
|
|
f'<text x="28" y="28" font-family="sans-serif" font-size="19" fill="#222">{label}: whole-surface defect clusters (XY plan)</text>',
|
|
'<g fill="#c8ced2" fill-opacity="0.24" stroke="#8b959b" stroke-opacity="0.12" stroke-width="0.3">']
|
|
for tri in all_region:
|
|
coords = [project(p) for p in tri.verts]
|
|
lines.append('<polygon points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in coords)
|
|
+ '"/>')
|
|
lines.append('</g><g fill="#d7191c" fill-opacity="0.78" stroke="#7f0000" stroke-width="0.4">')
|
|
for tri in selected:
|
|
coords = [project(p) for p in tri.verts]
|
|
lines.append('<polygon points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in coords)
|
|
+ '"/>')
|
|
lines.append('</g>')
|
|
for feature in route["edges"]:
|
|
route_coords = [project(p) for p in feature["points"]]
|
|
lines.append('<polyline points="' + " ".join(f"{x:.1f},{y:.1f}" for x, y in route_coords)
|
|
+ '" fill="none" stroke="#1769aa" stroke-width="2"/>')
|
|
lines.append(f'<text x="28" y="{int(height - 24)}" font-family="sans-serif" font-size="15" fill="#333">Blue: intended route · Red: flagged triangle cells · Gray: remaining ROI triangles · flagged triangles: {len(selected)}</text>')
|
|
lines.append('</svg>')
|
|
path.parent.mkdir(parents=True, exist_ok=True)
|
|
path.write_text("\n".join(lines) + "\n", encoding="utf-8")
|
|
|
|
|
|
def audit_case(label, mesh_path, layout_path, expect, focus, svg_dir, phase):
|
|
triangles, voxel_um, tile_count, lod_counts = read_mesh(mesh_path)
|
|
route = None
|
|
if layout_path != Path("-"):
|
|
layout = json.loads(layout_path.read_text(encoding="utf-8"))
|
|
route = make_route(layout, phase * voxel_um / 1_000_000.0)
|
|
mark_region(triangles, route)
|
|
region_indices = [i for i, tri in enumerate(triangles) if tri.in_region]
|
|
else:
|
|
region_indices = list(range(len(triangles)))
|
|
for tri in triangles:
|
|
tri.in_region = True
|
|
tri.reference = True
|
|
if len(region_indices) < 24:
|
|
raise ValueError(f"{label}: only {len(region_indices)} triangles intersect the ROI")
|
|
reference_indices = [i for i in region_indices if triangles[i].reference]
|
|
if len(reference_indices) < 24:
|
|
reference_indices = region_indices
|
|
for index in region_indices:
|
|
triangles[index].reference = True
|
|
edge_pairs = make_edges(triangles, region_indices)
|
|
reference_edges = [edge[2] for edge in edge_pairs
|
|
if triangles[edge[0]].reference and triangles[edge[1]].reference]
|
|
edge_cutoff = max(35.0, quantile(reference_edges, 0.995))
|
|
reference_set = set(reference_indices)
|
|
ref_roughness = {index: 0.0 for index in reference_indices}
|
|
for i, j, angle, length, _midpoint in edge_pairs:
|
|
if i in reference_set and j in reference_set:
|
|
ref_roughness[i] += angle * length * 0.5
|
|
ref_roughness[j] += angle * length * 0.5
|
|
ref_values = [ref_roughness[index] / max(triangles[index].area, 1.0e-9)
|
|
for index in reference_indices]
|
|
rough_cutoff = 10.0
|
|
ref_density = 0.0
|
|
ref_area = sum(triangles[index].area for index in reference_indices)
|
|
ref_severe = 0.0
|
|
for i, j, angle, _length, _midpoint in edge_pairs:
|
|
if i in reference_set and j in reference_set and angle >= edge_cutoff:
|
|
ref_severe += 1.0
|
|
ref_density = ref_severe / max(ref_area, 1.0e-9)
|
|
if focus == "floor":
|
|
considered = [i for i in region_indices if route is not None
|
|
and triangles[i].route_s >= 1.5
|
|
and triangles[i].route_s <= 10.0
|
|
and triangles[i].route_distance <= 1.25
|
|
and abs(triangles[i].center[2] - triangles[i].floor_z) <= 0.24
|
|
and triangles[i].smooth[2] >= 0.55]
|
|
else:
|
|
considered = region_indices
|
|
if len(considered) < 24:
|
|
raise ValueError(f"{label}: focus={focus} has only {len(considered)} triangles")
|
|
scores = score_surface(triangles, considered, edge_pairs, edge_cutoff,
|
|
rough_cutoff, ref_density)
|
|
floor_profile = floor_profile_audit(triangles, considered, region_indices, route)
|
|
wall_band = wall_band_measure(edge_pairs, route, edge_cutoff)
|
|
radial_wall = radial_wall_measure(edge_pairs, route, edge_cutoff)
|
|
found = bool(scores["clusters"]) or floor_profile["detected"]
|
|
passed = expect == "measure" or (expect == "fail" and found) or (expect == "pass" and not found)
|
|
if svg_dir and route is not None:
|
|
svg_for(svg_dir / f"flagged-regions-{label}.svg", triangles, route,
|
|
scores["hot_triangles"], label)
|
|
cluster_summary = scores["clusters"][:8]
|
|
result = {
|
|
"label": label,
|
|
"mesh": str(mesh_path),
|
|
"layout": None if layout_path == Path("-") else str(layout_path),
|
|
"focus": focus,
|
|
"grid_phase_h": phase,
|
|
"grid_phase_m": phase * voxel_um / 1_000_000.0,
|
|
"expected": expect,
|
|
"observed": "fail" if found else "pass",
|
|
"status": "PASS" if passed else "FAIL",
|
|
"voxel_um": voxel_um,
|
|
"lod0_tiles": tile_count,
|
|
"lod_tile_counts": lod_counts,
|
|
"triangles_loaded_lod0": len(triangles),
|
|
"triangles_in_roi": len(region_indices),
|
|
"triangles_audited": len(considered),
|
|
"adjacent_pairs": len(edge_pairs),
|
|
"reference_triangles": len(reference_indices),
|
|
"reference_edge_p995_deg": quantile(reference_edges, 0.995),
|
|
"reference_p995_normal_variation_deg_per_m": quantile(ref_values, 0.995),
|
|
"thresholds": scores["thresholds"],
|
|
"flagged_cells": len(scores["hot_cells"]),
|
|
"flagged_triangles": len(scores["hot_triangles"]),
|
|
"floor_profile": floor_profile,
|
|
"wall_band_height_histogram": wall_band,
|
|
"radial_wall_histogram": radial_wall,
|
|
"clusters": cluster_summary,
|
|
}
|
|
print(f"[WholeSurfaceAudit] label={label} focus={focus} phase_h={phase:+.1f} "
|
|
f"edges={len(route['edges']) if route else 0} rooms={len(route['rooms']) if route else 0} "
|
|
f"lod0_tiles={tile_count} "
|
|
f"triangles_roi={len(region_indices)} audited={len(considered)} "
|
|
f"adjacent_pairs={len(edge_pairs)} reference={len(reference_indices)} "
|
|
f"edge_cutoff_deg={edge_cutoff:.2f} "
|
|
f"roughness_cutoff_deg_m={scores['thresholds']['roughness_deg_per_m']:.2f} "
|
|
f"hot_cells={len(scores['hot_cells'])} flagged_triangles={len(scores['hot_triangles'])} "
|
|
f"clusters={len(scores['clusters'])} observed={result['observed'].upper()} "
|
|
f"expected={expect.upper()} status={result['status']}")
|
|
print(f"[WholeSurfaceFloorProfile] label={label} triangles={floor_profile['triangles']} "
|
|
f"bins={floor_profile['bins']} lateral_bands={floor_profile['lateral_bands']} "
|
|
f"wall_triangles={floor_profile.get('wall_triangles', 0)} "
|
|
f"wall_bins={floor_profile.get('wall_bins', 0)} "
|
|
f"ripple_amplitude_m={floor_profile['ripple_amplitude_m']:.4f} "
|
|
f"max_adjacent_step_m={floor_profile['max_adjacent_step_m']:.4f} "
|
|
f"observed={'FAIL' if floor_profile['detected'] else ('PASS' if floor_profile['present'] else 'NOT_PRESENT')}")
|
|
if wall_band:
|
|
print(f"[WholeSurfaceWallBandMeasure] label={label} phase_h={phase:+.1f} "
|
|
f"voxel_m={voxel_um / 1_000_000.0:.3f} recarve_edge_height_m=0.250 "
|
|
f"height_histogram={json.dumps(wall_band, sort_keys=True, separators=(',', ':'))}")
|
|
print(f"[WholeSurfaceRadialWallMeasure] label={label} phase_h={phase:+.1f} "
|
|
f"route_s_m=-5..3 height_above_floor_m=-1..0.75 "
|
|
f"radial_histogram={json.dumps(radial_wall, sort_keys=True, separators=(',', ':'))}")
|
|
for index, cluster in enumerate(cluster_summary):
|
|
print(f"[WholeSurfaceCluster] label={label} index={index} cells={cluster['cells']} "
|
|
f"area_m2={cluster['area_m2']:.3f} severe_edges={cluster['severe_edges']:.1f} "
|
|
f"slivers={cluster['sliver_triangles']} bounds_m={cluster['bounds_m']} "
|
|
f"reasons={','.join(cluster['reasons'])}")
|
|
hot_cells = scores["cells"]
|
|
for key in sorted(scores["hot_cells"], key=lambda item: (
|
|
hot_cells[item]["severe"], hot_cells[item]["variation"]), reverse=True)[:12]:
|
|
cell = hot_cells[key]
|
|
local = route_nearest(tuple((axis + 0.5) * CELL_M for axis in key), route) if route else None
|
|
print(f"[WholeSurfaceHotCell] label={label} cell_m=({key[0]},{key[1]},{key[2]}) "
|
|
f"route_s_m={local[1] if local else 0.0:.2f} radial_m={local[0] if local else 0.0:.2f} "
|
|
f"height_above_floor_m={(key[2] + 0.5 - local[2]) if local else 0.0:.2f} "
|
|
f"area_m2={cell['area']:.3f} roughness_deg_m={cell['roughness_deg_per_m']:.2f} "
|
|
f"severe_edges={cell['severe']:.1f} edge_density_m2={cell['severe_edge_density_m2']:.2f} "
|
|
f"slivers={cell['slivers']} max_jump_deg={cell['max_jump']:.2f} "
|
|
f"reasons={','.join(cell['reasons'])}")
|
|
return result
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument("--case", nargs=5, action="append", required=True,
|
|
metavar=("LABEL", "MESH", "LAYOUT_OR_DASH", "EXPECT", "FOCUS"))
|
|
parser.add_argument("--json", type=Path, required=True)
|
|
parser.add_argument("--svg-dir", type=Path)
|
|
parser.add_argument("--phase", type=float, choices=(-0.5, 0.0, 0.5), default=0.0,
|
|
help="vertical analysis phase in multiples of the base voxel spacing")
|
|
args = parser.parse_args()
|
|
results = []
|
|
for label, mesh, layout, expect, focus in args.case:
|
|
if expect not in ("pass", "fail", "measure") or focus not in ("whole", "floor"):
|
|
parser.error("EXPECT must be pass/fail/measure and FOCUS must be whole/floor")
|
|
results.append(audit_case(label, Path(mesh), Path(layout), expect, focus,
|
|
args.svg_dir, args.phase))
|
|
args.json.parent.mkdir(parents=True, exist_ok=True)
|
|
args.json.write_text(json.dumps({"audit": "whole-surface-l54-v1", "cases": results},
|
|
indent=2) + "\n", encoding="utf-8")
|
|
failures = [result for result in results if result["status"] != "PASS"]
|
|
print(f"[WholeSurfaceAuditSummary] cases={len(results)} failed={len(failures)} "
|
|
f"json={args.json} status={'FAIL' if failures else 'PASS'}")
|
|
return 1 if failures else 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|