From 4f65c2ce4f735eec0e73156bfcb7551b99148179 Mon Sep 17 00:00:00 2001 From: Martino Ferrari Date: Fri, 3 Jul 2026 11:32:31 +0200 Subject: [PATCH] Replace cosine-compression geometry model with a full pinhole CameraModel GeometryCalibration now performs true perspective forward/inverse projection (with genuine keystoning) around a shared CameraModel, paired with a CameraModelTolerance that will drive ModalFitter's per-field fixed/refined behavior in a later task. Co-Authored-By: Claude Sonnet 5 --- he11lib/__init__.py | 4 +- he11lib/geometry.py | 262 ++++++++++++++++++++++++++++++++++------- tests/test_geometry.py | 182 +++++++++++++++++++--------- 3 files changed, 352 insertions(+), 96 deletions(-) diff --git a/he11lib/__init__.py b/he11lib/__init__.py index 790d428..b5ec294 100644 --- a/he11lib/__init__.py +++ b/he11lib/__init__.py @@ -6,7 +6,7 @@ See docs/ for the full API and design documentation. from .data import MeasurementPlane, ReconstructionResult from .deconvolution import DiffusionDeconvolver from .fitting import ModalFitter, generate_mode_shells -from .geometry import GeometryCalibration +from .geometry import CameraModel, CameraModelTolerance, GeometryCalibration from .modes import LGBasis from .noise import NoiseEstimator from .phase_retrieval import PhaseRetrievalResult, PhaseRetriever, propagate_angular_spectrum @@ -21,6 +21,8 @@ __all__ = [ "DiffusionDeconvolver", "ModalFitter", "generate_mode_shells", + "CameraModel", + "CameraModelTolerance", "GeometryCalibration", "LGBasis", "NoiseEstimator", diff --git a/he11lib/geometry.py b/he11lib/geometry.py index 0b423e5..edd25b4 100644 --- a/he11lib/geometry.py +++ b/he11lib/geometry.py @@ -1,64 +1,244 @@ -"""Camera geometry correction: pixel-to-physical scale and viewing angle. +"""Camera geometry: a shared pinhole camera model and pixel<->physical mapping. -Converts a MeasurementPlane's pixel grid into physical (x, y) coordinates in -the beam's transverse plane, compensating for an oblique camera viewing -angle (which compresses the image along the tilt axis by cos(angle)). -Known calibration values (on the MeasurementPlane) are used directly; when a -value is unknown, an override must be supplied (e.g. by ModalFitter while -exploring it as a free parameter). +Models the camera as a full pinhole camera (3D position + yaw/pitch/roll +orientation + focal length + principal point) shared across all measurement +planes in one reconstruction. Every nominal value on `CameraModel` is paired +with a `CameraModelTolerance` entry that determines whether `ModalFitter` +holds it fixed (tolerance == 0) or refines it within a bound +(tolerance > 0) -- `CameraModel` alone is never trusted as exact. + +Coordinate conventions +---------------------- +World frame: `x` increases along the pixel-column direction, `y` increases +along the pixel-row direction, `z` is distance from the output window along +the beam axis (target planes live at `z = const > 0`). + +Camera frame: `X_cam` = right (pixel-column direction), `Y_cam` = down +(pixel-row direction), `Z_cam` = boresight (depth). At +`orientation_deg == (0, 0, 0)`, the camera frame is axis-aligned with the +world frame, so the boresight points along `+z` -- normal to every +`z = const` target plane, with no in-plane rotation. + +`orientation_deg = (yaw, pitch, roll)` composes as +`R = R_yaw(about Y) @ R_pitch(about X) @ R_roll(about Z)`, applied to the +camera axes to obtain their world-frame directions. """ from __future__ import annotations +from dataclasses import dataclass, fields +from typing import Sequence + import numpy as np -from .data import MeasurementPlane +CAMERA_FIELD_NAMES: tuple[str, ...] = ( + "focal_length_px", + "position_x", + "position_y", + "position_z", + "yaw_deg", + "pitch_deg", + "roll_deg", + "principal_point_x", + "principal_point_y", +) + + +@dataclass +class CameraModel: + """Nominal pinhole camera parameters, shared across all measurement planes. + + Never trusted as exact by itself -- pair with a `CameraModelTolerance` + to express how much each field may be refined during fitting. + + Parameters + ---------- + focal_length_px : focal length, in pixel units. + position : (x, y, z) camera position in the world (beam-axis) frame, + in meters. z=0 is the output window. + orientation_deg : (yaw, pitch, roll), in degrees. All-zero means the + boresight is normal to every z=const target plane, no in-plane + rotation (see module docstring for the full convention). + principal_point : (px, px) offset of the principal point from the frame + center, in pixels. + """ + + focal_length_px: float + position: tuple[float, float, float] + orientation_deg: tuple[float, float, float] + principal_point: tuple[float, float] = (0.0, 0.0) + + +@dataclass +class CameraModelTolerance: + """+/- bound (same units as `CameraModel`) within which each field is refined. + + `0` holds the paired `CameraModel` field fixed at its nominal value; + `> 0` bounds it to `[nominal - tolerance, nominal + tolerance]` during + fitting. All fields must be `>= 0`. + """ + + focal_length_px: float + position: tuple[float, float, float] + orientation_deg: tuple[float, float, float] + principal_point: tuple[float, float] = (0.0, 0.0) + + def __post_init__(self) -> None: + for f in fields(self): + value = getattr(self, f.name) + components = value if isinstance(value, tuple) else (value,) + for component in components: + if component < 0: + raise ValueError( + f"CameraModelTolerance.{f.name} must be >= 0, got {value}" + ) + + +def camera_to_values(camera: CameraModel) -> list[float]: + """Flatten a `CameraModel` into the 9 scalars named by `CAMERA_FIELD_NAMES`.""" + return [ + camera.focal_length_px, + camera.position[0], + camera.position[1], + camera.position[2], + camera.orientation_deg[0], + camera.orientation_deg[1], + camera.orientation_deg[2], + camera.principal_point[0], + camera.principal_point[1], + ] + + +def tolerance_to_values(tolerance: CameraModelTolerance) -> list[float]: + """Flatten a `CameraModelTolerance` into the 9 scalars named by `CAMERA_FIELD_NAMES`.""" + return [ + tolerance.focal_length_px, + tolerance.position[0], + tolerance.position[1], + tolerance.position[2], + tolerance.orientation_deg[0], + tolerance.orientation_deg[1], + tolerance.orientation_deg[2], + tolerance.principal_point[0], + tolerance.principal_point[1], + ] + + +def camera_from_values(values: Sequence[float]) -> CameraModel: + """Inverse of `camera_to_values`: rebuild a `CameraModel` from 9 scalars.""" + return CameraModel( + focal_length_px=values[0], + position=(values[1], values[2], values[3]), + orientation_deg=(values[4], values[5], values[6]), + principal_point=(values[7], values[8]), + ) + + +def _rotation_matrix(yaw_deg: float, pitch_deg: float, roll_deg: float) -> np.ndarray: + """3x3 rotation matrix mapping camera-frame axes to world-frame directions.""" + yaw = np.deg2rad(yaw_deg) + pitch = np.deg2rad(pitch_deg) + roll = np.deg2rad(roll_deg) + + cy, sy = np.cos(yaw), np.sin(yaw) + cx, sx = np.cos(pitch), np.sin(pitch) + cz, sz = np.cos(roll), np.sin(roll) + + r_yaw = np.array([[cy, 0.0, sy], [0.0, 1.0, 0.0], [-sy, 0.0, cy]]) + r_pitch = np.array([[1.0, 0.0, 0.0], [0.0, cx, -sx], [0.0, sx, cx]]) + r_roll = np.array([[cz, -sz, 0.0], [sz, cz, 0.0], [0.0, 0.0, 1.0]]) + + return r_yaw @ r_pitch @ r_roll class GeometryCalibration: - """Resolves pixel scale / viewing angle and builds a physical coordinate grid.""" + """Resolves the pixel<->physical mapping for a shared pinhole `CameraModel`.""" - def __init__(self, plane: MeasurementPlane): - self.plane = plane + def __init__(self, camera: CameraModel): + self.camera = camera + self._rotation = _rotation_matrix(*camera.orientation_deg) - @property - def pixel_scale_known(self) -> bool: - return self.plane.pixel_scale is not None + def pixel_coordinates( + self, x: np.ndarray, y: np.ndarray, z: float + ) -> tuple[np.ndarray, np.ndarray]: + """Forward pinhole projection: physical (x, y) at depth z -> centered pixel (row, col).""" + px, py, pz = self.camera.position + dx = x - px + dy = y - py + dz = z - pz - @property - def viewing_angle_known(self) -> bool: - return self.plane.viewing_angle_deg is not None + r = self._rotation + xc = r[0, 0] * dx + r[1, 0] * dy + r[2, 0] * dz + yc = r[0, 1] * dx + r[1, 1] * dy + r[2, 1] * dz + zc = r[0, 2] * dx + r[1, 2] * dy + r[2, 2] * dz + + if np.any(zc <= 0): + raise ValueError( + f"One or more target points are behind or edge-on to the " + f"camera at z={z}; cannot project." + ) + + f = self.camera.focal_length_px + cx, cy = self.camera.principal_point + col = f * xc / zc + cx + row = f * yc / zc + cy + return row, col def physical_coordinates( - self, - pixel_scale: float | None = None, - viewing_angle_deg: float | None = None, + self, image_shape: tuple[int, int], z: float ) -> tuple[np.ndarray, np.ndarray]: - """Physical (x, y) grid matching the plane's flux array shape. + """Inverse pinhole projection: pixel grid at depth z -> physical (x, y). - Known values on the MeasurementPlane take precedence; overrides are - only used to fill in values that are not known/calibrated. + Casts a ray from the camera through each pixel and intersects it + with the world plane z=const. Raises ValueError if the target + plane is edge-on to (parallel to) the view direction or behind the + camera for this pose. """ - scale = self.plane.pixel_scale if self.pixel_scale_known else pixel_scale - angle_deg = ( - self.plane.viewing_angle_deg if self.viewing_angle_known else viewing_angle_deg - ) - - if scale is None: - raise ValueError( - "pixel_scale is not known for this MeasurementPlane and no override was given" - ) - if angle_deg is None: - raise ValueError( - "viewing_angle_deg is not known for this MeasurementPlane and no override was given" - ) - - rows, cols = self.plane.flux.shape + rows, cols = image_shape row_idx = np.arange(rows) - rows // 2 col_idx = np.arange(cols) - cols // 2 col_grid, row_grid = np.meshgrid(col_idx, row_idx) - cos_angle = np.cos(np.deg2rad(angle_deg)) - x = col_grid * scale / cos_angle - y = row_grid * scale + f = self.camera.focal_length_px + cx, cy = self.camera.principal_point + dir_cam_x = (col_grid - cx) / f + dir_cam_y = (row_grid - cy) / f + dir_cam_z = np.ones_like(dir_cam_x) + + r = self._rotation + dir_world_x = r[0, 0] * dir_cam_x + r[0, 1] * dir_cam_y + r[0, 2] * dir_cam_z + dir_world_y = r[1, 0] * dir_cam_x + r[1, 1] * dir_cam_y + r[1, 2] * dir_cam_z + dir_world_z = r[2, 0] * dir_cam_x + r[2, 1] * dir_cam_y + r[2, 2] * dir_cam_z + + if np.any(np.abs(dir_world_z) < 1e-12): + raise ValueError( + f"Camera pose is edge-on to the target plane z={z}; no " + "valid ray-plane intersection." + ) + + px, py, pz = self.camera.position + t = (z - pz) / dir_world_z + if np.any(t <= 0): + raise ValueError( + f"Target plane z={z} is behind the camera for this pose; " + "no valid ray-plane intersection." + ) + + x = px + t * dir_world_x + y = py + t * dir_world_y return x, y + + def effective_pixel_scale(self, image_shape: tuple[int, int], z: float) -> float: + """Isotropic finite-difference approximation of the local pixel scale. + + `DiffusionDeconvolver` assumes one isotropic pixel-space blur + kernel; this is only exact for an on-axis, zero-orientation + camera, and an approximation whenever the true projection is + keystoned. + """ + rows, cols = image_shape + x, y = self.physical_coordinates(image_shape, z) + mid_row, mid_col = rows // 2, cols // 2 + dx = abs(x[mid_row, mid_col + 1] - x[mid_row, mid_col]) + dy = abs(y[mid_row + 1, mid_col] - y[mid_row, mid_col]) + return float((dx + dy) / 2) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 6fdd399..6179841 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -1,77 +1,151 @@ import numpy as np import pytest -from he11lib.data import MeasurementPlane -from he11lib.geometry import GeometryCalibration +from he11lib.geometry import CameraModel, CameraModelTolerance, GeometryCalibration -def test_pixel_scale_known_reflects_plane(): - plane_known = MeasurementPlane(flux=np.ones((5, 5)), z=0.3, pixel_scale=1e-4) - plane_unknown = MeasurementPlane(flux=np.ones((5, 5)), z=0.3) - - assert GeometryCalibration(plane_known).pixel_scale_known is True - assert GeometryCalibration(plane_unknown).pixel_scale_known is False +def test_camera_model_tolerance_accepts_zero_and_positive(): + CameraModelTolerance( + focal_length_px=0.0, + position=(0.0, 0.0, 0.0), + orientation_deg=(1.0, 2.0, 3.0), + principal_point=(0.5, 0.5), + ) # should not raise -def test_viewing_angle_known_reflects_plane(): - plane_known = MeasurementPlane(flux=np.ones((5, 5)), z=0.3, viewing_angle_deg=10.0) - plane_unknown = MeasurementPlane(flux=np.ones((5, 5)), z=0.3) - - assert GeometryCalibration(plane_known).viewing_angle_known is True - assert GeometryCalibration(plane_unknown).viewing_angle_known is False +def test_camera_model_tolerance_rejects_negative_scalar_field(): + with pytest.raises(ValueError, match="focal_length_px"): + CameraModelTolerance( + focal_length_px=-1.0, + position=(0.0, 0.0, 0.0), + orientation_deg=(0.0, 0.0, 0.0), + ) -def test_physical_coordinates_uses_known_calibration(): - plane = MeasurementPlane( - flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=0.0 +def test_camera_model_tolerance_rejects_negative_tuple_component(): + with pytest.raises(ValueError, match="position"): + CameraModelTolerance( + focal_length_px=1.0, + position=(0.0, -0.5, 0.0), + orientation_deg=(0.0, 0.0, 0.0), + ) + + +def make_on_axis_camera(focal_length_px=2000.0, camera_z=-2.0): + return CameraModel( + focal_length_px=focal_length_px, + position=(0.0, 0.0, camera_z), + orientation_deg=(0.0, 0.0, 0.0), ) - calib = GeometryCalibration(plane) - x, y = calib.physical_coordinates() - - row_idx = np.arange(5) - 2 - col_idx = np.arange(5) - 2 - expected_x = col_idx * 2e-4 - expected_y = row_idx * 2e-4 - np.testing.assert_allclose(x[2, :], expected_x) - np.testing.assert_allclose(y[:, 2], expected_y) -def test_physical_coordinates_compresses_x_for_viewing_angle(): - plane = MeasurementPlane( - flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=60.0 +def make_tilted_camera(): + return CameraModel( + focal_length_px=2000.0, + position=(0.05, -0.03, -2.0), + orientation_deg=(8.0, -5.0, 3.0), ) - calib = GeometryCalibration(plane) - x, y = calib.physical_coordinates() - - col_idx = np.arange(5) - 2 - expected_x = col_idx * 2e-4 / np.cos(np.deg2rad(60.0)) - np.testing.assert_allclose(x[2, :], expected_x) -def test_physical_coordinates_raises_without_calibration_or_override(): - plane = MeasurementPlane(flux=np.ones((5, 5)), z=0.3) - calib = GeometryCalibration(plane) +@pytest.mark.parametrize( + "camera", + [make_on_axis_camera(), make_tilted_camera()], + ids=["on_axis", "tilted_off_center"], +) +@pytest.mark.parametrize("z", [0.3, 0.5, 0.8]) +def test_projection_round_trip_recovers_pixel_grid(camera, z): + image_shape = (41, 41) + calib = GeometryCalibration(camera) - with pytest.raises(ValueError, match="pixel_scale"): - calib.physical_coordinates() + x, y = calib.physical_coordinates(image_shape, z) + row, col = calib.pixel_coordinates(x, y, z) + + rows, cols = image_shape + row_idx = np.arange(rows) - rows // 2 + col_idx = np.arange(cols) - cols // 2 + expected_col, expected_row = np.meshgrid(col_idx, row_idx) + + np.testing.assert_allclose(row, expected_row, atol=1e-6) + np.testing.assert_allclose(col, expected_col, atol=1e-6) -def test_physical_coordinates_accepts_override_for_unknown_values(): - plane = MeasurementPlane(flux=np.ones((5, 5)), z=0.3) - calib = GeometryCalibration(plane) +def test_keystone_regression_uniform_for_on_axis_camera(): + # A camera with zero orientation, centered on the beam axis, produces + # uniform pixel spacing for evenly spaced physical points (no keystoning). + camera = make_on_axis_camera() + calib = GeometryCalibration(camera) + z = 0.5 - x, y = calib.physical_coordinates(pixel_scale=1e-4, viewing_angle_deg=0.0) - col_idx = np.arange(5) - 2 - np.testing.assert_allclose(x[2, :], col_idx * 1e-4) + xs = np.array([-0.02, -0.01, 0.0, 0.01, 0.02]) + ys = np.zeros_like(xs) + _, col = calib.pixel_coordinates(xs, ys, z) + + spacings = np.diff(col) + np.testing.assert_allclose(spacings, spacings[0], rtol=1e-6) -def test_known_calibration_takes_precedence_over_override(): - plane = MeasurementPlane( - flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=0.0 +def test_keystone_regression_nonuniform_for_tilted_camera(): + # A tilted/off-axis camera produces non-uniform pixel spacing for the + # same evenly spaced physical points -- genuine keystoning. + camera = make_tilted_camera() + calib = GeometryCalibration(camera) + z = 0.5 + + xs = np.array([-0.02, -0.01, 0.0, 0.01, 0.02]) + ys = np.zeros_like(xs) + _, col = calib.pixel_coordinates(xs, ys, z) + + spacings = np.diff(col) + assert not np.allclose(spacings, spacings[0], rtol=1e-3) + + +def test_pixel_coordinates_raises_when_point_behind_camera(): + camera = CameraModel( + focal_length_px=2000.0, + position=(0.0, 0.0, 10.0), + orientation_deg=(0.0, 0.0, 0.0), ) - calib = GeometryCalibration(plane) + calib = GeometryCalibration(camera) - # override should be ignored since plane already specifies calibration - x, _ = calib.physical_coordinates(pixel_scale=999.0, viewing_angle_deg=45.0) - col_idx = np.arange(5) - 2 - np.testing.assert_allclose(x[2, :], col_idx * 2e-4) + with pytest.raises(ValueError): + calib.pixel_coordinates(np.array([0.0]), np.array([0.0]), z=0.5) + + +def test_physical_coordinates_raises_when_plane_behind_camera(): + # Camera sits downstream of the target plane and looks further + # downstream (boresight = +z world) -- the z=0.5 plane is behind it. + camera = CameraModel( + focal_length_px=2000.0, + position=(0.0, 0.0, 10.0), + orientation_deg=(0.0, 0.0, 0.0), + ) + calib = GeometryCalibration(camera) + + with pytest.raises(ValueError): + calib.physical_coordinates((21, 21), z=0.5) + + +def test_physical_coordinates_raises_when_edge_on(): + # Pitch=90 deg points the boresight along world -y, making the + # z=const target plane edge-on (parallel to the view direction). + camera = CameraModel( + focal_length_px=2000.0, + position=(0.0, 0.0, -2.0), + orientation_deg=(0.0, 90.0, 0.0), + ) + calib = GeometryCalibration(camera) + + with pytest.raises(ValueError): + calib.physical_coordinates((41, 41), z=0.5) + + +def test_effective_pixel_scale_matches_on_axis_focal_length(): + focal_length_px = 2000.0 + camera_z = -2.0 + z = 0.5 + camera = make_on_axis_camera(focal_length_px=focal_length_px, camera_z=camera_z) + calib = GeometryCalibration(camera) + + scale = calib.effective_pixel_scale((41, 41), z) + expected = (z - camera_z) / focal_length_px + assert scale == pytest.approx(expected, rel=1e-6)