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/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_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..3024d1e0a 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,66 @@ def capture(self) -> Image.Image: 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 - # Apply camera-specific crop and rotation - raw_capture = self.profile.crop_and_rotate(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) + + 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 + # 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 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 +173,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()) - - # 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 + self.shared_state.set_cam_raw(cropped.copy()) - # rescale to 8 bit - raw_capture = ( - raw_capture - * 255 - / (2**self.profile.bit_depth - self.profile.bias_offset - 1) - ) + return Image.fromarray(self._to_8bit(cropped)).resize((512, 512)) - # clip to avoid <0 or >255 values - raw_capture = np.clip(raw_capture.astype(np.int32), 0, 255).astype(np.uint8) + 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))) - # 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 +345,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 +363,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..9c680e196 100644 --- a/python/PiFinder/state.py +++ b/python/PiFinder/state.py @@ -309,7 +309,13 @@ 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 + # 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 +386,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 @@ -563,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 new file mode 100644 index 000000000..a79d67b92 --- /dev/null +++ b/python/tests/test_full_frame_solve.py @@ -0,0 +1,193 @@ +"""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 + +# 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") + +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