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 <noreply@anthropic.com>

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",

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
+Models the camera as a full pinhole camera (3D position + yaw/pitch/roll
-the beam's transverse plane, compensating for an oblique camera viewing
+orientation + focal length + principal point) shared across all measurement
-angle (which compresses the image along the tilt axis by cos(angle)).
+planes in one reconstruction. Every nominal value on `CameraModel` is paired
-Known calibration values (on the MeasurementPlane) are used directly; when a
+with a `CameraModelTolerance` entry that determines whether `ModalFitter`
-value is unknown, an override must be supplied (e.g. by ModalFitter while
+holds it fixed (tolerance == 0) or refines it within a bound
-exploring it as a free parameter).
+(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):
+    def __init__(self, camera: CameraModel):
-        self.plane = plane
+        self.camera = camera
+        self._rotation = _rotation_matrix(*camera.orientation_deg)
 
-    @property
+    def pixel_coordinates(
-    def pixel_scale_known(self) -> bool:
+        self, x: np.ndarray, y: np.ndarray, z: float
-        return self.plane.pixel_scale is not None
+    ) -> 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
+        r = self._rotation
-    def viewing_angle_known(self) -> bool:
+        xc = r[0, 0] * dx + r[1, 0] * dy + r[2, 0] * dz
-        return self.plane.viewing_angle_deg is not None
+        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,
+        self, image_shape: tuple[int, int], z: float
-        pixel_scale: float | None = None,
-        viewing_angle_deg: float | None = None,
     ) -> 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
+        Casts a ray from the camera through each pixel and intersects it
-        only used to fill in values that are not known/calibrated.
+        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
+        rows, cols = image_shape
-        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
         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))
+        f = self.camera.focal_length_px
-        x = col_grid * scale / cos_angle
+        cx, cy = self.camera.principal_point
-        y = row_grid * scale
+        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)

tests/test_geometry.py

@@ -1,77 +1,151 @@
 import numpy as np
 import pytest
 
-from he11lib.data import MeasurementPlane
+from he11lib.geometry import CameraModel, CameraModelTolerance, GeometryCalibration
-from he11lib.geometry import GeometryCalibration
 
 
-def test_pixel_scale_known_reflects_plane():
+def test_camera_model_tolerance_accepts_zero_and_positive():
-    plane_known = MeasurementPlane(flux=np.ones((5, 5)), z=0.3, pixel_scale=1e-4)
+    CameraModelTolerance(
-    plane_unknown = MeasurementPlane(flux=np.ones((5, 5)), z=0.3)
+        focal_length_px=0.0,
-
+        position=(0.0, 0.0, 0.0),
-    assert GeometryCalibration(plane_known).pixel_scale_known is True
+        orientation_deg=(1.0, 2.0, 3.0),
-    assert GeometryCalibration(plane_unknown).pixel_scale_known is False
+        principal_point=(0.5, 0.5),
+    )  # should not raise
 
 
-def test_viewing_angle_known_reflects_plane():
+def test_camera_model_tolerance_rejects_negative_scalar_field():
-    plane_known = MeasurementPlane(flux=np.ones((5, 5)), z=0.3, viewing_angle_deg=10.0)
+    with pytest.raises(ValueError, match="focal_length_px"):
-    plane_unknown = MeasurementPlane(flux=np.ones((5, 5)), z=0.3)
+        CameraModelTolerance(
-
+            focal_length_px=-1.0,
-    assert GeometryCalibration(plane_known).viewing_angle_known is True
+            position=(0.0, 0.0, 0.0),
-    assert GeometryCalibration(plane_unknown).viewing_angle_known is False
+            orientation_deg=(0.0, 0.0, 0.0),
+        )
 
 
-def test_physical_coordinates_uses_known_calibration():
+def test_camera_model_tolerance_rejects_negative_tuple_component():
-    plane = MeasurementPlane(
+    with pytest.raises(ValueError, match="position"):
-        flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=0.0
+        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():
+def make_tilted_camera():
-    plane = MeasurementPlane(
+    return CameraModel(
-        flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=60.0
+        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():
+@pytest.mark.parametrize(
-    plane = MeasurementPlane(flux=np.ones((5, 5)), z=0.3)
+    "camera",
-    calib = GeometryCalibration(plane)
+    [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"):
+    x, y = calib.physical_coordinates(image_shape, z)
-        calib.physical_coordinates()
+    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():
+def test_keystone_regression_uniform_for_on_axis_camera():
-    plane = MeasurementPlane(flux=np.ones((5, 5)), z=0.3)
+    # A camera with zero orientation, centered on the beam axis, produces
-    calib = GeometryCalibration(plane)
+    # 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)
+    xs = np.array([-0.02, -0.01, 0.0, 0.01, 0.02])
-    col_idx = np.arange(5) - 2
+    ys = np.zeros_like(xs)
-    np.testing.assert_allclose(x[2, :], col_idx * 1e-4)
+    _, 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():
+def test_keystone_regression_nonuniform_for_tilted_camera():
-    plane = MeasurementPlane(
+    # A tilted/off-axis camera produces non-uniform pixel spacing for the
-        flux=np.ones((5, 5)), z=0.3, pixel_scale=2e-4, viewing_angle_deg=0.0
+    # 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
+    with pytest.raises(ValueError):
-    x, _ = calib.physical_coordinates(pixel_scale=999.0, viewing_angle_deg=45.0)
+        calib.pixel_coordinates(np.array([0.0]), np.array([0.0]), z=0.5)
-    col_idx = np.arange(5) - 2
+
-    np.testing.assert_allclose(x[2, :], col_idx * 2e-4)
+
+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)