diff --git a/default_config.json b/default_config.json index 7684bff0f..d2ed312c1 100644 --- a/default_config.json +++ b/default_config.json @@ -10,6 +10,9 @@ "screen_direction": "right", "mount_type": "Alt/Az", "solver_debug": 0, + "solver_shadow_detect": false, + "solver_sep_fallback": false, + "solver_sep_sigma": 4.0, "sleep_timeout": "30s", "screen_off_timeout": "Off", "chart_radec": "Off", diff --git a/python/PiFinder/camera_interface.py b/python/PiFinder/camera_interface.py index ac80cdc37..d64a435b4 100644 --- a/python/PiFinder/camera_interface.py +++ b/python/PiFinder/camera_interface.py @@ -150,6 +150,9 @@ class CameraInterface: _camera_started = False _save_next_to = None # Filename to save next capture to (None = don't save) + # Publish the uncropped raw frame for the SEP full-frame detection path. + # Set from config at loop start; False keeps captures copy-free. + _publish_solver_raw = False _auto_exposure_enabled = False _auto_exposure_mode = "pid" # "pid" or "snr" _auto_exposure_pid: Optional[ExposurePIDController] = None @@ -300,6 +303,17 @@ def get_image_loop( # Store shared_state for access by capture() methods self.shared_state = shared_state + # SEP full-frame detection path (shadow logging / fallback solve) + # needs the uncropped raw published per frame. Read once at start; + # changing these keys requires an app restart. Off by default: + # no per-frame copy is made unless the path is enabled. + self._publish_solver_raw = bool( + cfg.get_option("solver_shadow_detect") + or cfg.get_option("solver_sep_fallback") + ) + if self._publish_solver_raw: + logger.info("Publishing full-frame solver_raw (SEP path enabled)") + # Store camera type in shared state for SQM calibration camera_type_str = self.get_cam_type() # e.g., "PI imx296", "PI hq" if " " in camera_type_str: diff --git a/python/PiFinder/camera_pi.py b/python/PiFinder/camera_pi.py index 4364bcd2f..0cec8c9e6 100644 --- a/python/PiFinder/camera_pi.py +++ b/python/PiFinder/camera_pi.py @@ -109,6 +109,16 @@ def capture(self) -> Image.Image: _request.release() + # Uncropped frame for the SEP full-frame detection path (shadow / + # fallback). Same orientation convention as the cropped frame: + # profile rotation applied, crop skipped. The reference is free -- + # the manager proxy pickles (copies) on set_solver_raw below. + solver_full = None + if getattr(self, "_publish_solver_raw", False): + solver_full = raw_capture + if self.profile.rotation_90 != 0: + solver_full = np.rot90(solver_full, self.profile.rotation_90) + # Apply camera-specific crop and rotation raw_capture = self.profile.crop_and_rotate(raw_capture) @@ -134,6 +144,15 @@ 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()) + if solver_full is not None: + self.shared_state.set_solver_raw( + { + "frame": solver_full, + "timestamp": time.time(), + "exposure_us": metadata.get("ExposureTime"), + "gain": metadata.get("AnalogueGain"), + } + ) # covert to 32 bit int to avoid overflow raw_capture = raw_capture.astype(np.float32) diff --git a/python/PiFinder/sep_detect.py b/python/PiFinder/sep_detect.py new file mode 100644 index 000000000..adedadbcd --- /dev/null +++ b/python/PiFinder/sep_detect.py @@ -0,0 +1,295 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +SEP (Source Extractor) star detection on the full-sensor RAW frame. + +The production detector (cedar-detect) works on the processed 8-bit +512x512 solver frame, where the 12->8-bit stretch has already crushed +faint stars into a couple of levels under a bright (light-polluted) +sky background. This module detects in the 12-bit domain instead, on +the *uncropped* sensor frame: + +1. 2x2 mean binning, for SNR (x2) and PSF energy concentration; on + Bayer mosaics it also removes the RGGB modulation. +2. Estimate and subtract a mesh background (``sep.Background``) -- this + removes light-pollution gradients and cloud glow, which is exactly + the failure mode of a global threshold under a Seoul sky. +3. Extract sources against the local background RMS with a small + matched filter, then rank by flux. + +Returned centroids are in FULL-frame pixel coordinates (y, x), ready +for the solver-frame mapping in ``solver_frame_map``. + +``sep`` is an optional dependency: importing this module never fails, +and ``detect_stars`` returns None when sep is unavailable. +""" + +import logging +import time +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +logger = logging.getLogger("Solver.SepDetect") + +_sep = None +_sep_import_failed = False + +# 3x3 gaussian-ish matched filter (SExtractor convention) -- correlates +# neighbouring pixels so PSF-shaped bumps beat single-pixel noise. +MATCHED_FILTER = np.array( + [[1.0, 2.0, 1.0], [2.0, 4.0, 2.0], [1.0, 2.0, 1.0]], dtype=np.float32 +) + + +def _sep_module(): + """Import sep lazily; remember a failure so we log it only once.""" + global _sep, _sep_import_failed + if _sep is None and not _sep_import_failed: + try: + import sep + + _sep = sep + except ImportError: + _sep_import_failed = True + logger.warning("sep not installed; SEP detection disabled") + return _sep + + +@dataclass +class SepDetection: + """Result of one SEP extraction, in full-frame pixel coordinates.""" + + centroids: np.ndarray # (N, 2) float (y, x), flux-descending + fluxes: np.ndarray # (N,) float, same order + background_median: float # binned-domain ADU + background_rms: float # binned-domain ADU + elapsed_ms: float + # Otherwise-keepable detections removed by the warm-pixel map. High + # values on an empty sky are expected (the map is doing its job). + masked_count: int = 0 + + +def warm_pixel_excess(frame: np.ndarray) -> np.ndarray: + """Per-pixel excess over the median of the 4 distance-2 neighbours + (the same-Bayer-channel positions on a colour sensor; on a mono + sensor simply a sparse neighbourhood -- valid either way). + + A warm/hot pixel is a single-pixel spike, so its excess is its full + amplitude; extended structure (sky gradient, cloud, defocused + star) raises the neighbours too and mostly cancels. A tightly focused + star also shows excess -- which is why map *building* additionally + requires recurrence at a fixed position across frames (stars move + with the sky, warm pixels don't). + """ + arr = np.asarray(frame, dtype=np.float32) + h, w = arr.shape + p = np.pad(arr, 2, mode="edge") + neighbours = np.stack( + [ + p[0:h, 2 : w + 2], # same channel, y-2 + p[4 : h + 4, 2 : w + 2], # y+2 + p[2 : h + 2, 0:w], # x-2 + p[2 : h + 2, 4 : w + 4], # x+2 + ] + ) + return arr - np.median(neighbours, axis=0) + + +def build_warm_pixel_map( + frames, + min_excess_adu: float = 45.0, + min_recurrence: float = 0.7, +) -> np.ndarray: + """Warm-pixel positions recurring across frames, as (N, 2) int (y, x). + + Args: + frames: Iterable of 2D raw arrays, all the same shape and in the + same orientation the map will be applied in (solver_raw + orientation: profile rot90 applied, no crop -- stage dumps are + already in this orientation). + min_excess_adu: Same-channel neighbour excess for a candidate. + min_recurrence: Fraction of frames a position must be a candidate + in. Static defects recur near 1.0; stars drift out within one + frame interval, single-frame noise almost never repeats. + """ + counts: Optional[np.ndarray] = None + n_frames = 0 + for frame in frames: + candidate = warm_pixel_excess(frame) > min_excess_adu + if counts is None: + counts = np.zeros(candidate.shape, dtype=np.uint16) + counts += candidate + n_frames += 1 + if counts is None or n_frames == 0: + return np.empty((0, 2), dtype=np.int32) + needed = max(1, int(np.ceil(min_recurrence * n_frames))) + ys, xs = np.nonzero(counts >= needed) + return np.column_stack((ys, xs)).astype(np.int32) + + +def bin2x2(frame: np.ndarray) -> np.ndarray: + """Mean-bin each 2x2 block; trims odd edges.""" + arr = np.asarray(frame) + h, w = arr.shape[0] // 2 * 2, arr.shape[1] // 2 * 2 + arr = arr[:h, :w].astype(np.float32) + return ( + arr[0::2, 0::2] + arr[0::2, 1::2] + arr[1::2, 0::2] + arr[1::2, 1::2] + ) * 0.25 + + +def detect_stars( + raw_frame: np.ndarray, + sigma: float = 3.5, + minarea: int = 3, + max_stars: int = 48, + edge_margin_px: int = 48, + saturation_level: Optional[float] = None, + warm_pixel_map: Optional[np.ndarray] = None, + warm_pixel_radius_px: float = 4.0, + max_semimajor_px: float = 2.0, + max_npix: int = 40, + cluster_radius_px: float = 50.0, + cluster_max_neighbors: int = 1, +) -> Optional[SepDetection]: + """ + Detect stars on a raw sensor frame (uint16 mosaic, any shape). + + Field lesson (Seoul, 2026-07-28 night): the uncropped frame's + vignetted borders under cloud glow produce dozens of spurious + extractions -- on a saturated-interior frame ALL "detections" sat at + the frame edge with junk fluxes (huge blob, zeros, negatives). Hence + the three quality filters here: an edge margin, a positive-flux + requirement, and a saturation guard that reports an honest zero when + the sky has burned the interior flat. + + Args: + raw_frame: 2D raw sensor array (Bayer mosaic or mono). + sigma: Extraction threshold in units of the local background RMS. + minarea: Minimum connected pixels above threshold. + max_stars: Keep at most this many, brightest (by flux) first. + edge_margin_px: Drop detections within this many full-res pixels + of the frame border (vignette / background-mesh edge zone). + saturation_level: Sensor full scale (e.g. 4095 for 12-bit). When + given and the binned interior median is at it, return zero + detections instead of edge noise. + warm_pixel_map: (N, 2) int (y, x) static-defect positions from + ``build_warm_pixel_map``, same orientation as ``raw_frame``. + Detections within ``warm_pixel_radius_px`` of a mapped position + are dropped and counted in ``masked_count``. + warm_pixel_radius_px: Match radius in full-res pixels (binning + quantises centroids to a 2 px grid, so keep this >= 4). + max_semimajor_px: Reject sources with a larger fitted semi-major + axis (binned px). Cloud-edge texture is extended; real stars + measured a <= 0.86 (p95), junk past 1.5 and NaN (degenerate + fits, also rejected). Default 2.0 leaves defocus headroom. + 2026-07-28 night corpus. + max_npix: Reject sources covering more binned pixels (stars p95 + 10, cloud blobs to 188 -- same corpus; 40 = defocus headroom). + cluster_radius_px / cluster_max_neighbors: Drop detections with + more than ``cluster_max_neighbors`` others within the radius + (full-res px). SEP deblends a bright cloud edge into tight + clumps; measured real (tetra3-matched) stars had ZERO + neighbours within 50 px in every case, junk up to 4. + + Returns: + SepDetection with centroids in full-frame (y, x) pixels, or None + if sep is unavailable or the frame is unusable. + """ + sep = _sep_module() + if sep is None: + return None + arr = np.asarray(raw_frame) + if arr.ndim != 2 or arr.shape[0] < 8 or arr.shape[1] < 8: + return None + + t0 = time.perf_counter() + binned = bin2x2(arr) + # sep requires C-contiguous native-endian float32 + data = np.ascontiguousarray(binned, dtype=np.float32) + bkg = sep.Background(data, bw=32, bh=32) + + def _empty() -> SepDetection: + return SepDetection( + centroids=np.empty((0, 2)), + fluxes=np.empty(0), + background_median=float(bkg.globalback), + background_rms=float(bkg.globalrms), + elapsed_ms=(time.perf_counter() - t0) * 1000.0, + ) + + if saturation_level is not None: + h2, w2 = data.shape + interior = data[h2 // 4 : -h2 // 4 or None, w2 // 4 : -w2 // 4 or None] + if np.median(interior) >= 0.98 * saturation_level: + return _empty() + + data_sub = data - bkg.back() + objects = sep.extract( + data_sub, + thresh=sigma, + err=bkg.rms(), + filter_kernel=MATCHED_FILTER, + minarea=minarea, + ) + + # A binned pixel (i, j) covers full-res pixels (2i, 2i+1) x (2j, 2j+1), + # so its centre sits at 2*coord + 0.5 in full-frame coordinates. + full_y = np.asarray(objects["y"]) * 2.0 + 0.5 + full_x = np.asarray(objects["x"]) * 2.0 + 0.5 + fluxes = np.asarray(objects["flux"], dtype=np.float64) + semimajor = np.asarray(objects["a"], dtype=np.float64) + npix = np.asarray(objects["npix"], dtype=np.int64) + + h, w = arr.shape + keep = ( + (fluxes > 0) + & (full_y >= edge_margin_px) + & (full_y < h - edge_margin_px) + & (full_x >= edge_margin_px) + & (full_x < w - edge_margin_px) + # Point-source shape gate: extended or degenerate (NaN) fits are + # cloud texture, not stars (thresholds measured -- see docstring). + & np.isfinite(semimajor) + & (semimajor <= max_semimajor_px) + & (npix <= max_npix) + ) + + # Warm-pixel map: drop otherwise-keepable detections sitting on a known + # static defect (before the top-N cap, so defects can't crowd out stars). + masked_count = 0 + if warm_pixel_map is not None and len(warm_pixel_map) and keep.any(): + wp = np.asarray(warm_pixel_map, dtype=np.float64) + d2 = ( + (full_y[:, None] - wp[None, :, 0]) ** 2 + + (full_x[:, None] - wp[None, :, 1]) ** 2 + ).min(axis=1) + warm = d2 <= warm_pixel_radius_px**2 + masked_count = int((keep & warm).sum()) + keep &= ~warm + + full_y, full_x, fluxes = full_y[keep], full_x[keep], fluxes[keep] + + # Cluster gate: SEP deblends bright cloud edges into tight clumps of + # "sources"; real stars at this plate scale are isolated (measured 0 + # neighbours within 50 px on every tetra3-matched star). + if len(full_y) > 1: + d2 = (full_y[:, None] - full_y[None, :]) ** 2 + ( + full_x[:, None] - full_x[None, :] + ) ** 2 + neighbours = (d2 <= cluster_radius_px**2).sum(axis=1) - 1 + isolated = neighbours <= cluster_max_neighbors + full_y, full_x, fluxes = full_y[isolated], full_x[isolated], fluxes[isolated] + + elapsed_ms = (time.perf_counter() - t0) * 1000.0 + + order = np.argsort(fluxes)[::-1][:max_stars] + return SepDetection( + centroids=np.column_stack((full_y[order], full_x[order])), + fluxes=fluxes[order], + background_median=float(bkg.globalback), + background_rms=float(bkg.globalrms), + elapsed_ms=elapsed_ms, + masked_count=masked_count, + ) diff --git a/python/PiFinder/sep_shadow.py b/python/PiFinder/sep_shadow.py new file mode 100644 index 000000000..29653a620 --- /dev/null +++ b/python/PiFinder/sep_shadow.py @@ -0,0 +1,331 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +Shadow / fallback runner for the SEP full-frame detection path. + +Purpose: evaluate the SEP candidate detector (12-bit uncropped +detection, ``sep_detect``) against the production cedar-detect path in a +single field session: + +* **Shadow**: on every solve attempt, also run SEP on the uncropped raw + frame and append one CSV row comparing both detectors. Zero effect on + the production solve. +* **Fallback** (opt-in on top of shadow data): when the production + solve fails and SEP found enough stars, attempt a real solve from the + SEP centroids in the rotated full frame -- the solution feeds the + normal pointing chain, so tracking works from it. Guarded so an + in-progress alignment never runs through this path (its y/x_target + would be in full-frame space). + +All entry points are defensive: any exception is logged and swallowed, +so the experiment can never take down the production solver. + +Config keys (restart to apply): ``solver_shadow_detect``, +``solver_sep_fallback``, ``solver_sep_sigma``. + +CSV: ``solver_shadow_log.csv`` in PiFinder_data (one row per solve +attempt while shadow logging is enabled). +""" + +import csv +import logging +import time +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from PiFinder import sep_detect, utils +from PiFinder import solver_frame_map as sfm +from PiFinder.sep_detect import SepDetection +from PiFinder.sqm.camera_profiles import get_camera_profile + +logger = logging.getLogger("Solver.SepShadow") + +CSV_FIELDS = [ + "timestamp", + "exposure_us", + "gain", + "cedar_centroids", + "matches", + "solved", + "sep_centroids", + "sep_top_flux", + "sep_bkg", + "sep_rms", + "sep_ms", + "fallback_used", + "fallback_rmse", + "sep_masked", +] + +# A solver_raw older than this no longer matches the attempt being logged +# (camera wedged, SEP path disabled mid-run); skip rather than mislabel. +MAX_FRAME_AGE_S = 15.0 + +# Warm-pixel map: (N, 2) int (y, x) in solver_raw orientation, built by +# ``python -m PiFinder.sep_warm_map`` from raw-frame corpora. Optional -- +# missing file just means no masking. Regenerate when the sensor ages or +# after long temperature shifts (warm pixels grow over both). +WARM_MAP_PATH = utils.data_dir / "sep_warm_pixels.npy" + + +@dataclass +class SepRun: + """One SEP pass over the freshest full-frame raw.""" + + detection: SepDetection + frame_hw: tuple + exposure_us: Optional[float] + gain: Optional[float] + + +class SepShadowRunner: + def __init__( + self, + shadow_enabled: bool, + fallback_enabled: bool, + sigma: float, + rotation_deg: float, + crop_width_px: int, + # 5, paired with sigma ~4: the 2026-07-28 live-sky sweep showed + # half the genuine rescues carry only 5-7 detections at that + # threshold (a gate of 8 was calibrated for a lower sigma's + # junk-inflated counts). No observed solve had fewer than 5. + min_fallback_stars: int = 5, + saturation_level: Optional[float] = None, + csv_path=None, + warm_pixel_map: Optional[np.ndarray] = None, + ): + self.shadow_enabled = shadow_enabled + self.fallback_enabled = fallback_enabled + self.sigma = sigma + self.rotation_deg = rotation_deg + self.crop_width_px = crop_width_px + self.min_fallback_stars = min_fallback_stars + self.saturation_level = saturation_level + self.csv_path = csv_path or (utils.data_dir / "solver_shadow_log.csv") + self.warm_pixel_map = warm_pixel_map + # Fallback backoff state (see fallback_should_attempt): a fallback + # solve on unsolvable input burns up to solve_timeout (1 s) of solver + # CPU per attempt -- indoors/under cloud that is every attempt. + self._attempt_counter = 0 + self._fallback_fail_streak = 0 + self._fallback_skip_until = 0 + self._last_failed_sep_count: Optional[int] = None + logger.info( + "SEP shadow runner: shadow=%s fallback=%s sigma=%.1f " + "rotation=%.0f° crop_width=%dpx warm_pixels=%d log=%s", + shadow_enabled, + fallback_enabled, + sigma, + rotation_deg, + crop_width_px, + 0 if warm_pixel_map is None else len(warm_pixel_map), + self.csv_path, + ) + + @classmethod + def create_if_enabled(cls, cfg, camera_type: Optional[str]): + """Build a runner from config, or None when the path is disabled + or the camera profile (crop geometry) is not resolvable yet.""" + try: + shadow = bool(cfg.get_option("solver_shadow_detect")) + fallback = bool(cfg.get_option("solver_sep_fallback")) + if not (shadow or fallback): + return None + if not camera_type: + return None + profile = get_camera_profile(camera_type) + crop_width = int( + profile.raw_size[0] - profile.crop_x[0] - profile.crop_x[1] + ) + rotation = sfm.stage5_rotation_deg( + cfg.get_option("screen_direction"), + cfg.get_option("camera_rotation"), + ) + sigma = float(cfg.get_option("solver_sep_sigma") or 4.0) + warm_map = None + try: + if WARM_MAP_PATH.exists(): + warm_map = np.asarray(np.load(WARM_MAP_PATH), dtype=np.int32) + logger.info( + "Loaded %d warm pixels from %s", len(warm_map), WARM_MAP_PATH + ) + except Exception: + logger.exception("Warm-pixel map load failed; continuing unmasked") + warm_map = None + return cls( + shadow_enabled=shadow, + fallback_enabled=fallback, + sigma=sigma, + rotation_deg=rotation, + crop_width_px=crop_width, + saturation_level=float(2**profile.bit_depth - 1), + warm_pixel_map=warm_map, + ) + except Exception: + logger.exception("SEP shadow runner init failed; disabled") + return None + + def fallback_should_attempt(self, sep_count: int) -> bool: + """Backoff gate for the fallback solve. + + A failed fallback solve costs up to solve_timeout (1 s) of solver + CPU. When the scene is persistently unsolvable (indoors, thick + cloud) the SEP count passes the star gate every attempt and that + cost recurs forever. After each consecutive failure we skip the + next ``min(2**streak, 8)`` attempts -- but re-arm IMMEDIATELY when + the SEP count rises to 1.5x the last failed attempt, which is what + a cloud gap opening on real stars looks like. Rescue solves in a + star window are therefore not delayed (2026-07-27 field: counts + jumped from <=5 masked to ~30 when stars appeared). + """ + if self._fallback_fail_streak == 0: + return True + if ( + self._last_failed_sep_count is not None + and sep_count >= 1.5 * self._last_failed_sep_count + ): + return True + return self._attempt_counter >= self._fallback_skip_until + + def record_fallback_result(self, solved: bool, sep_count: int) -> None: + if solved: + self._fallback_fail_streak = 0 + self._last_failed_sep_count = None + return + self._fallback_fail_streak += 1 + self._last_failed_sep_count = sep_count + self._fallback_skip_until = self._attempt_counter + min( + 2**self._fallback_fail_streak, 8 + ) + + def note_solved(self) -> None: + """A production solve succeeded: the sky is workable, so the next + cedar failure deserves an immediate fallback try again.""" + self._fallback_fail_streak = 0 + self._last_failed_sep_count = None + + def detect(self, shared_state) -> Optional[SepRun]: + """Run SEP on the freshest published full-frame raw, or None.""" + self._attempt_counter += 1 + try: + entry = shared_state.solver_raw() + if not entry or "frame" not in entry: + return None + if time.time() - float(entry.get("timestamp") or 0) > MAX_FRAME_AGE_S: + return None + frame = np.asarray(entry["frame"]) + detection = sep_detect.detect_stars( + frame, + sigma=self.sigma, + saturation_level=self.saturation_level, + warm_pixel_map=self.warm_pixel_map, + ) + if detection is None: + return None + return SepRun( + detection=detection, + frame_hw=(frame.shape[0], frame.shape[1]), + exposure_us=entry.get("exposure_us"), + gain=entry.get("gain"), + ) + except Exception: + logger.exception("SEP shadow detect failed") + return None + + def solve(self, t3, run: SepRun, shared_state) -> Optional[dict]: + """Solve from SEP centroids in the rotated full frame. + + Rotation, canvas size, fov and target_pixel are all mapped so the + resulting RA/Dec/Roll -- and the aligned pointing at target_pixel + -- carry the exact semantics of the production 512-frame solve + (see solver_frame_map). + """ + try: + cents, canvas = sfm.rotate_centroids( + run.detection.centroids, run.frame_hw, self.rotation_deg + ) + target_pixel = sfm.map_target_pixel_to_frame( + shared_state.target_pixel(), canvas, self.crop_width_px + ) + fov = sfm.fov_estimate_deg(canvas[1], self.crop_width_px) + return t3.solve_from_centroids( + cents, + canvas, + fov_estimate=fov, + fov_max_error=fov / 3.0, + match_max_error=0.005, + return_matches=True, + target_pixel=target_pixel, + solve_timeout=1000, + ) + except Exception: + logger.exception("SEP fallback solve failed") + return None + + def _rotate_csv_on_schema_change(self) -> None: + """Sideline a CSV written with an older field list, once. + + Mixed-width rows break offline analysis; the sidelined file keeps + its data under ``.old``. + """ + if getattr(self, "_csv_schema_checked", False): + return + self._csv_schema_checked = True + if not self.csv_path.exists(): + return + with open(self.csv_path, newline="") as f: + header = f.readline().strip() + if header != ",".join(CSV_FIELDS): + old = self.csv_path.with_suffix(self.csv_path.suffix + ".old") + self.csv_path.replace(old) + logger.info("Shadow CSV schema changed; previous file moved to %s", old) + + def log_attempt( + self, + exposure_us, + gain, + cedar_count: int, + matches, + solved: bool, + run: Optional[SepRun], + fallback_used: bool = False, + fallback_rmse=None, + ) -> None: + """Append one attempt-comparison row; never raises.""" + if not self.shadow_enabled: + return + try: + row = { + "timestamp": f"{time.time():.3f}", + "exposure_us": exposure_us, + "gain": gain, + "cedar_centroids": cedar_count, + "matches": matches, + "solved": int(bool(solved)), + "sep_centroids": len(run.detection.centroids) if run else "", + "sep_top_flux": ( + f"{run.detection.fluxes[0]:.0f}" + if run and len(run.detection.fluxes) + else "" + ), + "sep_bkg": f"{run.detection.background_median:.1f}" if run else "", + "sep_rms": f"{run.detection.background_rms:.2f}" if run else "", + "sep_ms": f"{run.detection.elapsed_ms:.0f}" if run else "", + "fallback_used": int(bool(fallback_used)), + "fallback_rmse": ( + f"{fallback_rmse:.2f}" if fallback_rmse is not None else "" + ), + "sep_masked": run.detection.masked_count if run else "", + } + self._rotate_csv_on_schema_change() + write_header = not self.csv_path.exists() + with open(self.csv_path, "a", newline="") as f: + writer = csv.DictWriter(f, fieldnames=CSV_FIELDS) + if write_header: + writer.writeheader() + writer.writerow(row) + except Exception: + logger.exception("SEP shadow CSV append failed") diff --git a/python/PiFinder/sep_warm_map.py b/python/PiFinder/sep_warm_map.py new file mode 100644 index 000000000..12f1d6a7f --- /dev/null +++ b/python/PiFinder/sep_warm_map.py @@ -0,0 +1,139 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +Build the SEP warm-pixel map from stage-dump corpora. + +Warm pixels are static single-pixel sensor defects that dominate SEP's +detection counts on empty sky (field bench 2026-07-28: 19 positions +accounted for 55% of all detections across a night of frames). This +tool finds them by same-channel neighbour excess recurring at a fixed +position across many frames, and writes the map ``sep_shadow`` loads at +startup. + +Usage (device, any mix of directories holding raw-frame dumps):: + + python -m PiFinder.sep_warm_map ~/PiFinder_data/captures + python -m PiFinder.sep_warm_map ~/frames --dry-run + +Inputs are 16-bit raw rasters found anywhere under the given dirs, in +solver_raw orientation (profile rot90 applied): + +* ``*_raw_full.png`` -- uncropped frames, used as-is. +* ``*_raw_cropped.png`` -- crop-window frames; positions are offset back + into full-frame coordinates via the camera profile. Only valid for + profiles without a rotation (crop happens before rot90), so cropped + input is refused when ``rotation_90 != 0``. + +The two corpora complement each other: cropped dumps exist on every +solve-failure streak, full-frame dumps cover the vignette border region +outside the crop window. + +Feed DARK corpora (night / lens-cap / dim indoor). Warm pixels are only +detectable against a low background; bright twilight frames dilute the +per-group recurrence fraction and drop legitimate positions (observed +2026-07-28: adding a bright evening's dumps shrank the map 57 -> 40 and +lost census-validated defects). +""" + +import argparse +import logging +from collections import defaultdict +from pathlib import Path + +import numpy as np +from PIL import Image + +from PiFinder import utils +from PiFinder.sep_detect import build_warm_pixel_map +from PiFinder.sqm.camera_profiles import get_camera_profile + +logger = logging.getLogger("Solver.SepWarmMap") + +DEFAULT_OUT = utils.data_dir / "sep_warm_pixels.npy" + + +def _load_frames(dump_dirs): + """Yield (kind, array) for every raw stage raster under the dirs.""" + for root in dump_dirs: + for pattern, kind in ( + ("*_raw_full.png", "full"), + ("*_raw_cropped.png", "cropped"), + ): + for path in sorted(Path(root).rglob(pattern)): + yield kind, np.asarray(Image.open(path)).astype(np.uint16) + + +def build_map( + dump_dirs, + camera_type: str = "imx462", + min_excess_adu: float = 45.0, + min_recurrence: float = 0.7, +) -> np.ndarray: + """Aggregate all dumps into one full-frame warm-pixel map.""" + profile = get_camera_profile(camera_type) + # Group frames by kind+shape: recurrence is only meaningful within a + # group of identically-framed captures. + groups = defaultdict(list) + for kind, frame in _load_frames(dump_dirs): + if frame.ndim != 2: + continue + if kind == "cropped" and profile.rotation_90 != 0: + raise SystemExit( + f"cropped dumps can't be mapped for {camera_type} " + "(crop precedes rot90); use raw_full dumps" + ) + groups[(kind, frame.shape)].append(frame) + + positions = [] + for (kind, shape), frames in sorted(groups.items()): + pts = build_warm_pixel_map( + frames, min_excess_adu=min_excess_adu, min_recurrence=min_recurrence + ) + if kind == "cropped": + pts = pts + np.array([profile.crop_y[0], profile.crop_x[0]], dtype=np.int32) + logger.info( + "%s %s: %d frames -> %d warm pixels", kind, shape, len(frames), len(pts) + ) + positions.append(pts) + + if not positions: + return np.empty((0, 2), dtype=np.int32) + merged = np.vstack(positions) + # Dedupe positions found by both corpora (2 px grid, well inside the + # 4 px match radius detect_stars uses). + merged = np.unique(merged // 2, axis=0) * 2 + return merged.astype(np.int32) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__.split("\n")[1]) + parser.add_argument("dump_dirs", nargs="+", help="dirs containing stages_* dumps") + parser.add_argument("--camera", default="imx462", help="camera profile name") + parser.add_argument("--out", type=Path, default=DEFAULT_OUT) + # Defaults validated on the 2026-07-27 night corpus: a 55-position map + # covered all 19 recurring SEP-census cells and cut empty-sky counts + # 25.0 -> 5.5 mean, while masking only ~0.14% of the frame area. + parser.add_argument("--min-excess-adu", type=float, default=45.0) + parser.add_argument("--min-recurrence", type=float, default=0.7) + parser.add_argument( + "--dry-run", action="store_true", help="report only, write nothing" + ) + args = parser.parse_args() + logging.basicConfig(level=logging.INFO, format="%(message)s") + + warm = build_map( + args.dump_dirs, + camera_type=args.camera, + min_excess_adu=args.min_excess_adu, + min_recurrence=args.min_recurrence, + ) + print(f"{len(warm)} warm pixels total") + if args.dry_run: + return + args.out.parent.mkdir(parents=True, exist_ok=True) + np.save(args.out, warm) + print(f"wrote {args.out}") + + +if __name__ == "__main__": + main() diff --git a/python/PiFinder/solver.py b/python/PiFinder/solver.py index f6fc17d88..ffb2d47ab 100644 --- a/python/PiFinder/solver.py +++ b/python/PiFinder/solver.py @@ -24,9 +24,11 @@ from multiprocessing import shared_memory import grpc +from PiFinder import config as config_mod from PiFinder import state_utils from PiFinder import utils from PiFinder import timez +from PiFinder.sep_shadow import SepShadowRunner from PiFinder.sqm import SQM as SQMCalculator from PiFinder.sqm.wings import WingEstimator @@ -861,6 +863,16 @@ def solver( sqm_radiometer = RadiometerAccumulator() last_stellar_diagnostic = 0.0 + # SEP shadow/fallback runner (full-frame detection path). Needs the + # camera type for crop geometry, which the camera process publishes + # after startup -- so creation is retried in the loop until it works. + _sep_cfg = config_mod.Config() + sep_shadow = None + sep_shadow_wanted = bool( + _sep_cfg.get_option("solver_shadow_detect") + or _sep_cfg.get_option("solver_sep_fallback") + ) + while True: logger.info("Starting Solver Loop") # Try to start cedar detect server, fall back to tetra3 centroider if unavailable @@ -1010,6 +1022,64 @@ def solver( **_solver_args, ) + # SEP full-frame path: shadow-detect on every attempt; + # optionally rescue a failed production solve from the + # SEP centroids (sep_shadow module docstring). + if sep_shadow is None and sep_shadow_wanted: + sep_shadow = SepShadowRunner.create_if_enabled( + _sep_cfg, shared_state.camera_type() + ) + sep_run = None + sep_fallback_used = False + if sep_shadow is not None: + sep_run = sep_shadow.detect(shared_state) + if ( + sep_run is not None + and sep_shadow.fallback_enabled + and (not solution or solution.get("RA") is None) + # An in-progress alignment must resolve through + # the production frame: this path cannot answer + # the alignment coordinate, so skip the rescue + # rather than let a success silently consume the + # pending alignment request. + and align_ra == 0 + and align_dec == 0 + and len(sep_run.detection.centroids) + >= sep_shadow.min_fallback_stars + # Backoff: persistently unsolvable scenes + # (indoors, thick cloud) otherwise burn up to + # solve_timeout per attempt, starving the whole + # solver loop. Re-arms instantly on a SEP count + # jump (cloud gap opening on stars). + and sep_shadow.fallback_should_attempt( + len(sep_run.detection.centroids) + ) + ): + fb_solution = sep_shadow.solve(t3, sep_run, shared_state) + sep_shadow.record_fallback_result( + bool(fb_solution and fb_solution.get("RA") is not None), + len(sep_run.detection.centroids), + ) + if fb_solution and fb_solution.get("RA") is not None: + # Per-centroid outputs are in full-frame + # coordinates; strip them so SQM photometry + # (which reads the 512 frame) never mixes + # coordinate spaces. + fb_solution.pop("matched_centroids", None) + fb_solution.pop("matched_stars", None) + # catID parallels the stripped arrays; keep the + # matched-* trio consistent on the message. + fb_solution.pop("matched_catID", None) + solution = fb_solution + sep_fallback_used = True + logger.debug( + "SEP fallback solve SUCCESS - %d SEP " + "centroids (cedar saw %d), RMSE %.1f", + len(sep_run.detection.centroids), + len(centroids), + solution.get("RMSE") or -1.0, + ) + if "matched_centroids" in solution: if sqm_calculator is None: sqm_calculator = create_sqm_calculator(shared_state) @@ -1060,6 +1130,11 @@ def solver( if solution and solution.get("RA") is not None: last_solve_success = last_solve_attempt + if sep_shadow is not None and not sep_fallback_used: + # Production solve: sky is workable, clear the + # fallback backoff so the next cedar failure gets + # an immediate rescue try. + sep_shadow.note_solved() solve_result = _build_successful_solve( solution=solution, last_image_metadata=last_image_metadata, @@ -1113,6 +1188,22 @@ def solver( t_extract_ms=t_extract, ) ) + + if sep_shadow is not None: + sep_shadow.log_attempt( + exposure_us=last_image_metadata.get("exposure_time"), + gain=last_image_metadata.get("gain"), + cedar_count=len(centroids), + matches=solution.get("Matches") if solution else None, + solved=bool(solution and solution.get("RA") is not None), + run=sep_run, + fallback_used=sep_fallback_used, + fallback_rmse=( + solution.get("RMSE") + if sep_fallback_used and solution + else None + ), + ) except Exception as e: logger.error( f"Exception during solve attempt: {e.__class__.__name__}: {str(e)}" diff --git a/python/PiFinder/solver_frame_map.py b/python/PiFinder/solver_frame_map.py new file mode 100644 index 000000000..e0198e033 --- /dev/null +++ b/python/PiFinder/solver_frame_map.py @@ -0,0 +1,129 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +Coordinate mapping between the production solver frame and the +full-sensor frame used by SEP detection. + +The production solver frame is the 512x512 processed image: the sensor +frame cropped to a centred square, resized, then rotated per +screen_direction / camera_rotation (camera_interface stage 5). Plate +solutions -- including ``target_pixel``, the persisted eyepiece +alignment point -- live in that space. + +The SEP path detects on the *uncropped* raw frame. To keep RA/Dec/Roll +and the alignment semantics identical, its solve runs in a frame with +the SAME rotation applied ("the rotated full frame"). Because the crop +is centred and the resize is isotropic, mapping the alignment point +between the two rotated frames reduces to a scale about the frame +centre -- the rotations cancel (proof: both rotations are about their +frame centres and the optical centre coincides with both). + +Rotation conventions are pinned by tests against PIL's ``Image.rotate`` +(counterclockwise, expand=False), which is what stage 5 uses. +""" + +import math +from typing import Tuple + +import numpy as np + +# The production solve calls tetra3 with (512, 512) and fov_estimate 12.0 +# (solver.py). That makes the plate scale 12 deg across the cropped +# square, whatever the crop width in sensor pixels. +SOLVER_FRAME_PX = 512 +SOLVER_FOV_DEG = 12.0 + + +def stage5_rotation_deg(screen_direction, camera_rotation) -> float: + """The rotation camera_interface stage 5 applies, in PIL CCW degrees. + + Reads ``SCREEN_ROTATE_AMOUNTS`` from camera_interface (the single + source of the per-variant rotation) via a lazy import so this module + stays cheap to import in tests and offline tools. + """ + if camera_rotation is not None: + return (-int(camera_rotation)) % 360 + from PiFinder.camera_interface import SCREEN_ROTATE_AMOUNTS + + return float(SCREEN_ROTATE_AMOUNTS.get(screen_direction, 270)) + + +def rotate_centroids( + centroids: np.ndarray, frame_hw: Tuple[int, int], angle_deg: float +) -> Tuple[np.ndarray, Tuple[int, int]]: + """ + Rotate (y, x) centroids the way stage 5 rotates the image (CCW). + + Quarter turns use the exact integer mapping with the canvas dims + swapped (np.rot90-style, no pixels lost). Other angles rotate about + the canvas centre with the canvas size unchanged (PIL expand=False + semantics; verified empirically against Image.rotate). + + Returns (rotated centroids, rotated canvas (h, w)). + """ + cents = np.asarray(centroids, dtype=np.float64).reshape(-1, 2) + h, w = frame_hw + angle = angle_deg % 360 + if angle == 0: + return cents.copy(), (h, w) + + if angle % 90 == 0: + out = cents.copy() + hh, ww = h, w + for _ in range(int(angle // 90) % 4): + y, x = out[:, 0], out[:, 1] + out = np.column_stack((ww - 1 - x, y)) + hh, ww = ww, hh + return out, (hh, ww) + + cy, cx = (h - 1) / 2.0, (w - 1) / 2.0 + rad = math.radians(angle) + cos_a, sin_a = math.cos(rad), math.sin(rad) + dy = cents[:, 0] - cy + dx = cents[:, 1] - cx + new_x = cx + dx * cos_a + dy * sin_a + new_y = cy - dx * sin_a + dy * cos_a + return np.column_stack((new_y, new_x)), (h, w) + + +def map_target_pixel_to_frame( + target_pixel_yx, frame_hw: Tuple[int, int], crop_width_px: int +) -> Tuple[float, float]: + """ + Map ``target_pixel`` (stored in rotated-512 space) into a rotated + full-frame canvas of size ``frame_hw``. + + Scale about the centre by crop_width/512; the rotations cancel (see + module docstring). ``crop_width_px`` is the cropped square's width + in sensor pixels (e.g. 980 for imx462). + """ + scale = crop_width_px / float(SOLVER_FRAME_PX) + c512 = (SOLVER_FRAME_PX - 1) / 2.0 + cy, cx = (frame_hw[0] - 1) / 2.0, (frame_hw[1] - 1) / 2.0 + ty, tx = float(target_pixel_yx[0]), float(target_pixel_yx[1]) + return (cy + (ty - c512) * scale, cx + (tx - c512) * scale) + + +def map_frame_pixel_to_target( + pixel_yx, frame_hw: Tuple[int, int], crop_width_px: int +) -> Tuple[float, float]: + """Inverse of :func:`map_target_pixel_to_frame`: a pixel in the rotated + full-frame canvas back into rotated-512 ``target_pixel`` space. + + Used by the SEP-path alignment: tetra3 returns the alignment target's + y/x in the full-frame canvas, but the production chain stores and + consumes target pixels in 512 space, so the result must come back + through the same centre-scale relation (same proof as the forward + mapping -- centre-symmetric crop, isotropic resize). + """ + scale = SOLVER_FRAME_PX / float(crop_width_px) + c512 = (SOLVER_FRAME_PX - 1) / 2.0 + cy, cx = (frame_hw[0] - 1) / 2.0, (frame_hw[1] - 1) / 2.0 + py, px = float(pixel_yx[0]), float(pixel_yx[1]) + return (c512 + (py - cy) * scale, c512 + (px - cx) * scale) + + +def fov_estimate_deg(frame_width_px: int, crop_width_px: int) -> float: + """FOV across ``frame_width_px`` sensor pixels, from the production + calibration of SOLVER_FOV_DEG across the cropped square.""" + return SOLVER_FOV_DEG * frame_width_px / float(crop_width_px) diff --git a/python/PiFinder/state.py b/python/PiFinder/state.py index a88880303..0c89affb5 100644 --- a/python/PiFinder/state.py +++ b/python/PiFinder/state.py @@ -309,6 +309,10 @@ 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 + # Uncropped raw sensor frame for the SEP full-frame detection path + # (dict with "frame" uint16 ndarray + capture metadata). Only + # published while solver_shadow_detect / solver_sep_fallback is on. + self.__solver_raw = None self.__sqm_radiometer_sample = None # Are we prepared to do alt/az math # We need gps lock and datetime @@ -571,6 +575,13 @@ def cam_raw(self): def set_cam_raw(self, v): self.__cam_raw = v + def solver_raw(self): + """Uncropped raw frame entry for the SEP full-frame path, or None.""" + return self.__solver_raw + + def set_solver_raw(self, v): + self.__solver_raw = v + def sqm_radiometer_sample(self): return self.__sqm_radiometer_sample diff --git a/python/pyproject.toml b/python/pyproject.toml index ad0f02e62..c1bd5311b 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -141,6 +141,7 @@ module = [ 'picamera2', 'bottle', 'libinput', + 'sep', ] ignore_missing_imports = true ignore_errors = true diff --git a/python/requirements.txt b/python/requirements.txt index 796057bdf..41b21b41e 100644 --- a/python/requirements.txt +++ b/python/requirements.txt @@ -11,6 +11,7 @@ json5==0.9.25 luma.oled==3.12.0 luma.lcd==2.11.0 numpy==1.26.4 +sep==1.4.1 numpy-quaternion==2023.0.4 pam==0.2.0 pandas==2.0.3 diff --git a/python/tests/test_sep_detect.py b/python/tests/test_sep_detect.py new file mode 100644 index 000000000..1601abeca --- /dev/null +++ b/python/tests/test_sep_detect.py @@ -0,0 +1,216 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +Unit tests for the SEP full-frame detection path and its coordinate +mapping into the production solver frame. + +The rotation conventions are pinned against PIL's Image.rotate (what +camera_interface stage 5 actually uses) so the SEP solve produces the +same Roll / target-pixel semantics as the production path. +""" + +import numpy as np +import pytest +from PIL import Image + +from PiFinder import solver_frame_map as sfm + +sep = pytest.importorskip("sep") + +from PiFinder import sep_detect # noqa: E402 + + +def _synthetic_frame(stars, shape=(540, 960), bg=1200.0, peak=400.0): + """Raw-like uint16 mosaic with a gradient background and gaussian stars. + + The checkerboard gain below is a worst-case robustness input: a + mono sensor has no phase response, so field frames are easier than + this fixture. + """ + h, w = shape + yy, xx = np.mgrid[0:h, 0:w] + frame = bg + 300.0 * (xx / w) + np.random.default_rng(3).normal(0, 8, (h, w)) + bayer = np.ones((h, w)) + bayer[::2, ::2] = 1.1 + bayer[1::2, 1::2] = 0.92 + frame *= bayer + for sy, sx in stars: + frame += peak * np.exp(-((xx - sx) ** 2 + (yy - sy) ** 2) / (2 * 1.5**2)) + return np.clip(frame, 0, 4095).astype(np.uint16) + + +@pytest.mark.unit +class TestSepDetect: + def test_detects_planted_stars_in_full_frame_coords(self): + stars = [(100, 200), (300, 700), (450, 120), (250, 480)] + frame = _synthetic_frame(stars) + result = sep_detect.detect_stars(frame, sigma=4.0) + assert result is not None + assert len(result.centroids) >= len(stars) + for sy, sx in stars: + d = np.hypot(result.centroids[:, 0] - sy, result.centroids[:, 1] - sx) + # full-frame coordinates: within 2 px of the planted position + assert d.min() < 2.0 + # flux-descending order + assert np.all(np.diff(result.fluxes) <= 0) + + def test_max_stars_cap(self): + stars = [(50 + 40 * i, 60 + 70 * (i % 12)) for i in range(30)] + frame = _synthetic_frame(stars) + result = sep_detect.detect_stars(frame, sigma=4.0, max_stars=10) + assert result is not None + assert len(result.centroids) <= 10 + + def test_unusable_frame_returns_none(self): + assert sep_detect.detect_stars(np.zeros((4, 4), dtype=np.uint16)) is None + assert sep_detect.detect_stars(np.zeros((10, 10, 3), dtype=np.uint16)) is None + + def test_bin2x2_geometry(self): + arr = np.arange(16, dtype=np.uint16).reshape(4, 4) + binned = sep_detect.bin2x2(arr) + assert binned.shape == (2, 2) + assert binned[0, 0] == pytest.approx((0 + 1 + 4 + 5) / 4) + + def test_edge_margin_drops_border_detections(self): + """Vignetted-border artifacts are excluded (field lesson: on a + saturated-interior frame every 'detection' hugged the frame edge).""" + stars = [(20, 300), (300, 20), (270, 480)] # two in the border zone + frame = _synthetic_frame(stars) + result = sep_detect.detect_stars(frame, sigma=4.0, edge_margin_px=48) + assert result is not None + for y, x in result.centroids: + assert 48 <= y < 540 - 48 + assert 48 <= x < 960 - 48 + # the interior star survives + d = np.hypot(result.centroids[:, 0] - 270, result.centroids[:, 1] - 480) + assert d.min() < 2.0 + + def test_saturated_interior_returns_zero_detections(self): + frame = np.full((540, 960), 4095, dtype=np.uint16) + # borders darker (vignette) so naive extraction would find edges + frame[:40, :] = 2000 + frame[-40:, :] = 2000 + result = sep_detect.detect_stars(frame, sigma=3.5, saturation_level=4095) + assert result is not None + assert len(result.centroids) == 0 + + +@pytest.mark.unit +class TestWarmPixelMap: + """Static single-pixel defects dominated SEP counts on empty sky + (2026-07-28 bench); the map removes them without touching stars.""" + + def test_excess_isolates_single_pixel_spike(self): + frame = np.full((64, 64), 1000.0, dtype=np.float32) + frame[20, 30] += 80.0 # warm pixel + excess = sep_detect.warm_pixel_excess(frame) + assert excess[20, 30] == pytest.approx(80.0) + # neighbours of the spike are not implicated + assert abs(excess[22, 30]) < 1.0 + assert abs(excess[20, 32]) < 1.0 + + def test_build_map_keeps_static_defects_drops_moving_star(self): + rng = np.random.default_rng(7) + warm = [(20, 30), (100, 200)] + frames = [] + for i in range(6): + f = 1000.0 + rng.normal(0, 5, (128, 256)) + for wy, wx in warm: + f[wy, wx] += 60.0 + f[50, 40 + 20 * i] += 300.0 # bright star drifting with the sky + frames.append(f.astype(np.uint16)) + pts = sep_detect.build_warm_pixel_map(frames, min_excess_adu=25.0) + assert {tuple(p) for p in pts} == set(warm) + + def test_build_map_empty_input(self): + assert len(sep_detect.build_warm_pixel_map([])) == 0 + + def test_detect_stars_masks_mapped_position_and_counts(self): + stars = [(100, 200), (300, 700), (450, 120), (250, 480)] + frame = _synthetic_frame(stars) + warm_map = np.array([[300, 700]]) # mask one planted "star" + result = sep_detect.detect_stars(frame, sigma=4.0, warm_pixel_map=warm_map) + assert result is not None + assert result.masked_count >= 1 + d = np.hypot(result.centroids[:, 0] - 300, result.centroids[:, 1] - 700) + assert len(d) == 0 or d.min() > 4.0 + # the unmasked stars all survive + for sy, sx in [(100, 200), (450, 120), (250, 480)]: + d = np.hypot(result.centroids[:, 0] - sy, result.centroids[:, 1] - sx) + assert d.min() < 2.0 + + def test_detect_stars_no_map_reports_zero_masked(self): + frame = _synthetic_frame([(100, 200)]) + result = sep_detect.detect_stars(frame, sigma=4.0) + assert result is not None + assert result.masked_count == 0 + + +@pytest.mark.unit +class TestRotationConvention: + """rotate_centroids must match what PIL does to the image.""" + + @pytest.mark.parametrize("angle", [0, 90, 180, 270]) + def test_quarter_turns_match_pil_on_square(self, angle): + h = w = 64 + y0, x0 = 10, 45 + img = np.zeros((h, w), dtype=np.uint8) + img[y0, x0] = 255 + rotated = np.asarray(Image.fromarray(img).rotate(angle)) + expect = np.unravel_index(rotated.argmax(), rotated.shape) + got, (nh, nw) = sfm.rotate_centroids( + np.array([[y0, x0]], dtype=float), (h, w), angle + ) + assert (round(got[0, 0]), round(got[0, 1])) == expect + assert (nh, nw) == (h, w) + + def test_quarter_turn_swaps_rect_canvas(self): + got, (nh, nw) = sfm.rotate_centroids(np.array([[0.0, 0.0]]), (1080, 1920), 90) + assert (nh, nw) == (1920, 1080) + # top-left pixel goes to bottom-left under CCW + assert got[0, 0] == pytest.approx(1919.0) + assert got[0, 1] == pytest.approx(0.0) + + def test_arbitrary_angle_matches_pil_blob(self): + h = w = 101 + y0, x0 = 30, 70 + img = np.zeros((h, w), dtype=np.float32) + yy, xx = np.mgrid[0:h, 0:w] + img += 255 * np.exp(-((xx - x0) ** 2 + (yy - y0) ** 2) / (2 * 2.0**2)) + rotated = np.asarray( + Image.fromarray(img.astype(np.uint8)).rotate(30, resample=Image.BILINEAR) + ) + py, px = np.unravel_index(rotated.argmax(), rotated.shape) + got, _ = sfm.rotate_centroids(np.array([[y0, x0]], dtype=float), (h, w), 30) + assert got[0, 0] == pytest.approx(py, abs=1.5) + assert got[0, 1] == pytest.approx(px, abs=1.5) + + +@pytest.mark.unit +class TestTargetPixelMapping: + def test_center_is_invariant(self): + # (255.5, 255.5) is the 512-frame centre -> maps to the canvas centre + y, x = sfm.map_target_pixel_to_frame((255.5, 255.5), (1920, 1080), 980) + assert y == pytest.approx((1920 - 1) / 2) + assert x == pytest.approx((1080 - 1) / 2) + + def test_offset_scales_by_crop_ratio(self): + # 100 px right of centre in 512-space = 100 * 980/512 sensor px + y, x = sfm.map_target_pixel_to_frame((255.5, 355.5), (1080, 1920), 980) + assert y == pytest.approx((1080 - 1) / 2) + assert x == pytest.approx((1920 - 1) / 2 + 100 * 980 / 512) + + def test_fov_scales_with_width(self): + assert sfm.fov_estimate_deg(980, 980) == pytest.approx(12.0) + assert sfm.fov_estimate_deg(1920, 980) == pytest.approx(23.51, abs=0.01) + + def test_stage5_rotation_matches_camera_interface_rules(self): + from PiFinder.camera_interface import SCREEN_ROTATE_AMOUNTS + + # camera_rotation overrides the screen_direction map entirely + assert sfm.stage5_rotation_deg("right", 45) == 315.0 + assert sfm.stage5_rotation_deg(None, 0) == 0.0 + # every mapped variant, plus the documented 270 fallback + for direction, rotation in SCREEN_ROTATE_AMOUNTS.items(): + assert sfm.stage5_rotation_deg(direction, None) == float(rotation) + assert sfm.stage5_rotation_deg("no_such_variant", None) == 270.0 diff --git a/python/tests/test_sep_fullframe_solve.py b/python/tests/test_sep_fullframe_solve.py new file mode 100644 index 000000000..246d4ad61 --- /dev/null +++ b/python/tests/test_sep_fullframe_solve.py @@ -0,0 +1,124 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +End-to-end equivalence of the SEP full-frame solve path. + +Projects real catalog stars (tetra3's own star table) onto a synthetic +full-sensor frame, then solves the same sky twice: + +* production path: centred crop -> 512 resize -> stage-5 rotation +* SEP path: full-frame centroids -> rotate_centroids -> + map_target_pixel_to_frame -> fov_estimate_deg + +and asserts both return the same Roll and the same aligned pointing at +the (off-centre) target pixel. This is the proof that swapping the +detection frame cannot disturb tracking or push-to. +""" + +import numpy as np +import pytest + +from PiFinder import solver_frame_map as sfm +from PiFinder import utils + +tetra3 = pytest.importorskip("tetra3") + +# imx462 geometry (the fielded camera) +FULL_H, FULL_W = 1080, 1920 +CROP_Y0, CROP_X0 = 50, 470 +CROP_W = 980 +PLATE_SCALE = sfm.SOLVER_FOV_DEG / CROP_W # deg per full-res px +TARGET_512 = (300.0, 340.0) # off-centre alignment point, rotated-512 space + + +def _project_stars(t3, ra0_deg, dec0_deg, max_stars=80): + """Gnomonic projection of catalog stars onto the full frame (y, x).""" + st = t3.star_table # columns: ra, dec (rad), x, y, z, mag + ra0, dec0 = np.deg2rad(ra0_deg), np.deg2rad(dec0_deg) + cosd = np.sin(dec0) * np.sin(st[:, 1]) + np.cos(dec0) * np.cos(st[:, 1]) * np.cos( + st[:, 0] - ra0 + ) + sel = np.where(cosd > np.cos(np.deg2rad(13)))[0] + sel = sel[np.argsort(st[sel, -1])][:max_stars] + ra, dec = st[sel, 0], st[sel, 1] + cosc = np.sin(dec0) * np.sin(dec) + np.cos(dec0) * np.cos(dec) * np.cos(ra - ra0) + x_ang = np.rad2deg(np.cos(dec) * np.sin(ra - ra0) / cosc) + y_ang = np.rad2deg( + (np.cos(dec0) * np.sin(dec) - np.sin(dec0) * np.cos(dec) * np.cos(ra - ra0)) + / cosc + ) + cy, cx = (FULL_H - 1) / 2, (FULL_W - 1) / 2 + x_px = cx - x_ang / PLATE_SCALE + y_px = cy - y_ang / PLATE_SCALE + inside = (x_px >= 0) & (x_px < FULL_W) & (y_px >= 0) & (y_px < FULL_H) + return np.column_stack((y_px[inside], x_px[inside])) + + +@pytest.mark.unit +def test_sep_fullframe_solve_matches_production_pointing(): + db_path = utils.tetra3_dir / "data" / "default_database.npz" + if not db_path.exists(): + pytest.skip("tetra3 default database not present (submodule not populated)") + t3 = tetra3.Tetra3(str(db_path)) + cents_full = _project_stars(t3, ra0_deg=84.0, dec0_deg=0.0) + assert len(cents_full) > 30 + + # target_pixel as persisted: rotated-512 space (stage-5 rotation 90) + tp512r, _ = sfm.rotate_centroids(np.array([TARGET_512]), (512, 512), 90) + tp512r = tuple(tp512r[0]) + + # --- production path + yc = cents_full[:, 0] - CROP_Y0 + xc = cents_full[:, 1] - CROP_X0 + ok = (yc >= 0) & (yc < CROP_W) & (xc >= 0) & (xc < CROP_W) + c512 = np.column_stack((yc[ok], xc[ok])) * (sfm.SOLVER_FRAME_PX / CROP_W) + c512r, canvas512 = sfm.rotate_centroids(c512, (512, 512), 90) + sol_prod = t3.solve_from_centroids( + c512r, + canvas512, + fov_estimate=12.0, + fov_max_error=4.0, + match_max_error=0.005, + target_pixel=tp512r, + solve_timeout=1000, + ) + assert sol_prod.get("RA") is not None + + # --- SEP full-frame path (exactly what sep_shadow.solve does) + cfull_r, canvas_full = sfm.rotate_centroids(cents_full, (FULL_H, FULL_W), 90) + tp_full = sfm.map_target_pixel_to_frame(tp512r, canvas_full, CROP_W) + fov = sfm.fov_estimate_deg(canvas_full[1], CROP_W) + sol_sep = t3.solve_from_centroids( + cfull_r, + canvas_full, + fov_estimate=fov, + fov_max_error=fov / 3, + match_max_error=0.005, + target_pixel=tp_full, + solve_timeout=1000, + ) + assert sol_sep.get("RA") is not None + + # The wider frame must see MORE of the sky's stars + assert sol_sep["Matches"] > sol_prod["Matches"] + + # Same camera pointing and identical Roll convention + assert abs(sol_prod["RA"] - sol_sep["RA"]) * 3600 < 120 + assert abs(sol_prod["Dec"] - sol_sep["Dec"]) * 3600 < 120 + droll = abs((sol_prod["Roll"] - sol_sep["Roll"] + 180) % 360 - 180) + assert droll < 0.05 + + # Aligned pointing at the target pixel agrees to within a fit residual + assert abs(sol_prod["RA_target"] - sol_sep["RA_target"]) * 3600 < 60 + assert abs(sol_prod["Dec_target"] - sol_sep["Dec_target"]) * 3600 < 60 + + +@pytest.mark.unit +def test_target_pixel_mapping_round_trip(): + """map_frame_pixel_to_target inverts map_target_pixel_to_frame exactly.""" + canvas = (1920, 1080) # rotated full frame (imx462, 90 deg) + for tp in [(256.0, 256.0), (100.5, 400.25), (0.0, 511.0)]: + full = sfm.map_target_pixel_to_frame(tp, canvas, 980) + back = sfm.map_frame_pixel_to_target(full, canvas, 980) + assert abs(back[0] - tp[0]) < 1e-9 + assert abs(back[1] - tp[1]) < 1e-9 diff --git a/python/tests/test_sep_shadow.py b/python/tests/test_sep_shadow.py new file mode 100644 index 000000000..51a4fbea9 --- /dev/null +++ b/python/tests/test_sep_shadow.py @@ -0,0 +1,87 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +""" +Unit tests for the SEP fallback backoff (SepShadowRunner). + +A failed fallback solve costs up to solve_timeout (1 s) of solver CPU; +on persistently unsolvable scenes (indoors, thick cloud) that recurs on +every attempt. The backoff skips a growing number of attempts after +consecutive failures, but must re-arm immediately when the SEP count +jumps -- a cloud gap opening on real stars must not wait out a backoff +window. +""" + +import pytest + +from PiFinder.sep_shadow import SepShadowRunner + + +def _runner(tmp_path): + return SepShadowRunner( + shadow_enabled=False, + fallback_enabled=True, + sigma=4.0, + rotation_deg=90.0, + crop_width_px=980, + csv_path=tmp_path / "shadow.csv", + ) + + +def _tick(runner, n=1): + """Advance the per-attempt counter the way detect() does.""" + runner._attempt_counter += n + + +@pytest.mark.unit +class TestFallbackBackoff: + def test_first_attempt_always_allowed(self, tmp_path): + runner = _runner(tmp_path) + _tick(runner) + assert runner.fallback_should_attempt(28) is True + + def test_failures_open_growing_skip_windows(self, tmp_path): + runner = _runner(tmp_path) + _tick(runner) + runner.record_fallback_result(False, 28) + # streak 1 -> skip 2 attempts + _tick(runner) + assert runner.fallback_should_attempt(28) is False + _tick(runner) + assert runner.fallback_should_attempt(28) is True + runner.record_fallback_result(False, 28) + # streak 2 -> skip 4 attempts + _tick(runner, 3) + assert runner.fallback_should_attempt(28) is False + _tick(runner) + assert runner.fallback_should_attempt(28) is True + + def test_skip_window_caps_at_eight_attempts(self, tmp_path): + runner = _runner(tmp_path) + for _ in range(10): # streak far past the cap + _tick(runner) + runner.record_fallback_result(False, 28) + _tick(runner, 8) + assert runner.fallback_should_attempt(28) is True + + def test_sep_count_jump_rearms_immediately(self, tmp_path): + """Cloud gap opens on stars: masked count jumps 5 -> 30. The + rescue solve must run right away, not wait out the window.""" + runner = _runner(tmp_path) + _tick(runner) + runner.record_fallback_result(False, 20) + _tick(runner) + assert runner.fallback_should_attempt(20) is False + assert runner.fallback_should_attempt(30) is True # >= 1.5x + + def test_success_and_production_solve_clear_the_streak(self, tmp_path): + runner = _runner(tmp_path) + _tick(runner) + runner.record_fallback_result(False, 28) + runner.record_fallback_result(True, 28) + _tick(runner) + assert runner.fallback_should_attempt(28) is True + + runner.record_fallback_result(False, 28) + runner.note_solved() + _tick(runner) + assert runner.fallback_should_attempt(28) is True