#!/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)