diff --git a/CHANGELOG.md b/CHANGELOG.md index 4ac48aa..d5a9706 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,38 @@ match the `VERSION` file and `v*` git tags. ## [Unreleased] +### Changed + +- **`--auto` no longer applies one preset's settings to every target.** A setting only some target + presets define was blended over those presets alone, so e.g. the galaxy-only coarse chroma pass, + the nebula anisotropic diffusion and `pre_gradient_removal` ran on star fields and clusters too. + Presets that don't list a setting now vote for its default. **Expect different `--auto` output** + on most targets, and re-save any `_config.toml` you rely on. + +### Fixed + +- Darks whose exposure differs from the lights are now scaled on the default parallel path (pool + workers never received `dark_exptime` and subtracted the dark unscaled, silently). +- Hierarchical runs: sessions with darks no longer lose their `info.json` metadata (Bayer pattern, + WCS, GPS, target name) to a variable-shadowing bug. +- `--photometry` and `--color-calibrate` work again: the stacked cube's 3-axis WCS made Gaia + projection fail silently, and the colour-calibration scale was always exactly 1. +- With `--auto` on 15+ frames, `--stack-method median`/`linear_fit`/`ivw`/`wavelet` are honoured + instead of silently becoming an unrejected patch-weighted mean. +- Seeded affine star matching failed for any shift over ~5 px (sign error), leaving frames on + translation-only registration whenever the blind matcher also failed. +- Resuming from a phase-1 or phase-2 checkpoint now runs target inference and the `--auto` advisor. +- Blind PSF estimation (used by `--deconvolve` under `--auto`) no longer flattens the PSF; the + broken "blind RL refinement" was removed. +- `--transient-detect`: S_corr has unit variance when the two epochs' flux scales differ (the ZOGY + kernel flux powers were swapped), so thresholds mean what they say. +- Origin `DATE-OBS` is local time: the `TIMEZONE` keyword is now applied (and kept in the stacked + header), fixing the `--fix-atmospheric-dispersion` zenith angle and light-curve MJD/airmass. + A below-horizon zenith angle now declines instead of correcting the wrong way. +- `--stream`: a pixel seeded by a single frame no longer freezes on that frame's value. +- `--use-gpu`: the next target in the same process no longer reuses the previous target's GPU-side + calibration masters. + ## [2.3.0] - 2026-09-23 ### Added diff --git a/ext/astro_native/Cargo.lock b/ext/astro_native/Cargo.lock index 031c263..cd78254 100644 --- a/ext/astro_native/Cargo.lock +++ b/ext/astro_native/Cargo.lock @@ -49,7 +49,7 @@ checksum = "fb5dfbc6d8d2675589ccbe4d0fd61df2419075625f8c1a62325e718e2b0049f9" [[package]] name = "astro_native" -version = "0.36.0" +version = "0.37.0" dependencies = [ "numpy", "pyo3", diff --git a/ext/astro_native/Cargo.toml b/ext/astro_native/Cargo.toml index a40dd11..4c0e272 100644 --- a/ext/astro_native/Cargo.toml +++ b/ext/astro_native/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "astro_native" -version = "0.36.0" +version = "0.37.0" edition = "2021" description = "Native (Rust) hot-path kernels for OriginStack: stacking combine, etc." diff --git a/ext/astro_native/pyproject.toml b/ext/astro_native/pyproject.toml index 24850f1..14cd8ba 100644 --- a/ext/astro_native/pyproject.toml +++ b/ext/astro_native/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "astro_native" -version = "0.36.0" +version = "0.37.0" description = "Native Rust hot-path kernels for OriginStack" requires-python = ">=3.10" classifiers = ["Programming Language :: Rust"] diff --git a/ext/astro_native/src/lib.rs b/ext/astro_native/src/lib.rs index abde232..759d824 100644 --- a/ext/astro_native/src/lib.rs +++ b/ext/astro_native/src/lib.rs @@ -569,25 +569,34 @@ fn burnin_seed_pixel( /// borders uncovered). Running the accept-test against that fabricated state /// would reject the first real sample forever; instead initialize directly /// from it. +/// +/// Below `ONLINE_CLIP_MIN_SAMPLES` accepted samples there is no spread to +/// clip against (one sample has M2 = 0, so every later sample failed the +/// 1e-6 test and the pixel froze on a single frame), so samples are accepted +/// unconditionally until the state holds that many. #[inline] fn fold_pixel(mean: f64, m2: f64, n_acc: f64, x: f64, sigma: f64) -> (f64, f64, f64, bool) { if n_acc <= 0.0 { return (x, 0.0, 1.0, true); } - let var_est = m2 / n_acc; - let std_est = var_est.max(1e-12).sqrt(); - if (x - mean).abs() <= sigma * std_est { - let n_acc_new = n_acc + 1.0; - let delta = x - mean; - let new_mean = mean + delta / n_acc_new; - let delta2 = x - new_mean; - let new_m2 = m2 + delta * delta2; - (new_mean, new_m2, n_acc_new, true) - } else { - (mean, m2, n_acc, false) + if n_acc >= ONLINE_CLIP_MIN_SAMPLES { + let std_est = (m2 / n_acc).max(1e-12).sqrt(); + if (x - mean).abs() > sigma * std_est { + return (mean, m2, n_acc, false); + } } + let n_acc_new = n_acc + 1.0; + let delta = x - mean; + let new_mean = mean + delta / n_acc_new; + let delta2 = x - new_mean; + let new_m2 = m2 + delta * delta2; + (new_mean, new_m2, n_acc_new, true) } +/// Accepted samples a pixel needs before the online sigma-clip starts +/// rejecting. Mirrored by `ONLINE_CLIP_MIN_SAMPLES` in src/stacking.py. +const ONLINE_CLIP_MIN_SAMPLES: f64 = 3.0; + /// Per-pixel online sigma-clip: a MAD-rejected burn-in window (first `k = /// min(burn_in, n)` samples) seeds a running (mean, M2) Welford state; each /// remaining sample is tested against that running estimate before being diff --git a/src/auto_settings.py b/src/auto_settings.py index bd62412..8a426a7 100644 --- a/src/auto_settings.py +++ b/src/auto_settings.py @@ -503,13 +503,12 @@ def _apply_dynamic_settings( ) -> List[str]: """Continuous replacement for _apply_target_settings(): every parameter is blended across all 8 presets weighted by _blend_weights, instead of - looking up one bucket's fixed table. Numeric parameters get a weighted - average (renormalized over just the presets that define that - parameter); boolean/string parameters -- which can't be fractionally - blended -- take the value from whichever contributing preset has the - single highest weight (nearest-neighbor-in-signal-space), so the - decision still shifts continuously with the signals even though the - final choice at any instant is binary. + looking up one bucket's fixed table. A preset that does not list a + parameter contributes the baseline value args already holds. Numeric + parameters get a weighted average; boolean/string parameters -- which + can't be fractionally blended -- take the value with the largest total + weight, so the decision still shifts continuously with the signals even + though the final choice at any instant is binary. """ changes: List[str] = [] _explicit = getattr(args, '_explicit_cli_dests', set()) @@ -532,18 +531,33 @@ def _set(attr: str, val: object) -> None: seen.add(attr) all_attrs.append(attr) + # A preset that does not list an attr means "leave it at the baseline", + # not "abstain": it votes for the value args held before this blend, with + # its own weight. Renormalizing over only the listing presets applied a + # one-preset value at full strength to every target -- galaxy's coarse + # chroma pass (chroma_nr_large_sigma 0 -> 50) on a 97%-star-field blend. + # A baseline of None (attr absent from args) has nothing to vote for, so + # that attr falls back to the listing presets alone. for attr in all_attrs: - contributors = [(t, val) for t, rows in _TARGET_SETTINGS.items() - for a, val in rows if a == attr] - w_sum = sum(weights.get(t, 0.0) for t, _ in contributors) + listed = {t: val for t, rows in _TARGET_SETTINGS.items() + for a, val in rows if a == attr} + sample_val = next(iter(listed.values())) + baseline = getattr(args, attr, None) + votes = [(weights.get(t, 0.0), listed[t] if t in listed else baseline) + for t in _TARGET_SETTINGS + if t in listed or baseline is not None] + w_sum = sum(w for w, _ in votes) if w_sum <= 0: continue - sample_val = contributors[0][1] if isinstance(sample_val, (bool, str)): - best_t, best_val = max(contributors, key=lambda tv: weights.get(tv[0], 0.0)) - _set(attr, best_val) + # Largest total weight per value, not the single heaviest preset: + # six presets that leave a flag off outweigh one that turns it on. + totals: Dict[object, float] = {} + for w, val in votes: + totals[val] = totals.get(val, 0.0) + w + _set(attr, max(totals, key=totals.get)) else: - blended = sum(weights.get(t, 0.0) * val for t, val in contributors) / w_sum + blended = sum(w * float(val) for w, val in votes) / w_sum if isinstance(sample_val, int): blended = int(round(blended)) else: diff --git a/src/cli.py b/src/cli.py index 14cb170..04e0441 100644 --- a/src/cli.py +++ b/src/cli.py @@ -6,7 +6,6 @@ import logging import math import os -import re import sys import tempfile import threading @@ -27,6 +26,7 @@ from src.models import Config, ProcessingStats from src.pipeline import stack_target from src.utils import ( + TZ_OFFSET_RE, disable_astropy_network, format_time, print_header, @@ -676,10 +676,6 @@ def _run_combined_sessions(subdirs: list, output: str, args: argparse.Namespace) _ROTATION_SPLIT_THRESHOLD_DEG = 3.0 -# A FITS TIMEZONE value this code can append to DATE-OBS: '-0700', '+05:30'. -# Anything else ('PDT', 'US/Pacific', '-07:00 (PDT)') would produce an -# unparseable timestamp and silently disable the rotation prediction. -_TZ_OFFSET_RE = re.compile(r'[+-]\d{2}:?\d{2}') def _predict_rotation_spread(subdirs: List[str]) -> Optional[float]: @@ -728,7 +724,7 @@ def _give_up(d, why): h0, h1 = lights[0].header or {}, lights[-1].header or {} offset = str(h0.get('TIMEZONE', '') or '').strip() - if offset and not _TZ_OFFSET_RE.fullmatch(offset): + if offset and not TZ_OFFSET_RE.fullmatch(offset): return _give_up(d, f"TIMEZONE {offset!r} is not a +/-HHMM offset") ra, dec = math.degrees(si.ra_rad), math.degrees(si.dec_rad) @@ -1050,9 +1046,11 @@ def process_directory(directory: str, output: str, args: argparse.Namespace): safe_print(f" Bias: pedestal={b_med:.1f} ADU " f"noise={b_std:.1f} ADU -> {b_quality}") if masters.get('dark') is not None: - d = masters['dark'] - dark_med = float(np.median(d)) - dark_peak = float(d.max()) + # not `d`: that is the enclosing loop's target directory, + # read again below for args._input_directory + dk = masters['dark'] + dark_med = float(np.median(dk)) + dark_peak = float(dk.max()) dark_et = masters.get('dark_exptime') dark_hdr = frames['dark'][0].header if frames.get('dark') else {} dark_temp_c = dark_hdr.get('CCD-TEMP') diff --git a/src/color_calibrate.py b/src/color_calibrate.py index 16a101b..26bd371 100644 --- a/src/color_calibrate.py +++ b/src/color_calibrate.py @@ -47,7 +47,7 @@ def _field_radius_deg(header) -> float: """Estimate the field-of-view radius in degrees from WCS.""" try: - wcs = WCS(header) + wcs = WCS(header).celestial # (3, H, W) cube header -> 2-axis naxis1 = int(header.get("NAXIS1", 0)) naxis2 = int(header.get("NAXIS2", 0)) if naxis1 == 0 or naxis2 == 0: @@ -236,6 +236,21 @@ def _synthetic_channel_flux_batch(teff_k: np.ndarray, channel_response=None return fluxes[0], fluxes[1], fluxes[2] +def _channel_correction(meas: np.ndarray, expected: np.ndarray) -> float: + """Multiplicative factor that brings a channel's measured star fluxes to + the catalogue's expected ones: median(expected / measured). + + ``apply_photometric_calibration`` multiplies the image by this, so a + channel reading too bright gets a factor below 1. Units (ADU vs. the + catalogue's 10^(-0.4 m)) are common to all three channels and cancel in + the callers' mean normalisation. This used to divide each star's ratio by + the channel's own median ratio before taking the median -- always 1.0 -- + so --color-calibrate never changed the image. + """ + ratio = np.maximum(expected, 1e-30) / np.maximum(meas, 1e-30) + return float(np.median(ratio)) + + def fit_channel_scales_spcc(img: np.ndarray, header, catalog, channel_response=None, verbose: bool = False) -> Tuple[float, float, float]: @@ -313,14 +328,9 @@ def fit_channel_scales_spcc(img: np.ndarray, header, catalog, flux_r_expected[fallback_mask] = 10.0 ** (-0.4 * fallback_r[fallback_mask]) flux_b_expected[fallback_mask] = 10.0 ** (-0.4 * fallback_b[fallback_mask]) - def _robust_ratio(meas: np.ndarray, expected: np.ndarray) -> float: - ratio = meas / np.maximum(expected, 1e-30) - ratio = ratio / np.median(ratio) - return float(np.median(ratio)) - - scale_r = _robust_ratio(fluxes[:, 0], flux_r_expected) - scale_g = _robust_ratio(fluxes[:, 1], flux_g_expected) - scale_b = _robust_ratio(fluxes[:, 2], flux_b_expected) + scale_r = _channel_correction(fluxes[:, 0], flux_r_expected) + scale_g = _channel_correction(fluxes[:, 1], flux_g_expected) + scale_b = _channel_correction(fluxes[:, 2], flux_b_expected) mean_scale = (scale_r + scale_g + scale_b) / 3.0 if mean_scale > 0: @@ -403,15 +413,10 @@ def fit_channel_scales(img: np.ndarray, header, flux_b_expected = 10.0 ** (-0.4 * k) # Compute per-star measured ratios vs expected ratios - # scale_R = median(flux_R_measured / flux_R_expected) normalised to G - def _robust_ratio(meas: np.ndarray, expected: np.ndarray) -> float: - ratio = meas / np.maximum(expected, 1e-30) - ratio = ratio / np.median(ratio) # normalise so green ≈ 1 - return float(np.median(ratio)) - - scale_r = _robust_ratio(fluxes[:, 0], flux_r_expected) - scale_g = _robust_ratio(fluxes[:, 1], flux_g_expected) - scale_b = _robust_ratio(fluxes[:, 2], flux_b_expected) + # Per-channel correction = median(expected / measured); mean-normalised below. + scale_r = _channel_correction(fluxes[:, 0], flux_r_expected) + scale_g = _channel_correction(fluxes[:, 1], flux_g_expected) + scale_b = _channel_correction(fluxes[:, 2], flux_b_expected) # Normalise so that mean(scale) = 1 (preserve overall brightness) mean_scale = (scale_r + scale_g + scale_b) / 3.0 diff --git a/src/difference_imaging.py b/src/difference_imaging.py index 5681702..08f5092 100644 --- a/src/difference_imaging.py +++ b/src/difference_imaging.py @@ -276,9 +276,13 @@ def zogy(new: np.ndarray, ref: np.ndarray, del d_hat, pd_hat # --- Noise terms for S_corr (eq. 26-30) --- - # Matched-filter kernels for each input image. - kn_hat = fr * fn ** 2 * np.conj(pn_hat) * abs_pr2 / denom - kr_hat = fn * fr ** 2 * np.conj(pr_hat) * abs_pn2 / denom + # Matched-filter kernels for each input image: S = kn*N - kr*R, which is + # F_D * D (x) P_D expanded (eq. 28-29). The flux powers were once swapped, + # invisible at flux_new == flux_ref == 1 but mis-scaling the variance + # whenever the epochs differ (S_corr std 0.4-0.65 on pure noise at a + # flux ratio of 3 or 0.3, so a "5 sigma" cut really sat at 8-12 sigma). + kn_hat = fn * fr ** 2 * np.conj(pn_hat) * abs_pr2 / denom + kr_hat = fr * fn ** 2 * np.conj(pr_hat) * abs_pn2 / denom del abs_pn2, abs_pr2, denom, pn_hat, pr_hat # Variance maps live on the padded grid too. The padding is sky, so it diff --git a/src/frame_processor.py b/src/frame_processor.py index 40a061a..0dbb33c 100644 --- a/src/frame_processor.py +++ b/src/frame_processor.py @@ -216,13 +216,16 @@ def _process_single_frame(path: str, header: dict, masters: Dict[str, Optional[n _use_gpu_calib = False if _gpu_ctx.active and data.ndim == 2: _shape = data.shape - if _shape not in _gpu_calib_cache: - with _gpu_calib_lock: - if _shape not in _gpu_calib_cache: - _ensure_gpu_masters(masters, _gpu_ctx) - _gpu_calib_cache[_shape] = _probe_gpu_calibration( - _shape[0], _shape[1], _gpu_ctx, masters) - _use_gpu_calib = _gpu_calib_cache.get(_shape, False) + with _gpu_calib_lock: + # Every frame, not only on a new shape: the next target in the same + # process (hierarchical run, second desktop-app job) has the same + # sensor size but different masters. The check is an identity + # comparison; a change re-uploads and clears the probe cache. + _ensure_gpu_masters(masters, _gpu_ctx) + if _shape not in _gpu_calib_cache: + _gpu_calib_cache[_shape] = _probe_gpu_calibration( + _shape[0], _shape[1], _gpu_ctx, masters) + _use_gpu_calib = _gpu_calib_cache[_shape] # Calibration — preserve negative noise through bias/dark subtraction, # clip only once after all steps to avoid cumulative truncation of shadow detail @@ -524,26 +527,38 @@ def _process_single_frame(path: str, header: dict, masters: Dict[str, Optional[n _gpu_calib_cache: Dict[Tuple[int, int], bool] = {} _gpu_calib_lock = threading.Lock() _gpu_masters: Dict[str, Any] = {} # GPU-resident copies of master arrays -_gpu_masters_sig: Optional[Tuple] = None # signature to detect master changes +_gpu_masters_sig: Optional[Tuple] = None # the host arrays they were uploaded from + +_GPU_MASTER_KEYS = ('bias', 'dark', 'flat', '_flat_norm', 'hot_pixel_map') def _masters_sig(masters: Dict) -> Tuple: - """Lightweight signature so we can detect when masters change between sessions.""" - def _s(key): - arr = masters.get(key) - return (arr.shape, arr.dtype.str) if isinstance(arr, np.ndarray) else None - return (_s('dark'), _s('flat'), _s('bias')) + """The host master arrays themselves, compared by identity. + + A (shape, dtype) signature matched every target from the same sensor, so + the next target in the process was calibrated with the previous target's + GPU-side dark/flat -- even one with no dark at all. Holding the arrays + (not their ``id()``) means a freed array's address can't be reused into a + false match. + """ + return tuple(masters.get(k) for k in _GPU_MASTER_KEYS) + + +def _same_masters(a: Optional[Tuple], b: Optional[Tuple]) -> bool: + return a is not None and b is not None and all(x is y for x, y in zip(a, b)) def _ensure_gpu_masters(masters: Dict, gpu) -> None: - """Upload master calibration arrays to GPU once per session; re-uploads on change.""" + """Upload master calibration arrays to GPU; re-uploads (and drops the + per-shape calibration probe results) when the host masters change.""" global _gpu_masters, _gpu_masters_sig sig = _masters_sig(masters) - if sig == _gpu_masters_sig and _gpu_masters: + if _same_masters(sig, _gpu_masters_sig) and _gpu_masters: return + _gpu_calib_cache.clear() xp = gpu.xp result: Dict[str, Any] = {} - for key in ('bias', 'dark', 'flat', '_flat_norm', 'hot_pixel_map'): + for key in _GPU_MASTER_KEYS: arr = masters.get(key) if isinstance(arr, np.ndarray): result[key] = xp.asarray(arr) @@ -691,9 +706,37 @@ def _banding_cfg(args) -> Optional[tuple]: float(getattr(args, 'banding_sigma', 3.0))) +def _share_masters(masters: Dict[str, Any]) -> Tuple[list, Dict[str, tuple], Dict[str, Any]]: + """Split *masters* for pool workers: numpy arrays go into shared memory + (returned as the blocks to unlink plus name -> (shm_name, dtype, shape) + specs), everything else -- scalars such as ``dark_exptime`` -- is + returned as a plain dict passed through the initializer. + + The scalars used to be dropped, so every worker saw no ``dark_exptime`` + and subtracted a dark of a different exposure unscaled. + """ + shm_blocks: list = [] + shm_specs: Dict[str, tuple] = {} + scalars: Dict[str, Any] = {} + for name, arr in masters.items(): + if arr is None or name.startswith('_shm_'): + continue + if not isinstance(arr, np.ndarray): + scalars[name] = arr + continue + arr_c = np.ascontiguousarray(arr) + shm = SharedMemory(create=True, size=arr_c.nbytes) + shm_arr = np.ndarray(arr_c.shape, dtype=arr_c.dtype, buffer=shm.buf) + shm_arr[:] = arr_c + shm_blocks.append(shm) + shm_specs[name] = (shm.name, arr_c.dtype.str, arr_c.shape) + return shm_blocks, shm_specs, scalars + + def _init_worker_shm(shm_specs: Dict[str, tuple], trail_reject: bool = False, banding: Optional[tuple] = None, - session_cfa: Optional[dict] = None) -> None: + session_cfa: Optional[dict] = None, + scalar_masters: Optional[Dict[str, Any]] = None) -> None: """Initializer for pool workers — attach to shared-memory calibration arrays. *shm_specs* maps master name → (shm_name, dtype_str, shape). Workers @@ -706,7 +749,7 @@ def _init_worker_shm(shm_specs: Dict[str, tuple], trail_reject: bool = False, _worker_trail_reject = bool(trail_reject) _worker_banding = banding set_session_cfa(session_cfa) - _worker_masters = {} + _worker_masters = dict(scalar_masters or {}) for name, (shm_name, dtype_str, shape) in shm_specs.items(): shm = SharedMemory(name=shm_name, create=False) arr = np.ndarray(shape, dtype=np.dtype(dtype_str), buffer=shm.buf) @@ -1000,17 +1043,7 @@ def _accum(timings: Optional[dict]) -> None: # Share calibration arrays via shared memory — zero disk I/O, one copy # in RAM shared across all workers (read-only view per worker process). - shm_blocks: list = [] - shm_specs: Dict[str, tuple] = {} - for name, arr in masters.items(): - if arr is None or name.startswith('_shm_') or not isinstance(arr, np.ndarray): - continue - arr_c = np.ascontiguousarray(arr) - shm = SharedMemory(create=True, size=arr_c.nbytes) - shm_arr = np.ndarray(arr_c.shape, dtype=arr_c.dtype, buffer=shm.buf) - shm_arr[:] = arr_c - shm_blocks.append(shm) - shm_specs[name] = (shm.name, arr_c.dtype.str, arr_c.shape) + shm_blocks, shm_specs, scalar_masters = _share_masters(masters) _ca = getattr(args, 'ca_correction', False) _cr = getattr(args, 'cosmic_ray_rejection', False) @@ -1032,7 +1065,7 @@ def _accum(timings: Optional[dict]) -> None: with ProcessPoolExecutor(max_workers=workers, mp_context=mp_context(), initializer=_init_worker_shm, initargs=(shm_specs, _tr, _banding_cfg(args), - _session_cfa)) as pool: + _session_cfa, scalar_masters)) as pool: futures = {pool.submit(_parallel_frame_worker, t): t[1] for t in tasks} _wv = _get_ui_events() _wv_done = 0 @@ -1442,17 +1475,7 @@ def _worker_count_cap(n_workers: int) -> int: safe_print(f" Reloading {n} accepted frames ({workers} workers, " f"quality analysis skipped)...") - shm_blocks: list = [] - shm_specs: Dict[str, tuple] = {} - for name, arr in masters.items(): - if arr is None or name.startswith('_shm_') or not isinstance(arr, np.ndarray): - continue - arr_c = np.ascontiguousarray(arr) - shm = SharedMemory(create=True, size=arr_c.nbytes) - shm_arr = np.ndarray(arr_c.shape, dtype=arr_c.dtype, buffer=shm.buf) - shm_arr[:] = arr_c - shm_blocks.append(shm) - shm_specs[name] = (shm.name, arr_c.dtype.str, arr_c.shape) + shm_blocks, shm_specs, scalar_masters = _share_masters(masters) tasks = [(final[i].path, final_indices[i], args.debayer_method, args.white_balance, mm_rgb_path, mm_lum_path, rgb_shape, lum_shape, @@ -1465,7 +1488,7 @@ def _worker_count_cap(n_workers: int) -> int: with ProcessPoolExecutor(max_workers=workers, mp_context=mp_context(), initializer=_init_worker_shm, initargs=(shm_specs, _tr, _banding_cfg(args), - _session_cfa)) as pool: + _session_cfa, scalar_masters)) as pool: futures = {pool.submit(_parallel_frame_worker, t): t[1] for t in tasks} for future in tqdm(as_completed(futures), total=n, desc=" Reloading", unit="frame", diff --git a/src/io_fits.py b/src/io_fits.py index a6075a3..4618eb7 100644 --- a/src/io_fits.py +++ b/src/io_fits.py @@ -461,7 +461,9 @@ def populate_fits_header(header: fits.Header, frames: List[FrameInfo], # image, not a raw CFA mosaic. Viewers such as Siril use BAYERPAT to # detect raw frames and will attempt to debayer the already-processed # image if the keyword is present, producing garbage. - copy_keys = ['TELESCOP', 'INSTRUME', 'OBSERVER', 'OBJECT', 'DATE-OBS', + # TIMEZONE travels with DATE-OBS: Origin stamps DATE-OBS in local time + # and readers (utils.obs_time_utc_iso) need the offset to get UTC. + copy_keys = ['TELESCOP', 'INSTRUME', 'OBSERVER', 'OBJECT', 'DATE-OBS', 'TIMEZONE', 'EXPTIME', 'CCD-TEMP', 'GAIN', 'OFFSET', 'XBINNING', 'YBINNING', 'XPIXSZ', 'YPIXSZ', 'FOCALLEN', 'APTDIA'] for key in copy_keys: diff --git a/src/models.py b/src/models.py index 361f52a..b1e4fd9 100644 --- a/src/models.py +++ b/src/models.py @@ -71,9 +71,6 @@ class Config: # (swept 1.0-2.0 on real data: 1.0-1.25 matches the old RBF's # large-scale flatness; larger trades flatness for smoothness) - # Blind PSF estimation - BLIND_PSF_ITERATIONS = 8 # RL iterations for PSF update - # Total Variation deconvolution TV_LAMBDA = 0.02 # TV regularisation weight TV_ITERATIONS = 50 # Gradient descent steps diff --git a/src/observing_geometry.py b/src/observing_geometry.py index d1460c9..26cde34 100644 --- a/src/observing_geometry.py +++ b/src/observing_geometry.py @@ -114,9 +114,14 @@ def airmass(ra_deg: float, dec_deg: float, lat_deg: float, lon_deg: float, def zenith_angle_deg(ra_deg: float, dec_deg: float, lat_deg: float, lon_deg: float, height_m: float, when_iso: str) -> Optional[float]: - """Zenith angle (90 - altitude), or None.""" + """Zenith angle (90 - altitude), or None -- also None when the target is + at or below the horizon, which only happens with a wrong time or site + (tan(z) changes sign past 90 deg, so a dispersion correction built on it + would shift the channels the wrong way).""" aa = altaz(ra_deg, dec_deg, lat_deg, lon_deg, height_m, when_iso) - return None if aa is None else 90.0 - aa[0] + if aa is None or aa[0] <= 0.0: + return None + return 90.0 - aa[0] def parallactic_angle_deg(ra_deg: float, dec_deg: float, lat_deg: float, diff --git a/src/photometry.py b/src/photometry.py index 4b14e42..217276b 100644 --- a/src/photometry.py +++ b/src/photometry.py @@ -41,7 +41,7 @@ import numpy as np from src.photometry_core import _field_centre_and_radius, _id_str, _pixel_coords, aperture_photometry_batch, row_nanmax -from src.utils import header_get_first +from src.utils import header_get_first, obs_time_utc_iso _log = logging.getLogger("originstack") @@ -67,12 +67,9 @@ def _airmass(header, session_info, ra_deg, dec_deg, when=None): """ if session_info is None or not getattr(session_info, "has_gps", False): return None - time_str = str(when) if when else None - if time_str is None: - time_str = header_get_first(header, ("DATE-OBS", "DATE_OBS", "DATEOBS"), - cast=str) - if time_str is None: - time_str = getattr(session_info, "date_time", None) + # obs_time_utc_iso applies the header's TIMEZONE: Origin DATE-OBS is local. + time_str = str(when) if when else obs_time_utc_iso( + header, fallback=getattr(session_info, "date_time", None)) if not time_str: return None from src.observing_geometry import airmass as _airmass_of diff --git a/src/photometry_core.py b/src/photometry_core.py index 8512432..a221fb4 100644 --- a/src/photometry_core.py +++ b/src/photometry_core.py @@ -62,7 +62,9 @@ def _pixel_coords(table, header) -> Optional[np.ndarray]: except Exception: return None try: - wcs = WCS(header) + # .celestial: the stacked product is a (3, H, W) cube, so a bare + # WCS(header) is 3-axis and the 2-argument all_world2pix below raises. + wcs = WCS(header).celestial if "ra" in table.colnames: ra = np.array(table["ra"], dtype=float) dec = np.array(table["dec"], dtype=float) diff --git a/src/photometry_timeseries.py b/src/photometry_timeseries.py index cedb303..3ed87a2 100644 --- a/src/photometry_timeseries.py +++ b/src/photometry_timeseries.py @@ -30,7 +30,7 @@ from src.photometry import _GAIA_BAND_FOR_CHANNEL, _airmass, _read_gain, match_gaia_field from src.photometry_core import _id_str, aperture_photometry_batch, row_nanmax -from src.utils import header_get_first +from src.utils import header_get_first, obs_time_utc_iso _log = logging.getLogger("originstack") _CH = ("r", "g", "b") @@ -66,13 +66,16 @@ def _cropped_session_wcs_header(session_info, shape_hw, left, top): def _frame_time_iso(frame, session_info, j, n): - """ISO UTC timestamp for sub *j*: the frame's own DATE-OBS if present, - else interpolated from the session start + total duration.""" + """Naive-UTC ISO timestamp for sub *j*: the frame's own DATE-OBS (with + its TIMEZONE applied -- Origin stamps local time) if present, else + interpolated from the session start + total duration. Always offset-free: + astropy Time cannot parse info.json's '-0700', which blanked every MJD.""" hdr = getattr(frame, "header", {}) or {} - own = header_get_first(hdr, ("DATE-OBS", "DATE_OBS", "DATEOBS"), cast=str) + own = obs_time_utc_iso(hdr) if own: return own - start = getattr(session_info, "date_time", None) if session_info else None + start = obs_time_utc_iso(fallback=getattr(session_info, "date_time", None) + if session_info else None) dur_ms = getattr(session_info, "total_duration_ms", None) if session_info else None if start and dur_ms and n > 1: try: diff --git a/src/pipeline.py b/src/pipeline.py index 8d141df..ba13297 100644 --- a/src/pipeline.py +++ b/src/pipeline.py @@ -34,7 +34,7 @@ from src.postprocess import postprocess_stack from src.registration import run_registration_phase, select_reference_frame from src.stacking import run_stacking_phase -from src.utils import format_time, get_memory_usage_mb, header_get_first, print_phase, safe_print +from src.utils import format_time, get_memory_usage_mb, obs_time_utc_iso, print_phase, safe_print _log = logging.getLogger("originstack") @@ -172,6 +172,47 @@ def _export_frame_jpegs(final: List[FrameInfo], final_indices: List[int], safe_print(f" Exported {len(final)} frame JPEGs → {export_dir}") +def _infer_target_and_advise(final: List[FrameInfo], args, directory: str, + use_simbad: bool) -> None: + """Target inference, the originvision session prior, then the + auto-advisor: the sequence every path into Phases 2-4 must run, fresh or + resumed from any checkpoint phase. + + Inference always runs (the result is written to the FITS header). With + --originvision + --auto, a fast 3-frame sample + (``sample_session_priors``, not the full-session scorer) feeds the same + prior_type/prior_confidence boost as the metadata inference, taking over + only when it is more confident -- e.g. a dense star field where metadata + gave nothing but originvision still recognizes the galaxy -- and stashes + a defect flag for ``_apply_quality_settings``' defensive nudge. + """ + from src.target_inference import infer_target_from_metadata + _si = getattr(args, '_session_info', None) + name, ttype, conf, src = infer_target_from_metadata( + directory, final, use_simbad=use_simbad, + session_name=_si.object_name if _si else None) + args._inferred_target = name + args._inferred_type = ttype + args._inferred_confidence = conf + args._inferred_source = src + if name and ttype and ttype != 'unknown': + safe_print(f"\n Target: {name} [{ttype.replace('_', ' ').title()}] " + f"conf={conf:.0%} source={src}") + + prior_type, prior_conf = ttype, conf + if getattr(args, 'originvision', False) and getattr(args, 'auto', False): + from src.originvision import map_originvision_category, sample_session_priors + result = sample_session_priors(final, args) + if result: + args._originvision_defect_flagged = result.get('defect_flagged', False) + mapped_type = map_originvision_category(result.get('category')) + mapped_conf = result.get('category_confidence', 0.0) + if mapped_type and mapped_conf > prior_conf: + prior_type, prior_conf = mapped_type, mapped_conf + + _run_auto_advisor(final, args, prior_type=prior_type, prior_confidence=prior_conf) + + def _run_auto_advisor(final: List[FrameInfo], args, prior_type: Optional[str] = None, prior_confidence: float = 0.0) -> None: @@ -461,17 +502,7 @@ def stack_target(frames: List[FrameInfo], output_path: str, args: argparse.Names # runs with bare CLI defaults — reintroducing colour casts. # Recover the metadata prior (folder/header/session) locally so the # classifier matches a fresh run; skip Simbad to avoid a network hit. - if getattr(args, 'auto', False): - from src.target_inference import infer_target_from_metadata - _r_name, _r_type, _r_conf, _r_src = infer_target_from_metadata( - _directory, final, use_simbad=False, - session_name=_session_info.object_name if _session_info else None) - if _r_name and _r_type and _r_type != 'unknown': - safe_print(f"\n Target: {_r_name} " - f"[{_r_type.replace('_', ' ').title()}] " - f"conf={_r_conf:.0%} source={_r_src}") - _run_auto_advisor(final, args, - prior_type=_r_type, prior_confidence=_r_conf) + _infer_target_and_advise(final, args, _directory, use_simbad=False) if getattr(args, 'photometry_timeseries', False): safe_print("\n NOTE: --photometry-timeseries needs the registered " @@ -561,6 +592,10 @@ def stack_target(frames: List[FrameInfo], output_path: str, args: argparse.Names rgb_shape=rgb_shape, lum_shape=lum_shape) stats.quality_time = 0.0 + # Same sequence as the fresh path and the phase-3 resume: + # without it a run resumed at phase 1 or 2 registered, stacked + # and post-processed with bare CLI defaults (no --auto presets). + _infer_target_and_advise(final, args, _directory, use_simbad=False) else: print_phase(1, "Processing & Quality Analysis") phase_start = time.time() @@ -623,52 +658,10 @@ def stack_target(frames: List[FrameInfo], output_path: str, args: argparse.Names # Save checkpoint after phase 1 save_checkpoint(output_path, phase=1, lights=lights, final=final, stats=stats) - # Target inference from metadata (always runs; result written to header) - from src.target_inference import infer_target_from_metadata - _si = getattr(args, '_session_info', None) - _inferred_name, _inferred_type, _inferred_conf, _inferred_src = \ - infer_target_from_metadata( - _directory, - final, - use_simbad=not getattr(args, 'offline', False), - session_name=_si.object_name if _si else None, - ) - args._inferred_target = _inferred_name - args._inferred_type = _inferred_type - args._inferred_confidence = _inferred_conf - args._inferred_source = _inferred_src - if _inferred_name and _inferred_type and _inferred_type != 'unknown': - safe_print( - f"\n Target: {_inferred_name} " - f"[{_inferred_type.replace('_', ' ').title()}] " - f"conf={_inferred_conf:.0%} source={_inferred_src}" - ) - - # originvision session-level signal (--originvision + --auto): a - # fast sample (src/originvision.py::sample_session_priors, NOT - # the full-session score_lights_with_originvision below) feeds - # the same prior_type/prior_confidence boost mechanism the - # metadata/SIMBAD inference above already uses -- only - # takes over when it's more confident than that prior (e.g. - # a dense star field where metadata gave nothing but - # originvision still recognizes the galaxy). Also stashes a - # defect flag for _apply_quality_settings' defensive nudge. - _prior_type, _prior_conf = _inferred_type, _inferred_conf - if getattr(args, 'originvision', False) and getattr(args, 'auto', False): - from src.originvision import map_originvision_category, sample_session_priors - _originvision_result = sample_session_priors(final, args) - if _originvision_result: - args._originvision_defect_flagged = _originvision_result.get( - 'defect_flagged', False) - _mapped_type = map_originvision_category(_originvision_result.get('category')) - _mapped_conf = _originvision_result.get('category_confidence', 0.0) - if _mapped_type and _mapped_conf > _prior_conf: - _prior_type, _prior_conf = _mapped_type, _mapped_conf - - # Heuristic auto-advisor - _run_auto_advisor(final, args, - prior_type=_prior_type, - prior_confidence=_prior_conf) + # Target inference + originvision prior + auto-advisor. + _infer_target_and_advise( + final, args, _directory, + use_simbad=not getattr(args, 'offline', False)) if not final: print('\n ERROR: No accepted frames after checkpoint restore!') @@ -1198,8 +1191,7 @@ def _psf_fallback_suffix() -> str: and getattr(_si_disp, 'has_gps', False) and getattr(_si_disp, 'has_wcs', False)): from src.observing_geometry import zenith_angle_deg - _t = header_get_first(hdu.header, ('DATE-OBS', 'DATE_OBS', 'DATEOBS'), - cast=str) or _si_disp.date_time + _t = obs_time_utc_iso(hdu.header, fallback=_si_disp.date_time) if _t: _za = zenith_angle_deg( math.degrees(_si_disp.ra_rad), math.degrees(_si_disp.dec_rad), diff --git a/src/psf_deconvolution.py b/src/psf_deconvolution.py index d5d3677..00463b6 100644 --- a/src/psf_deconvolution.py +++ b/src/psf_deconvolution.py @@ -228,8 +228,7 @@ def make_synthetic_psf(fwhm: float, psf_size: int = None, model: str = 'gaussian def estimate_psf_blind(img: np.ndarray, star_positions, - psf_size: int = None, - iterations: int = None) -> Tuple[Optional[np.ndarray], float]: + psf_size: int = None) -> Tuple[Optional[np.ndarray], float]: """Empirical (model-free) PSF estimation by stacking bright star cutouts. Instead of fitting a parametric model, this function extracts cutouts @@ -238,17 +237,18 @@ def estimate_psf_blind(img: np.ndarray, star_positions, *actual* on-sky PSF shape, including any asymmetry, coma, or tracking errors that a Gaussian/Moffat model cannot capture. - The ``iterations`` parameter optionally runs Richardson-Lucy blind - deconvolution update steps on top of the median-stack PSF to sharpen - the estimate. Each step refines the PSF by computing: - - psf_new = psf * correlate(img / (psf * img), img_reversed) + There used to be an optional "blind RL refinement" on top (on by + default, 8 iterations). It was not an RL PSF update: it correlated the + whole image with the ratio map and blended a 30% patch of that in each + step, flattening the PSF into a broad plateau (peak 0.068 -> 0.011, 93% + -> 22% of the energy within 3 px on a sigma=1.5 synthetic field). A + correct blind update needs an alternating object estimate; the median + stack is already the measured PSF, so the refinement was removed. Args: img: Float32 stacked image (H, W, 3). star_positions: Source table from detect_stars_auto. psf_size: Output kernel side length (default Config.RL_PSF_SIZE). - iterations: RL blind update iterations (0 = median stack only). Returns: (psf_kernel, fwhm_pixels) — (None, 0.0) if estimation fails. @@ -259,8 +259,6 @@ def estimate_psf_blind(img: np.ndarray, star_positions, if psf_size is None: psf_size = Config.RL_PSF_SIZE - if iterations is None: - iterations = Config.BLIND_PSF_ITERATIONS H, W = img.shape[:2] lum = (img if img.ndim == 2 @@ -318,35 +316,9 @@ def estimate_psf_blind(img: np.ndarray, star_positions, psf = np.maximum(psf, 0.0) psf /= psf.sum() - # Optional blind RL refinement (pure scipy.signal, no skimage involved) - if iterations > 0: - from scipy.signal import fftconvolve - lum_pos = lum - lum.min() + 1e-6 - psf_est = psf.copy() - for _ in range(iterations): - conv = fftconvolve(lum_pos, psf_est, mode='same') - ratio = lum_pos / (conv + 1e-12) - update = fftconvolve(ratio, psf_est[::-1, ::-1], mode='same') - # PSF update: cross-correlate ratio with lum - psf_update = fftconvolve(ratio[::-1, ::-1], lum_pos, mode='same') - # Trim to psf_size - cy, cx = np.unravel_index(np.argmax(psf_update), psf_update.shape) - y0 = max(0, cy - half) - x0 = max(0, cx - half) - patch = psf_update[y0:y0 + psf_size, x0:x0 + psf_size] - if patch.shape == (psf_size, psf_size): - patch = np.maximum(patch, 0.0) - s = patch.sum() - if s > 1e-12: - psf_est = 0.7 * psf_est + 0.3 * patch / s - - psf = np.maximum(psf_est, 0.0) - psf /= psf.sum() - # Radial apodization — force the kernel to zero at its edge. # An empirical PSF from stacked star cutouts retains non-zero energy in the - # square corners (star wings, residual noise), and blind RL refinement can - # inject blocky off-centre structure. FFT deconvolution with such a + # square corners (star wings, residual noise). FFT deconvolution with such a # hard-edged square kernel produces square ringing ("boxes") around every # point source. A Tukey (flat-core, cosine-taper) radial window keeps the # PSF core/wings intact while tapering the outer edge smoothly to zero, diff --git a/src/registration.py b/src/registration.py index f47c51a..8bba8cd 100644 --- a/src/registration.py +++ b/src/registration.py @@ -73,11 +73,12 @@ def match_stars_affine(ref_positions: Optional[Any], img_positions: Optional[Any for s in img_positions[:max_stars]]) # Shift img points by initial estimate for better matching. - # Initial shift is the amount needed to move 'img' to align with 'ref'. - # So img_points = ref_points + shift. To find correspondence, we map - # img_points back to ref space: img_points - shift. + # initial_shift is calculate_shift's (dy, dx): the shift that moves 'img' + # onto 'ref' (ndimage.shift(img, s) aligns it), so ref_points = + # img_points + shift. Subtracting it (as this once did) predicted every + # star 2*|shift| away, so any seed over ~AFFINE_MATCH_RADIUS/2 px failed. shift_vec = np.array([initial_shift[1], initial_shift[0]]) # [sx, sy] - img_pts_shifted = img_pts - shift_vec + img_pts_shifted = img_pts + shift_vec tree = cKDTree(ref_pts) distances, indices = tree.query(img_pts_shifted, k=1) diff --git a/src/stacking.py b/src/stacking.py index 37c555f..2da939b 100644 --- a/src/stacking.py +++ b/src/stacking.py @@ -927,6 +927,11 @@ def _process_tile(coords): return result +# Accepted samples a pixel needs before online_sigma_clip_fold_frame starts +# rejecting. Mirrored by ONLINE_CLIP_MIN_SAMPLES in ext/astro_native/src/lib.rs. +ONLINE_CLIP_MIN_SAMPLES = 3.0 + + def online_sigma_clip_seed_burnin(burn_stack: np.ndarray, coverage: np.ndarray, sigma: float = 3.0 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, int]: @@ -1028,10 +1033,14 @@ def online_sigma_clip_fold_frame(mean: np.ndarray, m2: np.ndarray, n_acc: np.nda # zero-variance estimate, which would reject it (and every sample after # it) forever. unseeded = (n_acc <= 0.0) & covered + # Below ONLINE_CLIP_MIN_SAMPLES there is no spread to clip against (one + # sample has m2 = 0, so every later sample failed the test and the pixel + # froze on a single frame): accept unconditionally until then. + warming = n_acc < ONLINE_CLIP_MIN_SAMPLES var_est = m2 / np.maximum(n_acc, 1.0) std_est = np.sqrt(np.maximum(var_est, 1e-12)) - accept = covered & (unseeded | (np.abs(x - mean) <= sigma * std_est)) + accept = covered & (unseeded | warming | (np.abs(x - mean) <= sigma * std_est)) n_acc_new = np.where(unseeded, 1.0, n_acc + accept) delta = x - mean @@ -1989,6 +1998,11 @@ def _combine_one(key: Tuple[str, int]) -> Tuple[Tuple[str, int], np.ndarray]: return result +# Stack methods run_stacking_phase's patch-weighted combine reproduces: plain +# mean, or a rejection method whose mask it builds before weighting. +_PATCH_WEIGHTED_METHODS = ('mean', 'sigma_clip', 'winsorized', 'percentile', 'esd') + + def run_stacking_phase( final: List[FrameInfo], final_indices: List[int], @@ -2166,7 +2180,16 @@ def _align_one(j): # aligned space (not full-res maps): both combine paths sample them # bilinearly at full-frame coordinates, so N full-resolution weight # maps (~5 GB at 200+ frames) are never materialised. - if quality_maps is not None and len(quality_maps) == n_final: + # Only for methods the patch path reproduces (mean, or a rejection + # method it builds a mask for). median/linear_fit/ivw/wavelet used to + # fall in here too and silently became an unrejected weighted mean -- + # --auto turns patch weighting on at 15+ frames, so a requested + # median kept every cosmic ray and trail while the log said "median". + _patch_ok = args.stack_method in _PATCH_WEIGHTED_METHODS + if quality_maps is not None and len(quality_maps) == n_final and not _patch_ok: + safe_print(f" NOTE: patch weighting is not available for --stack-method " + f"{args.stack_method}; using {args.stack_method} without it") + if quality_maps is not None and len(quality_maps) == n_final and _patch_ok: print(f" Patch-weighted mean combine ({n_final} frames × {top},{bottom},{left},{right} crop)...") qgrids = np.ascontiguousarray(np.stack(quality_maps), dtype=np.float32) qgrid_geom = (float(H), float(W), float(top), float(left)) diff --git a/src/utils.py b/src/utils.py index 1888b54..6a087ee 100644 --- a/src/utils.py +++ b/src/utils.py @@ -3,6 +3,7 @@ import logging import os +import re import sys from typing import Optional @@ -314,6 +315,40 @@ def parse_timestamp(when: str): return dt +# A FITS TIMEZONE value that can be appended to DATE-OBS: '-0700', '+05:30'. +# Anything else ('PDT', 'US/Pacific') would make an unparseable timestamp. +TZ_OFFSET_RE = re.compile(r'[+-]\d{2}:?\d{2}') + +_DATE_OBS_KEYS = ('DATE-OBS', 'DATE_OBS', 'DATEOBS') + + +def obs_time_utc_iso(header=None, fallback=None) -> Optional[str]: + """Observation time as a naive-UTC ISO string, or None. + + Celestron Origin lights stamp ``DATE-OBS`` in *local* time with the + offset in a separate ``TIMEZONE`` keyword; ``parse_timestamp`` treats an + offset-less timestamp as UTC, so reading DATE-OBS alone put every + time-dependent result seven hours out (a zenith angle of 150 deg -- below + the horizon -- for a target at 72 deg). A valid TIMEZONE is appended when + DATE-OBS carries no offset of its own. ``fallback`` (e.g. an info.json + ``dateTime`` like ``2026-08-31T20:40:32-0700``) is used when the header + has no DATE-OBS. The result has no offset, so astropy ``Time`` parses it + too (it rejects ``-0700``). + """ + when = header_get_first(header, _DATE_OBS_KEYS, cast=str) if header is not None else None + if when: + when = when.strip() + tz = str(header.get('TIMEZONE', '') or '').strip() + has_time = len(when) > 10 # a bare date takes no offset + if (has_time and tz and TZ_OFFSET_RE.fullmatch(tz) + and not TZ_OFFSET_RE.search(when[10:]) and not when.endswith('Z')): + when += tz + else: + when = fallback + dt = parse_timestamp(when) if when else None + return dt.isoformat() if dt is not None else None + + def header_get_first(header, keys, cast=None, default=None): """First present, non-None value among ``keys`` in a FITS-header-like mapping (anything with ``.get``). With ``cast`` given, the value is run diff --git a/tests/test_auto_settings_blend.py b/tests/test_auto_settings_blend.py index 353ce4c..8e7bb49 100644 --- a/tests/test_auto_settings_blend.py +++ b/tests/test_auto_settings_blend.py @@ -158,3 +158,46 @@ def test_emission_nebula_is_the_dominant_weight(self): w = a._blend_weights(sig, prior_type='emission_nebula', prior_confidence=0.90) assert max(w, key=w.get) == 'emission_nebula' assert w['emission_nebula'] > 0.5 + + +def _parser_args(): + """Real parser defaults, so every preset attr has a baseline -- the + state _apply_dynamic_settings actually runs against.""" + from src.cli import build_parser + ns = build_parser().parse_args(['-d', 'x', '-o', 'y']) + ns._explicit_cli_dests = set() + return ns + + +class TestUnlistedAttrsKeepBaseline: + """A preset that does not list an attr votes for the baseline. Before + this, a value only one preset listed was applied at full strength to + every target (renormalized over the listing presets alone).""" + + def test_star_field_does_not_inherit_single_preset_settings(self): + sig = dict(a._TYPE_ANCHORS['star_field']) + sig['fwhm'] = 3.0 + base = _parser_args() + args = _parser_args() + a._apply_dynamic_settings(sig, a._blend_weights(sig), args) + # galaxy-only / nebula-only settings stay at (or within a few percent + # of the preset-vs-default gap from) their defaults + for attr in ('denoise_aniso', 'masked_correlation', 'pre_gradient_removal'): + assert getattr(args, attr) == getattr(base, attr), attr + assert args.chroma_nr_large_sigma < 2.0 # galaxy-only 50; was 50 + assert args.ghs_hp == pytest.approx(base.ghs_hp, abs=0.005) # planetary-only 0.92 + assert abs(args.aniso_iterations - base.aniso_iterations) <= 1 + + @pytest.mark.parametrize("ttype", list(a._TYPE_ANCHORS.keys())) + def test_anchor_still_gets_its_own_values_with_real_baselines(self, ttype): + sig = dict(a._TYPE_ANCHORS[ttype]) + sig['fwhm'] = 3.0 + args = _parser_args() + a._apply_dynamic_settings(sig, a._blend_weights(sig), args) + for attr, expected in a._TARGET_SETTINGS.get(ttype, []): + got = getattr(args, attr) + if isinstance(expected, bool): + assert got == expected, f"{ttype}.{attr}: got {got}, want {expected}" + else: + tol = max(0.15 * abs(expected), 0.1) + assert abs(got - expected) < tol, f"{ttype}.{attr}: got {got}, want ~{expected}" diff --git a/tests/test_blind_psf.py b/tests/test_blind_psf.py new file mode 100644 index 0000000..8ad0ed6 --- /dev/null +++ b/tests/test_blind_psf.py @@ -0,0 +1,31 @@ +"""estimate_psf_blind must return the stars' own PSF. Its old default +'blind RL refinement' blended a patch of a whole-image correlation into the +estimate every iteration, flattening a sigma=1.5 px Gaussian PSF to a broad +plateau (peak 0.068 -> 0.011) -- and --auto enables the blind PSF for most +--deconvolve runs.""" +import numpy as np +from astropy.table import Table + +from src.psf_deconvolution import estimate_psf_blind + + +def test_blind_psf_matches_the_star_profile(): + rng = np.random.default_rng(7) + h, w, sigma = 256, 256, 1.5 + yy, xx = np.mgrid[:h, :w] + xs, ys = rng.uniform(20, w - 20, 60), rng.uniform(20, h - 20, 60) + fluxes = rng.uniform(2000, 6000, 60) + lum = np.full((h, w), 100.0) + for x, y, f in zip(xs, ys, fluxes): + lum += f / (2 * np.pi * sigma ** 2) * np.exp(-((xx - x) ** 2 + (yy - y) ** 2) / (2 * sigma ** 2)) + lum += rng.normal(0, 1.0, lum.shape) + img = np.repeat(lum[:, :, None], 3, axis=2).astype(np.float32) + stars = Table({'xcentroid': xs, 'ycentroid': ys, 'flux': fluxes}) + + psf, _ = estimate_psf_blind(img, stars) + assert psf is not None + c = psf.shape[0] // 2 + py, px = np.mgrid[:psf.shape[0], :psf.shape[1]] + within3 = psf[(py - c) ** 2 + (px - c) ** 2 <= 9].sum() / psf.sum() + assert within3 > 0.7 # ideal ~0.86 continuous, 0.81 sampled; was 0.22 + assert psf.max() > 0.8 / (2 * np.pi * sigma ** 2) # ideal peak 0.071; was 0.011 diff --git a/tests/test_cli_target_loop.py b/tests/test_cli_target_loop.py new file mode 100644 index 0000000..ac96006 --- /dev/null +++ b/tests/test_cli_target_loop.py @@ -0,0 +1,44 @@ +"""process_directory's per-target loop must hand each target its own +directory. The calibration-analysis block used to assign the master dark to +``d`` -- the loop variable holding the target directory -- so every session +with darks got ``args._input_directory`` = the dark array, and pipeline.py +silently fell back to the parent folder for info.json and target inference.""" +import argparse +import os +import tempfile +from unittest import mock + +import numpy as np + +from src import cli + + +def test_each_target_gets_its_own_directory_when_darks_exist(): + seen = [] + + def fake_stack_target(frames, outp, args, masters, stats): + seen.append(args._input_directory) + return None + + light = mock.Mock(header={'NAXIS1': 8, 'NAXIS2': 8}) + dark = mock.Mock(header={'NAXIS1': 8, 'NAXIS2': 8, 'EXPTIME': 10.0}) + with tempfile.TemporaryDirectory() as root: + dirs = [os.path.join(root, f"session{i}") for i in range(2)] + for p in dirs: + os.mkdir(p) + args = argparse.Namespace( + skip_step=[], hierarchical=True, mosaic=False, combine_sessions=False, + dry_run=False, health_check=False, preset=None, verbose=False, + stack_method='auto', _explicit_cli_dests=set()) + with mock.patch.object(cli, 'discover_frames', return_value={ + 'light': [light], 'dark': [dark], 'flat': [], 'bias': []}), \ + mock.patch.object(cli, '_load_calibration_dir', return_value={ + 'dark': [], 'flat': [], 'bias': []}), \ + mock.patch.object(cli, 'group_lights_by_filter', + side_effect=lambda lights: {'L': lights}), \ + mock.patch.object(cli, '_build_masters', side_effect=lambda f, s, a: { + 'dark': np.full((8, 8), 50.0, np.float32), 'dark_exptime': 10.0}), \ + mock.patch.object(cli, 'stack_target', side_effect=fake_stack_target): + cli.process_directory(root, os.path.join(root, 'out.fits'), args) + + assert seen == sorted(dirs) diff --git a/tests/test_difference_imaging.py b/tests/test_difference_imaging.py index a785cc1..a4f4c68 100644 --- a/tests/test_difference_imaging.py +++ b/tests/test_difference_imaging.py @@ -701,3 +701,21 @@ def test_matches_the_reference_greedy_algorithm_on_ties_and_nans(self): got = [(t.y, t.x) for t in detect_transients(s, 4.0, sep, cap)] self.assertEqual(got, expected, f"trial {trial}") + + +@pytest.mark.parametrize('flux_new', [0.3, 1.0, 3.0]) +def test_score_corr_is_unit_variance_on_pure_noise_at_any_flux_ratio(flux_new): + """S_corr is a significance only if it has unit variance on noise. The + matched-filter kernels once had their flux powers swapped -- invisible at + flux_new == 1, but S_corr std was 0.40-0.65 at ratios of 3 and 0.3, so a + '5 sigma' threshold really sat at 8-12 sigma and real transients were + missed.""" + rng = np.random.default_rng(11) + shape = (256, 256) + new = rng.normal(0.0, 5.0, shape) + ref = rng.normal(0.0, 3.0, shape) + res = zogy(new, ref, _gaussian_psf(25, 3.0), _gaussian_psf(25, 4.0), + sigma_new=5.0, sigma_ref=3.0, flux_new=flux_new, flux_ref=1.0, + var_new=np.full(shape, 25.0), var_ref=np.full(shape, 9.0)) + interior = res.score_corr[32:-32, 32:-32] + assert float(np.std(interior)) == pytest.approx(1.0, abs=0.08) diff --git a/tests/test_gpu_masters_cache.py b/tests/test_gpu_masters_cache.py new file mode 100644 index 0000000..7f7f064 --- /dev/null +++ b/tests/test_gpu_masters_cache.py @@ -0,0 +1,54 @@ +"""With --use-gpu, the GPU-side copies of the calibration masters were uploaded +only the first time a frame shape was seen and keyed on (shape, dtype), so the +next target in the same process -- a hierarchical run, or a second desktop-app +job -- with the same sensor was calibrated with the previous target's dark and +flat, even when it had no dark at all. A fake GPU (numpy as xp) drives the real +_process_single_frame GPU branch.""" +import numpy as np +import pytest +from astropy.io import fits + +import src.frame_processor as fp + + +class _FakeGpu: + active = True + xp = np + + def is_oom(self, exc): + return False + + def free_pool(self): + pass + + def to_host(self, arr): + return arr + + +@pytest.fixture +def fake_gpu(monkeypatch): + monkeypatch.setattr(fp, 'get_gpu', lambda: _FakeGpu()) + monkeypatch.setattr(fp, '_probe_gpu_calibration', lambda h, w, g, m: True) + monkeypatch.setattr(fp, '_gpu_calib_cache', {}) + monkeypatch.setattr(fp, '_gpu_masters', {}) + monkeypatch.setattr(fp, '_gpu_masters_sig', None) + + +def test_next_target_does_not_reuse_previous_gpu_masters(tmp_path, fake_gpu): + h = w = 64 + light = (1000 + np.random.default_rng(1).normal(0, 5, (h, w))).astype(np.float32) + path = str(tmp_path / 'l.fits') + fits.writeto(path, light, fits.Header({'EXPTIME': 10.0, 'BAYERPAT': 'RGGB'})) + + with_dark = {'bias': None, 'dark': np.full((h, w), 300.0, np.float32), 'flat': None, + 'dark_exptime': 10.0} + no_dark = {'bias': None, 'dark': None, 'flat': None, 'dark_exptime': None} + other_dark = {'bias': None, 'dark': np.full((h, w), 100.0, np.float32), 'flat': None, + 'dark_exptime': 10.0} + + means = [float(fp._process_single_frame(path, {}, m, 'malvar', 'none', + skip_quality=True)['rgb'].mean()) + for m in (with_dark, no_dark, other_dark)] + assert means[0] == pytest.approx(700.0, abs=5.0) + assert means[1] == pytest.approx(1000.0, abs=5.0) # was 700: target A's dark + assert means[2] == pytest.approx(900.0, abs=5.0) # same shape/dtype as A's dark diff --git a/tests/test_match_stars_affine_seed.py b/tests/test_match_stars_affine_seed.py new file mode 100644 index 0000000..3022de8 --- /dev/null +++ b/tests/test_match_stars_affine_seed.py @@ -0,0 +1,44 @@ +"""match_stars_affine's seed is calculate_shift's (dy, dx) -- the shift that +moves the frame onto the reference, i.e. ref = img + shift. It used to +subtract the seed, predicting every star 2*|shift| away, so any seed larger +than ~AFFINE_MATCH_RADIUS/2 failed and the frame dropped to translation-only +(field-rotation smear) whenever the blind matcher also failed.""" +import numpy as np +import pytest + +from src.registration import calculate_shift, match_stars_affine + + +def _table(xy): + return [{'xcentroid': float(x), 'ycentroid': float(y)} for x, y in xy] + + +def test_seed_from_calculate_shift_recovers_rotation_and_shift(): + rng = np.random.default_rng(1) + ref_xy = np.column_stack([rng.uniform(40, 560, 60), rng.uniform(40, 360, 60)]) + ang = np.radians(0.5) + rot = np.array([[np.cos(ang), -np.sin(ang)], [np.sin(ang), np.cos(ang)]]) + img_xy = ref_xy @ rot.T + [-22.0, 37.0] # frame drifted 37 px down, 22 left + # The seed calculate_shift would return: (dy, dx) with ref = img + seed. + seed = tuple((ref_xy - img_xy).mean(axis=0)[::-1]) + tf = match_stars_affine(_table(ref_xy), _table(img_xy), initial_shift=seed) + assert tf is not None + mapped = img_xy @ tf.params[:2, :2].T + tf.params[:2, 2] + np.testing.assert_allclose(mapped, ref_xy, atol=1e-6) + + +def test_seed_convention_matches_calculate_shift(): + """Pins the convention the fix relies on, on rendered images.""" + h, w = 200, 300 + yy, xx = np.mgrid[:h, :w] + ref_xy = np.array([[60.0, 50.0], [200.0, 80.0], [120.0, 150.0], [250.0, 170.0]]) + + def render(xy): + im = np.full((h, w), 10.0, np.float32) + for x, y in xy: + im += 500 * np.exp(-((xx - x) ** 2 + (yy - y) ** 2) / 4.5) + return im + + dy, dx = calculate_shift(render(ref_xy), render(ref_xy + [-12.0, 18.0]), verbose=False) + assert dy == pytest.approx(-18.0, abs=0.1) + assert dx == pytest.approx(12.0, abs=0.1) diff --git a/tests/test_obs_time_utc.py b/tests/test_obs_time_utc.py new file mode 100644 index 0000000..883582e --- /dev/null +++ b/tests/test_obs_time_utc.py @@ -0,0 +1,64 @@ +"""Origin lights stamp DATE-OBS in local time with the offset in a separate +TIMEZONE keyword. Reading DATE-OBS alone (parse_timestamp treats an +offset-less time as UTC) put every time-dependent result seven hours out: +a --fix-atmospheric-dispersion zenith angle of 150 deg (below the horizon, +so the correction went the wrong way) and light-curve MJDs/airmasses off.""" +import types + +import pytest + +from src.observing_geometry import zenith_angle_deg +from src.photometry_timeseries import _frame_time_iso, _to_mjd +from src.utils import obs_time_utc_iso + + +@pytest.mark.parametrize('hdr, expected', [ + ({'DATE-OBS': '2026-08-31T20:40:32', 'TIMEZONE': '-0700'}, '2026-09-01T03:40:32'), + ({'DATE-OBS': '2026-08-31T20:40:32', 'TIMEZONE': '-07:00'}, '2026-09-01T03:40:32'), + ({'DATE-OBS': '2026-08-31T20:40:32-0700', 'TIMEZONE': '-0700'}, '2026-09-01T03:40:32'), # not twice + ({'DATE-OBS': '2026-08-31T20:40:32Z', 'TIMEZONE': '-0700'}, '2026-08-31T20:40:32'), + ({'DATE-OBS': '2026-08-31T20:40:32', 'TIMEZONE': 'PDT'}, '2026-08-31T20:40:32'), # unusable + ({'DATE-OBS': '2026-08-31T20:40:32'}, '2026-08-31T20:40:32'), + ({'DATE-OBS': '2026-08-31', 'TIMEZONE': '-0700'}, '2026-08-31T00:00:00'), +]) +def test_obs_time_applies_timezone(hdr, expected): + assert obs_time_utc_iso(hdr) == expected + + +def test_obs_time_fallback_is_normalised_for_astropy(): + iso = obs_time_utc_iso({}, fallback='2026-08-31T20:40:32-0700') + assert iso == '2026-09-01T03:40:32' + assert _to_mjd(iso) == pytest.approx(61284.15315, abs=1e-4) + + +def test_frame_time_uses_timezone_and_session_fallback(): + frame = types.SimpleNamespace(header={'DATE-OBS': '2026-08-31T20:40:32', 'TIMEZONE': '-0700'}) + assert _frame_time_iso(frame, None, 0, 1) == '2026-09-01T03:40:32' + # No per-frame DATE-OBS: interpolate from info.json's offset-bearing start. + # This used to hand '...-0700' to astropy Time, fail, and blank every MJD. + si = types.SimpleNamespace(date_time='2026-08-31T20:40:32-0700', total_duration_ms=600_000) + iso = _frame_time_iso(types.SimpleNamespace(header={}), si, 1, 2) + assert iso.startswith('2026-09-01T03:50:32') + assert _to_mjd(iso) == _to_mjd(iso) # not NaN + + +def test_zenith_angle_is_none_below_the_horizon(): + # Antares from 34N, 118W: up in the evening (local), down 7 h later. + ra, dec, lat, lon = 247.35, -26.43, 34.0, -118.0 + assert zenith_angle_deg(ra, dec, lat, lon, 0.0, '2026-06-15T05:30:00') is not None + assert zenith_angle_deg(ra, dec, lat, lon, 0.0, '2026-06-15T15:30:00') is None + + +def test_stacked_header_keeps_timezone(): + from astropy.io import fits + + from src.cli import parse_args + from src.io_fits import populate_fits_header + from src.models import FrameInfo, ProcessingStats + f = FrameInfo(path='a.fits', type='light', + header={'DATE-OBS': '2026-08-31T20:40:32', 'TIMEZONE': '-0700', 'EXPTIME': 10.0}) + hdr = fits.Header() + populate_fits_header(hdr, [f], ProcessingStats(), + parse_args(['-d', 'x', '-o', 'y.fits', '--no-auto']), (3, 8, 8), + [(0.0, 0.0)], {'bias': None, 'dark': None, 'flat': None}) + assert obs_time_utc_iso(hdr) == '2026-09-01T03:40:32' diff --git a/tests/test_online_sigma_clip_warmup.py b/tests/test_online_sigma_clip_warmup.py new file mode 100644 index 0000000..3da1819 --- /dev/null +++ b/tests/test_online_sigma_clip_warmup.py @@ -0,0 +1,63 @@ +"""The streaming (--stream) sigma-clip accumulator must not freeze a pixel on +its first sample. With one accepted sample, M2 = 0 and the spread estimate +floored at 1e-6, so every later sample was rejected forever: a border pixel +only one burn-in frame covered stayed at that single noisy value for the whole +stream. Samples are now accepted unconditionally until a pixel holds +ONLINE_CLIP_MIN_SAMPLES of them. Checked on the numpy path and, when built, +the native kernel (which must agree).""" +import numpy as np +import pytest + +from src import stacking as st + + +def _paths(): + yield 'numpy' + if st.HAS_NATIVE: + yield 'native' + + +def _run(path, fn, *a, **kw): + had = st.HAS_NATIVE + st.HAS_NATIVE = had and path == 'native' + try: + return fn(*a, **kw) + finally: + st.HAS_NATIVE = had + + +@pytest.mark.parametrize('path', list(_paths())) +def test_pixel_seeded_by_one_frame_keeps_accumulating(path): + rng = np.random.default_rng(0) + k, h, w, c = 10, 8, 8, 3 + burn = (100 + rng.normal(0, 5, (k, h, w, c))).astype(np.float32) + cov = np.ones((k, h, w), np.float32) + cov[1:, 0, 0] = 0.0 # pixel (0,0): only the first burn-in frame + mean, m2, n_acc, _ = _run(path, st.online_sigma_clip_seed_burnin, burn, cov, sigma=3.0) + assert n_acc[0, 0, 0] == 1.0 + + full = np.ones((h, w), np.float32) + for _ in range(50): + frame = (100 + rng.normal(0, 5, (h, w, c))).astype(np.float32) + _run(path, st.online_sigma_clip_fold_frame, mean, m2, n_acc, frame, full, sigma=3.0) + # Was stuck at 1 for good. Not ~51: three samples give a noisy spread, so + # some early rejections are expected (40 seeds: mean 47, worst 36 of 51 + # at a 3-sample warm-up, vs 59 of 60 for a burn-in-seeded pixel). + assert n_acc[0, 0, 0] > 15 + assert abs(mean[0, 0, 0] - 100.0) < 3.0 + + +@pytest.mark.parametrize('path', list(_paths())) +def test_outliers_are_still_rejected_once_warmed_up(path): + rng = np.random.default_rng(1) + k, h, w, c = 10, 8, 8, 3 + burn = (100 + rng.normal(0, 5, (k, h, w, c))).astype(np.float32) + mean, m2, n_acc, _ = _run(path, st.online_sigma_clip_seed_burnin, burn, + np.ones((k, h, w), np.float32), sigma=3.0) + frame = (100 + rng.normal(0, 5, (h, w, c))).astype(np.float32) + frame[4, 4, :] = 5000.0 # cosmic ray + before = n_acc[4, 4, 0] + _run(path, st.online_sigma_clip_fold_frame, mean, m2, n_acc, frame, + np.ones((h, w), np.float32), sigma=3.0) + assert n_acc[4, 4, 0] == before + assert abs(mean[4, 4, 0] - 100.0) < 10.0 diff --git a/tests/test_photometry.py b/tests/test_photometry.py index a58108c..b67d112 100644 --- a/tests/test_photometry.py +++ b/tests/test_photometry.py @@ -268,3 +268,24 @@ def test_bright_unsaturated_cluster_not_flagged(tmp_path): n_sat = sum(int(r["saturated"]) for r in rows) # At most the single literal-max-pixel star; the bright cluster is fine. assert n_sat <= 1 + + +def test_pixel_coords_on_the_stacked_cube_header(): + """The stacked product is a (3, H, W) cube: its header is NAXIS=3, and a + bare WCS(header) is 3-axis, so _pixel_coords raised inside its broad + except and returned None -- --photometry then skipped every run with + 'could not project Gaia onto the image'.""" + from astropy.table import Table + + from src.color_calibrate import _field_radius_deg + from src.photometry_core import _pixel_coords + + hdr = fits.PrimaryHDU(np.zeros((3, 100, 120), np.float32)).header + hdr.update(CTYPE1='RA---TAN', CTYPE2='DEC--TAN', CRVAL1=10.0, CRVAL2=20.0, + CRPIX1=60.5, CRPIX2=50.5, CD1_1=-1e-3, CD1_2=0.0, CD2_1=0.0, CD2_2=1e-3) + assert hdr['NAXIS'] == 3 + xy = _pixel_coords(Table({'ra': [10.0], 'dec': [20.0]}), hdr) + assert xy is not None + np.testing.assert_allclose(xy[0], [59.5, 49.5], atol=1e-6) # CRPIX, 0-based + # corners are ~0.078 deg from centre; the old fallback returned 0.5 + assert 0.05 < _field_radius_deg(hdr) < 0.15 diff --git a/tests/test_resume_runs_advisor.py b/tests/test_resume_runs_advisor.py new file mode 100644 index 0000000..e096181 --- /dev/null +++ b/tests/test_resume_runs_advisor.py @@ -0,0 +1,48 @@ +"""Every way into Phases 2-4 must run target inference and the --auto advisor. +Only the phase-3 resume and the fresh path did: a run resumed from a phase-1 +or phase-2 checkpoint registered, stacked and post-processed with bare CLI +defaults (no galaxy mode, SCNR or deconvolution presets).""" +import json +import os +import tempfile +from unittest import mock + +import pytest + +from src import pipeline +from src.checkpoint import _checkpoint_dir, _ckpt_json_path +from src.models import FrameInfo, ProcessingStats +from tests.test_e2e import _create_synthetic_dataset, _make_minimal_args + + +def _lights(paths): + return [FrameInfo(path=p, type='light', + header={'BAYERPAT': 'RGGB', 'EXPTIME': 120.0, + 'NAXIS1': paths['W'], 'NAXIS2': paths['H']}) + for p in paths['light']] + + +@pytest.mark.parametrize('phase', [1, 2, 3]) +def test_resume_from_any_phase_runs_the_advisor(phase): + with tempfile.TemporaryDirectory() as tmp: + paths = _create_synthetic_dataset(tmp) + out = os.path.join(tmp, 'stacked.fits') + masters = {'dark': None, 'flat': None, 'bias': None, 'dark_exptime': None} + # First run leaves a phase-3 checkpoint behind. + first = _make_minimal_args(no_resume=False, keep_checkpoint=True) + assert pipeline.stack_target(_lights(paths), out, first, masters, ProcessingStats()) + ckpt = _ckpt_json_path(_checkpoint_dir(out)) + with open(ckpt) as fh: + state = json.load(fh) + state['phase'] = phase + with open(ckpt, 'w') as fh: + json.dump(state, fh) + + calls = [] + again = _make_minimal_args(no_resume=False, keep_checkpoint=True, auto=True) + with mock.patch.object(pipeline, '_run_auto_advisor', + side_effect=lambda final, args, **kw: calls.append(len(final))): + assert pipeline.stack_target(_lights(paths), out, again, masters, + ProcessingStats()) + assert calls == [len(paths['light'])] + assert hasattr(again, '_inferred_type') diff --git a/tests/test_spcc.py b/tests/test_spcc.py index 5c3786d..8f6638d 100644 --- a/tests/test_spcc.py +++ b/tests/test_spcc.py @@ -15,6 +15,7 @@ from unittest.mock import patch import numpy as np +import pytest from src.color_calibrate import ( _blackbody_spectrum, @@ -143,3 +144,43 @@ def test_falls_back_to_colorindex_without_teff_column(self): patch('src.color_calibrate._aperture_flux', return_value=fluxes): scales = fit_channel_scales_spcc(img, header, catalog, verbose=False) assert all(np.isfinite(s) for s in scales) + + +class TestScalesActuallyCorrectTheCast: + """Both fitters used to divide each star's ratio by its own channel median + before taking the median, which is always 1 -- so the scales were always + (1, 1, 1) and --color-calibrate never changed the image. Here the red + channel reads 1.6x too bright: the fitted scales must pull it back.""" + + def _run(self, fitter, with_teff): + from src.color_calibrate import fit_channel_scales + rng = np.random.default_rng(3) + n = 40 + g = rng.uniform(9.0, 12.0, n) + bp_rp = rng.uniform(0.3, 1.5, n) + cols = {'ra': np.zeros(n), 'dec': np.zeros(n), 'phot_g_mean_mag': g, + 'phot_bp_mean_mag': g + 0.5 * bp_rp, 'phot_rp_mean_mag': g - 0.5 * bp_rp} + if with_teff: + cols['teff_gspphot'] = np.full(n, np.nan) # colour-index fallback for all + catalog = _FakeTable(cols) + # Measured = the fitter's own expected fluxes x a per-channel gain. + from src.color_calibrate import _bp_rp_to_bv + bv = _bp_rp_to_bv(bp_rp) + expected = np.column_stack([10 ** (-0.4 * (g - 0.5 * bv)), 10 ** (-0.4 * g), + 10 ** (-0.4 * (g + bv))]) + fluxes = expected * np.array([1.6, 1.0, 1.0]) * 1e7 + fn = fit_channel_scales_spcc if fitter == 'spcc' else fit_channel_scales + with patch('src.color_calibrate._pixel_coords', + return_value=np.zeros((n, 2))), \ + patch('src.color_calibrate._aperture_flux', return_value=fluxes): + return fn(np.ones((64, 64, 3), np.float32), {}, catalog) + + def test_colorindex_scales_undo_a_red_cast(self): + r, g, b = self._run('colorindex', with_teff=False) + assert r / g == pytest.approx(1 / 1.6, rel=1e-6) + assert b / g == pytest.approx(1.0, rel=1e-6) + + def test_spcc_scales_undo_a_red_cast(self): + r, g, b = self._run('spcc', with_teff=True) + assert r / g == pytest.approx(1 / 1.6, rel=1e-6) + assert b / g == pytest.approx(1.0, rel=1e-6) diff --git a/tests/test_stacking_patch_methods.py b/tests/test_stacking_patch_methods.py new file mode 100644 index 0000000..bed2c7c --- /dev/null +++ b/tests/test_stacking_patch_methods.py @@ -0,0 +1,40 @@ +"""The patch-weighted combine only reproduces mean and the rejection methods +it builds a mask for. It used to run for every --stack-method whenever +quality maps existed (--auto turns patch weighting on at 15+ frames), so a +requested median/linear_fit/ivw/wavelet silently became an unrejected +weighted mean: a single cosmic ray survived at ~1/N of its amplitude.""" +import numpy as np +import pytest + +from src import stacking as st +from src.cli import parse_args +from src.models import FrameInfo, ProcessingStats + + +def _stack(method, capsys): + n, h, w, c = 12, 48, 48, 3 + args = parse_args(['-d', 'x', '-o', 'y.fits', '--stack-method', method, '--no-auto']) + rng = np.random.default_rng(0) + mem = (100 + rng.normal(0, 1, (n, h, w, c))).astype(np.float32) + mem[3, 24, 24, :] = 10000.0 # one cosmic ray in one frame + final = [FrameInfo(path=f'f{i}.fits', type='light', header={}, + metrics={'score': 1.0, 'noise': 1.0, 'fwhm': 3.0}) for i in range(n)] + qmaps = [np.ones((4, 4), np.float32) for _ in range(n)] + _, fits_stacked, top, _, left, _ = st.run_stacking_phase( + final, list(range(n)), mem, [(0.0, 0.0)] * n, [None] * n, h, w, c, args, + ProcessingStats(), quality_maps=qmaps) + return fits_stacked[24 - top, 24 - left], capsys.readouterr().out + + +@pytest.mark.parametrize('method', ['median', 'linear_fit']) +def test_non_patch_methods_keep_their_own_rejection(method, capsys): + px, out = _stack(method, capsys) + assert np.all(np.abs(px - 100.0) < 5.0), px # was ~925 + assert 'patch weighting is not available' in out + + +@pytest.mark.parametrize('method', ['sigma_clip', 'percentile']) +def test_rejection_methods_still_use_the_patch_path(method, capsys): + px, out = _stack(method, capsys) + assert np.all(np.abs(px - 100.0) < 5.0), px + assert 'Patch-weighted mean combine' in out diff --git a/tests/test_worker_masters.py b/tests/test_worker_masters.py new file mode 100644 index 0000000..27d0e5d --- /dev/null +++ b/tests/test_worker_masters.py @@ -0,0 +1,37 @@ +"""Pool workers must see the same masters as the sequential path, scalars +included. Only numpy arrays used to be copied into shared memory, so +``dark_exptime`` never reached a worker and a dark of a different exposure +was subtracted unscaled on the default (ProcessPool) path, silently.""" +import numpy as np +from astropy.io import fits + +import src.frame_processor as fp + + +def test_scalar_masters_reach_worker_and_scale_the_dark(tmp_path): + H = W = 64 + rng = np.random.default_rng(0) + # sky 1000 + 30 s of dark current (300); the master dark is 10 s (100) + light = (1300 + rng.normal(0, 5, (H, W))).astype(np.float32) + path = str(tmp_path / 'l.fits') + fits.writeto(path, light, fits.Header({'EXPTIME': 30.0, 'BAYERPAT': 'RGGB'})) + masters = {'bias': None, 'dark': np.full((H, W), 100.0, np.float32), + 'flat': None, 'hot_pixel_map': None, 'dark_exptime': 10.0} + + blocks, specs, scalars = fp._share_masters(masters) + saved = fp._worker_masters + try: + assert scalars == {'dark_exptime': 10.0} + fp._init_worker_shm(specs, scalar_masters=scalars) + assert fp._worker_masters['dark_exptime'] == 10.0 + pooled = fp._process_single_frame(path, {}, fp._worker_masters, 'malvar', 'none', + skip_quality=True) + seq = fp._process_single_frame(path, {'EXPTIME': 30.0}, masters, 'malvar', 'none', + skip_quality=True) + assert pooled['rgb'][..., 1].mean() == np.float32(seq['rgb'][..., 1].mean()) + assert abs(float(pooled['rgb'][..., 1].mean()) - 1000.0) < 5.0 + finally: + fp._worker_masters = saved + for shm in blocks: + shm.close() + shm.unlink()