#!/usr/bin/env python3 """ ee_validate_mesh.py Validate an EEAE spherical triangle mesh sidecar produced by ee_mesh_build.py. This script is intentionally diagnostic-only. It does not modify the mesh and it does not attempt reverse-time deformation. It checks whether the mesh is a reasonable closed spherical triangulation and emits a richer diagnostics report that the EEAE can display in the Mesh / Solver cabinet. Primary outputs: 1. build1.mesh.validation.json 2. build1.mesh.validation.hotspots.geojson Validation categories: - schema and basic counts - vertex coordinate/unit-vector sanity - triangle index validity and degeneracy - edge manifold / Euler characteristic checks - adjacency consistency - edge-length, triangle-area, angle, and aspect-ratio distributions - Voronoi dual consistency - cell/material/age/young-crust summaries - zipper-anchor snap-distance summaries - hotspot GeoJSON for long edges, skinny triangles, and zipper snap outliers """ from __future__ import annotations import argparse import collections import datetime as _dt import json import math import statistics from pathlib import Path from typing import Any, Iterable, Optional DEFAULT_RADIUS_KM = 6371.0 SCHEMA_VERSION = "ee-mesh-validation-v1" def read_json(path: Path) -> Any: with path.open("r", encoding="utf-8") as f: return json.load(f) def write_json(path: Path, payload: Any, *, pretty: bool = False) -> None: path.parent.mkdir(parents=True, exist_ok=True) with path.open("w", encoding="utf-8") as f: if pretty: json.dump(payload, f, ensure_ascii=False, indent=2, allow_nan=False) else: json.dump(payload, f, ensure_ascii=False, separators=(",", ":"), allow_nan=False) f.write("\n") def is_finite_number(value: Any) -> bool: try: n = float(value) except (TypeError, ValueError): return False return math.isfinite(n) def normalize_lon(lon: float) -> float: out = ((float(lon) + 180.0) % 360.0) - 180.0 if out == -180.0 and lon > 0: return 180.0 return out def clamp(x: float, lo: float, hi: float) -> float: return max(lo, min(hi, x)) def vector_norm(v: list[float]) -> float: return math.sqrt(sum(float(x) * float(x) for x in v)) def normalize_vector(v: list[float]) -> list[float]: n = vector_norm(v) if not math.isfinite(n) or n <= 0: return [1.0, 0.0, 0.0] return [float(v[0]) / n, float(v[1]) / n, float(v[2]) / n] def dot(a: list[float], b: list[float]) -> float: return float(a[0]) * float(b[0]) + float(a[1]) * float(b[1]) + float(a[2]) * float(b[2]) def cross(a: list[float], b: list[float]) -> list[float]: return [ float(a[1]) * float(b[2]) - float(a[2]) * float(b[1]), float(a[2]) * float(b[0]) - float(a[0]) * float(b[2]), float(a[0]) * float(b[1]) - float(a[1]) * float(b[0]), ] def sub(a: list[float], b: list[float]) -> list[float]: return [float(a[0]) - float(b[0]), float(a[1]) - float(b[1]), float(a[2]) - float(b[2])] def scale(a: list[float], s: float) -> list[float]: return [float(a[0]) * s, float(a[1]) * s, float(a[2]) * s] def add(a: list[float], b: list[float]) -> list[float]: return [float(a[0]) + float(b[0]), float(a[1]) + float(b[1]), float(a[2]) + float(b[2])] def lonlat_to_unit(lon: float, lat: float) -> list[float]: lam = math.radians(float(lon)) phi = math.radians(float(lat)) cp = math.cos(phi) return normalize_vector([cp * math.cos(lam), cp * math.sin(lam), math.sin(phi)]) def unit_to_lonlat(unit: list[float]) -> list[float]: u = normalize_vector(unit) lon = normalize_lon(math.degrees(math.atan2(u[1], u[0]))) lat = math.degrees(math.asin(clamp(u[2], -1.0, 1.0))) return [round(lon, 8), round(lat, 8)] def central_angle_rad(a: list[float], b: list[float]) -> float: c = cross(a, b) return math.atan2(vector_norm(c), clamp(dot(a, b), -1.0, 1.0)) def edge_length_km(a: list[float], b: list[float], radius_km: float) -> float: return central_angle_rad(a, b) * radius_km def spherical_area_sr(a: list[float], b: list[float], c: list[float]) -> float: det = dot(a, cross(b, c)) denom = 1.0 + dot(a, b) + dot(b, c) + dot(c, a) area = 2.0 * math.atan2(abs(det), denom) if area < 0: area += 4.0 * math.pi return float(area) def tangent_angle_deg(center: list[float], a: list[float], b: list[float]) -> float: # Project both neighboring points into the tangent plane at center, then # measure the angle between those tangent vectors. ta = sub(a, scale(center, dot(center, a))) tb = sub(b, scale(center, dot(center, b))) nta = vector_norm(ta) ntb = vector_norm(tb) if nta <= 1e-15 or ntb <= 1e-15: return 0.0 ta = [ta[0] / nta, ta[1] / nta, ta[2] / nta] tb = [tb[0] / ntb, tb[1] / ntb, tb[2] / ntb] return math.degrees(math.acos(clamp(dot(ta, tb), -1.0, 1.0))) def triangle_angles_deg(a: list[float], b: list[float], c: list[float]) -> list[float]: return [ tangent_angle_deg(a, b, c), tangent_angle_deg(b, c, a), tangent_angle_deg(c, a, b), ] def slerp(a: list[float], b: list[float], t: float) -> list[float]: omega = central_angle_rad(a, b) if omega < 1e-12: return normalize_vector(add(scale(a, 1.0 - t), scale(b, t))) so = math.sin(omega) return normalize_vector(add( scale(a, math.sin((1.0 - t) * omega) / so), scale(b, math.sin(t * omega) / so), )) def interpolate_edge_lonlat(a: list[float], b: list[float], max_segment_km: float, radius_km: float) -> list[list[float]]: distance_km = edge_length_km(a, b, radius_km) steps = max(1, int(math.ceil(distance_km / max(1.0, float(max_segment_km))))) return [unit_to_lonlat(slerp(a, b, i / steps)) for i in range(steps + 1)] def split_at_antimeridian(path: list[list[float]]) -> list[list[list[float]]]: if len(path) < 2: return [] out: list[list[list[float]]] = [] current = [path[0]] for point in path[1:]: prev = current[-1] if abs(float(point[0]) - float(prev[0])) > 180.0: if len(current) >= 2: out.append(current) current = [point] else: current.append(point) if len(current) >= 2: out.append(current) return out def percentile(values: list[float], q: float) -> Optional[float]: if not values: return None if len(values) == 1: return float(values[0]) sorted_values = sorted(values) pos = (len(sorted_values) - 1) * float(q) lower = int(math.floor(pos)) upper = int(math.ceil(pos)) if lower == upper: return float(sorted_values[lower]) frac = pos - lower return float(sorted_values[lower] * (1.0 - frac) + sorted_values[upper] * frac) def stats(values: list[float], digits: int = 4) -> dict[str, Any]: finite = [float(v) for v in values if math.isfinite(float(v))] if not finite: return {"count": 0} return { "count": len(finite), "min": round(min(finite), digits), "p05": round(percentile(finite, 0.05), digits), "median": round(statistics.median(finite), digits), "mean": round(statistics.fmean(finite), digits), "p95": round(percentile(finite, 0.95), digits), "max": round(max(finite), digits), } def count_by(items: Iterable[Any]) -> dict[str, int]: counter = collections.Counter(str(item) for item in items if item is not None) return dict(sorted(counter.items(), key=lambda kv: (-kv[1], kv[0]))) def tri_indices(triangle: dict[str, Any], vertex_id_to_index: dict[str, int]) -> list[int]: if isinstance(triangle.get("vi"), list): out = [] for value in triangle["vi"]: try: out.append(int(value)) except (TypeError, ValueError): pass return out if isinstance(triangle.get("v"), list): out = [] for value in triangle["v"]: idx = vertex_id_to_index.get(str(value)) if idx is not None: out.append(idx) return out return [] def edge_key(a: int, b: int) -> tuple[int, int]: return (a, b) if a < b else (b, a) def vertex_record(raw: dict[str, Any]) -> dict[str, Any]: lon = float(raw.get("lon")) if is_finite_number(raw.get("lon")) else None lat = float(raw.get("lat")) if is_finite_number(raw.get("lat")) else None unit_raw = raw.get("unit") unit: list[float] if isinstance(unit_raw, list) and len(unit_raw) >= 3 and all(is_finite_number(x) for x in unit_raw[:3]): unit = [float(unit_raw[0]), float(unit_raw[1]), float(unit_raw[2])] elif lon is not None and lat is not None: unit = lonlat_to_unit(lon, lat) else: unit = [1.0, 0.0, 0.0] unit_norm = vector_norm(unit) unit = normalize_vector(unit) if lon is None or lat is None: lon, lat = unit_to_lonlat(unit) return { "id": str(raw.get("id", "")), "lon": lon, "lat": lat, "unit": unit, "rawUnitNorm": unit_norm, "unitNormError": abs(unit_norm - 1.0), "material": raw.get("material"), "cellType": raw.get("cellType"), "ageMa": raw.get("ageMa"), "groupId": raw.get("groupId"), "isMor": bool(raw.get("isMor")), "youngContinentalWeight": float(raw.get("youngContinentalWeight", 0.0) or 0.0), "seedKinds": raw.get("seedKinds") if isinstance(raw.get("seedKinds"), dict) else {}, } def add_warning(warnings: list[dict[str, Any]], severity: str, code: str, message: str, **extra: Any) -> None: warnings.append({"severity": severity, "code": code, "message": message, **extra}) def validation_status(warnings: list[dict[str, Any]]) -> str: severities = {w.get("severity") for w in warnings} if "fail" in severities: return "fail" if "warning" in severities: return "warning" return "pass" def make_edge_feature(mesh_id: str, source: str, a: list[float], b: list[float], props: dict[str, Any], radius_km: float, max_segment_km: float) -> list[dict[str, Any]]: path = interpolate_edge_lonlat(a, b, max_segment_km, radius_km) features = [] for idx, split_path in enumerate(split_at_antimeridian(path)): features.append({ "type": "Feature", "properties": { "schemaVersion": "ee-mesh-validation-hotspot-v1", "sourceMeshId": mesh_id, "kind": source, "splitIndex": idx, **props, }, "geometry": {"type": "LineString", "coordinates": split_path}, }) return features def build_hotspot_geojson( *, mesh_id: str, vertices: list[dict[str, Any]], long_edges: list[dict[str, Any]], skinny_triangles: list[dict[str, Any]], zipper_outliers: list[dict[str, Any]], radius_km: float, max_segment_km: float, max_features: int, ) -> dict[str, Any]: features: list[dict[str, Any]] = [] for item in sorted(long_edges, key=lambda x: x.get("lengthKm", 0), reverse=True)[:max_features]: a_idx, b_idx = item["edge"] a = vertices[a_idx]["unit"] b = vertices[b_idx]["unit"] features.extend(make_edge_feature( mesh_id, "long_edge", a, b, { "edge": [a_idx, b_idx], "lengthKm": round(float(item.get("lengthKm", 0)), 4), "triangleIds": item.get("triangleIds", []), }, radius_km, max_segment_km, )) # Draw the three edges of each worst skinny triangle. This avoids polygon # fill artifacts and handles antimeridian splitting edge-by-edge. skinny_budget = max(0, max_features - len(features)) for item in sorted(skinny_triangles, key=lambda x: x.get("minAngleDeg", 999))[:skinny_budget]: vi = item.get("vi", []) if len(vi) != 3: continue for a_idx, b_idx in [(vi[0], vi[1]), (vi[1], vi[2]), (vi[2], vi[0])]: a = vertices[a_idx]["unit"] b = vertices[b_idx]["unit"] features.extend(make_edge_feature( mesh_id, "skinny_triangle_edge", a, b, { "triangleId": item.get("triangleId"), "minAngleDeg": round(float(item.get("minAngleDeg", 0)), 4), "aspectRatio": round(float(item.get("aspectRatio", 0)), 4), }, radius_km, max_segment_km, )) point_budget = max(0, max_features - len(features)) for item in sorted(zipper_outliers, key=lambda x: x.get("snapDistanceKm", 0), reverse=True)[:point_budget]: v_idx = item.get("meshVertexIndex") if not isinstance(v_idx, int) or v_idx < 0 or v_idx >= len(vertices): continue features.append({ "type": "Feature", "properties": { "schemaVersion": "ee-mesh-validation-hotspot-v1", "sourceMeshId": mesh_id, "kind": "zipper_snap_outlier", "zipperControlId": item.get("zipperControlId"), "anchor": item.get("anchor"), "snapDistanceKm": round(float(item.get("snapDistanceKm", 0)), 4), "meshVertexId": item.get("meshVertexId"), }, "geometry": {"type": "Point", "coordinates": [vertices[v_idx]["lon"], vertices[v_idx]["lat"]]}, }) return { "type": "FeatureCollection", "name": "ee_mesh_validation_hotspots", "properties": { "schemaVersion": "ee-mesh-validation-hotspots-v1", "sourceMeshId": mesh_id, "featureCount": len(features), "maxFeatures": max_features, }, "features": features[:max_features], } def validate_mesh(mesh: dict[str, Any], args: argparse.Namespace) -> tuple[dict[str, Any], dict[str, Any]]: mesh_id = str(mesh.get("meshId") or "unknown_mesh") params = mesh.get("parameters") if isinstance(mesh.get("parameters"), dict) else {} counts_in = mesh.get("counts") if isinstance(mesh.get("counts"), dict) else {} radius_km = float(args.radius_km or params.get("radius_km") or DEFAULT_RADIUS_KM) target_spacing_km = float(args.target_spacing_km or params.get("target_spacing_km") or 180.0) long_edge_threshold_km = float(args.long_edge_threshold_km or max(target_spacing_km * 4.0, target_spacing_km + 300.0)) small_angle_threshold_deg = float(args.small_angle_threshold_deg) aspect_ratio_threshold = float(args.aspect_ratio_threshold) zipper_snap_warning_km = float(args.zipper_snap_warning_km) unit_norm_tolerance = float(args.unit_norm_tolerance) warnings: list[dict[str, Any]] = [] raw_vertices = mesh.get("vertices") raw_triangles = mesh.get("triangles") if not isinstance(raw_vertices, list): raw_vertices = [] add_warning(warnings, "fail", "missing_vertices", "Mesh JSON does not contain a vertices array.") if not isinstance(raw_triangles, list): raw_triangles = [] add_warning(warnings, "fail", "missing_triangles", "Mesh JSON does not contain a triangles array.") vertices = [vertex_record(v if isinstance(v, dict) else {}) for v in raw_vertices] vertex_id_to_index = {v["id"]: idx for idx, v in enumerate(vertices) if v.get("id")} unit_errors = [float(v["unitNormError"]) for v in vertices] bad_unit_count = sum(1 for err in unit_errors if err > unit_norm_tolerance) invalid_lonlat_count = sum( 1 for v in vertices if not (math.isfinite(float(v["lon"])) and math.isfinite(float(v["lat"])) and abs(float(v["lat"])) <= 90.000001) ) if bad_unit_count: add_warning( warnings, "warning", "unit_norm_error", f"{bad_unit_count} vertices have unit-vector norm error above tolerance.", tolerance=unit_norm_tolerance, ) if invalid_lonlat_count: add_warning( warnings, "fail", "invalid_lonlat", f"{invalid_lonlat_count} vertices have invalid lon/lat coordinates.", ) valid_triangles: list[dict[str, Any]] = [] invalid_triangle_count = 0 duplicate_triangle_counter: collections.Counter[tuple[int, int, int]] = collections.Counter() edge_counts: collections.Counter[tuple[int, int]] = collections.Counter() edge_triangle_ids: collections.defaultdict[tuple[int, int], list[str]] = collections.defaultdict(list) edge_lengths: list[float] = [] triangle_areas_sr: list[float] = [] triangle_areas_km2: list[float] = [] min_angles: list[float] = [] max_angles: list[float] = [] aspect_ratios: list[float] = [] long_edges: list[dict[str, Any]] = [] skinny_triangles: list[dict[str, Any]] = [] degenerate_triangles: list[str] = [] cell_types: list[str] = [] age_values: list[float] = [] young_weights: list[float] = [float(v["youngContinentalWeight"]) for v in vertices] for idx, raw_tri in enumerate(raw_triangles): tri = raw_tri if isinstance(raw_tri, dict) else {} tid = str(tri.get("id") or f"tri_{idx:06d}") vi = tri_indices(tri, vertex_id_to_index) if len(vi) != 3 or len(set(vi)) != 3 or any(i < 0 or i >= len(vertices) for i in vi): invalid_triangle_count += 1 continue tri_key = tuple(sorted(vi)) duplicate_triangle_counter[tri_key] += 1 a, b, c = [vertices[i]["unit"] for i in vi] lengths = [ edge_length_km(a, b, radius_km), edge_length_km(b, c, radius_km), edge_length_km(c, a, radius_km), ] area_sr = spherical_area_sr(a, b, c) area_km2 = area_sr * radius_km * radius_km angles = triangle_angles_deg(a, b, c) min_angle = min(angles) max_angle = max(angles) shortest = min(lengths) longest = max(lengths) aspect = (longest / shortest) if shortest > 1e-12 else float("inf") edge_pairs = [(vi[0], vi[1]), (vi[1], vi[2]), (vi[2], vi[0])] for (edge_a, edge_b), length in zip(edge_pairs, lengths): key = edge_key(edge_a, edge_b) edge_counts[key] += 1 edge_triangle_ids[key].append(tid) edge_lengths.append(length) if length > long_edge_threshold_km: long_edges.append({ "edge": [key[0], key[1]], "lengthKm": length, "triangleIds": edge_triangle_ids[key], }) triangle_areas_sr.append(area_sr) triangle_areas_km2.append(area_km2) min_angles.append(min_angle) max_angles.append(max_angle) aspect_ratios.append(aspect) cell_types.append(str(tri.get("cellType") or "unknown")) if is_finite_number(tri.get("ageMa")): age_values.append(float(tri.get("ageMa"))) if area_sr <= args.degenerate_area_sr or min_angle <= 1e-8 or not math.isfinite(aspect): degenerate_triangles.append(tid) if min_angle < small_angle_threshold_deg or aspect > aspect_ratio_threshold: skinny_triangles.append({ "triangleId": tid, "vi": vi, "minAngleDeg": min_angle, "maxAngleDeg": max_angle, "aspectRatio": aspect, "areaKm2": area_km2, }) valid_triangles.append({ "id": tid, "vi": vi, "adjacent": tri.get("adjacent") if isinstance(tri.get("adjacent"), list) else [], }) duplicate_triangle_count = sum(count - 1 for count in duplicate_triangle_counter.values() if count > 1) if invalid_triangle_count: add_warning(warnings, "fail", "invalid_triangles", f"{invalid_triangle_count} triangles have invalid vertex indices or repeated vertices.") if duplicate_triangle_count: add_warning(warnings, "fail", "duplicate_triangles", f"{duplicate_triangle_count} duplicate triangle definitions were detected.") if degenerate_triangles: add_warning(warnings, "fail", "degenerate_triangles", f"{len(degenerate_triangles)} triangles are degenerate or nearly zero-area.") edge_count = len(edge_counts) boundary_edges = [edge for edge, count in edge_counts.items() if count == 1] nonmanifold_edges = [edge for edge, count in edge_counts.items() if count > 2] euler_characteristic = len(vertices) - edge_count + len(valid_triangles) if boundary_edges: add_warning(warnings, "fail", "boundary_edges", f"{len(boundary_edges)} edges have only one incident triangle; a closed spherical mesh should have zero.") if nonmanifold_edges: add_warning(warnings, "fail", "nonmanifold_edges", f"{len(nonmanifold_edges)} edges have more than two incident triangles.") if euler_characteristic != 2: add_warning(warnings, "warning", "euler_characteristic", f"Euler characteristic is {euler_characteristic}, expected 2 for a closed triangulated sphere.") if long_edges: add_warning(warnings, "warning", "long_edges", f"{len(long_edges)} triangle edges exceed {long_edge_threshold_km:.1f} km.", thresholdKm=long_edge_threshold_km) if skinny_triangles: add_warning(warnings, "warning", "skinny_triangles", f"{len(skinny_triangles)} triangles violate min-angle/aspect-ratio thresholds.", smallAngleThresholdDeg=small_angle_threshold_deg, aspectRatioThreshold=aspect_ratio_threshold) # Adjacency checks. tri_by_id = {tri["id"]: tri for tri in valid_triangles} adjacency_ref_count = 0 invalid_adjacency_refs = 0 nonreciprocal_adjacency = 0 adjacency_without_shared_edge = 0 adjacency_size_counter: collections.Counter[int] = collections.Counter() for tri in valid_triangles: adj = [str(x) for x in tri.get("adjacent", [])] adjacency_size_counter[len(adj)] += 1 for neighbor_id in adj: adjacency_ref_count += 1 neighbor = tri_by_id.get(neighbor_id) if neighbor is None: invalid_adjacency_refs += 1 continue if tri["id"] not in [str(x) for x in neighbor.get("adjacent", [])]: nonreciprocal_adjacency += 1 if len(set(tri["vi"]).intersection(neighbor["vi"])) != 2: adjacency_without_shared_edge += 1 if invalid_adjacency_refs: add_warning(warnings, "fail", "invalid_adjacency_refs", f"{invalid_adjacency_refs} adjacency references point to missing triangles.") if nonreciprocal_adjacency: add_warning(warnings, "warning", "nonreciprocal_adjacency", f"{nonreciprocal_adjacency} adjacency references are not reciprocal.") if adjacency_without_shared_edge: add_warning(warnings, "warning", "adjacency_without_shared_edge", f"{adjacency_without_shared_edge} adjacency references do not share an edge.") # Voronoi checks. voronoi = mesh.get("voronoi") if isinstance(mesh.get("voronoi"), dict) else {} voronoi_status = str(voronoi.get("status") or "missing") voronoi_vertices = voronoi.get("vertices") if isinstance(voronoi.get("vertices"), list) else [] voronoi_regions = voronoi.get("regions") if isinstance(voronoi.get("regions"), list) else [] voronoi_region_sizes = [len(region) for region in voronoi_regions if isinstance(region, list)] invalid_voronoi_region_indices = 0 for region in voronoi_regions: if not isinstance(region, list): invalid_voronoi_region_indices += 1 continue for idx in region: try: vi = int(idx) except (TypeError, ValueError): invalid_voronoi_region_indices += 1 continue if vi < 0 or vi >= len(voronoi_vertices): invalid_voronoi_region_indices += 1 if voronoi_status != "ok": add_warning(warnings, "warning", "voronoi_status", f"Voronoi status is {voronoi_status!r}, expected 'ok'.") if voronoi_regions and len(voronoi_regions) != len(vertices): add_warning(warnings, "warning", "voronoi_region_count", f"Voronoi region count is {len(voronoi_regions)}, expected vertex count {len(vertices)}.") if invalid_voronoi_region_indices: add_warning(warnings, "fail", "invalid_voronoi_region_indices", f"{invalid_voronoi_region_indices} invalid Voronoi region references were detected.") # Zipper snap distances and other constraint summaries. anchors = mesh.get("constraintAnchors") if isinstance(mesh.get("constraintAnchors"), dict) else {} zipper_controls = anchors.get("zipperControls") if isinstance(anchors.get("zipperControls"), list) else [] zipper_snap_distances: list[float] = [] zipper_outliers: list[dict[str, Any]] = [] for zipper in zipper_controls: if not isinstance(zipper, dict): continue for anchor_name in ("fromIsochron", "toMor", "oppositeIsochron"): anchor = zipper.get(anchor_name) if not isinstance(anchor, dict): continue if not is_finite_number(anchor.get("snapDistanceKm")): continue distance = float(anchor.get("snapDistanceKm")) zipper_snap_distances.append(distance) if distance > zipper_snap_warning_km: mesh_vertex_id = str(anchor.get("nearestMeshVertexId") or "") zipper_outliers.append({ "zipperControlId": zipper.get("id"), "anchor": anchor_name, "snapDistanceKm": distance, "meshVertexId": mesh_vertex_id, "meshVertexIndex": vertex_id_to_index.get(mesh_vertex_id, -1), }) if zipper_outliers: add_warning(warnings, "warning", "zipper_snap_outliers", f"{len(zipper_outliers)} zipper anchor snaps exceed {zipper_snap_warning_km:.2f} km.", thresholdKm=zipper_snap_warning_km) young_vertex_count = sum(1 for w in young_weights if w > 0) strong_young_vertex_count = sum(1 for w in young_weights if w >= 0.5) mor_vertex_count = sum(1 for v in vertices if v["isMor"]) source_anchor_vertex_count = sum(1 for v in vertices if v["seedKinds"]) validation = { "schemaVersion": SCHEMA_VERSION, "generatedAtUtc": _dt.datetime.now(_dt.UTC).replace(microsecond=0).isoformat().replace("+00:00", "Z"), "sourceMeshId": mesh_id, "sourceMeshSchemaVersion": mesh.get("schemaVersion"), "status": validation_status(warnings), "parameters": { "radiusKm": radius_km, "targetSpacingKm": target_spacing_km, "longEdgeThresholdKm": long_edge_threshold_km, "smallAngleThresholdDeg": small_angle_threshold_deg, "aspectRatioThreshold": aspect_ratio_threshold, "zipperSnapWarningKm": zipper_snap_warning_km, "unitNormTolerance": unit_norm_tolerance, }, "counts": { "vertices": len(vertices), "triangles": len(raw_triangles), "validTriangles": len(valid_triangles), "invalidTriangles": invalid_triangle_count, "edges": edge_count, "boundaryEdges": len(boundary_edges), "nonmanifoldEdges": len(nonmanifold_edges), "duplicateTriangles": duplicate_triangle_count, "degenerateTriangles": len(degenerate_triangles), "longEdges": len(long_edges), "skinnyTriangles": len(skinny_triangles), "eulerCharacteristic": euler_characteristic, }, "sourceCounts": counts_in, "topology": { "closedSphereExpectedEulerCharacteristic": 2, "eulerCharacteristic": euler_characteristic, "edgeIncidence": { "incidence1": len(boundary_edges), "incidence2": sum(1 for count in edge_counts.values() if count == 2), "incidenceGt2": len(nonmanifold_edges), }, "adjacency": { "references": adjacency_ref_count, "invalidReferences": invalid_adjacency_refs, "nonReciprocalReferences": nonreciprocal_adjacency, "withoutSharedEdge": adjacency_without_shared_edge, "adjacencySizeHistogram": {str(k): v for k, v in sorted(adjacency_size_counter.items())}, }, }, "quality": { "edgeLengthKm": stats(edge_lengths), "triangleAreaSteradians": stats(triangle_areas_sr, digits=8), "triangleAreaKm2": stats(triangle_areas_km2), "minAngleDeg": stats(min_angles), "maxAngleDeg": stats(max_angles), "aspectRatio": stats(aspect_ratios), "unitNormError": stats(unit_errors, digits=10), }, "classification": { "cellTypes": count_by(cell_types), "triangleAgeMa": stats(age_values) if age_values else {"count": 0}, "materials": count_by(v.get("material") for v in vertices), "groups": { "vertexGroupCount": len({v.get("groupId") for v in vertices if v.get("groupId")}), }, "sourceAnchors": { "verticesWithSeedKinds": source_anchor_vertex_count, "morVertices": mor_vertex_count, }, }, "youngContinentalCrust": { "vertexCountWithWeightGt0": young_vertex_count, "vertexCountWithWeightGte0_5": strong_young_vertex_count, "weight": stats(young_weights, digits=6), }, "voronoi": { "status": voronoi_status, "vertexCount": len(voronoi_vertices), "regionCount": len(voronoi_regions), "regionSize": stats([float(x) for x in voronoi_region_sizes]), "invalidRegionReferences": invalid_voronoi_region_indices, }, "constraints": { "zipperControls": len(zipper_controls), "zipperAnchorSnapDistanceKm": stats(zipper_snap_distances), "zipperSnapOutlierCount": len(zipper_outliers), "groupAssociations": len(anchors.get("groupAssociations") or []) if isinstance(anchors.get("groupAssociations"), list) else 0, "morAssignments": len(anchors.get("morAssignments") or []) if isinstance(anchors.get("morAssignments"), list) else 0, }, "warnings": warnings, "hotspotSummary": { "longEdgesConsidered": len(long_edges), "skinnyTrianglesConsidered": len(skinny_triangles), "zipperSnapOutliersConsidered": len(zipper_outliers), }, } hotspots = build_hotspot_geojson( mesh_id=mesh_id, vertices=vertices, long_edges=long_edges, skinny_triangles=skinny_triangles, zipper_outliers=zipper_outliers, radius_km=radius_km, max_segment_km=float(args.hotspot_segment_km), max_features=int(args.max_hotspot_features), ) validation["hotspotSummary"]["featureCount"] = len(hotspots.get("features", [])) return validation, hotspots def build_arg_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(description="Validate an EEAE mesh sidecar.") parser.add_argument("mesh_json", type=Path, help="Path to build1.mesh.v1.json or another ee-mesh-v1 sidecar.") parser.add_argument("--out", type=Path, required=True, help="Output validation JSON path.") parser.add_argument("--hotspots-out", type=Path, required=True, help="Output validation hotspot GeoJSON path.") parser.add_argument("--radius-km", type=float, default=None, help="Override radius for geometric measurements.") parser.add_argument("--target-spacing-km", type=float, default=None, help="Override target spacing for threshold derivation.") parser.add_argument("--long-edge-threshold-km", type=float, default=None, help="Warn when an edge exceeds this length.") parser.add_argument("--small-angle-threshold-deg", type=float, default=8.0, help="Warn when a triangle min angle is below this threshold.") parser.add_argument("--aspect-ratio-threshold", type=float, default=8.0, help="Warn when triangle longest/shortest edge exceeds this threshold.") parser.add_argument("--zipper-snap-warning-km", type=float, default=1.0, help="Warn when zipper anchor snap distance exceeds this value.") parser.add_argument("--unit-norm-tolerance", type=float, default=1e-6, help="Warn when raw unit vector norm differs from 1 by more than this.") parser.add_argument("--degenerate-area-sr", type=float, default=1e-14, help="Fail triangles below this spherical area threshold.") parser.add_argument("--hotspot-segment-km", type=float, default=80.0, help="Maximum display segment length for hotspot line interpolation.") parser.add_argument("--max-hotspot-features", type=int, default=750, help="Maximum hotspot features to emit.") parser.add_argument("--pretty", action="store_true", help="Pretty-print output JSON.") return parser def main(argv: Optional[list[str]] = None) -> int: parser = build_arg_parser() args = parser.parse_args(argv) mesh = read_json(args.mesh_json) validation, hotspots = validate_mesh(mesh, args) write_json(args.out, validation, pretty=args.pretty) write_json(args.hotspots_out, hotspots, pretty=False) print(json.dumps({ "validationOut": str(args.out), "hotspotsOut": str(args.hotspots_out), "sourceMeshId": validation.get("sourceMeshId"), "status": validation.get("status"), "warnings": len(validation.get("warnings", [])), "vertices": validation.get("counts", {}).get("vertices"), "triangles": validation.get("counts", {}).get("triangles"), "edges": validation.get("counts", {}).get("edges"), "boundaryEdges": validation.get("counts", {}).get("boundaryEdges"), "nonmanifoldEdges": validation.get("counts", {}).get("nonmanifoldEdges"), "longEdges": validation.get("counts", {}).get("longEdges"), "skinnyTriangles": validation.get("counts", {}).get("skinnyTriangles"), "hotspotFeatures": validation.get("hotspotSummary", {}).get("featureCount"), }, indent=2)) return 0 if validation.get("status") != "fail" or args.pretty else 0 if __name__ == "__main__": raise SystemExit(main())