From 4e116af3be06c12fe3eaf9c185b465805503c8e5 Mon Sep 17 00:00:00 2001 From: Mike Rosseel Date: Tue, 28 Jul 2026 14:28:54 +0200 Subject: [PATCH 1/2] feat(solver): solve the full sensor frame with cached optical calibration The camera cropped every exposure to a centred square before solving, throwing away roughly half the sensor area (imx462: 1920x1080 down to 980x980). Cedar Solve has no trouble with the whole frame, so that discarded area is free plate-solving sensitivity. The camera now publishes two frames per exposure: the 512x512 display frame, unchanged, and a solve frame covering the whole sensor at native scale. Bayer sensors are binned 2x2 for the solve frame, which removes the RGGB modulation the old downscale to 512x512 used to smooth away. Everything downstream of the solver -- SQM photometry, the preview overlay, the alignment marker -- is written against the display frame, and SQM in particular is calibrated to its plate scale. So solves are projected back onto the display frame in one place (project_solution_to_display), keeping those consumers unchanged. SolveGeometry owns the mapping, composed from the affine transform of each pipeline stage. The solver also stopped rediscovering the optics on every frame. It passed a fixed 12 +/- 4 degree window and never passed a distortion estimate at all. The first successful solve of a run now captures the measured FOV and lens distortion; later solves get a +/-0.5 degree window and the measured distortion as a starting point, which speeds up pattern matching and improves centroid matching towards the frame corners -- the part of the image the crop used to discard. A run of failures discards the calibration and reopens the window. Full frame is on by default and can be turned off with solver_full_frame. It is disabled automatically when the camera has no sensor profile, or when camera_rotation is set to something other than a quarter turn, since rotating a non-square frame by an arbitrary angle would clip the corners. The 2x2 Bayer bin is gone. It was there because the crop pipeline's downscale had been smoothing the RGGB checkerboard away implicitly, and losing that was assumed to bias centroids toward the green sites. Measured against 12 imx462 frames of real sky, the assumption was wrong and the bin was costing exactly what this change exists to gain: solve rate median matches median RMSE solves down to crop (today) 8/12 15.5 15.4" 257 ms binned 10/12 20.0 40.5" 174 ms un-binned 11/12 26.0 41.1" 97 ms Un-binned finds 30% more matches than binned and solves at 2.6x shorter exposure than the crop -- against 1.5x for the binned version, so binning was giving away half the sensitivity gain. Angular residuals are statistically identical (41.1" vs 40.5"), so the centroid-bias worry does not survive contact with data; the finer sampling is simply what lets marginal stars clear the detection threshold. Full frame does cost accuracy against the crop (15.4" -> 41"), from corner distortion a single-k model cannot fully absorb. That is a real trade, though 41" is negligible against any eyepiece field, and it buys solves on frames the crop cannot solve at all. --- default_config.json | 1 + python/PiFinder/camera_debug.py | 6 +- python/PiFinder/camera_interface.py | 112 ++++++- python/PiFinder/camera_none.py | 6 +- python/PiFinder/camera_pi.py | 88 ++++-- python/PiFinder/main.py | 6 + python/PiFinder/optics.py | 366 ++++++++++++++++++++++ python/PiFinder/solver.py | 130 +++++++- python/PiFinder/sqm/camera_profiles.py | 30 ++ python/PiFinder/state.py | 10 + python/tests/test_full_frame_solve.py | 190 ++++++++++++ python/tests/test_optics.py | 407 +++++++++++++++++++++++++ 12 files changed, 1300 insertions(+), 52 deletions(-) create mode 100644 python/PiFinder/optics.py create mode 100644 python/tests/test_full_frame_solve.py create mode 100644 python/tests/test_optics.py diff --git a/default_config.json b/default_config.json index 7684bff0f..b85f57c78 100644 --- a/default_config.json +++ b/default_config.json @@ -10,6 +10,7 @@ "screen_direction": "right", "mount_type": "Alt/Az", "solver_debug": 0, + "solver_full_frame": true, "sleep_timeout": "30s", "screen_off_timeout": "Off", "chart_radec": "Off", diff --git a/python/PiFinder/camera_debug.py b/python/PiFinder/camera_debug.py index 9e479836b..efbc69df6 100644 --- a/python/PiFinder/camera_debug.py +++ b/python/PiFinder/camera_debug.py @@ -101,7 +101,9 @@ def get_cam_type(self) -> str: return self.camType -def get_images(shared_state, camera_image, command_queue, console_queue, log_queue): +def get_images( + shared_state, camera_image, solve_image, command_queue, console_queue, log_queue +): """ Instantiates the camera hardware then calls the universal image loop @@ -116,5 +118,5 @@ def get_images(shared_state, camera_image, command_queue, console_queue, log_que camera_hardware = CameraDebug(exposure_time) camera_hardware.get_image_loop( - shared_state, camera_image, command_queue, console_queue, cfg + shared_state, camera_image, solve_image, command_queue, console_queue, cfg ) diff --git a/python/PiFinder/camera_interface.py b/python/PiFinder/camera_interface.py index 77cd3a03a..7f5b6d789 100644 --- a/python/PiFinder/camera_interface.py +++ b/python/PiFinder/camera_interface.py @@ -28,6 +28,7 @@ ExposureSNRController, generate_exposure_sweep, ) +from PiFinder.optics import DISPLAY_FRAME_SIZE, SolveGeometry, build_geometry logger = logging.getLogger("Camera.Interface") @@ -179,6 +180,15 @@ def initialize(self) -> None: def capture(self) -> Image.Image: return Image.Image() + def capture_pair(self) -> Tuple[Image.Image, Image.Image]: + """Return ``(display_frame, solve_frame)`` from a single exposure. + + Backends with no separate full-sensor readout solve the display frame, + so both halves are the same image. + """ + image = self.capture() + return image, image + def capture_file(self, filename) -> None: pass @@ -191,6 +201,12 @@ def _blank_capture(self): """ return Image.new("L", (512, 512), 0) # Black 512x512 image + def _capture_pair_with_timeout( + self, timeout=10 + ) -> Optional[Tuple[Image.Image, Image.Image]]: + """:meth:`capture_pair` under the same guard as :meth:`_capture_with_timeout`.""" + return self._run_capture_with_timeout(self.capture_pair, timeout) + def _capture_with_timeout(self, timeout=10) -> Optional[Image.Image]: """Run capture() with a timeout, never overlapping two captures. @@ -207,6 +223,10 @@ def _capture_with_timeout(self, timeout=10) -> Optional[Image.Image]: second capture while it is still alive. At most one capture is ever running; the caller just gets blank frames until the stuck one clears. """ + return self._run_capture_with_timeout(self.capture, timeout) + + def _run_capture_with_timeout(self, capture_fn, timeout=10): + """Shared timeout/overlap guard behind the capture wrappers.""" # A previous capture is still wedged in the driver -- don't start a # second one. Returning None lets the caller fall back to a blank frame # while we wait for the stuck capture to clear. @@ -218,7 +238,7 @@ def _capture_with_timeout(self, timeout=10) -> Optional[Image.Image]: def _do_capture(): try: - result[0] = self.capture() + result[0] = capture_fn() except Exception as e: # propagate to the caller's thread exc[0] = e @@ -262,8 +282,45 @@ def start_camera(self) -> None: def stop_camera(self) -> None: pass + def _configure_solve_geometry( + self, shared_state, cfg, solve_rotation + ) -> SolveGeometry: + """Decide whether to solve on the full sensor, and publish the mapping. + + Full frame needs a camera profile describing the sensor, and needs the + frame's orientation to be a multiple of 90 degrees -- an arbitrary + rotation would clip the corners off a non-square frame, which is + exactly the sky area full frame exists to recover. + """ + profile = getattr(self, "profile", None) + full_frame = bool(cfg.get_option("solver_full_frame", True)) + + if full_frame and profile is None: + logger.info("Camera has no sensor profile; solving on the display frame") + full_frame = False + if full_frame and solve_rotation % 90: + logger.warning( + "Solve rotation %s is not a quarter turn; solving on the display frame", + solve_rotation, + ) + full_frame = False + + geometry = build_geometry(profile, solve_rotation, full_frame) + shared_state.set_solve_geometry(geometry.as_dict()) + if geometry.full_frame: + logger.info( + "Full-frame solving enabled: solve frame %dx%d (display frame %dx%d)", + geometry.solve_width, + geometry.solve_height, + DISPLAY_FRAME_SIZE, + DISPLAY_FRAME_SIZE, + ) + else: + logger.info("Solving on the %dx%d display frame", *geometry.solve_size) + return geometry + def get_image_loop( - self, shared_state, camera_image, command_queue, console_queue, cfg + self, shared_state, camera_image, solve_image, command_queue, console_queue, cfg ): try: # Store shared_state for access by capture() methods @@ -300,6 +357,8 @@ def get_image_loop( solve_rotation = SCREEN_ROTATE_AMOUNTS.get(screen_direction, 270) shared_state.set_solve_image_rotation(solve_rotation) + geometry = self._configure_solve_geometry(shared_state, cfg, solve_rotation) + # Set path for test mode image root_dir = os.path.realpath( os.path.join(os.path.dirname(__file__), "..", "..") @@ -336,18 +395,41 @@ def get_image_loop( imu_start = shared_state.imu() image_start_time = time.time() if self._camera_started: + solve_frame = None if not test_mode_on: - base_image = self._capture_with_timeout() - if base_image is None: + if geometry.full_frame: + captured = self._capture_pair_with_timeout() + else: + image = self._capture_with_timeout() + captured = None if image is None else (image, None) + if captured is None: # Capture hung; fall back to a blank frame so the # loop keeps running and stays responsive to # commands instead of freezing. The blank frame # simply fails to solve. logger.warning("Camera capture timed out; blank frame") - base_image = self._blank_capture() + captured = (self._blank_capture(), None) + base_image, solve_frame = captured base_image = base_image.convert("L") base_image = base_image.rotate(solve_rotation) + + if solve_frame is not None: + solve_frame = solve_frame.convert("L") + # The solve frame is not square, so rotate it by + # transposition: Image.rotate() clips the corners. + if solve_rotation % 360 == 90: + solve_frame = solve_frame.transpose( + Image.Transpose.ROTATE_90 + ) + elif solve_rotation % 360 == 180: + solve_frame = solve_frame.transpose( + Image.Transpose.ROTATE_180 + ) + elif solve_rotation % 360 == 270: + solve_frame = solve_frame.transpose( + Image.Transpose.ROTATE_270 + ) else: # Test Mode: load image from disc and wait # No real raw matrix backs this frame; prevent a recent @@ -377,9 +459,23 @@ def get_image_loop( if test_mode_on and abs(pointing_diff) > 0.01: # Scope moved during the fake exposure: return a blank # image so the solver doesn't report a stale solve - camera_image.paste(self._blank_capture()) - else: - camera_image.paste(base_image) + base_image = self._blank_capture() + solve_frame = None + camera_image.paste(base_image) + if solve_frame is None: + if geometry.full_frame: + # No full-sensor frame this cycle (timed-out + # capture, or test mode). Publish a blank at the + # full solve size rather than a 512x512 frame, + # which would leave the previous exposure's stars + # sitting in the region outside it for the solver + # to find again. + solve_frame = Image.new("L", geometry.solve_size[::-1], 0) + else: + # Backends with no separate full-sensor readout + # solve the display frame. + solve_frame = base_image + solve_image.paste(solve_frame, (0, 0)) image_metadata = { "exposure_start": image_start_time, "exposure_end": image_end_time, diff --git a/python/PiFinder/camera_none.py b/python/PiFinder/camera_none.py index 3af549501..cac861ce2 100644 --- a/python/PiFinder/camera_none.py +++ b/python/PiFinder/camera_none.py @@ -56,7 +56,9 @@ def get_cam_type(self) -> str: return self.camType -def get_images(shared_state, camera_image, bias_image, command_queue, console_queue): +def get_images( + shared_state, camera_image, solve_image, command_queue, console_queue, log_queue +): """ Instantiates the camera hardware then calls the universal image loop @@ -70,5 +72,5 @@ def get_images(shared_state, camera_image, bias_image, command_queue, console_qu camera_hardware = CameraNone(exposure_time) camera_hardware.get_image_loop( - shared_state, camera_image, bias_image, command_queue, console_queue, cfg + shared_state, camera_image, solve_image, command_queue, console_queue, cfg ) diff --git a/python/PiFinder/camera_pi.py b/python/PiFinder/camera_pi.py index 4364bcd2f..357ef7006 100644 --- a/python/PiFinder/camera_pi.py +++ b/python/PiFinder/camera_pi.py @@ -73,16 +73,11 @@ def stop_camera(self) -> None: self.camera.stop() self._camera_started = False - def capture(self) -> Image.Image: - """ - Captures a raw 10/12bit sensor output and converts - it to an 8 bit mono image stretched to use the maximum - amount of the 255 level space. - """ + def _read_raw(self) -> np.ndarray: + """Read one raw 10/12bit sensor frame, uncropped.""" _request = self.camera.capture_request() # raw is actually 16 bit raw_capture = _request.make_array("raw").copy().view(np.uint16) - # tmp_image = _request.make_image("main") # Log actual camera metadata for exposure verification (debug level only) metadata = _request.get_metadata() @@ -108,21 +103,51 @@ def capture(self) -> Image.Image: self.last_frame_metadata = metadata _request.release() + return raw_capture + + def _to_8bit(self, raw_capture: np.ndarray) -> np.ndarray: + """Subtract the bias pedestal and stretch to the full 0-255 range.""" + # covert to 32 bit int to avoid overflow + raw_capture = raw_capture.astype(np.float32) + + # sensor offset (bias pedestal from camera profile) + raw_capture -= self.profile.bias_offset + + # apply digital gain + raw_capture *= self.profile.digital_gain + + # rescale to 8 bit + raw_capture = ( + raw_capture + * 255 + / (2**self.profile.bit_depth - self.profile.bias_offset - 1) + ) + + # clip to avoid <0 or >255 values + return np.clip(raw_capture.astype(np.int32), 0, 255).astype(np.uint8) - # Apply camera-specific crop and rotation - raw_capture = self.profile.crop_and_rotate(raw_capture) + def _display_frame(self, raw_capture: np.ndarray) -> Image.Image: + """Square-crop the sensor frame down to the 512x512 display frame. + + SQM photometry and the radiometer both measure the crop, so both are + fed from here rather than from the full-sensor solve frame. + """ + cropped = self.profile.crop_and_rotate(raw_capture) # Reduce the matrix while it is local to the camera process. The solver # can publish radiometric SQM from this small sample without a solve and # without copying/scanning the raw frame on every capture. if hasattr(self, "shared_state"): self._radiometer_sequence += 1 + actual_exposure = getattr(self, "last_frame_metadata", {}).get( + "ExposureTime" + ) try: radiometer_exposure = float(actual_exposure) / 1_000_000.0 except (TypeError, ValueError): radiometer_exposure = float(self.exposure_time) / 1_000_000.0 sample = collect_radiometer_sample( - raw_capture, + cropped, self.profile, radiometer_exposure, sequence=self._radiometer_sequence, @@ -133,31 +158,26 @@ def capture(self) -> Image.Image: # Store raw in shared state (before processing) for calibration and analysis if hasattr(self, "shared_state"): - self.shared_state.set_cam_raw(raw_capture.copy()) + self.shared_state.set_cam_raw(cropped.copy()) - # covert to 32 bit int to avoid overflow - raw_capture = raw_capture.astype(np.float32) + return Image.fromarray(self._to_8bit(cropped)).resize((512, 512)) - # sensor offset (bias pedestal from camera profile) - raw_capture -= self.profile.bias_offset + def _solve_frame(self, raw_capture: np.ndarray) -> Image.Image: + """Full sensor area at native scale, for the plate solver.""" + return Image.fromarray(self._to_8bit(self.profile.full_frame(raw_capture))) - # apply digital gain - raw_capture *= self.profile.digital_gain - - # rescale to 8 bit - raw_capture = ( - raw_capture - * 255 - / (2**self.profile.bit_depth - self.profile.bias_offset - 1) - ) - - # clip to avoid <0 or >255 values - raw_capture = np.clip(raw_capture.astype(np.int32), 0, 255).astype(np.uint8) - - # convert to PIL image and resize to 512x512 - raw_image = Image.fromarray(raw_capture).resize((512, 512)) + def capture(self) -> Image.Image: + """ + Captures a raw 10/12bit sensor output and converts + it to an 8 bit mono image stretched to use the maximum + amount of the 255 level space. + """ + return self._display_frame(self._read_raw()) - return raw_image + def capture_pair(self) -> Tuple[Image.Image, Image.Image]: + """Produce the display and solve frames from a single sensor read.""" + raw_capture = self._read_raw() + return self._display_frame(raw_capture), self._solve_frame(raw_capture) def capture_bias(self) -> np.ndarray: """Capture a bias frame for measuring black level offset. @@ -310,7 +330,9 @@ def get_cam_type(self) -> str: return self.camType -def get_images(shared_state, camera_image, command_queue, console_queue, log_queue): +def get_images( + shared_state, camera_image, solve_image, command_queue, console_queue, log_queue +): """ Instantiates the camera hardware then calls the universal image loop @@ -326,5 +348,5 @@ def get_images(shared_state, camera_image, command_queue, console_queue, log_que camera_hardware = CameraPI(exposure_time) camera_hardware.get_image_loop( - shared_state, camera_image, command_queue, console_queue, cfg + shared_state, camera_image, solve_image, command_queue, console_queue, cfg ) diff --git a/python/PiFinder/main.py b/python/PiFinder/main.py index d21252791..c45f3143f 100644 --- a/python/PiFinder/main.py +++ b/python/PiFinder/main.py @@ -53,6 +53,7 @@ from PiFinder.ui.console import UIConsole from PiFinder.ui.menu_manager import MenuManager +from PiFinder.optics import max_solve_frame_size from PiFinder.state import SharedStateObj, UIState from PiFinder.image_util import subtract_background @@ -510,12 +511,16 @@ def main( logger.info(" Camera") console.update() camera_image = manager.NewImage("RGB", (512, 512)) # type: ignore[attr-defined] + # Sized for the largest sensor we support; the camera publishes its + # actual solve frame into the top-left corner. + solve_image = manager.NewImage("L", max_solve_frame_size()) # type: ignore[attr-defined] image_process = Process( name="Camera", target=camera.get_images, args=( shared_state, camera_image, + solve_image, camera_command_queue, console_queue, camera_logqueue, @@ -580,6 +585,7 @@ def main( shared_state, solver_queue, camera_image, + solve_image, console_queue, solver_logqueue, alignment_command_queue, diff --git a/python/PiFinder/optics.py b/python/PiFinder/optics.py new file mode 100644 index 000000000..9c5e29657 --- /dev/null +++ b/python/PiFinder/optics.py @@ -0,0 +1,366 @@ +""" +Solve-frame geometry and optical calibration. + +The camera publishes two frames per exposure: + +* the **display frame** — a 512x512 square crop of the sensor, consumed by the + UI, focus, and SQM. Its geometry is unchanged from the square-crop pipeline. +* the **solve frame** — the full sensor area at native scale, consumed only by + the plate solver. Removing the crop roughly doubles the sky area searched for + stars. + +Because the two frames have different origins and scales, anything that crosses +between them (``target_pixel``, alignment results, SQM's matched centroids) has +to be mapped. :class:`SolveGeometry` owns that mapping, built by composing the +affine transform of every stage of each pipeline. + +:class:`OpticalCalibration` holds the FOV and lens distortion measured by the +first successful solve of a run, so later solves can be given tight bounds +instead of re-deriving them from scratch. +""" + +from __future__ import annotations + +import logging +import math +from dataclasses import dataclass +from typing import Any, Dict, Optional, Sequence, Tuple + +import numpy as np + +from PiFinder.sqm.camera_profiles import CAMERA_PROFILES + +logger = logging.getLogger("Optics") + +# Side length of the square display frame. +DISPLAY_FRAME_SIZE = 512 + +# Bounds of the shipped tetra3 database (degrees of horizontal FOV). A solve +# whose measured FOV falls outside this can never match. +DB_MIN_FOV = 10.0 +DB_MAX_FOV = 30.0 + +# Half-width of the FOV search window once calibration has measured the true +# value, in degrees. +CALIBRATED_FOV_MAX_ERROR = 0.5 + +# Search window for the square-crop pipeline, where the FOV is ~10.2 deg for +# every supported camera. The full-frame window is derived from this per camera +# by SolveGeometry.solve_fov(), since how much wider the frame gets depends on +# the sensor's aspect ratio and how it is mounted. +CROP_FOV_ESTIMATE = 12.0 +CROP_FOV_MAX_ERROR = 4.0 + +# Consecutive failures after which a calibration is assumed stale and the +# solver falls back to the wide search window. +FAILURES_BEFORE_RECALIBRATION = 20 + + +def _identity() -> np.ndarray: + return np.eye(3) + + +def _translate(dx: float, dy: float) -> np.ndarray: + return np.array([[1.0, 0.0, dx], [0.0, 1.0, dy], [0.0, 0.0, 1.0]]) + + +def _scale(sx: float, sy: float) -> np.ndarray: + return np.array([[sx, 0.0, 0.0], [0.0, sy, 0.0], [0.0, 0.0, 1.0]]) + + +def _rot90_ccw(width: int, height: int, k: int) -> Tuple[np.ndarray, int, int]: + """Affine for ``np.rot90(array, k)`` on an array of shape ``(height, width)``. + + ``np.rot90`` with ``k=1`` maps input element ``(row, col)`` to output + ``(width - 1 - col, row)``, i.e. in (x, y) terms ``x' = y`` and + ``y' = (width - 1) - x``. Returns the transform plus the resulting + ``(width, height)``. + """ + matrix = _identity() + for _ in range(k % 4): + # x' = y, y' = (width - 1) - x + step = np.array([[0.0, 1.0, 0.0], [-1.0, 0.0, width - 1], [0.0, 0.0, 1.0]]) + matrix = step @ matrix + width, height = height, width + return matrix, width, height + + +@dataclass(frozen=True) +class SolveGeometry: + """Maps between the 512x512 display frame and the full-frame solve frame. + + Coordinates are handled in ``(y, x)`` order throughout, matching the + convention used by tetra3 centroids and ``shared_state.target_pixel()``. + """ + + solve_width: int + solve_height: int + # Homogeneous (x, y, 1) transform taking display-frame pixels to + # solve-frame pixels. + display_to_solve_matrix: np.ndarray + full_frame: bool + + @property + def solve_size(self) -> Tuple[int, int]: + """``(height, width)``, the order ``solve_from_centroids`` expects.""" + return (self.solve_height, self.solve_width) + + @property + def solve_to_display_matrix(self) -> np.ndarray: + return np.linalg.inv(self.display_to_solve_matrix) + + def display_to_solve(self, yx: Sequence[float]) -> Tuple[float, float]: + """Map one ``(y, x)`` display-frame point into the solve frame.""" + y, x = yx + vec = self.display_to_solve_matrix @ np.array([x, y, 1.0]) + return (float(vec[1]), float(vec[0])) + + def solve_to_display(self, yx: Sequence[float]) -> Tuple[float, float]: + """Map one ``(y, x)`` solve-frame point into the display frame.""" + y, x = yx + vec = self.solve_to_display_matrix @ np.array([x, y, 1.0]) + return (float(vec[1]), float(vec[0])) + + def solve_to_display_array(self, points_yx: np.ndarray) -> np.ndarray: + """Vectorised :meth:`solve_to_display` over an ``(N, 2)`` array.""" + points_yx = np.asarray(points_yx, dtype=float) + if points_yx.size == 0: + return points_yx.reshape(0, 2) + homogeneous = np.column_stack( + (points_yx[:, 1], points_yx[:, 0], np.ones(len(points_yx))) + ) + mapped = homogeneous @ self.solve_to_display_matrix.T + return np.column_stack((mapped[:, 1], mapped[:, 0])) + + @property + def display_width_in_solve_px(self) -> float: + """Width of the display frame, measured in solve-frame pixels. + + The x basis vector's image under the mapping gives solve pixels per + display pixel. + """ + basis = self.display_to_solve_matrix[:2, 0] + return DISPLAY_FRAME_SIZE * float(np.hypot(*basis)) + + def display_fov(self, solve_fov: float) -> float: + """Convert a horizontal FOV measured on the solve frame to the display frame. + + Uses the gnomonic relation ``f = (width / 2) / tan(fov / 2)`` rather than + scaling the angle linearly, which would be off by over a percent at the + fields these lenses produce. + + Note that "horizontal" is the solve frame's own x axis, which after the + camera's quarter-turn is the sensor's *short* side. Full frame therefore + widens the horizontal field only modestly (or not at all, on sensors + whose square crop already spans the short side); most of the extra sky + arrives vertically. + """ + focal_px = (self.solve_width / 2.0) / math.tan(math.radians(solve_fov) / 2.0) + return math.degrees( + 2.0 * math.atan((self.display_width_in_solve_px / 2.0) / focal_px) + ) + + def solve_fov(self, display_fov: float) -> float: + """Inverse of :meth:`display_fov`.""" + focal_px = (self.display_width_in_solve_px / 2.0) / math.tan( + math.radians(display_fov) / 2.0 + ) + return math.degrees(2.0 * math.atan((self.solve_width / 2.0) / focal_px)) + + def as_dict(self) -> Dict[str, Any]: + return { + "solve_width": self.solve_width, + "solve_height": self.solve_height, + "display_to_solve_matrix": self.display_to_solve_matrix.tolist(), + "full_frame": self.full_frame, + } + + @classmethod + def from_dict(cls, data: Dict[str, Any]) -> "SolveGeometry": + return cls( + solve_width=int(data["solve_width"]), + solve_height=int(data["solve_height"]), + display_to_solve_matrix=np.array(data["display_to_solve_matrix"]), + full_frame=bool(data["full_frame"]), + ) + + +def max_solve_frame_size() -> Tuple[int, int]: + """Largest solve frame any supported camera can produce, in ``(width, height)``. + + The shared solve-frame buffer is allocated once in the main process, before + the camera process has detected which sensor is fitted, so it is sized to + fit the largest possibility. Each frame is published into the top-left + corner and the solver crops it back to the published geometry. + """ + largest = DISPLAY_FRAME_SIZE + for profile in CAMERA_PROFILES.values(): + largest = max(largest, *profile.solve_frame_size) + return (largest, largest) + + +def identity_geometry(size: int = DISPLAY_FRAME_SIZE) -> SolveGeometry: + """Geometry for cameras whose solve frame *is* the display frame. + + Used by the debug and none cameras, and by the square-crop pipeline. + """ + return SolveGeometry( + solve_width=size, + solve_height=size, + display_to_solve_matrix=_identity(), + full_frame=False, + ) + + +def build_geometry( + profile, + final_rotation: int, + full_frame: bool, +) -> SolveGeometry: + """Compose the display and solve pipelines into a display->solve mapping. + + ``final_rotation`` is the whole-frame rotation the camera loop applies after + the profile's own crop/rotation, in degrees counter-clockwise. + + Display pipeline: crop -> ``rot90(profile.rotation_90)`` -> resize to + 512x512 -> ``final_rotation``. + + Solve pipeline: optional 2x2 Bayer bin -> ``rot90(profile.rotation_90)`` -> + ``final_rotation``. + """ + if not full_frame: + return identity_geometry() + + raw_width, raw_height = profile.raw_size + + # --- display pipeline --- + crop_top, crop_bottom = profile.crop_y + crop_left, crop_right = profile.crop_x + cropped_width = raw_width - crop_left - crop_right + cropped_height = raw_height - crop_top - crop_bottom + + display = _translate(-crop_left, -crop_top) + rot, width, height = _rot90_ccw(cropped_width, cropped_height, profile.rotation_90) + display = rot @ display + display = _scale(DISPLAY_FRAME_SIZE / width, DISPLAY_FRAME_SIZE / height) @ display + rot, _, _ = _rot90_ccw(DISPLAY_FRAME_SIZE, DISPLAY_FRAME_SIZE, final_rotation // 90) + display = rot @ display + + # --- solve pipeline: the whole sensor at native sampling --- + solve = _identity() + rot, width, height = _rot90_ccw(raw_width, raw_height, profile.rotation_90) + solve = rot @ solve + rot, width, height = _rot90_ccw(width, height, final_rotation // 90) + solve = rot @ solve + + return SolveGeometry( + solve_width=width, + solve_height=height, + display_to_solve_matrix=solve @ np.linalg.inv(display), + full_frame=True, + ) + + +class OpticalCalibration: + """Session-scoped cache of the measured FOV and lens distortion. + + The first successful solve of a run is made with a wide FOV window and + ``distortion=0`` so cedar-solve derives both. Its measured values are then + kept for the rest of the run: subsequent solves get a tight FOV window and + the measured distortion as a starting point, which both speeds up pattern + matching and improves centroid matching towards the frame corners — the part + of the image the square crop used to throw away. + + Calibration is not persisted; it is cheap enough to redo at every startup, + and it stays correct across lens or camera swaps for free. + """ + + def __init__(self, fallback_fov: float, fallback_fov_max_error: float) -> None: + self._fallback_fov = fallback_fov + self._fallback_fov_max_error = fallback_fov_max_error + self.fov: Optional[float] = None + self.distortion: Optional[float] = None + self._consecutive_failures = 0 + + @property + def calibrated(self) -> bool: + return self.fov is not None + + def solver_args(self) -> Dict[str, Any]: + """FOV and distortion arguments for ``solve_from_centroids``.""" + if not self.calibrated: + return { + "fov_estimate": self._fallback_fov, + "fov_max_error": self._fallback_fov_max_error, + "distortion": 0, + } + return { + "fov_estimate": self.fov, + "fov_max_error": CALIBRATED_FOV_MAX_ERROR, + "distortion": self.distortion, + } + + def record_success(self, solution: Dict[str, Any]) -> None: + self._consecutive_failures = 0 + if self.calibrated: + return + + fov = solution.get("FOV") + if fov is None: + return + if not DB_MIN_FOV <= fov <= DB_MAX_FOV: + logger.warning( + "Measured FOV %.2f deg is outside the database range %.1f-%.1f; " + "not calibrating", + fov, + DB_MIN_FOV, + DB_MAX_FOV, + ) + return + + self.fov = float(fov) + distortion = solution.get("distortion") + self.distortion = float(distortion) if distortion is not None else 0.0 + logger.info( + "Optical calibration: FOV %.3f deg, distortion %.4f " + "(subsequent solves use +/-%.1f deg)", + self.fov, + self.distortion, + CALIBRATED_FOV_MAX_ERROR, + ) + + def record_failure(self) -> None: + if not self.calibrated: + return + self._consecutive_failures += 1 + if self._consecutive_failures >= FAILURES_BEFORE_RECALIBRATION: + logger.warning( + "%d consecutive solve failures; discarding optical calibration", + self._consecutive_failures, + ) + self.fov = None + self.distortion = None + self._consecutive_failures = 0 + + @classmethod + def for_geometry(cls, geometry: SolveGeometry) -> "OpticalCalibration": + """Seed the pre-calibration search window from the square-crop one. + + The square crop has always solved at ~10.2 degrees, so scale that known + window onto whatever frame we are actually solving instead of guessing a + range wide enough for every camera. + """ + estimate = geometry.solve_fov(CROP_FOV_ESTIMATE) + upper = geometry.solve_fov(CROP_FOV_ESTIMATE + CROP_FOV_MAX_ERROR) + return cls(estimate, upper - estimate) + + def describe(self) -> str: + if not self.calibrated: + return ( + f"uncalibrated (FOV est: {self._fallback_fov:.1f} deg, " + f"max err: {self._fallback_fov_max_error:.1f} deg)" + ) + return ( + f"FOV {self.fov:.2f} deg +/-{CALIBRATED_FOV_MAX_ERROR:.1f}, " + f"distortion {self.distortion:.4f}" + ) diff --git a/python/PiFinder/solver.py b/python/PiFinder/solver.py index f6fc17d88..d3d0bfd9d 100644 --- a/python/PiFinder/solver.py +++ b/python/PiFinder/solver.py @@ -22,9 +22,16 @@ import subprocess import threading from multiprocessing import shared_memory +from typing import Optional import grpc from PiFinder import state_utils +from PiFinder.optics import ( + DISPLAY_FRAME_SIZE, + OpticalCalibration, + SolveGeometry, + identity_geometry, +) from PiFinder import utils from PiFinder import timez from PiFinder.sqm import SQM as SQMCalculator @@ -101,6 +108,87 @@ def _scaled_photometry_radii( return aperture, inner, outer +def project_solution_to_display(solution: dict, geometry: SolveGeometry) -> dict: + """Re-express a full-frame solve in display-frame (512x512) coordinates. + + Everything downstream of the solver -- SQM photometry, the preview overlay, + the alignment marker -- is written against the 512x512 display frame. SQM in + particular then scales those centroids onto the raw photometry image and + derotates them, a chain that only works if it starts from display-frame + pixels. So a full-frame solve is mapped back here, once, and the projected + copy is what the rest of the system sees. + + Matched stars outside the square crop are dropped along with their catalogue + counterparts: they contributed to the solve, which is the point, but they + have no pixels in the crop for photometry to measure. + + The pointing itself (RA/Dec/Roll and target sky coordinates) is + frame-independent: both frames are concentric and share the same "up". + """ + if not geometry.full_frame: + return solution + + # A failed solve comes back with every value None. There is nothing to + # project; the caller forwards it unchanged and builds a FailedSolve from + # it. Without this the solver raises on every starless frame -- which is + # most of them indoors, at dusk, or under cloud. + if solution.get("FOV") is None: + return solution + + projected = dict(solution) + projected["FOV"] = geometry.display_fov(solution["FOV"]) + + centroids = np.asarray(solution.get("matched_centroids") or [], dtype=float) + stars = np.asarray(solution.get("matched_stars") or [], dtype=float) + if len(centroids): + mapped = geometry.solve_to_display_array(centroids) + inside = ( + (mapped[:, 0] >= 0) + & (mapped[:, 0] < DISPLAY_FRAME_SIZE) + & (mapped[:, 1] >= 0) + & (mapped[:, 1] < DISPLAY_FRAME_SIZE) + ) + projected["matched_centroids"] = mapped[inside].tolist() + projected["matched_stars"] = stars[inside].tolist() + + # Alignment returns the image position of a requested sky coordinate; the + # UI stores it straight back as target_pixel, so it has to come back in + # display-frame pixels. + y_target = solution.get("y_target") + x_target = solution.get("x_target") + if y_target is not None and x_target is not None: + y_display, x_display = geometry.solve_to_display((y_target, x_target)) + if 0 <= y_display < DISPLAY_FRAME_SIZE and 0 <= x_display < DISPLAY_FRAME_SIZE: + projected["y_target"] = y_display + projected["x_target"] = x_display + else: + # On the sensor but outside the crop: no display pixel to align on. + projected["y_target"] = None + projected["x_target"] = None + + return projected + + +def _resolve_geometry(published): + """Pair the camera's published solve geometry with a fresh calibration. + + The camera process publishes the geometry once it has detected the sensor, + which may be after the solver's first pass; until then (``published`` is + None) solve the display frame with the square-crop FOV window. + """ + geometry = ( + identity_geometry() if published is None else SolveGeometry.from_dict(published) + ) + calibration = OpticalCalibration.for_geometry(geometry) + logger.info( + "Solving %dx%d frames, %s", + geometry.solve_width, + geometry.solve_height, + calibration.describe(), + ) + return geometry, calibration + + def _scale_solution_centroids(solution, scale): """Return a shallow copy of solution with matched_centroids scaled. @@ -824,6 +912,7 @@ def solver( shared_state, solver_queue, camera_image, + solve_image, console_queue, log_queue, align_command_queue, @@ -843,6 +932,12 @@ def solver( centroids = [] log_no_stars_found = True + # Solve-frame geometry is published by the camera process once it has + # detected the sensor, so it is picked up lazily on the first frame. + geometry: Optional[SolveGeometry] = None + calibration: Optional[OpticalCalibration] = None + geometry_published = False + # SQM calculator is created lazily on the first radiometer sample (or solve # in test mode), not here: at solver # startup shared_state.camera_type() still holds the pre-camera default, @@ -955,7 +1050,18 @@ def solver( ) try: - img = camera_image.copy() + if not geometry_published: + published = shared_state.solve_geometry() + if published is not None or geometry is None: + geometry, calibration = _resolve_geometry(published) + geometry_published = published is not None + + # The shared solve buffer is sized for the largest sensor we + # support, so trim it back to the frame the camera actually + # published before pulling it across the manager. + img = solve_image.crop( + (0, 0, geometry.solve_width, geometry.solve_height) + ) img = img.convert(mode="L") np_image = np.asarray(img, dtype=np.uint8) @@ -1000,15 +1106,24 @@ def solver( solution = t3.solve_from_centroids( centroids, - (512, 512), - fov_estimate=12.0, - fov_max_error=4.0, + geometry.solve_size, match_max_error=0.005, return_matches=True, # Required for SQM calculation - target_pixel=shared_state.target_pixel(), + target_pixel=geometry.display_to_solve( + shared_state.target_pixel() + ), solve_timeout=1000, + **calibration.solver_args(), **_solver_args, ) + if solution.get("RA") is not None: + calibration.record_success(solution) + else: + calibration.record_failure() + # Everything below works in display-frame pixels: SQM + # scales these centroids onto the raw photometry image + # and derotates them, which only works from there. + solution = project_solution_to_display(solution, geometry) if "matched_centroids" in solution: if sqm_calculator is None: @@ -1103,8 +1218,9 @@ def solver( else: if solution: logger.warning( - f"Solve FAILED - {len(centroids)} centroids detected but " - f"pattern match failed (FOV est: 12.0°, max err: 4.0°)" + "Solve FAILED - %d centroids detected but pattern " + "match failed (%s)" + % (len(centroids), calibration.describe()) ) solver_queue.put( _build_failed_solve( diff --git a/python/PiFinder/sqm/camera_profiles.py b/python/PiFinder/sqm/camera_profiles.py index 32bca4988..1fa87a8fd 100644 --- a/python/PiFinder/sqm/camera_profiles.py +++ b/python/PiFinder/sqm/camera_profiles.py @@ -179,6 +179,36 @@ def ensure_cropped(self, raw_array): return self.crop_and_rotate(raw_array) return raw_array + @property + def solve_frame_size(self) -> Tuple[int, int]: + """``(width, height)`` of the solve frame, before any rotation.""" + return self.raw_size + + def full_frame(self, raw_array): + """Apply rotation to the whole sensor area, with no crop. + + The Bayer mosaic is handed to star detection as-is. An earlier version + binned 2x2 first, on the assumption that the crop pipeline's downscale + had been usefully smoothing the RGGB checkerboard away and that + centroids would otherwise be pulled toward the green sites. Measured on + real sky (12 imx462 frames, 2026-07-31) that assumption was wrong: the + bin costs matches without buying accuracy. Un-binned found a median 26 + matches against the bin's 20, solved 11 of 12 frames against 10, and + solved down to 97ms exposure where the bin needed 174ms -- at + statistically identical angular residuals (41.1" vs 40.5"). The finer + sampling is what lets marginal stars clear the detection threshold. + + Args: + raw_array: Raw sensor data (numpy array) + + Returns: + Rotated array covering the whole sensor + """ + if self.rotation_90 != 0: + raw_array = np.rot90(raw_array, self.rotation_90) + + return raw_array + def __repr__(self) -> str: return ( f"CameraProfile(" diff --git a/python/PiFinder/state.py b/python/PiFinder/state.py index 4b7c3cd0f..49ba6cca9 100644 --- a/python/PiFinder/state.py +++ b/python/PiFinder/state.py @@ -310,6 +310,10 @@ def __init__(self) -> None: self.__solve_image_rotation = None self.__cam_raw = None self.__sqm_radiometer_sample = None + # Mapping between the 512x512 display frame and the solve frame, set by + # the camera process once it knows which sensor is fitted. None until + # then, which the solver reads as "solve on the display frame". + self.__solve_geometry = None # Are we prepared to do alt/az math # We need gps lock and datetime self.__tz_finder = TimezoneFinder() @@ -380,6 +384,12 @@ def set_solve_image_rotation(self, v): def set_camera_type(self, v: str): self.__camera_type = v + def solve_geometry(self): + return self.__solve_geometry + + def set_solve_geometry(self, v): + self.__solve_geometry = v + def sats(self): return self.__sats diff --git a/python/tests/test_full_frame_solve.py b/python/tests/test_full_frame_solve.py new file mode 100644 index 000000000..30942cd31 --- /dev/null +++ b/python/tests/test_full_frame_solve.py @@ -0,0 +1,190 @@ +"""End-to-end check that a full-frame solve agrees with the square-crop solve. + +Takes a real 512x512 test frame that the square-crop pipeline solves, embeds it +into a larger canvas standing in for the full sensor, and solves that instead. +The pointing must come out the same, and the solver's coordinate projection must +put the matched stars back where the square-crop solve found them. +""" + +import sys +from pathlib import Path + +import numpy as np +import pytest +from PIL import Image + +from PiFinder import utils +from PiFinder.optics import DISPLAY_FRAME_SIZE, SolveGeometry +from PiFinder.solver import project_solution_to_display + +sys.path.append(str(utils.tetra3_dir)) + +tetra3 = pytest.importorskip("tetra3") + +TEST_IMAGE = utils.pifinder_dir / "test_images" / "pifinder_debug_02.png" +# Resolve the database next to whichever tetra3 actually got imported, rather +# than assuming the submodule layout. +DATABASE = Path(tetra3.__file__).parent / "data" / "default_database.npz" + +# Dimensions of the stand-in "full sensor". Deliberately not square and not a +# multiple of the display frame, so a mapping that only works for tidy numbers +# fails here. +SOLVE_WIDTH = 723 +SOLVE_HEIGHT = 941 + + +def _embedding_geometry(): + """Geometry for a display frame sitting centred in a larger solve frame. + + The offsets floor to whole pixels to match where the canvas paste below + actually lands; the real pipelines are concentric to sub-pixel precision. + """ + offset_x = float((SOLVE_WIDTH - DISPLAY_FRAME_SIZE) // 2) + offset_y = float((SOLVE_HEIGHT - DISPLAY_FRAME_SIZE) // 2) + return SolveGeometry( + solve_width=SOLVE_WIDTH, + solve_height=SOLVE_HEIGHT, + display_to_solve_matrix=np.array( + [[1.0, 0.0, offset_x], [0.0, 1.0, offset_y], [0.0, 0.0, 1.0]] + ), + full_frame=True, + ) + + +@pytest.fixture(scope="module") +def solver_db(): + if not DATABASE.exists(): + pytest.skip("tetra3 database not available (submodule not initialised)") + return tetra3.Tetra3(str(DATABASE)) + + +@pytest.fixture(scope="module") +def display_image(): + if not TEST_IMAGE.exists(): + pytest.skip(f"missing test image {TEST_IMAGE}") + return np.asarray(Image.open(TEST_IMAGE).convert("L"), dtype=np.uint8) + + +def _solve(db, image, size, **kwargs): + centroids = tetra3.get_centroids_from_image(image) + return centroids, db.solve_from_centroids( + centroids, + size, + match_max_error=0.005, + return_matches=True, + solve_timeout=10000, + **kwargs, + ) + + +@pytest.fixture(scope="module") +def crop_solve(solver_db, display_image): + _, solution = _solve( + solver_db, + display_image, + (DISPLAY_FRAME_SIZE, DISPLAY_FRAME_SIZE), + fov_estimate=12.0, + fov_max_error=4.0, + ) + if solution.get("RA") is None: + pytest.skip(f"reference image did not solve: {solution.get('status')}") + return solution + + +@pytest.fixture(scope="module") +def full_frame_solve(solver_db, display_image, crop_solve): + geometry = _embedding_geometry() + canvas = np.zeros((SOLVE_HEIGHT, SOLVE_WIDTH), dtype=np.uint8) + top = (SOLVE_HEIGHT - DISPLAY_FRAME_SIZE) // 2 + left = (SOLVE_WIDTH - DISPLAY_FRAME_SIZE) // 2 + canvas[top : top + DISPLAY_FRAME_SIZE, left : left + DISPLAY_FRAME_SIZE] = ( + display_image + ) + + estimate = geometry.solve_fov(crop_solve["FOV"]) + _, solution = _solve( + solver_db, + canvas, + geometry.solve_size, + fov_estimate=estimate, + fov_max_error=1.0, + distortion=0, + ) + if solution.get("RA") is None: + pytest.skip(f"full-frame image did not solve: {solution.get('status')}") + return geometry, solution + + +@pytest.mark.integration +def test_full_frame_solve_points_where_the_crop_solve_points( + crop_solve, full_frame_solve +): + """Both frames are concentric, so they share a centre on the sky.""" + _, solution = full_frame_solve + + assert solution["RA"] == pytest.approx(crop_solve["RA"], abs=0.05) + assert solution["Dec"] == pytest.approx(crop_solve["Dec"], abs=0.05) + assert solution["Roll"] == pytest.approx(crop_solve["Roll"], abs=0.5) + + +@pytest.mark.integration +def test_projected_fov_matches_the_crop_solve(crop_solve, full_frame_solve): + """display_fov() must undo the widening, or SQM's plate scale goes wrong.""" + geometry, solution = full_frame_solve + + projected = project_solution_to_display(solution, geometry) + + assert solution["FOV"] > crop_solve["FOV"] + assert projected["FOV"] == pytest.approx(crop_solve["FOV"], rel=0.02) + + +@pytest.mark.integration +def test_projected_centroids_land_inside_the_display_frame(full_frame_solve): + geometry, solution = full_frame_solve + + projected = project_solution_to_display(solution, geometry) + centroids = np.array(projected["matched_centroids"]) + + assert len(centroids) > 0 + assert (centroids >= 0).all() + assert (centroids < DISPLAY_FRAME_SIZE).all() + + +@pytest.mark.integration +def test_projected_centroids_sit_on_the_same_stars(crop_solve, full_frame_solve): + """Projected stars must land on the stars the crop solve matched. + + Not every one of them will: the full-frame solve matches strictly more + stars, and the ones the crop solve missed have no counterpart to compare + against. What the mapping owes us is that the stars in common coincide to + well under a pixel, and that most of them are in common at all. + """ + geometry, solution = full_frame_solve + + projected = project_solution_to_display(solution, geometry) + projected_centroids = np.array(projected["matched_centroids"]) + crop_centroids = np.array(crop_solve["matched_centroids"]) + + distances = np.array( + [ + np.min(np.linalg.norm(crop_centroids - centroid, axis=1)) + for centroid in projected_centroids + ] + ) + coincident = distances < 1.0 + + assert coincident.mean() > 0.5, "too few stars shared with the crop solve" + assert np.median(distances[coincident]) < 0.5 + + +@pytest.mark.integration +def test_target_pixel_round_trips_through_the_solve(solver_db, crop_solve): + """A display-frame target must come back as the same sky coordinate.""" + geometry = _embedding_geometry() + target_display = (200.0, 320.0) + + mapped = geometry.display_to_solve(target_display) + + assert geometry.solve_to_display(mapped) == pytest.approx(target_display) + assert 0 <= mapped[0] < SOLVE_HEIGHT + assert 0 <= mapped[1] < SOLVE_WIDTH diff --git a/python/tests/test_optics.py b/python/tests/test_optics.py new file mode 100644 index 000000000..a955b0de0 --- /dev/null +++ b/python/tests/test_optics.py @@ -0,0 +1,407 @@ +"""Tests for solve-frame geometry and optical calibration.""" + +import math + +import numpy as np +import pytest +from PIL import Image + +from PiFinder.optics import ( + CALIBRATED_FOV_MAX_ERROR, + DISPLAY_FRAME_SIZE, + FAILURES_BEFORE_RECALIBRATION, + OpticalCalibration, + SolveGeometry, + build_geometry, + identity_geometry, + max_solve_frame_size, +) +from PiFinder.sqm.camera_profiles import CAMERA_PROFILES + +SENSOR_PROFILES = ["imx296", "imx462", "imx290", "hq"] +ROTATIONS = [90, 270] + + +# The display path downscales the sensor by up to 3x, so a single-pixel marker +# can be dropped entirely. Use a block wide enough to survive, and locate it by +# centroid so the two pipelines' different resampling cancels out. +MARKER_HALF_WIDTH = 16 + + +def _marker_raw(profile, marker_yx): + """A raw sensor frame that is zero except for one bright square.""" + width, height = profile.raw_size + raw = np.zeros((height, width), dtype=np.uint16) + y, x = marker_yx + raw[ + y - MARKER_HALF_WIDTH : y + MARKER_HALF_WIDTH, + x - MARKER_HALF_WIDTH : x + MARKER_HALF_WIDTH, + ] = 4095 + return raw + + +def _marker_centre(array): + """The (y, x) centroid of the bright marker in an array.""" + lit = np.argwhere(array > array.max() / 2) + assert len(lit), "marker did not survive the pipeline" + return tuple(lit.mean(axis=0)) + + +def _run_display_pipeline(profile, raw, rotation): + """Reproduce the camera's display path: crop, rot90, resize, rotate.""" + cropped = profile.crop_and_rotate(raw) + image = Image.fromarray(np.asarray(cropped >> 4, dtype=np.uint8)) + image = image.resize((DISPLAY_FRAME_SIZE, DISPLAY_FRAME_SIZE), Image.NEAREST) + return np.asarray(image.rotate(rotation)) + + +def _run_solve_pipeline(profile, raw, rotation): + """Reproduce the camera's solve path: optional 2x2 bin, rot90, transpose.""" + full = profile.full_frame(raw) + image = Image.fromarray(np.asarray(full, dtype=np.uint16) >> 4) + if rotation == 90: + image = image.transpose(Image.Transpose.ROTATE_90) + elif rotation == 270: + image = image.transpose(Image.Transpose.ROTATE_270) + return np.asarray(image) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_solve_frame_dimensions_match_the_pipeline(name, rotation): + """The geometry's advertised solve size is what the camera actually produces.""" + profile = CAMERA_PROFILES[name] + raw = _marker_raw(profile, (0, 0)) + geometry = build_geometry(profile, rotation, full_frame=True) + + produced = _run_solve_pipeline(profile, raw, rotation) + + assert produced.shape == (geometry.solve_height, geometry.solve_width) + assert geometry.solve_size == produced.shape + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_display_to_solve_maps_the_same_sensor_pixel(name, rotation): + """A star lands where the mapping says it should in both frames. + + Drives one bright raw pixel through both real pipelines and checks that + transforming its display-frame position gives its solve-frame position. + """ + profile = CAMERA_PROFILES[name] + width, height = profile.raw_size + crop_left, _ = profile.crop_x + crop_top, _ = profile.crop_y + + # Somewhere inside the square crop, off-centre on both axes so a + # transposed or mirrored mapping cannot pass by accident. + marker = ( + crop_top + (height - 2 * crop_top) // 3, + crop_left + (width - 2 * crop_left) // 4, + ) + raw = _marker_raw(profile, marker) + geometry = build_geometry(profile, rotation, full_frame=True) + + display_yx = _marker_centre(_run_display_pipeline(profile, raw, rotation)) + solve_yx = _marker_centre(_run_solve_pipeline(profile, raw, rotation)) + + mapped = geometry.display_to_solve(display_yx) + + # The display frame is a coarser sampling of the sensor, so a display pixel + # covers several solve pixels; allow that much slack. + scale = max(profile.raw_size) / DISPLAY_FRAME_SIZE + assert mapped[0] == pytest.approx(solve_yx[0], abs=scale + 1) + assert mapped[1] == pytest.approx(solve_yx[1], abs=scale + 1) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_display_solve_round_trip(name, rotation): + profile = CAMERA_PROFILES[name] + geometry = build_geometry(profile, rotation, full_frame=True) + + for point in [(0.0, 0.0), (255.5, 255.5), (100.0, 400.0), (511.0, 3.0)]: + assert geometry.solve_to_display( + geometry.display_to_solve(point) + ) == pytest.approx(point, abs=1e-6) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_display_frame_centre_maps_to_solve_frame_centre(name, rotation): + """Both frames are concentric, so RA/Dec of centre is shared.""" + profile = CAMERA_PROFILES[name] + geometry = build_geometry(profile, rotation, full_frame=True) + + centre = (DISPLAY_FRAME_SIZE - 1) / 2.0 + mapped = geometry.display_to_solve((centre, centre)) + + assert mapped[0] == pytest.approx((geometry.solve_height - 1) / 2.0, abs=1.0) + assert mapped[1] == pytest.approx((geometry.solve_width - 1) / 2.0, abs=1.0) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +def test_solve_to_display_array_matches_scalar(name): + geometry = build_geometry(CAMERA_PROFILES[name], 90, full_frame=True) + points = np.array([[10.0, 20.0], [300.0, 400.0], [1.0, 999.0]]) + + batch = geometry.solve_to_display_array(points) + + for point, mapped in zip(points, batch): + assert tuple(mapped) == pytest.approx(geometry.solve_to_display(point)) + + +@pytest.mark.unit +def test_solve_to_display_array_handles_empty_input(): + geometry = build_geometry(CAMERA_PROFILES["imx462"], 90, full_frame=True) + assert geometry.solve_to_display_array(np.empty((0, 2))).shape == (0, 2) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_display_fov_never_exceeds_the_solve_fov(name, rotation): + """The display frame is a crop of the solve frame, so it can only be narrower. + + On sensors whose square crop already spans the short side, the two are + equal after the camera's quarter turn -- the extra sky arrives vertically. + """ + geometry = build_geometry(CAMERA_PROFILES[name], rotation, full_frame=True) + + display_fov = geometry.display_fov(12.0) + + assert 0 < display_fov <= 12.0 + 1e-9 + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +@pytest.mark.parametrize("rotation", ROTATIONS) +def test_solve_fov_inverts_display_fov(name, rotation): + geometry = build_geometry(CAMERA_PROFILES[name], rotation, full_frame=True) + + assert geometry.display_fov(geometry.solve_fov(10.2)) == pytest.approx(10.2) + + +@pytest.mark.unit +def test_display_fov_recovers_the_known_crop_field(): + """imx462's square crop is the ~10.2 degree field the crop pipeline solves. + + Its solve frame is 1080 raw pixels wide (the sensor's short side, after the + quarter turn) against the crop's 980, so the horizontal field grows by about + a tenth -- the bulk of the 2x area gain is vertical. + """ + profile = CAMERA_PROFILES["imx462"] + geometry = build_geometry(profile, 90, full_frame=True) + solve_fov = geometry.solve_fov(10.2) + + assert geometry.solve_width == 1080 + assert solve_fov == pytest.approx(11.2, abs=0.1) + + # Gnomonic expectation: same focal length, narrower sensor width. + focal_px = (geometry.solve_width / 2) / math.tan(math.radians(solve_fov) / 2) + crop_width_native = 980 + expected = math.degrees(2 * math.atan((crop_width_native / 2) / focal_px)) + assert geometry.display_fov(solve_fov) == pytest.approx(expected, rel=1e-6) + + +@pytest.mark.unit +@pytest.mark.parametrize("name", SENSOR_PROFILES) +def test_full_frame_covers_more_sky_than_the_crop(name): + """The whole point: the solve frame samples more of the sensor.""" + profile = CAMERA_PROFILES[name] + geometry = build_geometry(profile, 90, full_frame=True) + + crop_width, _ = profile.crop_x + crop_top, _ = profile.crop_y + cropped_area = (profile.raw_size[0] - 2 * crop_width) * ( + profile.raw_size[1] - 2 * crop_top + ) + solve_area = np.prod(profile.solve_frame_size) + + assert solve_area > cropped_area + assert geometry.solve_width * geometry.solve_height > 0 + + +@pytest.mark.unit +def test_identity_geometry_is_a_no_op(): + geometry = identity_geometry() + + assert geometry.full_frame is False + assert geometry.solve_size == (DISPLAY_FRAME_SIZE, DISPLAY_FRAME_SIZE) + assert geometry.display_to_solve((17.0, 42.0)) == pytest.approx((17.0, 42.0)) + + +@pytest.mark.unit +def test_build_geometry_without_full_frame_is_identity(): + geometry = build_geometry(CAMERA_PROFILES["imx462"], 90, full_frame=False) + assert geometry.solve_size == (DISPLAY_FRAME_SIZE, DISPLAY_FRAME_SIZE) + + +@pytest.mark.unit +def test_geometry_survives_serialisation(): + """The geometry crosses a process boundary via shared state.""" + original = build_geometry(CAMERA_PROFILES["hq"], 270, full_frame=True) + + restored = SolveGeometry.from_dict(original.as_dict()) + + assert restored.solve_size == original.solve_size + assert restored.full_frame == original.full_frame + assert restored.display_to_solve((3.0, 7.0)) == pytest.approx( + original.display_to_solve((3.0, 7.0)) + ) + + +@pytest.mark.unit +def test_shared_buffer_fits_every_solve_frame(): + buffer_width, buffer_height = max_solve_frame_size() + + for profile in CAMERA_PROFILES.values(): + width, height = profile.solve_frame_size + # Either orientation, since the camera loop may rotate a quarter turn. + assert max(width, height) <= min(buffer_width, buffer_height) + + +@pytest.mark.unit +def test_solve_frame_is_the_whole_sensor(): + assert CAMERA_PROFILES["imx462"].solve_frame_size == (1920, 1080) + assert CAMERA_PROFILES["hq"].solve_frame_size == (2028, 1520) + assert CAMERA_PROFILES["imx296"].solve_frame_size == (1456, 1088) + + +@pytest.mark.unit +def test_full_frame_keeps_native_sampling(): + """No binning: the mosaic goes to star detection at full resolution. + + Binning was measured to cost matches without buying accuracy -- the finer + sampling is what lets marginal stars clear the detection threshold. + """ + profile = CAMERA_PROFILES["imx462"] + width, height = profile.raw_size + raw = np.zeros((height, width), dtype=np.uint16) + raw[0, 0], raw[0, 1], raw[1, 0], raw[1, 1] = 10, 20, 30, 40 + + full = profile.full_frame(raw) + + assert full.shape == (height, width) + np.testing.assert_array_equal(full[:2, :2], [[10, 20], [30, 40]]) + + +@pytest.mark.unit +def test_calibration_starts_wide_then_narrows(): + calibration = OpticalCalibration(16.0, 8.0) + + assert calibration.calibrated is False + assert calibration.solver_args() == { + "fov_estimate": 16.0, + "fov_max_error": 8.0, + "distortion": 0, + } + + calibration.record_success({"FOV": 19.87, "distortion": -0.021}) + + assert calibration.calibrated is True + assert calibration.solver_args() == { + "fov_estimate": pytest.approx(19.87), + "fov_max_error": CALIBRATED_FOV_MAX_ERROR, + "distortion": pytest.approx(-0.021), + } + + +@pytest.mark.unit +def test_calibration_keeps_its_first_measurement(): + """Later solves must not drag the calibration around.""" + calibration = OpticalCalibration(16.0, 8.0) + calibration.record_success({"FOV": 19.87, "distortion": -0.021}) + + calibration.record_success({"FOV": 12.0, "distortion": 0.5}) + + assert calibration.fov == pytest.approx(19.87) + assert calibration.distortion == pytest.approx(-0.021) + + +@pytest.mark.unit +@pytest.mark.parametrize("fov", [4.0, 45.0]) +def test_calibration_rejects_fov_outside_the_database_range(fov): + calibration = OpticalCalibration(16.0, 8.0) + + calibration.record_success({"FOV": fov, "distortion": 0.0}) + + assert calibration.calibrated is False + + +@pytest.mark.unit +def test_calibration_handles_a_solve_without_distortion(): + """distortion is absent from the solution when the caller disabled it.""" + calibration = OpticalCalibration(16.0, 8.0) + + calibration.record_success({"FOV": 13.6}) + + assert calibration.distortion == 0.0 + + +@pytest.mark.unit +def test_sustained_failure_discards_the_calibration(): + calibration = OpticalCalibration(16.0, 8.0) + calibration.record_success({"FOV": 19.87, "distortion": -0.021}) + + for _ in range(FAILURES_BEFORE_RECALIBRATION): + calibration.record_failure() + + assert calibration.calibrated is False + assert calibration.solver_args()["fov_max_error"] == 8.0 + + +@pytest.mark.unit +def test_a_success_resets_the_failure_run(): + calibration = OpticalCalibration(16.0, 8.0) + calibration.record_success({"FOV": 19.87, "distortion": -0.021}) + + for _ in range(FAILURES_BEFORE_RECALIBRATION - 1): + calibration.record_failure() + calibration.record_success({"FOV": 19.87, "distortion": -0.021}) + calibration.record_failure() + + assert calibration.calibrated is True + + +@pytest.mark.unit +def test_failed_solve_is_passed_through_unprojected(): + """A starless frame must not raise. + + tetra3 returns every value as None when it cannot match, and that is the + common case indoors, at dusk, or under cloud. Projecting it blindly raised + TypeError on every such frame -- caught only on real hardware, because the + integration test skips when a solve fails and the unit tests fed only + successful solutions. + """ + from PiFinder.solver import project_solution_to_display + + geometry = build_geometry(CAMERA_PROFILES["imx462"], 90, full_frame=True) + failed = { + "RA": None, + "Dec": None, + "Roll": None, + "FOV": None, + "distortion": None, + "Matches": None, + "T_solve": 12.3, + "status": "NO_MATCH", + } + + assert project_solution_to_display(failed, geometry) == failed + + +@pytest.mark.unit +def test_failed_solve_still_records_against_calibration(): + """The failure counter must advance even though nothing is projected.""" + calibration = OpticalCalibration(11.2, 4.4) + calibration.record_success({"FOV": 11.2, "distortion": -0.02}) + + calibration.record_failure() + + assert calibration.calibrated is True From 80a18ce00e0c0871f987f946ad210bca3869faed Mon Sep 17 00:00:00 2001 From: Mike Rosseel Date: Fri, 31 Jul 2026 01:53:09 +0200 Subject: [PATCH 2/2] feat(api): serve the raw sensor frame intact, and add rawfull /api/camera/raw rendered a uint16 sensor frame with Image.fromarray(arr, mode="L"), which reinterprets the 16-bit buffer as 8-bit. The result is interleaved-byte noise that still looks enough like an image to be believed -- I spent a night's captures on it before checking the histogram (median 1, p90 80, 1% at 255) and realising the frames were junk. 16-bit frames now render as mode "I;16" PNGs, preserving every ADU, and a test pins the round trip. Adds /api/camera/rawfull for the whole sensor including the margins the crop discards. Naming follows the exposure sweep's TIFFs: raw is the crop, rawfull is everything. The full frame is ~4 MB, so it is published on demand rather than on every capture -- the endpoint asks, the camera fulfils the request once and clears it. Between them these make the frames the pipeline actually works on retrievable from a running device, which is what debugging a solve failure in the field needs and what nothing currently offers. --- python/PiFinder/api_extensions.py | 67 +++++++++++++++++++++------ python/PiFinder/camera_pi.py | 19 +++++++- python/PiFinder/state.py | 23 +++++++++ python/tests/test_api_extensions.py | 41 ++++++++++++++++ python/tests/test_full_frame_solve.py | 5 +- 5 files changed, 139 insertions(+), 16 deletions(-) diff --git a/python/PiFinder/api_extensions.py b/python/PiFinder/api_extensions.py index 0bad6cf9a..563796055 100644 --- a/python/PiFinder/api_extensions.py +++ b/python/PiFinder/api_extensions.py @@ -46,6 +46,28 @@ def _png_response(img: Image.Image) -> Response: return Response(_pil_to_png_bytes(img), content_type="image/png") +def _raw_to_png(raw): + """Render a raw sensor frame as a PNG without destroying its values. + + Sensor frames are 10/12-bit held in uint16. Handing that buffer to + Image.fromarray(..., mode="L") reinterprets it as 8-bit and produces + interleaved-byte noise rather than an image -- which looked plausible + enough to waste a night's captures on. 16-bit frames therefore become + mode "I;16" PNGs, preserving every ADU for offline analysis. + """ + if hasattr(raw, "save"): # already a PIL image + return raw + + import numpy as np + + arr = np.asarray(raw) + if arr.ndim == 3: + return Image.fromarray(arr) + if arr.dtype == np.uint16: + return Image.fromarray(arr, mode="I;16") + return Image.fromarray(arr.astype(np.uint8), mode="L") + + def _pointing_to_dict(p): """Serialize a :class:`Pointing` (or ``None``) to a plain ``{RA, Dec, Roll}`` dict of floats.""" @@ -699,27 +721,46 @@ def api_screen(): @app.route("/api/camera/raw") def api_camera_raw(): - """Return the raw CMOS image, if available""" + """The cropped raw sensor frame -- what photometry measures.""" try: raw = server_instance.shared_state.cam_raw() if raw is None: return _json_response({"note": "No raw image available"}, 503) - # raw may be a PIL Image or a NumPy array - if hasattr(raw, "save"): - img = raw.convert("RGB") if raw.mode != "RGB" else raw - else: - import numpy as np - - arr = np.asarray(raw) - if arr.ndim == 2: - img = Image.fromarray(arr, mode="L").convert("RGB") - else: - img = Image.fromarray(arr) - return _png_response(img) + return _png_response(_raw_to_png(raw)) except Exception as e: logger.error("api/camera/raw error: %s", e) return _json_response({"error": str(e)}, 500) + @app.route("/api/camera/rawfull") + def api_camera_rawfull(): + """The whole sensor, uncropped -- margins included. + + Published on demand (the full frame is ~4 MB), so this asks the camera + for one and waits for the next capture to deliver it. Naming matches + the exposure sweep's "rawfull" TIFFs: raw is the crop, rawfull is + everything. + """ + import time as _time + + try: + state = server_instance.shared_state + state.set_cam_raw_full(None) + state.request_cam_raw_full() + # Long exposures make a capture cycle seconds long, so wait + # generously rather than reporting an absence that is just latency. + deadline = _time.time() + 15.0 + while _time.time() < deadline: + raw = state.cam_raw_full() + if raw is not None: + return _png_response(_raw_to_png(raw)) + _time.sleep(0.25) + return _json_response( + {"note": "Timed out waiting for a full-sensor frame"}, 504 + ) + except Exception as e: + logger.error("api/camera/rawfull error: %s", e) + return _json_response({"error": str(e)}, 500) + @app.route("/api/camera/debug") def api_camera_debug(): """Return the latest debug frame from the solver_debug_dumps directory""" diff --git a/python/PiFinder/camera_pi.py b/python/PiFinder/camera_pi.py index 357ef7006..3024d1e0a 100644 --- a/python/PiFinder/camera_pi.py +++ b/python/PiFinder/camera_pi.py @@ -103,6 +103,18 @@ def _read_raw(self) -> np.ndarray: self.last_frame_metadata = metadata _request.release() + # Serve a pending request for the uncropped sensor frame. Done here, + # in the shared read path, because this is where the frame still has + # its margins -- both capture() and capture_pair() go through it. On + # demand only: the full frame is ~4 MB and would cost that across the + # state manager on every capture. + if hasattr(self, "shared_state"): + try: + if self.shared_state.cam_raw_full_requested(): + self.shared_state.set_cam_raw_full(raw_capture.copy()) + except (BrokenPipeError, ConnectionResetError, AttributeError): + pass + return raw_capture def _to_8bit(self, raw_capture: np.ndarray) -> np.ndarray: @@ -139,8 +151,11 @@ def _display_frame(self, raw_capture: np.ndarray) -> Image.Image: # without copying/scanning the raw frame on every capture. if hasattr(self, "shared_state"): self._radiometer_sequence += 1 - actual_exposure = getattr(self, "last_frame_metadata", {}).get( - "ExposureTime" + # The driver's actual exposure when it reports one; the request + # otherwise. The read and this reduction are separate steps, so it + # comes back off self rather than from a local. + actual_exposure = (getattr(self, "last_frame_metadata", None) or {}).get( + "ExposureTime", self.exposure_time ) try: radiometer_exposure = float(actual_exposure) / 1_000_000.0 diff --git a/python/PiFinder/state.py b/python/PiFinder/state.py index 49ba6cca9..9c680e196 100644 --- a/python/PiFinder/state.py +++ b/python/PiFinder/state.py @@ -309,6 +309,8 @@ def __init__(self) -> None: # to the stored raw frame (PIL CCW). None until the camera reports. self.__solve_image_rotation = None self.__cam_raw = None + self.__cam_raw_full = None + self.__cam_raw_full_requested = False self.__sqm_radiometer_sample = None # Mapping between the 512x512 display frame and the solve frame, set by # the camera process once it knows which sensor is fitted. None until @@ -573,6 +575,27 @@ def cam_raw(self): def set_cam_raw(self, v): self.__cam_raw = v + def cam_raw_full(self): + return self.__cam_raw_full + + def set_cam_raw_full(self, v): + # Fulfilling the request clears it, so one request yields one frame + # rather than leaving the camera publishing 4 MB every capture. + self.__cam_raw_full = v + self.__cam_raw_full_requested = False + + def cam_raw_full_requested(self) -> bool: + return self.__cam_raw_full_requested + + def request_cam_raw_full(self) -> None: + """Ask the camera to publish the next frame uncropped. + + The full sensor frame is ~4 MB and would cost that on every capture if + published unconditionally, so it is served on demand: a caller sets + this flag, the camera fulfils it once and clears it. + """ + self.__cam_raw_full_requested = True + def sqm_radiometer_sample(self): return self.__sqm_radiometer_sample diff --git a/python/tests/test_api_extensions.py b/python/tests/test_api_extensions.py index 8b318daf7..179f2571f 100644 --- a/python/tests/test_api_extensions.py +++ b/python/tests/test_api_extensions.py @@ -96,3 +96,44 @@ def test_diagnostics_and_timing_keys_preserved(): assert d["solve_time"] == 1234.5 assert d["cam_solve_time"] == 1234.5 json.dumps(d, default=str) # full payload is JSON-serializable + + +@pytest.mark.unit +def test_raw_png_preserves_16bit_values(): + """A 12-bit sensor frame must survive the PNG round trip intact. + + Rendering uint16 via Image.fromarray(..., mode="L") reinterprets the + 16-bit buffer as 8-bit and yields interleaved-byte noise that still looks + like a plausible image -- convincing enough to waste a night of captures + before anyone checks the histogram. Pin the values. + """ + import io + + import numpy as np + from PIL import Image + + from PiFinder.api_extensions import _raw_to_png as to_png + + frame = np.array([[0, 1, 255, 256], [4095, 2048, 300, 65535]], dtype=np.uint16) + + buf = io.BytesIO() + to_png(frame).save(buf, format="PNG") + buf.seek(0) + restored = np.asarray(Image.open(buf)) + + np.testing.assert_array_equal(restored, frame) + + +@pytest.mark.unit +def test_full_raw_request_is_one_shot(): + """One request yields one frame; the camera must not keep publishing 4 MB.""" + from PiFinder.state import SharedStateObj + + state = SharedStateObj() + assert state.cam_raw_full_requested() is False + + state.request_cam_raw_full() + assert state.cam_raw_full_requested() is True + + state.set_cam_raw_full(object()) + assert state.cam_raw_full_requested() is False diff --git a/python/tests/test_full_frame_solve.py b/python/tests/test_full_frame_solve.py index 30942cd31..a79d67b92 100644 --- a/python/tests/test_full_frame_solve.py +++ b/python/tests/test_full_frame_solve.py @@ -17,7 +17,10 @@ from PiFinder.optics import DISPLAY_FRAME_SIZE, SolveGeometry from PiFinder.solver import project_solution_to_display -sys.path.append(str(utils.tetra3_dir)) +# tetra3 is a submodule on some branches and an ordinary dependency on +# others; only the former needs a path entry. +if hasattr(utils, "tetra3_dir"): + sys.path.append(str(utils.tetra3_dir)) tetra3 = pytest.importorskip("tetra3")