added optics test (undergoing)

This commit is contained in:
ykcheungab
2026-09-23 10:01:08 +08:00
parent f9662df08c
commit bc5a928852
11 changed files with 2099 additions and 9 deletions
+144 -7
View File
@@ -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,
+9
View File
@@ -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"
+24
View File
@@ -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
@@ -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:
Binary file not shown.
+265 -1
View File
@@ -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: <export revision dir>/refraction_report.json)",
)
optics_p.add_argument(
"--plot", type=Path, default=None,
help="Comparison PNG path (default: <export revision dir>/"
"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:
+5 -1
View File
@@ -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:
+794
View File
@@ -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 <revision_dir>"
)
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
+228
View File
@@ -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,
}
+308
View File
@@ -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)
+318
View File
@@ -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()