The shared boundary-cap path now emits a cap wherever the field is solid just below a strate's configured top (it previously required a solid-to-air crossing there) and checks a 5x5 perimeter grid before suppressing a fallback cap. Strict audit: 0 misses / 0 edge gaps, LOD0-8. farHorizonVersion 3 -> 4. SurfaceWorld LOD0-2, full detail and exact queries unchanged; 28/28 matrix plus 14/14 crossing; route parity and proofs pass. Ledges: the lower stratum's terrain at the FI/SW interface, seen through the FI void, sampled at 0.5/1/2 m at LOD1-3 (up to 1.5 m / 7.7 m ring steps). Fix deferred to the units/layers round. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01CSFYgsdW8ejZSyiMf4qRSX
173 lines
6.9 KiB
Python
173 lines
6.9 KiB
Python
#!/usr/bin/env python3
|
|
"""Strict L19 audit of FloatingIslands roof coverage at every LOD."""
|
|
import array
|
|
import collections
|
|
import math
|
|
import struct
|
|
import sys
|
|
|
|
|
|
def read_exact(stream, size):
|
|
data = stream.read(size)
|
|
if len(data) != size:
|
|
raise EOFError(f"wanted {size} bytes, got {len(data)}")
|
|
return data
|
|
|
|
|
|
def height_at(triangles, x, y, preferred_z):
|
|
best = None
|
|
for (x0, y0, z0), (x1, y1, z1), (x2, y2, z2) in triangles:
|
|
den = (y1-y2)*(x0-x2)+(x2-x1)*(y0-y2)
|
|
if abs(den) < 1.0e-10:
|
|
continue
|
|
a = ((y1-y2)*(x-x2)+(x2-x1)*(y-y2))/den
|
|
b = ((y2-y0)*(x-x2)+(x0-x2)*(y-y2))/den
|
|
c = 1.0-a-b
|
|
if min(a,b,c) >= -1.0e-5:
|
|
z = a*z0+b*z1+c*z2
|
|
if best is None or abs(z-preferred_z) < abs(best-preferred_z):
|
|
best = z
|
|
return best
|
|
|
|
|
|
def quantile(values, fraction):
|
|
if not values:
|
|
return math.nan
|
|
ordered = sorted(values)
|
|
return ordered[math.floor(fraction*(len(ordered)-1))]
|
|
|
|
|
|
if len(sys.argv) != 3:
|
|
raise SystemExit("usage: CeilingContinuity.py scene.vfmesh expected-cap-plane-metres")
|
|
scene_path = sys.argv[1]
|
|
root_z = float(sys.argv[2])
|
|
tiles = {}
|
|
voxel_size_m = None
|
|
with open(scene_path, "rb") as stream:
|
|
if read_exact(stream, 8) != b"VFMSH001":
|
|
raise SystemExit("invalid VFMSH001 scene")
|
|
version, voxel_um, tile_count = struct.unpack("<3I", read_exact(stream, 12))
|
|
camera = struct.unpack("<3d", read_exact(stream, 24))
|
|
yaw, pitch = struct.unpack("<2d", read_exact(stream, 16))
|
|
voxel_size_m = voxel_um/1_000_000.0
|
|
for _ in range(tile_count):
|
|
key = struct.unpack("<4i", read_exact(stream, 16))
|
|
surface_vertices, skirt_count, vertex_count, index_count, mask_count = struct.unpack(
|
|
"<5I", read_exact(stream, 20))
|
|
positions = array.array("f")
|
|
positions.frombytes(read_exact(stream, vertex_count*12))
|
|
normals = array.array("f")
|
|
normals.frombytes(read_exact(stream, vertex_count*12))
|
|
indices = array.array("I")
|
|
indices.frombytes(read_exact(stream, index_count*4))
|
|
masks = read_exact(stream, mask_count)
|
|
if sys.byteorder != "little":
|
|
positions.byteswap()
|
|
normals.byteswap()
|
|
indices.byteswap()
|
|
extent = 32*(1 << key[3])*voxel_size_m
|
|
z0, z1 = key[2]*extent, (key[2]+1)*extent
|
|
ceiling = []
|
|
for triangle in range(index_count//3):
|
|
if masks and masks[triangle] != 0:
|
|
continue
|
|
a,b,c = indices[triangle*3:triangle*3+3]
|
|
if max(a,b,c) >= surface_vertices:
|
|
continue
|
|
nz = (normals[a*3+2]+normals[b*3+2]+normals[c*3+2])/3.0
|
|
if nz > -0.65:
|
|
continue
|
|
points = tuple((positions[v*3],positions[v*3+1],positions[v*3+2])
|
|
for v in (a,b,c))
|
|
center_z = sum(point[2] for point in points)/3.0
|
|
if abs(center_z-root_z) <= 1.0:
|
|
ceiling.append(points)
|
|
tiles[key] = {
|
|
"extent": extent, "x0": key[0]*extent, "y0": key[1]*extent,
|
|
"x1": (key[0]+1)*extent, "y1": (key[1]+1)*extent,
|
|
"z0": z0, "z1": z1, "ceiling": ceiling,
|
|
}
|
|
|
|
by_lod = collections.defaultdict(dict)
|
|
for key, tile in tiles.items():
|
|
if tile["z0"] <= root_z+1.0 and tile["z1"] >= root_z-1.0:
|
|
by_lod[key[3]][key] = tile
|
|
|
|
print(f"scene={scene_path} version={version} voxel_um={voxel_um} tiles={tile_count} "
|
|
f"camera_m={camera} yaw={yaw} pitch={pitch} expected_cap_m={root_z:.9f}")
|
|
print("lod zband_tiles cap_tiles missing_tiles cap_triangles face_z_m_p05/median/p95 "
|
|
"tile_samples hit_samples max_abs_root_delta_cm shared_edges expected_edge_samples "
|
|
"paired_edge_samples edge_gaps edge_delta_cm_median/p95/max")
|
|
all_pass = True
|
|
all_misses = []
|
|
for lod in range(9):
|
|
target = by_lod[lod]
|
|
cap_tiles = {key: tile for key,tile in target.items() if tile["ceiling"]}
|
|
face_z = [point[2] for tile in cap_tiles.values() for tri in tile["ceiling"]
|
|
for point in tri]
|
|
deltas = []
|
|
sample_count = 0
|
|
hit_count = 0
|
|
for key,tile in target.items():
|
|
if not tile["ceiling"]:
|
|
continue
|
|
for iy in range(16):
|
|
for ix in range(16):
|
|
x = tile["x0"]+(ix+0.5)*tile["extent"]/16
|
|
y = tile["y0"]+(iy+0.5)*tile["extent"]/16
|
|
sample_count += 1
|
|
height = height_at(tile["ceiling"],x,y,root_z)
|
|
if height is not None:
|
|
hit_count += 1
|
|
deltas.append(abs(height-root_z)*100.0)
|
|
else:
|
|
all_misses.append((lod,key,ix,iy,x,y))
|
|
|
|
edge_deltas = []
|
|
edge_count = 0
|
|
edge_gaps = 0
|
|
edge_gap_samples = 0
|
|
for key,a in cap_tiles.items():
|
|
x,y,z,level = key
|
|
for neighbor_key,axis in (((x+1,y,z,level),"x"),((x,y+1,z,level),"y")):
|
|
b = cap_tiles.get(neighbor_key)
|
|
if b is None:
|
|
continue
|
|
edge_count += 1
|
|
epsilon = min(0.02,a["extent"]*1.0e-4)
|
|
matched = 0
|
|
for index in range(33):
|
|
t = (index+0.5)/33
|
|
if axis == "x":
|
|
fixed = a["x1"]
|
|
varying = a["y0"]+(a["y1"]-a["y0"])*t
|
|
za = height_at(a["ceiling"],fixed-epsilon,varying,root_z)
|
|
zb = height_at(b["ceiling"],fixed+epsilon,varying,root_z)
|
|
else:
|
|
fixed = a["y1"]
|
|
varying = a["x0"]+(a["x1"]-a["x0"])*t
|
|
za = height_at(a["ceiling"],varying,fixed-epsilon,root_z)
|
|
zb = height_at(b["ceiling"],varying,fixed+epsilon,root_z)
|
|
if za is not None and zb is not None:
|
|
edge_deltas.append(abs(za-zb)*100.0)
|
|
matched += 1
|
|
edge_gap_samples += 33-matched
|
|
if matched != 33:
|
|
edge_gaps += 1
|
|
missing = len(target)-len(cap_tiles)
|
|
max_root = max(deltas) if deltas else math.nan
|
|
print(f"{lod} {len(target)} {len(cap_tiles)} {missing} "
|
|
f"{sum(len(tile['ceiling']) for tile in cap_tiles.values())} "
|
|
f"{quantile(face_z,.05):.6f}/{quantile(face_z,.5):.6f}/{quantile(face_z,.95):.6f} "
|
|
f"{sample_count} {hit_count} {max_root:.3f} {edge_count} {edge_count*33} {len(edge_deltas)} "
|
|
f"{edge_gaps} {quantile(edge_deltas,.5):.3f}/{quantile(edge_deltas,.95):.3f}/"
|
|
f"{max(edge_deltas) if edge_deltas else math.nan:.3f}")
|
|
if target and (missing or hit_count != sample_count or edge_gap_samples):
|
|
all_pass = False
|
|
for lod,key,ix,iy,x,y in all_misses:
|
|
print(f" miss lod={lod} key={key} sample=({ix},{iy}) xy_m=({x:.6f},{y:.6f})")
|
|
print(f"[CeilingContinuityL19] total_misses={len(all_misses)} "
|
|
f"result={'PASS' if all_pass else 'FAIL'} criteria=LOD0-8 every 16x16 interior ray hits; "
|
|
"every shared same-LOD edge sample has paired hits")
|
|
raise SystemExit(0 if all_pass else 1)
|