diff --git a/blender/scripts/export.py b/blender/scripts/export.py index 13a886c..5e90948 100644 --- a/blender/scripts/export.py +++ b/blender/scripts/export.py @@ -21,6 +21,7 @@ from __future__ import annotations import argparse import json +import math import sys from pathlib import Path @@ -53,8 +54,8 @@ def _require_object(name: str, expected_type: str | None = None) -> bpy.types.Ob return obj -def _load_job(path: Path) -> tuple[str | None, dict, list[str] | None]: - """Load output path, objects map, and optional housing collections.""" +def _load_job(path: Path) -> tuple[str | None, dict, list[str] | None, dict | None]: + """Load output path, objects map, housing collections and the optics stack.""" data = json.loads(path.read_text(encoding="utf-8")) if not isinstance(data, dict): raise ValueError(f"Export job must be a JSON object: {path}") @@ -70,7 +71,13 @@ def _load_job(path: Path) -> tuple[str | None, dict, list[str] | None]: raise ValueError( f"Export job 'housing_collections' must be a list of names: {path}" ) - return (str(output) if output else None), objects, housing + optics = data.get("optics") + if optics is not None: + if not isinstance(optics, dict) or not isinstance(optics.get("layers"), list): + raise ValueError( + f"Export job 'optics' must be a mapping with 'layers': {path}" + ) + return (str(output) if output else None), objects, housing, optics def export_camera_intrinsics(camera_obj: bpy.types.Object, scene: bpy.types.Scene): @@ -255,6 +262,124 @@ def export_housing_in_sensor_frame( return len(part_names), int(len(vertices)), int(len(faces)) +def _mesh_world(obj: bpy.types.Object): + """Triangulated mesh plus its world matrix, with a cleanup callback.""" + return _triangulated_mesh(obj) + + +def _measure_layer_thickness( + outer_obj: bpy.types.Object, inner_obj: bpy.types.Object +) -> dict[str, float | int | None]: + """Thickness of one layer, by casting from the outer surface onto the inner. + + Each outer vertex is sent inward along its negated vertex normal and the + distance to the inner surface is recorded, which measures the real layer + rather than the nearest-neighbour gap between two point clouds (that + overestimates wherever the surfaces are oblique to each other). + """ + from mathutils.bvhtree import BVHTree + + inner_mesh, inner_world, inner_cleanup = _mesh_world(inner_obj) + try: + inner_vertices = [inner_world @ vertex.co for vertex in inner_mesh.vertices] + inner_faces = [tuple(poly.vertices) for poly in inner_mesh.polygons] + if not inner_faces: + return {"count": 0, "p05": None, "p50": None, "p95": None} + tree = BVHTree.FromPolygons(inner_vertices, inner_faces, all_triangles=True) + finally: + inner_cleanup() + + outer_mesh, outer_world, outer_cleanup = _mesh_world(outer_obj) + try: + rotation = outer_world.to_3x3() + distances: list[float] = [] + for vertex in outer_mesh.vertices: + origin = outer_world @ vertex.co + normal = rotation @ vertex.normal + if normal.length <= 0.0: + continue + normal.normalize() + hit = tree.ray_cast(origin - normal * 1e-4, -normal) + location = hit[0] + if location is None: + continue + distances.append(float((location - origin).length)) + finally: + outer_cleanup() + + if not distances: + return {"count": 0, "p05": None, "p50": None, "p95": None} + distances.sort() + + def pick(fraction: float) -> float: + position = (len(distances) - 1) * fraction + low = int(math.floor(position)) + high = int(math.ceil(position)) + if low == high: + return distances[low] + return distances[low] + (distances[high] - distances[low]) * (position - low) + + return { + "count": len(distances), + "p05": pick(0.05), + "p50": pick(0.50), + "p95": pick(0.95), + } + + +def export_optics_in_sensor_frame( + layers: list[dict], + t_sensor_from_world: Matrix, + filepath: Path, +) -> int: + """Export the refracting interfaces, in the sensor frame, one pair per layer.""" + payload: dict = {"layer_names": np.array([layer["name"] for layer in layers], dtype=np.str_)} + thickness_p50: list[float] = [] + ior_values: list[float] = [] + + for index, layer in enumerate(layers): + if layer.get("solid"): + targets = [("solid", _require_object(str(layer["solid"]), "MESH"))] + else: + targets = [ + ("outer", _require_object(str(layer["outer"]), "MESH")), + ("inner", _require_object(str(layer["inner"]), "MESH")), + ] + for prefix, obj in targets: + mesh, world, cleanup = _mesh_world(obj) + try: + transform = t_sensor_from_world @ world + payload[f"layer{index}_{prefix}_vertices"] = np.array( + [transform @ vertex.co for vertex in mesh.vertices], dtype=np.float32 + ) + payload[f"layer{index}_{prefix}_faces"] = np.asarray( + [tuple(poly.vertices) for poly in mesh.polygons], dtype=np.int32 + ) + finally: + cleanup() + if layer.get("solid"): + # Thickness of a solid: enter one wall and leave by the opposite one. + measured = _measure_layer_thickness(targets[0][1], targets[0][1]) + else: + measured = _measure_layer_thickness(targets[0][1], targets[1][1]) + declared = layer.get("thickness_mm") + thickness = float(declared) if declared is not None else (measured["p50"] or 0.0) + thickness_p50.append(thickness) + ior_values.append(float(layer["ior"])) + print( + f" optics layer {layer['name']!r}: ior={float(layer['ior']):g} " + f"thickness={thickness:.4f} mm " + f"({'declared' if declared is not None else 'measured p50'}), " + f"measured p05={measured['p05']} p95={measured['p95']}" + ) + + payload["ior"] = np.asarray(ior_values, dtype=np.float64) + payload["thickness_mm"] = np.asarray(thickness_p50, dtype=np.float64) + np.savez_compressed(filepath, **payload) + print(f"Optics exported to {filepath} ({len(layers)} layer(s))") + return len(layers) + + def export_camera_pose_sensor_frame( camera_obj: bpy.types.Object, t_sensor_from_world: Matrix, @@ -307,16 +432,17 @@ def parse_args(argv: list[str]) -> argparse.Namespace: def _resolve_config( args: argparse.Namespace, -) -> tuple[Path, str, str, str, list[str]]: +) -> tuple[Path, str, str, str, list[str], dict | None]: """Merge --job JSON with explicit flag overrides. Explicit flags win.""" job_output: str | None = None job_objects: dict = {} job_housing: list[str] | None = None + job_optics: dict | None = None if args.job is not None: job_path = args.job.resolve() if not job_path.is_file(): raise RuntimeError(f"Export job not found: {job_path}") - job_output, job_objects, job_housing = _load_job(job_path) + job_output, job_objects, job_housing, job_optics = _load_job(job_path) output_raw = args.output or job_output if not output_raw: @@ -342,13 +468,20 @@ def _resolve_config( ] elif job_housing is not None: housing = list(job_housing) - return Path(output_raw).resolve(), str(mesh), str(cam), str(origin), housing + return ( + Path(output_raw).resolve(), + str(mesh), + str(cam), + str(origin), + housing, + job_optics, + ) def main() -> int: args = parse_args(_argv_after_double_dash()) try: - output_dir, mesh_name, cam_name, origin_name, housing_collections = ( + output_dir, mesh_name, cam_name, origin_name, housing_collections, optics = ( _resolve_config(args) ) except Exception as exc: @@ -378,6 +511,10 @@ def main() -> int: export_housing_in_sensor_frame( _housing_objects(housing_collections), t_sensor_from_world, housing_path ) + if optics and optics.get("layers"): + export_optics_in_sensor_frame( + list(optics["layers"]), t_sensor_from_world, output_dir / "optics.npz" + ) K, width, height = export_camera_intrinsics(cam_obj, scene) np.savez( intrinsics_path, diff --git a/pyproject.toml b/pyproject.toml index 8161232..936b1ac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,15 @@ preview = [ "PyQt5>=5.15", "PyOpenGL>=3.1", ] +# The refraction check only needs arrays. +optics = [ + "numpy>=1.24", +] +# The side-by-side comparison figure. +viz = [ + "numpy>=1.24", + "matplotlib>=3.7", +] [project.scripts] orisys-design = "orisys_sensor_design.cli:main" diff --git a/sensors/finger-xynova/r01/design.yaml b/sensors/finger-xynova/r01/design.yaml index 8b18d46..466a576 100644 --- a/sensors/finger-xynova/r01/design.yaml +++ b/sensors/finger-xynova/r01/design.yaml @@ -63,6 +63,30 @@ printability: max_sliver_ratio: 4.0 # longest edge / minimum width min_patch_area_mm2: 1.0e-4 # area above which a polygon is a material feature +# Optical stack above the markers, ordered from the camera inward. +# The markers sit on the sensing surface underneath these layers, so their image +# is displaced relative to the straight-line pinhole model. +# +# ior values are LITERATURE ASSUMPTIONS, not measurements - that is what +# ior_source records, so an assumed index can never pass for a measured one. +# Replace them with measured values when available; nothing else needs to change. +optics: + marker_surface: sensing_surface + layers: + # Solids, not surface pairs: `shell` and `rubber-1` are each one closed mesh, + # so the ray enters and leaves the same object. + # Stack: camera > shell > rubber-1 > sensing_surface > rubber-2. + # rubber-2 is BEHIND the markers, so it does not refract camera light and is + # deliberately not a layer. + - name: plastic_shell + solid: shell + ior: 1.59 + ior_source: assumed_pc_literature + - name: clear_rubber + solid: rubber-1 + ior: 1.43 + ior_source: assumed_silicone_literature + collections: frames: FRAMES export_compute: EXPORT_COMPUTE diff --git a/sensors/finger-xynova/r01/profiles/compute.yaml b/sensors/finger-xynova/r01/profiles/compute.yaml index 3e1c3f0..5812710 100644 --- a/sensors/finger-xynova/r01/profiles/compute.yaml +++ b/sensors/finger-xynova/r01/profiles/compute.yaml @@ -11,6 +11,10 @@ artifacts: - housing.npz - camera_intrinsics.npz - world_to_cam.npy + # Refracting interfaces for the optical stack, so the SDK's rectification stage + # can account for refraction (workflow §7). Additive; absent when no stack is + # declared in design.yaml. + - optics.npz # Measured unroll maps stay in SDK (not overwritten by nominal export). runtime_owned_by_sdk: diff --git a/sensors/finger-xynova/r01/source/geometry.blend b/sensors/finger-xynova/r01/source/geometry.blend index 85cab04..01b3e2c 100644 Binary files a/sensors/finger-xynova/r01/source/geometry.blend and b/sensors/finger-xynova/r01/source/geometry.blend differ diff --git a/src/orisys_sensor_design/cli.py b/src/orisys_sensor_design/cli.py index 8abab8c..e73177d 100644 --- a/src/orisys_sensor_design/cli.py +++ b/src/orisys_sensor_design/cli.py @@ -212,7 +212,30 @@ def cmd_export( if housing_collections else None ) - job = build_export_job(paths, output_dir, profile, collections) + from orisys_sensor_design.optics import OpticsConfigError, resolve_optics + + try: + stack = resolve_optics(paths) + except OpticsConfigError as exc: + for message in str(exc).splitlines(): + print(f"ERROR: {message}", file=sys.stderr) + return 2 + optics_layers = ( + [ + { + "name": layer.name, + "solid": layer.solid, + "outer": layer.outer, + "inner": layer.inner, + "ior": layer.ior, + "thickness_mm": layer.thickness_mm, + } + for layer in stack.layers + ] + if stack is not None + else None + ) + job = build_export_job(paths, output_dir, profile, collections, optics_layers) job_path = write_export_job(job, export_revision_dir(root, paths, profile)) print(f"Wrote export job: {job_path}") @@ -542,6 +565,211 @@ def cmd_preview(revision_dir: Path, args: argparse.Namespace) -> int: return code +def cmd_optics(revision_dir: Path, args: argparse.Namespace) -> int: + """Quantify refraction through the declared optical stack. Advisory.""" + root = repo_root(revision_dir) + paths = RevisionPaths.from_revision_dir(revision_dir) + errors = validate_revision_metadata(paths) + if errors: + for message in errors: + print(f"ERROR: {message}", file=sys.stderr) + return 1 + + try: + import numpy as np + except ImportError: + print( + "ERROR: the optics check needs numpy.\n" + "Install it with: python -m pip install -e .[optics]", + file=sys.stderr, + ) + return 2 + + from orisys_sensor_design.optics import ( + OpticsConfigError, + build_report, + camera_rays, + load_interfaces, + resolve_optics, + sample_stack, + summarize, + write_report, + ) + from orisys_sensor_design.optics import REPORT_FILENAME as OPTICS_REPORT_FILENAME + + want_plot = not args.no_plot + if want_plot: + try: + import matplotlib # noqa: F401 + except ImportError: + print( + "ERROR: --plot needs matplotlib.\n" + "Install it with: python -m pip install -e .[viz]\n" + "Or re-run with --no-plot.", + file=sys.stderr, + ) + return 2 + + try: + stack = resolve_optics(paths) + except OpticsConfigError as exc: + for message in str(exc).splitlines(): + print(f"ERROR: {message}", file=sys.stderr) + return 2 + if stack is None: + print( + f"ERROR: {paths.design_yaml} declares no 'optics' block; " + "nothing to trace", + file=sys.stderr, + ) + return 2 + + export_dir = ( + args.output_dir.resolve() + if args.output_dir is not None + else export_output_dir(root, paths, args.profile) + ) + try: + layers, marker, measured_thickness, marker_at_last_exit = load_interfaces( + export_dir, stack + ) + except OpticsConfigError as exc: + print(f"ERROR: {exc}", file=sys.stderr) + return 3 + + with np.load(export_dir / "camera_intrinsics.npz", allow_pickle=False) as handle: + K = np.asarray(handle["K"], dtype=np.float64) + width = int(handle["width"]) + height = int(handle["height"]) + world_to_cam = np.load(export_dir / "world_to_cam.npy").astype(np.float64) + + origin, directions, pixel_u, pixel_v = camera_rays( + K, world_to_cam, width, height, args.samples + ) + + # Marker pitch on the surface: the printed page carries nx cells across. + texture = resolve_texture_settings(paths) + pitch_mm = texture.vector_width_mm / max(texture.nx, 1) + + samples = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=pitch_mm, + measured_thickness=measured_thickness, + marker_at_last_exit=marker_at_last_exit, + pixel_u=pixel_u, + pixel_v=pixel_v, + image_size=(width, height), + ) + result = summarize(samples) + report = build_report( + paths=paths, stack=stack, result=result, export_dir=export_dir + ) + out = ( + args.out + if args.out is not None + else export_dir.parent / OPTICS_REPORT_FILENAME + ) + # Write the report first: a plotting failure must never lose the measurement. + write_report(report, out) + + plot_path: Path | None = None + plot_target = ( + args.plot + if args.plot is not None + else export_dir.parent / "refraction_comparison.png" + ) + if want_plot: + from orisys_sensor_design.optics_plot import render_comparison + + caption = [ + f"revision {paths.root.name} | {len(stack.layers)} layer(s) -> " + f"{stack.marker_surface}", + ] + for layer in result["layers"]: + thickness = layer["thickness_mm"] + caption.append( + f" {layer['name']}: ior={layer['ior']:g} " + f"({layer['ior_source']}), " + + ("thickness n/a" if thickness is None else f"{thickness:.4f} mm") + ) + try: + plot_path = render_comparison( + samples=samples, + marker=marker, + K=K, + world_to_cam=world_to_cam, + out_path=plot_target, + sanity_mm=args.sanity_mm, + caption_lines=caption, + ) + except Exception as exc: # the report is already on disk + print( + f"ERROR: could not render the comparison figure: {exc}\n" + f"The report is still available at {out}", + file=sys.stderr, + ) + return 2 + report["plot"] = str(plot_path) + write_report(report, out) + + error_mm = result["error_mm"] + cells = result["error_cells"] + print(f"Refraction through {len(stack.layers)} layer(s) -> {stack.marker_surface}") + for layer in result["layers"]: + thickness = layer["thickness_mm"] + print( + f" {layer['name']}: ior={layer['ior']:g} ({layer['ior_source']}), " + f"thickness=" + + ("n/a" if thickness is None else f"{thickness:.4f} mm") + ) + print( + f" rays: {result['resolved']}/{result['samples']} resolved " + f"({result['statuses']})" + ) + print( + " marker displacement: " + + " ".join( + f"{key}={value:.4f} mm" + for key, value in error_mm.items() + if value is not None + ) + ) + print( + " in marker cells (pitch %.4f mm): " % pitch_mm + + " ".join( + f"{key}={value:.2f}" + for key, value in cells.items() + if value is not None + ) + ) + if result["fraction_over_half_cell"] is not None: + print( + f" {100.0 * result['fraction_over_half_cell']:.1f}% of the surface is " + "displaced by more than half a marker cell" + ) + print(f"Wrote {out}") + if plot_path is not None: + from orisys_sensor_design.optics_plot import figure_summary + + counts = figure_summary(samples, args.sanity_mm) + print(f"Wrote {plot_path}") + if counts["implausible"]: + print( + f"WARNING: {counts['implausible']} sample(s) exceed " + f"{args.sanity_mm:g} mm; they are marked separately in the figure " + "and the field is not trustworthy at those rays" + ) + print( + "NOTE: advisory only; the rectifying remap belongs to orisys-sdk " + "(workflow §7)." + ) + return 0 + + def cmd_manifest(revision_dir: Path, out: Path | None, profile: str) -> int: root = repo_root(revision_dir) paths = RevisionPaths.from_revision_dir(revision_dir) @@ -754,6 +982,40 @@ def main(argv: list[str] | None = None) -> None: help="Never split normals; average every face at a vertex", ) + optics_p = sub.add_parser( + "optics", + help="Quantify refraction through the declared optical stack (advisory)", + ) + optics_p.add_argument("revision_dir", type=Path) + optics_p.add_argument("--profile", default="compute", choices=["compute"]) + optics_p.add_argument( + "--output-dir", + type=Path, + default=None, + help="Export payload directory (default: the profile's output_dir)", + ) + optics_p.add_argument( + "--samples", type=int, default=96, + help="Ray grid resolution per axis (default: 96)", + ) + optics_p.add_argument( + "--out", type=Path, default=None, + help="Report path (default: /refraction_report.json)", + ) + optics_p.add_argument( + "--plot", type=Path, default=None, + help="Comparison PNG path (default: /" + "refraction_comparison.png)", + ) + optics_p.add_argument( + "--no-plot", dest="no_plot", action="store_true", + help="Skip the side-by-side figure and write only the JSON report", + ) + optics_p.add_argument( + "--sanity-mm", dest="sanity_mm", type=float, default=5.0, + help="Displacements above this are flagged as implausible (default: 5)", + ) + args = parser.parse_args(argv) if args.command == "validate": code = cmd_validate(args.revision_dir, args.blender) @@ -765,6 +1027,8 @@ def main(argv: list[str] | None = None) -> None: code = cmd_inspect(args.revision_dir, args) elif args.command == "preview": code = cmd_preview(args.revision_dir, args) + elif args.command == "optics": + code = cmd_optics(args.revision_dir, args) elif args.command == "manifest": code = cmd_manifest(args.revision_dir, args.out, args.profile) else: diff --git a/src/orisys_sensor_design/manifest.py b/src/orisys_sensor_design/manifest.py index 9667b10..9c88ff6 100644 --- a/src/orisys_sensor_design/manifest.py +++ b/src/orisys_sensor_design/manifest.py @@ -170,9 +170,10 @@ def build_export_job( output_dir: Path, profile: str, housing_collections: list[str] | None = None, + optics_layers: list[dict[str, Any]] | None = None, ) -> dict[str, Any]: """Build the JSON job dict that Blender's export.py consumes (stdlib JSON).""" - return { + job: dict[str, Any] = { "schema": EXPORT_JOB_SCHEMA, "output": str(output_dir.resolve()), "profile": profile, @@ -183,6 +184,9 @@ def build_export_job( else export_housing_collections(paths) ), } + if optics_layers: + job["optics"] = {"layers": list(optics_layers)} + return job def export_revision_dir(root: Path, paths: RevisionPaths, profile: str) -> Path: diff --git a/src/orisys_sensor_design/optics.py b/src/orisys_sensor_design/optics.py new file mode 100644 index 0000000..0acc41b --- /dev/null +++ b/src/orisys_sensor_design/optics.py @@ -0,0 +1,794 @@ +"""Refraction through the declared optical stack, and the error it causes. + +The markers sit on the sensing surface underneath a clear rubber layer and a +plastic shell. Light from a marker is bent at every interface, so the image the +camera forms is displaced relative to the straight-line (pinhole) model the +current pipeline uses. This module declares that stack, traces rays through it +with Snell's law, and reports how large the resulting error is - in millimetres +and in marker cells. + +Scope note: ``sensor-design-workflow.md`` §7 makes ``orisys-sdk`` the authority +for rectification and unroll. This module *quantifies* the error and exports the +interfaces; it does not rectify anything. + +Modelling assumptions, stated here because they are choices, not facts: + +* each layer is a separate part with **air on both sides**, so a ray goes + air -> outer -> [n] -> inner -> air for every layer. The measured ~2 mm gap + between the shell and the rubber is therefore air, which matches the geometry; +* layers are homogeneous, non-dispersive and non-absorbing: refraction only, no + scattering, absorption or Fresnel loss; +* the ray direction is normalised after transforming into the sensor frame, so a + uniformly scaled pose (``world_to_cam`` is not rigid here) still yields the + correct geometric ray. +""" + +from __future__ import annotations + +import json +import math +from dataclasses import dataclass +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np + +from orisys_sensor_design.manifest import RevisionPaths, load_yaml + +OPTICS_FILENAME = "optics.npz" +REPORT_SCHEMA = "orisys-sensor-design/refraction-report-v1" +REPORT_FILENAME = "refraction_report.json" + +#: Index of the medium outside every layer, and the value at which a trace must +#: reproduce the pinhole result exactly. +AIR_INDEX = 1.0 + +_BLOCK_KEYS = frozenset({"marker_surface", "layers"}) +_LAYER_KEYS = frozenset( + {"name", "solid", "outer", "inner", "ior", "ior_source", "thickness_mm"} +) +_MIN_IOR = 1.0 +_MAX_IOR = 3.0 + + +class OpticsConfigError(ValueError): + """Raised when the ``optics`` metadata or the exported interfaces are invalid.""" + + +@dataclass(frozen=True) +class OpticsLayer: + """One refracting layer. + + A layer is either a **closed solid** (``solid``) whose entry and exit are the + first and second hits on the same mesh, or a pair of separate surfaces + (``outer``/``inner``). The declared stack uses solids, because ``shell`` and + ``rubber-1`` are each a single closed mesh. + """ + + name: str + ior: float + ior_source: str + solid: str | None = None + outer: str | None = None + inner: str | None = None + thickness_mm: float | None = None + + @property + def is_solid(self) -> bool: + return self.solid is not None + + @property + def entry_name(self) -> str: + return str(self.solid if self.is_solid else self.outer) + + @property + def exit_name(self) -> str: + return str(self.solid if self.is_solid else self.inner) + + +@dataclass(frozen=True) +class OpticalStack: + """The stack from the camera inward, ending at the marker surface.""" + + marker_surface: str + layers: tuple[OpticsLayer, ...] + + @property + def total_thickness_mm(self) -> float: + return float(sum(layer.thickness_mm or 0.0 for layer in self.layers)) + + +@dataclass(frozen=True) +class Interface: + """A triangulated surface in the sensor frame.""" + + name: str + vertices: np.ndarray + faces: np.ndarray + + +@dataclass(frozen=True) +class TraceResult: + """Outcome of one traced ray.""" + + status: str # "ok" | "miss" | "tir" + point: np.ndarray | None = None + incidence_deg: float | None = None + + +# --------------------------------------------------------------------- config + + +def load_optics_block(paths: RevisionPaths) -> dict[str, Any] | None: + """Return the raw ``optics`` mapping from design.yaml (None when absent).""" + design = load_yaml(paths.design_yaml) + raw = design.get("optics") + if raw is None: + return None + if not isinstance(raw, dict): + raise OpticsConfigError( + f"{paths.design_yaml}: 'optics' must be a mapping, got {type(raw).__name__}" + ) + return raw + + +def validate_optics_block(block: dict[str, Any], where: str = "optics") -> list[str]: + """Return every problem with an ``optics`` block (empty when valid).""" + errors: list[str] = [] + unknown = sorted(set(block) - _BLOCK_KEYS) + if unknown: + errors.append( + f"{where}: unknown key(s) {', '.join(unknown)} " + f"(allowed: {', '.join(sorted(_BLOCK_KEYS))})" + ) + marker = block.get("marker_surface") + if not isinstance(marker, str) or not marker.strip(): + errors.append(f"{where}.marker_surface: must be a non-empty object name") + + layers = block.get("layers") + if not isinstance(layers, list) or not layers: + errors.append(f"{where}.layers: must be a non-empty list") + return errors + + seen: set[str] = set() + for index, layer in enumerate(layers): + label = f"{where}.layers[{index}]" + if not isinstance(layer, dict): + errors.append(f"{label}: must be a mapping") + continue + extra = sorted(set(layer) - _LAYER_KEYS) + if extra: + errors.append( + f"{label}: unknown key(s) {', '.join(extra)} " + f"(allowed: {', '.join(sorted(_LAYER_KEYS))})" + ) + name = layer.get("name") + if not isinstance(name, str) or not name.strip(): + errors.append(f"{label}.name: must be a non-empty string") + elif name in seen: + errors.append(f"{label}.name: duplicate layer name {name!r}") + else: + seen.add(name) + solid = layer.get("solid") + outer = layer.get("outer") + inner = layer.get("inner") + if solid is not None: + if not isinstance(solid, str) or not solid.strip(): + errors.append(f"{label}.solid: must be a non-empty object name") + if outer or inner: + errors.append( + f"{label}: declare either 'solid' or 'outer'/'inner', not both" + ) + else: + for key, value in (("outer", outer), ("inner", inner)): + if not isinstance(value, str) or not value.strip(): + errors.append( + f"{label}.{key}: required unless 'solid' is given" + ) + ior = layer.get("ior") + if not isinstance(ior, (int, float)) or isinstance(ior, bool): + errors.append(f"{label}.ior: must be a number") + elif not _MIN_IOR < float(ior) < _MAX_IOR: + errors.append( + f"{label}.ior: must be greater than {_MIN_IOR} and less than " + f"{_MAX_IOR} (got {ior!r})" + ) + source = layer.get("ior_source") + if not isinstance(source, str) or not source.strip(): + errors.append( + f"{label}.ior_source: required - record where the index came from " + "so an assumed value cannot pass for a measured one" + ) + thickness = layer.get("thickness_mm") + if thickness is not None: + if not isinstance(thickness, (int, float)) or isinstance(thickness, bool): + errors.append(f"{label}.thickness_mm: must be a number") + elif float(thickness) <= 0.0: + errors.append(f"{label}.thickness_mm: must be positive") + + # Consecutive layers must not share an interface. This model puts air on both + # sides of every layer, so naming one surface as both the exit of one layer and + # the entry of the next produces a zero-thickness air gap and a degenerate ray + # path (observed: every ray missed after that point). + for index in range(len(layers) - 1): + first, second = layers[index], layers[index + 1] + if not isinstance(first, dict) or not isinstance(second, dict): + continue + if first.get("inner") and first.get("inner") == second.get("outer"): + errors.append( + f"{where}.layers[{index}].inner and layers[{index + 1}].outer are both " + f"{first['inner']!r}; consecutive layers must not share an interface " + "- name each layer's own two surfaces (contiguous media are not modelled)" + ) + return errors + + +def resolve_optics(paths: RevisionPaths) -> OpticalStack | None: + """Validate and build the optical stack, or None when none is declared.""" + block = load_optics_block(paths) + if block is None: + return None + errors = validate_optics_block(block, f"{paths.design_yaml} optics") + if errors: + raise OpticsConfigError("\n".join(errors)) + layers = tuple( + OpticsLayer( + name=str(layer["name"]), + ior=float(layer["ior"]), + ior_source=str(layer["ior_source"]), + solid=( + str(layer["solid"]) if layer.get("solid") is not None else None + ), + outer=( + str(layer["outer"]) if layer.get("outer") is not None else None + ), + inner=( + str(layer["inner"]) if layer.get("inner") is not None else None + ), + thickness_mm=( + None if layer.get("thickness_mm") is None + else float(layer["thickness_mm"]) + ), + ) + for layer in block["layers"] + ) + return OpticalStack(marker_surface=str(block["marker_surface"]), layers=layers) + + +# ------------------------------------------------------------------- tracing + + +def refract( + direction: np.ndarray, + normal: np.ndarray, + n_from: float, + n_to: float, +) -> np.ndarray | None: + """Snell's law at one interface, or None under total internal reflection. + + ``normal`` is reoriented to oppose the incoming ray, so callers do not have + to guarantee a winding convention. + """ + direction = np.asarray(direction, dtype=np.float64) + normal = np.asarray(normal, dtype=np.float64) + cos_i = -float(np.dot(normal, direction)) + if cos_i < 0.0: + normal = -normal + cos_i = -cos_i + eta = float(n_from) / float(n_to) + k = 1.0 - eta * eta * (1.0 - cos_i * cos_i) + if k < 0.0: + return None + refracted = eta * direction + (eta * cos_i - math.sqrt(k)) * normal + length = float(np.linalg.norm(refracted)) + if length <= 1e-15: + return None + return refracted / length + + +def intersect_mesh( + origin: np.ndarray, + direction: np.ndarray, + interface: Interface, + *, + epsilon: float = 1e-9, +) -> tuple[float, np.ndarray, np.ndarray] | None: + """Closest ray/triangle hit on ``interface``: ``(t, point, normal)``. + + Möller–Trumbore, vectorised over the triangles. The returned normal is + flipped to oppose the incoming ray. + """ + vertices = interface.vertices + faces = interface.faces + if len(faces) == 0: + return None + v0 = vertices[faces[:, 0]] + v1 = vertices[faces[:, 1]] + v2 = vertices[faces[:, 2]] + edge1 = v1 - v0 + edge2 = v2 - v0 + direction = np.asarray(direction, dtype=np.float64) + origin = np.asarray(origin, dtype=np.float64) + + pvec = np.cross(direction, edge2) + det = np.einsum("ij,ij->i", edge1, pvec) + parallel = np.abs(det) <= 1e-15 + safe = np.where(parallel, 1.0, det) + inv = 1.0 / safe + tvec = origin - v0 + u = np.einsum("ij,ij->i", tvec, pvec) * inv + qvec = np.cross(tvec, edge1) + v = (qvec * direction).sum(axis=1) * inv + t = np.einsum("ij,ij->i", edge2, qvec) * inv + + hit = ( + ~parallel + & (u >= -1e-9) + & (v >= -1e-9) + & (u + v <= 1.0 + 1e-9) + & (t > epsilon) + ) + if not np.any(hit): + return None + indices = np.nonzero(hit)[0] + best = indices[int(np.argmin(t[indices]))] + point = origin + direction * float(t[best]) + normal = np.cross(edge1[best], edge2[best]) + length = float(np.linalg.norm(normal)) + if length <= 1e-15: + return None + normal = normal / length + if float(np.dot(normal, direction)) > 0.0: + normal = -normal + return float(t[best]), point, normal + + +def surfaces_coincide( + first: np.ndarray, second: np.ndarray, tolerance: float, min_fraction: float = 0.8 +) -> bool: + """Whether two vertex clouds describe the same surface, within ``tolerance``. + + The markers are printed on the inner face of the rubber, so the marker surface + and the last layer's inner interface are the same geometry. A ray that exits + the rubber is therefore already standing on the markers, and searching for a + further intersection finds nothing (observed: 9132 of 9216 rays "missed"). + """ + if len(first) == 0 or len(second) == 0: + return False + cell = max(float(tolerance), 1e-9) + keys = set(map(tuple, np.round(second / cell).astype(np.int64))) + # A vertex counts as coincident if any of its 27 neighbouring cells is occupied, + # so the comparison is robust to points landing either side of a cell boundary. + offsets = [ + (dx, dy, dz) + for dx in (-1, 0, 1) + for dy in (-1, 0, 1) + for dz in (-1, 0, 1) + ] + matched = 0 + for vertex in np.round(first / cell).astype(np.int64): + for offset in offsets: + if tuple(vertex + offset) in keys: + matched += 1 + break + return matched / len(first) >= min_fraction + + +def trace_ray( + origin: np.ndarray, + direction: np.ndarray, + layers: list[tuple[OpticsLayer, Interface, Interface]], + marker: Interface, + *, + marker_at_last_exit: bool = False, +) -> TraceResult: + """Trace one ray from the camera through the stack onto the marker surface. + + ``marker_at_last_exit`` is set when the markers sit on the final interface, in + which case the exit point *is* the marker hit. + """ + position = np.asarray(origin, dtype=np.float64) + heading = np.asarray(direction, dtype=np.float64) + heading = heading / float(np.linalg.norm(heading)) + incidence: float | None = None + + for index, (layer, outer, inner) in enumerate(layers): + entry = intersect_mesh(position, heading, outer) + if entry is None: + return TraceResult("miss") + _, position, normal = entry + if incidence is None: + incidence = math.degrees( + math.acos(min(1.0, abs(float(np.dot(normal, heading))))) + ) + heading = refract(heading, normal, AIR_INDEX, layer.ior) + if heading is None: + return TraceResult("tir") + + if layer.is_solid: + # Step just inside and take the next hit on the *same* mesh: that is + # the far wall of the solid, which is where the ray leaves it. + step = 1e-6 * max( + 1.0, float(np.linalg.norm(outer.vertices.max(axis=0) - outer.vertices.min(axis=0))) + ) + inner_hit = intersect_mesh(position + heading * step, heading, inner) + else: + inner_hit = intersect_mesh(position, heading, inner) + if inner_hit is None: + return TraceResult("miss") + _, position, normal = inner_hit + heading = refract(heading, normal, layer.ior, AIR_INDEX) + if heading is None: + return TraceResult("tir") + + if marker_at_last_exit: + return TraceResult("ok", point=position, incidence_deg=incidence) + hit = intersect_mesh(position, heading, marker) + if hit is None: + return TraceResult("miss") + _, point, _ = hit + return TraceResult("ok", point=point, incidence_deg=incidence) + + +# -------------------------------------------------------------- interfaces + + +def _require(data: dict[str, np.ndarray], keys: tuple[str, ...], where: str) -> None: + missing = [key for key in keys if key not in data] + if missing: + raise OpticsConfigError( + f"{where} is missing {', '.join(missing)} (found: {', '.join(sorted(data))})" + ) + + +def _interface(data: dict[str, np.ndarray], prefix: str, name: str, where: str) -> Interface: + _require(data, (f"{prefix}_vertices", f"{prefix}_faces"), where) + return Interface( + name=name, + vertices=np.asarray(data[f"{prefix}_vertices"], dtype=np.float64), + faces=np.asarray(data[f"{prefix}_faces"], dtype=np.int64), + ) + + +def load_interfaces( + export_dir: Path, stack: OpticalStack +) -> tuple[ + list[tuple[OpticsLayer, Interface, Interface]], + Interface, + list[float | None], + bool, +]: + """Layer interfaces, the marker surface, measured thicknesses, and coincidence.""" + export_dir = Path(export_dir) + optics_path = export_dir / OPTICS_FILENAME + if not optics_path.is_file(): + raise OpticsConfigError( + f"{optics_path} not found; run: orisys-design export " + ) + with np.load(optics_path, allow_pickle=False) as handle: + data = {key: handle[key] for key in handle.files} + _require(data, ("layer_names",), str(optics_path)) + available = [str(name) for name in data["layer_names"]] + + layers: list[tuple[OpticsLayer, Interface, Interface]] = [] + for index, layer in enumerate(stack.layers): + if layer.name not in available: + raise OpticsConfigError( + f"{optics_path} has no layer {layer.name!r} " + f"(found: {', '.join(available)})" + ) + if layer.is_solid: + # A solid provides both interfaces: entry and exit hit one mesh. + solid = _interface(data, f"layer{index}_solid", layer.entry_name, str(optics_path)) + layers.append((layer, solid, solid)) + else: + layers.append( + ( + layer, + _interface( + data, f"layer{index}_outer", layer.entry_name, str(optics_path) + ), + _interface( + data, f"layer{index}_inner", layer.exit_name, str(optics_path) + ), + ) + ) + + mesh_path = export_dir / "mesh.npz" + if not mesh_path.is_file(): + raise OpticsConfigError(f"{mesh_path} not found; run: orisys-design export") + with np.load(mesh_path, allow_pickle=False) as handle: + mesh = {key: handle[key] for key in handle.files} + marker = Interface( + name=stack.marker_surface, + vertices=np.asarray(mesh["vertices"], dtype=np.float64), + faces=np.asarray(mesh["faces"], dtype=np.int64), + ) + measured: list[float | None] = [] + thicknesses = data.get("thickness_mm") + for index in range(len(stack.layers)): + if thicknesses is not None and index < len(thicknesses): + value = float(thicknesses[index]) + measured.append(value if np.isfinite(value) else None) + else: + measured.append(None) + # Do the markers sit on the final interface? Then the exit point is the hit. + marker_at_last_exit = False + if layers: + last_inner = layers[-1][2] + # 0.1 mm expressed in the export frame's own units, taken from the payload + # rather than assumed, since the sensor frame is not metres. + span = float(np.linalg.norm(marker.vertices.max(axis=0) - marker.vertices.min(axis=0))) + tolerance = span * 1e-3 if span > 0 else 1e-6 + marker_at_last_exit = surfaces_coincide( + # Ask whether the *markers* lie on the inner surface, not whether the + # two clouds match each other: the inner shell has roughly twice as + # many vertices as the marker sheet, so the symmetric test fails. + marker.vertices, + last_inner.vertices, + tolerance, + ) + return layers, marker, measured, marker_at_last_exit + + +# -------------------------------------------------------------- analysis + + +def _percentile(values: list[float], fraction: float) -> float | None: + if not values: + return None + ordered = sorted(values) + if len(ordered) == 1: + return ordered[0] + position = (len(ordered) - 1) * fraction + low = math.floor(position) + high = math.ceil(position) + if low == high: + return ordered[int(position)] + return ordered[low] + (ordered[high] - ordered[low]) * (position - low) + + +def camera_rays( + K: np.ndarray, + world_to_cam: np.ndarray, + width: int, + height: int, + samples: int, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Camera origin, a ``samples x samples`` grid of unit ray directions, and pixels. + + Directions are built in camera space (Blender looks down ``-Z``) and rotated + into the sensor frame; normalising absorbs the uniform scale in the pose. The + pixel coordinates are returned so a caller can plot the field in image space. + """ + sensor_from_cam = np.linalg.inv(np.asarray(world_to_cam, dtype=np.float64)) + origin = sensor_from_cam @ np.array([0.0, 0.0, 0.0, 1.0]) + origin = origin[:3] + rotation = sensor_from_cam[:3, :3] + + us = (np.arange(samples) + 0.5) * (width / samples) + vs = (np.arange(samples) + 0.5) * (height / samples) + grid_u, grid_v = np.meshgrid(us, vs) + x = (grid_u - float(K[0, 2])) / float(K[0, 0]) + y = (grid_v - float(K[1, 2])) / float(K[1, 1]) + directions_cam = np.stack((x, y, -np.ones_like(x)), axis=-1).reshape(-1, 3) + directions = directions_cam @ rotation.T + norms = np.linalg.norm(directions, axis=1, keepdims=True) + directions = directions / np.where(norms > 0.0, norms, 1.0) + return origin, directions, grid_u.reshape(-1), grid_v.reshape(-1) + + +@dataclass(frozen=True) +class SampleField: + """Per-ray results over the pixel grid, for plotting and for summarising. + + Arrays are flat and row-major over a ``grid_shape`` grid; ``straight``, + ``refracted`` and ``error_mm`` hold ``NaN`` wherever a ray did not resolve, so + a caller can never mistake "no data" for "no error". + """ + + pixel_u: np.ndarray + pixel_v: np.ndarray + status: list[str] + straight: np.ndarray + refracted: np.ndarray + error_mm: np.ndarray + incidence_deg: np.ndarray + grid_shape: tuple[int, int] + image_size: tuple[int, int] + pitch_mm: float + layers: list[dict[str, Any]] + total_thickness_mm: float + + @property + def count(self) -> int: + return int(len(self.pixel_u)) + + @property + def resolved(self) -> int: + return int(np.count_nonzero(~np.isnan(self.error_mm))) + + def status_counts(self) -> dict[str, int]: + counts: dict[str, int] = {} + for status in self.status: + counts[status] = counts.get(status, 0) + 1 + return counts + + +def _percentile_of(values: np.ndarray, fraction: float) -> float | None: + """Percentile over the finite entries of an array, or None when there are none.""" + finite = values[np.isfinite(values)] + if finite.size == 0: + return None + return _percentile([float(value) for value in finite], fraction) + + +def sample_stack( + *, + stack: OpticalStack, + layers: list[tuple[OpticsLayer, Interface, Interface]], + marker: Interface, + origin: np.ndarray, + directions: np.ndarray, + pitch_mm: float, + measured_thickness: list[float | None] | None = None, + marker_at_last_exit: bool = False, + pixel_u: np.ndarray | None = None, + pixel_v: np.ndarray | None = None, + image_size: tuple[int, int] = (0, 0), +) -> SampleField: + """Trace every ray twice and keep the per-ray comparison.""" + count = len(directions) + straight_points = np.full((count, 3), np.nan, dtype=np.float64) + refracted_points = np.full((count, 3), np.nan, dtype=np.float64) + errors = np.full(count, np.nan, dtype=np.float64) + angles = np.full(count, np.nan, dtype=np.float64) + statuses: list[str] = [] + + for index, direction in enumerate(directions): + straight = intersect_mesh(origin, direction, marker) + result = trace_ray( + origin, + direction, + layers, + marker, + marker_at_last_exit=marker_at_last_exit, + ) + statuses.append(result.status) + if straight is None or result.status != "ok" or result.point is None: + continue + _, straight_point, _ = straight + straight_points[index] = straight_point + refracted_points[index] = result.point + errors[index] = float(np.linalg.norm(result.point - straight_point)) + if result.incidence_deg is not None: + angles[index] = result.incidence_deg + + resolved: list[float | None] = [] + for index, (layer, _, _) in enumerate(layers): + if layer.thickness_mm is not None: + resolved.append(layer.thickness_mm) + elif measured_thickness and index < len(measured_thickness): + resolved.append(measured_thickness[index]) + else: + resolved.append(None) + + if pixel_u is None: + pixel_u = np.arange(count, dtype=np.float64) + if pixel_v is None: + pixel_v = np.zeros(count, dtype=np.float64) + side = int(round(math.sqrt(count))) if count else 0 + grid_shape = (side, side) if side * side == count else (1, count) + + return SampleField( + pixel_u=np.asarray(pixel_u, dtype=np.float64), + pixel_v=np.asarray(pixel_v, dtype=np.float64), + status=statuses, + straight=straight_points, + refracted=refracted_points, + error_mm=errors, + incidence_deg=angles, + grid_shape=grid_shape, + image_size=(int(image_size[0]), int(image_size[1])), + pitch_mm=float(pitch_mm), + layers=[ + { + "name": layer.name, + "ior": layer.ior, + "ior_source": layer.ior_source, + "thickness_mm": resolved[index], + } + for index, (layer, _, _) in enumerate(layers) + ], + total_thickness_mm=float(sum(v for v in resolved if v is not None)), + ) + + +def summarize(samples: SampleField) -> dict[str, Any]: + """The aggregate report, computed from a :class:`SampleField`.""" + finite = np.isfinite(samples.error_mm) + errors_mm = samples.error_mm[finite] + pitch = samples.pitch_mm + cells = errors_mm / pitch if pitch > 0 else np.zeros_like(errors_mm) + + def summary(values: np.ndarray) -> dict[str, float | None]: + return { + "min": _percentile_of(values, 0.0), + "p05": _percentile_of(values, 0.05), + "p50": _percentile_of(values, 0.50), + "p95": _percentile_of(values, 0.95), + "max": _percentile_of(values, 1.0), + } + + return { + "samples": samples.count, + "resolved": samples.resolved, + "statuses": samples.status_counts(), + "pitch_mm": pitch, + "error_mm": summary(errors_mm), + "error_cells": summary(cells), + "incidence_deg": summary(samples.incidence_deg), + "fraction_over_half_cell": ( + float(np.count_nonzero(cells > 0.5) / cells.size) if cells.size else None + ), + "layers": samples.layers, + "total_thickness_mm": samples.total_thickness_mm, + } + + +def quantify( + *, + stack: OpticalStack, + layers: list[tuple[OpticsLayer, Interface, Interface]], + marker: Interface, + origin: np.ndarray, + directions: np.ndarray, + pitch_mm: float, + measured_thickness: list[float | None] | None = None, + marker_at_last_exit: bool = False, +) -> dict[str, Any]: + """Compare the refracted hit against the straight-line (pinhole) hit.""" + return summarize( + sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=pitch_mm, + measured_thickness=measured_thickness, + marker_at_last_exit=marker_at_last_exit, + ) + ) + + +def build_report( + *, + paths: RevisionPaths, + stack: OpticalStack, + result: dict[str, Any], + export_dir: Path, + blend_sha256: str | None = None, +) -> dict[str, Any]: + return { + "schema": REPORT_SCHEMA, + "generated_at": datetime.now(timezone.utc).replace(microsecond=0).isoformat(), + "revision": str(paths.root), + "export_dir": str(Path(export_dir)), + "marker_surface": stack.marker_surface, + "blend_sha256": blend_sha256, + "assumptions": [ + "each layer is modelled as a slab with air on both sides", + "layers are homogeneous, non-dispersive and non-absorbing", + "refraction only: no scattering, absorption or Fresnel loss", + ], + "result": result, + } + + +def write_report(report: dict[str, Any], path: Path) -> Path: + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8") + return path diff --git a/src/orisys_sensor_design/optics_plot.py b/src/orisys_sensor_design/optics_plot.py new file mode 100644 index 0000000..c5ab913 --- /dev/null +++ b/src/orisys_sensor_design/optics_plot.py @@ -0,0 +1,228 @@ +"""Two-panel refraction figure: marker grid beside the displacement field. + +Both panels share the camera image frame, so the left panel's marker grid and the +right panel's displacement colour map line up pixel for pixel. Plotting the left +panel on the surface instead would make the two impossible to compare by eye. + +The figure is deliberately honest about bad data: unresolved rays are drawn as a +mask and counted, and displacements beyond a sanity threshold get their own colour +and their own count. A broken trace must not be able to masquerade as a large but +plausible field. + +matplotlib is imported lazily and only by this module, so the core optics code and +the JSON report keep working without it. The ``Agg`` backend is selected before +pyplot is imported, so this works headlessly. +""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import numpy as np + +from orisys_sensor_design.optics import Interface, SampleField + +#: Displacements above this are not credible for a millimetre-scale stack; they +#: are drawn separately rather than dominating the colour scale. +DEFAULT_SANITY_MM = 5.0 + +_GRID_COLOR = "#3b6ea5" +_UNRESOLVED_COLOR = "#8a8a8a" +_IMPLAUSIBLE_COLOR = "#d62728" + + +def project_pinhole( + points: np.ndarray, K: np.ndarray, world_to_cam: np.ndarray +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Project sensor-frame points to pixels with the straight pinhole model. + + Returns ``(u, v, depth)``; ``depth`` is ``-z`` in camera space, matching the + SDK's ``project_world_to_image``. + """ + points = np.asarray(points, dtype=np.float64).reshape(-1, 3) + homogeneous = np.column_stack((points, np.ones(len(points), dtype=np.float64))) + camera = (np.asarray(world_to_cam, dtype=np.float64) @ homogeneous.T).T[:, :3] + depth = -camera[:, 2] + safe = np.where(np.abs(depth) > 1e-12, depth, np.nan) + u = float(K[0, 0]) * (camera[:, 0] / safe) + float(K[0, 2]) + v = float(K[1, 1]) * (camera[:, 1] / safe) + float(K[1, 2]) + return u, v, depth + + +def _grid_cells(samples: SampleField) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Reshape the field onto the sample grid for ``pcolormesh``. + + Returns ``(u_centres, v_centres, error_grid)``. + """ + rows, columns = samples.grid_shape + u = samples.pixel_u.reshape(rows, columns) + v = samples.pixel_v.reshape(rows, columns) + error = samples.error_mm.reshape(rows, columns) + return u, v, error + + +def render_comparison( + *, + samples: SampleField, + marker: Interface, + K: np.ndarray, + world_to_cam: np.ndarray, + out_path: Path, + title: str = "Refraction through the optical stack", + dpi: int = 150, + sanity_mm: float = DEFAULT_SANITY_MM, + caption_lines: list[str] | None = None, +) -> Path: + """Write the two-panel PNG and return its path.""" + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + out_path = Path(out_path) + out_path.parent.mkdir(parents=True, exist_ok=True) + + rows, columns = samples.grid_shape + u = samples.pixel_u.reshape(rows, columns) + v = samples.pixel_v.reshape(rows, columns) + error = samples.error_mm.reshape(rows, columns) + width, height = samples.image_size + if width <= 0 or height <= 0: + width = int(np.nanmax(u) * 2) or 1 + height = int(np.nanmax(v) * 2) or 1 + + finite = error[np.isfinite(error)] + implausible = int(np.count_nonzero(finite > sanity_mm)) + plausible = finite[finite <= sanity_mm] + vmax = float(plausible.max()) if plausible.size else 0.0 + # An all-zero field (an air-only stack) would collapse the colour range. + vmax = max(vmax, 1e-6) + + figure, axes = plt.subplots(1, 2, figsize=(13.0, 6.0), dpi=dpi) + figure.suptitle(title, fontsize=13) + + # ---------------------------------------------------------------- left panel + left = axes[0] + marker_u, marker_v, depth = project_pinhole(marker.vertices, K, world_to_cam) + visible = depth > 0 + left.scatter( + marker_u[visible], + marker_v[visible], + s=1.0, + c=_GRID_COLOR, + alpha=0.35, + linewidths=0, + label=f"surface samples ({int(np.count_nonzero(visible))})", + ) + unresolved = ~np.isfinite(error) + left.scatter( + u[unresolved], + v[unresolved], + s=4.0, + marker="s", + c=_UNRESOLVED_COLOR, + alpha=0.35, + linewidths=0, + label=f"unresolved rays ({int(np.count_nonzero(unresolved))})", + ) + left.set_title("Left: marker grid (pinhole model)") + left.legend(loc="upper right", fontsize=7, framealpha=0.85) + + # --------------------------------------------------------------- right panel + right = axes[1] + safe = np.where(np.isfinite(error) & (error <= sanity_mm), error, np.nan) + mesh = right.pcolormesh( + u, + v, + np.ma.masked_invalid(safe), + shading="nearest", + cmap="viridis", + vmin=0.0, + vmax=vmax, + ) + beyond = np.isfinite(error) & (error > sanity_mm) + if np.any(beyond): + right.scatter( + u[beyond], + v[beyond], + s=8.0, + marker="x", + c=_IMPLAUSIBLE_COLOR, + linewidths=0.8, + label=f"implausible > {sanity_mm:g} mm ({implausible})", + ) + right.legend(loc="upper right", fontsize=7, framealpha=0.85) + if not np.any(np.isfinite(error)): + right.text( + 0.5, + 0.5, + "no resolved rays", + transform=right.transAxes, + ha="center", + va="center", + fontsize=14, + color=_IMPLAUSIBLE_COLOR, + ) + right.set_title("Right: displacement magnitude (refracted vs pinhole)") + + bar = figure.colorbar(mesh, ax=right, fraction=0.046, pad=0.03) + bar.set_label("displacement (mm)") + ticks = list(bar.get_ticks()) + if samples.pitch_mm > 0: + # Pin the locator first; set_ticklabels alone is only valid on a fixed + # locator, and warns otherwise. + bar.set_ticks(ticks) + bar.set_ticklabels( + [f"{value:.3g}\n{value / samples.pitch_mm:.2f} cell" for value in ticks] + ) + + for axis in axes: + axis.set_xlim(0, width) + axis.set_ylim(height, 0) # image rows run downwards + axis.set_aspect("equal", adjustable="box") + axis.set_xlabel("image u (px)") + axes[0].set_ylabel("image v (px)") + + error_mm = samples.error_mm[np.isfinite(samples.error_mm)] + statuses = samples.status_counts() + lines = list(caption_lines or []) + lines.append( + f"pitch {samples.pitch_mm:.4f} mm | ray grid {rows}x{columns} | " + f"resolved {samples.resolved}/{samples.count} | " + + ", ".join(f"{key} {value}" for key, value in sorted(statuses.items())) + ) + if error_mm.size: + summary = "displacement mm: " + " ".join( + f"{label}={np.percentile(error_mm, q):.4g}" + for label, q in (("p05", 5), ("p50", 50), ("p95", 95)) + ) + summary += f" max={error_mm.max():.4g}" + if samples.pitch_mm > 0: + summary += ( + f" | cells p50={np.percentile(error_mm, 50) / samples.pitch_mm:.2f}" + ) + lines.append(summary) + if implausible: + lines.append( + f"WARNING: {implausible} sample(s) exceed {sanity_mm:g} mm and are " + "plotted separately - the trace is not trustworthy at those rays" + ) + figure.text(0.02, 0.015, "\n".join(line for line in lines if line), fontsize=7, va="bottom") + + figure.tight_layout(rect=(0, 0.07, 1, 0.95)) + figure.savefig(out_path) + plt.close(figure) + return out_path + + +def figure_summary(samples: SampleField, sanity_mm: float = DEFAULT_SANITY_MM) -> dict[str, Any]: + """Counts the caller prints alongside the figure.""" + finite = samples.error_mm[np.isfinite(samples.error_mm)] + return { + "resolved": samples.resolved, + "samples": samples.count, + "statuses": samples.status_counts(), + "implausible": int(np.count_nonzero(finite > sanity_mm)), + "sanity_mm": sanity_mm, + } diff --git a/tests/test_optics.py b/tests/test_optics.py new file mode 100644 index 0000000..59cc825 --- /dev/null +++ b/tests/test_optics.py @@ -0,0 +1,308 @@ +"""Tests for the optical stack: validation, Snell's law, and the slab oracle.""" + +from __future__ import annotations + +import math +from pathlib import Path + +import pytest + +np = pytest.importorskip("numpy") + +from orisys_sensor_design.optics import ( + AIR_INDEX, + Interface, + OpticsConfigError, + OpticsLayer, + OpticalStack, + intersect_mesh, + load_optics_block, + quantify, + refract, + resolve_optics, + surfaces_coincide, + trace_ray, + validate_optics_block, +) +from orisys_sensor_design.manifest import RevisionPaths, repo_root + +REPO = repo_root() + + +def _block(**overrides) -> dict: + base = { + "marker_surface": "sensing_surface", + "layers": [ + { + "name": "clear_rubber", + "outer": "rubber-2", + "inner": "rubber-1", + "ior": 1.43, + "ior_source": "assumed_silicone_literature", + } + ], + } + base.update(overrides) + return base + + +def _quad(z: float, half: float = 10.0, name: str = "q") -> Interface: + """A horizontal quad facing -Z, as two triangles.""" + vertices = np.array( + [ + [-half, -half, z], + [half, -half, z], + [half, half, z], + [-half, half, z], + ], + dtype=np.float64, + ) + faces = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64) + return Interface(name=name, vertices=vertices, faces=faces) + + +# ------------------------------------------------------------- validation + + +def test_validate_accepts_a_complete_block() -> None: + assert validate_optics_block(_block()) == [] + + +def test_validate_rejects_unknown_keys() -> None: + errors = validate_optics_block(_block(dpi=600)) + assert any("unknown key" in message for message in errors) + + +def test_validate_requires_ior_source() -> None: + block = _block() + del block["layers"][0]["ior_source"] + errors = validate_optics_block(block) + assert any("ior_source" in message for message in errors) + + +def test_validate_rejects_an_index_at_or_below_one() -> None: + for value in (1.0, 0.9): + block = _block() + block["layers"][0]["ior"] = value + errors = validate_optics_block(block) + assert any("ior" in message for message in errors), value + + +def test_validate_rejects_duplicate_layer_names() -> None: + block = _block() + block["layers"].append(dict(block["layers"][0])) + errors = validate_optics_block(block) + assert any("duplicate" in message for message in errors) + + +def test_validate_rejects_empty_layers() -> None: + assert any("non-empty" in m for m in validate_optics_block(_block(layers=[]))) + + +def test_validate_rejects_a_shared_interface_between_layers() -> None: + """The exact mistake that made every ray miss: rubber-2 as exit and entry.""" + block = _block( + layers=[ + { + "name": "plastic_shell", + "outer": "shell", + "inner": "rubber-2", + "ior": 1.59, + "ior_source": "assumed", + }, + { + "name": "clear_rubber", + "outer": "rubber-2", + "inner": "rubber-1", + "ior": 1.43, + "ior_source": "assumed", + }, + ] + ) + errors = validate_optics_block(block) + assert any("share an interface" in message for message in errors) + + +def test_resolve_optics_returns_none_when_absent(tmp_path: Path) -> None: + revision = tmp_path / "sensors" / "finger" / "r99" + revision.mkdir(parents=True) + (revision / "design.yaml").write_text("sensor: finger\nrevision: r99\n", encoding="utf-8") + paths = RevisionPaths.from_revision_dir(revision) + assert load_optics_block(paths) is None + assert resolve_optics(paths) is None + + +def test_r01_design_yaml_declares_a_valid_stack() -> None: + paths = RevisionPaths.from_revision_dir(REPO / "sensors" / "finger-xynova" / "r01") + stack = resolve_optics(paths) + assert stack is not None + assert stack.marker_surface == "sensing_surface" + # The stack is camera > shell > rubber-1 > sensing_surface > rubber-2, so the + # two refracting solids are shell and rubber-1; rubber-2 sits behind the + # markers and does not bend camera light. + assert [layer.name for layer in stack.layers] == ["plastic_shell", "clear_rubber"] + assert [layer.solid for layer in stack.layers] == ["shell", "rubber-1"] + assert all(layer.is_solid for layer in stack.layers) + assert all(layer.ior_source for layer in stack.layers), ( + "an assumed index must record its source" + ) + + +# ------------------------------------------------------------------- optics + + +def test_refract_obeys_snells_law() -> None: + normal = np.array([0.0, 0.0, 1.0]) + angle_in = math.radians(40.0) + direction = np.array([math.sin(angle_in), 0.0, -math.cos(angle_in)]) + out = refract(direction, normal, 1.0, 1.5) + assert out is not None + angle_out = math.acos(abs(float(out[2]))) + assert math.sin(angle_in) == pytest.approx(1.5 * math.sin(angle_out), rel=1e-9) + + +def test_refract_returns_none_under_total_internal_reflection() -> None: + normal = np.array([0.0, 0.0, 1.0]) + steep = math.radians(70.0) + direction = np.array([math.sin(steep), 0.0, math.cos(steep)]) + assert refract(direction, normal, 1.5, 1.0) is None + + +def test_intersect_mesh_finds_the_plane_and_reports_a_miss() -> None: + plane = _quad(0.0) + hit = intersect_mesh(np.array([0.0, 0.0, -5.0]), np.array([0.0, 0.0, 1.0]), plane) + assert hit is not None + distance, point, normal = hit + assert distance == pytest.approx(5.0) + assert point[2] == pytest.approx(0.0) + assert abs(float(normal[2])) == pytest.approx(1.0) + miss = intersect_mesh(np.array([0.0, 0.0, -5.0]), np.array([0.0, 1.0, 0.0]), plane) + assert miss is None + + +def test_surfaces_coincide_is_asymmetric_by_design() -> None: + dense = np.array([[x, 0.0, 0.0] for x in np.linspace(0.0, 1.0, 20)]) + sparse = dense[::4] + # Every sparse point lies on the dense set, but not the reverse. + assert surfaces_coincide(sparse, dense, 1e-6) is True + assert surfaces_coincide(dense, sparse, 1e-6) is False + + +def test_slab_displacement_matches_the_closed_form() -> None: + """The tracer against an independent analytic oracle, not against itself. + + A plane-parallel slab of index ``n`` and thickness ``t`` shifts a ray of + incidence ``theta1`` *along the interface* by ``t * (tan(theta1) - tan(theta2))`` + with ``sin(theta2) = sin(theta1)/n``. That in-plane shift is the one that moves + a marker's apparent position; the perpendicular gap between the incident and + exit rays is smaller by a factor of ``cos(theta1)``, and conflating the two + understates the error by 15% at 30 degrees. + """ + thickness = 2.0 + index = 1.5 + layer = OpticsLayer( + name="slab", outer="top", inner="bottom", ior=index, ior_source="test" + ) + layers = [(layer, _quad(0.0, name="top"), _quad(-thickness, name="bottom"))] + # Wide enough that a 40 degree ray still lands on it after 15 units of travel. + marker = _quad(-10.0, half=60.0, name="marker") + + for theta_deg in (10.0, 20.0, 30.0, 40.0): + theta = math.radians(theta_deg) + direction = np.array([math.sin(theta), 0.0, -math.cos(theta)]) + origin = np.array([0.0, 0.0, 5.0]) + + straight = intersect_mesh(origin, direction, marker) + traced = trace_ray(origin, direction, layers, marker) + assert straight is not None and traced.status == "ok", theta_deg + + _, straight_point, _ = straight + measured = float(np.linalg.norm(traced.point - straight_point)) + theta_out = math.asin(math.sin(theta) / index) + expected = thickness * (math.tan(theta) - math.tan(theta_out)) + assert measured == pytest.approx(expected, rel=1e-3), theta_deg + + +def test_index_one_reproduces_the_pinhole_result() -> None: + """With no refraction the trace must agree with the straight ray exactly.""" + layer = OpticsLayer(name="air", outer="top", inner="bottom", ior=1.0, ior_source="test") + layers = [(layer, _quad(0.0, name="top"), _quad(-2.0, name="bottom"))] + marker = _quad(-10.0, name="marker") + origin = np.array([0.0, 0.0, 5.0]) + direction = np.array([0.3, 0.1, -1.0]) + direction = direction / np.linalg.norm(direction) + + straight = intersect_mesh(origin, direction, marker) + traced = trace_ray(origin, direction, layers, marker) + assert straight is not None and traced.status == "ok" + _, straight_point, _ = straight + assert np.allclose(traced.point, straight_point, atol=1e-9) + + +def test_marker_at_last_exit_uses_the_exit_point() -> None: + """Markers printed on the inner face: the exit point is the marker hit.""" + layer = OpticsLayer(name="rubber", outer="top", inner="bottom", ior=1.43, ior_source="t") + layers = [(layer, _quad(0.0, name="top"), _quad(-1.0, name="bottom"))] + marker = _quad(-1.0, name="marker") + origin = np.array([0.0, 0.0, 5.0]) + direction = np.array([0.1, 0.0, -1.0]) + direction = direction / np.linalg.norm(direction) + + without = trace_ray(origin, direction, layers, marker) + with_flag = trace_ray(origin, direction, layers, marker, marker_at_last_exit=True) + assert with_flag.status == "ok" + assert with_flag.point[2] == pytest.approx(-1.0, abs=1e-6) + # Without the flag the ray sits on the marker already and finds nothing. + assert without.status == "miss" + + +def test_quantify_reports_zero_error_for_an_air_only_stack() -> None: + layer = OpticsLayer(name="air", outer="top", inner="bottom", ior=1.0, ior_source="t") + layers = [(layer, _quad(0.0, name="top"), _quad(-2.0, name="bottom"))] + marker = _quad(-10.0, name="marker") + stack = OpticalStack(marker_surface="marker", layers=(layer,)) + + directions = np.array([[x, 0.0, -1.0] for x in (-0.2, 0.0, 0.2)]) + directions = directions / np.linalg.norm(directions, axis=1, keepdims=True) + result = quantify( + stack=stack, + layers=layers, + marker=marker, + origin=np.array([0.0, 0.0, 5.0]), + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + ) + assert result["statuses"]["ok"] == 3 + assert result["error_mm"]["max"] == pytest.approx(0.0, abs=1e-9) + assert result["error_cells"]["max"] == pytest.approx(0.0, abs=1e-9) + assert result["layers"][0]["thickness_mm"] == pytest.approx(2.0) + + +def test_quantify_reports_measured_thickness_when_not_declared() -> None: + layer = OpticsLayer(name="slab", outer="top", inner="bottom", ior=1.5, ior_source="t") + layers = [(layer, _quad(0.0, name="top"), _quad(-2.0, name="bottom"))] + marker = _quad(-10.0, name="marker") + stack = OpticalStack(marker_surface="marker", layers=(layer,)) + result = quantify( + stack=stack, + layers=layers, + marker=marker, + origin=np.array([0.0, 0.0, 5.0]), + directions=np.array([[0.0, 0.0, -1.0]]), + pitch_mm=0.25, + measured_thickness=[1.75], + ) + assert result["layers"][0]["thickness_mm"] == pytest.approx(1.75) + assert result["total_thickness_mm"] == pytest.approx(1.75) + + +def test_load_optics_block_rejects_a_non_mapping(tmp_path: Path) -> None: + revision = tmp_path / "sensors" / "finger" / "r99" + revision.mkdir(parents=True) + (revision / "design.yaml").write_text( + "sensor: finger\nrevision: r99\noptics: 3\n", encoding="utf-8" + ) + paths = RevisionPaths.from_revision_dir(revision) + with pytest.raises(OpticsConfigError, match="must be a mapping"): + load_optics_block(paths) diff --git a/tests/test_optics_plot.py b/tests/test_optics_plot.py new file mode 100644 index 0000000..2487b8a --- /dev/null +++ b/tests/test_optics_plot.py @@ -0,0 +1,318 @@ +"""Tests for the two-panel refraction comparison figure and the sample field.""" + +from __future__ import annotations + +import math +import sys +from pathlib import Path + +import pytest + +np = pytest.importorskip("numpy") + +from orisys_sensor_design.optics import ( + AIR_INDEX, + Interface, + OpticsLayer, + OpticalStack, + camera_rays, + sample_stack, + summarize, + quantify, +) +from orisys_sensor_design.optics_plot import ( + DEFAULT_SANITY_MM, + figure_summary, + project_pinhole, + render_comparison, +) + +REPO = Path(__file__).resolve().parents[1] +XYNOVA_EXPORT = REPO / "build" / "export" / "finger-xynova" / "r01" / "compute" + + +def _quad(z: float, half: float = 10.0, name: str = "q") -> Interface: + vertices = np.array( + [[-half, -half, z], [half, -half, z], [half, half, z], [-half, half, z]], + dtype=np.float64, + ) + faces = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64) + return Interface(name=name, vertices=vertices, faces=faces) + + +def _air_stack(): + layer = OpticsLayer(name="air", outer="top", inner="bottom", ior=1.0, ior_source="t") + layers = [(layer, _quad(0.0, name="top"), _quad(-2.0, name="bottom"))] + marker = _quad(-10.0, half=60.0, name="marker") + return OpticalStack(marker_surface="marker", layers=(layer,)), layers, marker + + +def _rays(samples: int = 4) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + K = np.array([[100.0, 0.0, 50.0], [0.0, 100.0, 50.0], [0.0, 0.0, 1.0]]) + # Place the camera 5 units in front of the stack: at the sensor origin it + # would sit exactly on the first interface and every ray would miss. + world_to_cam = np.eye(4) + world_to_cam[2, 3] = -5.0 + return camera_rays(K, world_to_cam, 100, 100, samples) + + +# ----------------------------------------------------------- sample field + + +def test_sample_stack_returns_arrays_with_nan_for_unresolved() -> None: + stack, layers, marker = _air_stack() + origin, directions, pixel_u, pixel_v = _rays(4) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + pixel_u=pixel_u, + pixel_v=pixel_v, + image_size=(100, 100), + ) + assert field.count == 16 + assert field.error_mm.shape == (16,) + assert field.straight.shape == (16, 3) + # A plane-facing camera resolves every ray here. + assert field.resolved == 16 + assert np.all(np.isfinite(field.error_mm)) + assert field.grid_shape == (4, 4) + assert field.image_size == (100, 100) + assert np.allclose(field.error_mm, 0.0, atol=1e-9) + + +def test_sample_stack_marks_misses_as_nan() -> None: + stack, layers, marker = _air_stack() + # Rays pointing away from every surface must not become zeros. + directions = np.tile(np.array([0.0, 0.0, 1.0]), (5, 1)) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=np.array([0.0, 0.0, 5.0]), + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + ) + assert field.resolved == 0 + assert np.all(np.isnan(field.error_mm)) + assert field.status_counts().get("miss") == 5 + + +def test_summarize_agrees_with_quantify() -> None: + """The refactor must not change the JSON report's numbers.""" + stack, layers, marker = _air_stack() + origin, directions, pixel_u, pixel_v = _rays(4) + kwargs = dict( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + ) + summary = summarize(sample_stack(**kwargs)) + legacy = quantify(**kwargs) + assert summary == legacy + assert summary["statuses"]["ok"] == 16 + + +def test_camera_rays_pixel_grid_matches_the_direction_count() -> None: + origin, directions, pixel_u, pixel_v = _rays(5) + assert directions.shape == (25, 3) + assert pixel_u.shape == (25,) + assert pixel_v.shape == (25,) + assert pixel_u.min() > 0.0 and pixel_u.max() < 100.0 + assert origin.shape == (3,) + + +# ----------------------------------------------------------------- plotting + + +def test_project_pinhole_places_a_centred_point_at_the_principal_point() -> None: + K = np.array([[100.0, 0.0, 50.0], [0.0, 100.0, 50.0], [0.0, 0.0, 1.0]]) + u, v, depth = project_pinhole(np.array([[0.0, 0.0, -5.0]]), K, np.eye(4)) + assert u[0] == pytest.approx(50.0) + assert v[0] == pytest.approx(50.0) + assert depth[0] == pytest.approx(5.0) + + +matplotlib = pytest.importorskip("matplotlib") + + +def _render(tmp_path: Path, field, marker, *, sanity_mm: float = DEFAULT_SANITY_MM): + K = np.array([[100.0, 0.0, 50.0], [0.0, 100.0, 50.0], [0.0, 0.0, 1.0]]) + return render_comparison( + samples=field, + marker=marker, + K=K, + world_to_cam=np.eye(4), + out_path=tmp_path / "figure.png", + sanity_mm=sanity_mm, + ) + + +def test_render_comparison_writes_a_png(tmp_path: Path) -> None: + stack, layers, marker = _air_stack() + origin, directions, pixel_u, pixel_v = _rays(4) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + pixel_u=pixel_u, + pixel_v=pixel_v, + image_size=(100, 100), + ) + path = _render(tmp_path, field, marker) + assert path.is_file() and path.stat().st_size > 0 + assert path.read_bytes()[:8] == b"\x89PNG\r\n\x1a\n" + + +def test_render_comparison_survives_an_all_unresolved_field(tmp_path: Path) -> None: + stack, layers, marker = _air_stack() + directions = np.tile(np.array([0.0, 0.0, 1.0]), (4, 1)) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=np.array([0.0, 0.0, 5.0]), + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + pixel_u=np.array([1.0, 2.0, 3.0, 4.0]), + pixel_v=np.array([1.0, 2.0, 3.0, 4.0]), + image_size=(100, 100), + ) + path = _render(tmp_path, field, marker) + assert path.is_file() and path.stat().st_size > 0 + + +def test_render_comparison_handles_an_all_zero_field(tmp_path: Path) -> None: + """An air-only stack has zero displacement; the colour range must not collapse.""" + stack, layers, marker = _air_stack() + origin, directions, pixel_u, pixel_v = _rays(4) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + pixel_u=pixel_u, + pixel_v=pixel_v, + image_size=(100, 100), + ) + assert np.allclose(field.error_mm, 0.0, atol=1e-9) + assert _render(tmp_path, field, marker).is_file() + + +def test_figure_summary_counts_implausible_samples() -> None: + stack, layers, marker = _air_stack() + origin, directions, pixel_u, pixel_v = _rays(4) + field = sample_stack( + stack=stack, + layers=layers, + marker=marker, + origin=origin, + directions=directions, + pitch_mm=0.25, + measured_thickness=[2.0], + pixel_u=pixel_u, + pixel_v=pixel_v, + image_size=(100, 100), + ) + counts = figure_summary(field, sanity_mm=5.0) + assert counts["resolved"] == 16 + assert counts["implausible"] == 0 # air only: no displacement at all + assert counts["sanity_mm"] == 5.0 + + +# ---------------------------------------------------------------------- CLI + + +@pytest.mark.skipif( + not (XYNOVA_EXPORT / "mesh.npz").is_file(), reason="export not built" +) +def test_cli_writes_both_the_report_and_the_figure(tmp_path: Path) -> None: + from orisys_sensor_design import cli + + revision = REPO / "sensors" / "finger-xynova" / "r01" + args = type( + "Args", + (), + { + "profile": "compute", + "output_dir": XYNOVA_EXPORT, + "samples": 8, + "out": tmp_path / "report.json", + "plot": tmp_path / "figure.png", + "no_plot": False, + "sanity_mm": 5.0, + }, + )() + assert cli.cmd_optics(revision, args) == 0 + assert (tmp_path / "report.json").is_file() + assert (tmp_path / "figure.png").is_file() + + +@pytest.mark.skipif( + not (XYNOVA_EXPORT / "mesh.npz").is_file(), reason="export not built" +) +def test_cli_no_plot_skips_the_figure(tmp_path: Path) -> None: + from orisys_sensor_design import cli + + revision = REPO / "sensors" / "finger-xynova" / "r01" + args = type( + "Args", + (), + { + "profile": "compute", + "output_dir": XYNOVA_EXPORT, + "samples": 8, + "out": tmp_path / "report.json", + "plot": tmp_path / "figure.png", + "no_plot": True, + "sanity_mm": 5.0, + }, + )() + assert cli.cmd_optics(revision, args) == 0 + assert (tmp_path / "report.json").is_file() + assert not (tmp_path / "figure.png").exists() + + +@pytest.mark.skipif( + not (XYNOVA_EXPORT / "mesh.npz").is_file(), reason="export not built" +) +def test_cli_without_matplotlib_exits_two_without_tracing( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """A missing plotting dependency must fail fast, before the ray trace.""" + from orisys_sensor_design import cli + + monkeypatch.setitem(sys.modules, "matplotlib", None) + revision = REPO / "sensors" / "finger-xynova" / "r01" + args = type( + "Args", + (), + { + "profile": "compute", + "output_dir": XYNOVA_EXPORT, + "samples": 8, + "out": tmp_path / "report.json", + "plot": tmp_path / "figure.png", + "no_plot": False, + "sanity_mm": 5.0, + }, + )() + assert cli.cmd_optics(revision, args) == 2 + assert not (tmp_path / "report.json").exists()