cosmic op• Data kinds: images → images
• Call: import fullseye as fs; fs.ledger.cosmic_ray_reject_stack(frames, kappa=5.0, min_frames=3, read_sigma=None, gain=1.0) (to call the implementation directly, import astrostack; astrostack.cosmic_ray_reject_stack(frames, kappa=5.0, min_frames=3, read_sigma=None, gain=1.0); from the registry, opsastrostack.get("cosmic_ray_reject_stack"))
• Return value through the ledger: fullseye.ledger.cosmic_ray_reject_stack(...) returns only the declared out type images (the underlying function also returns auxiliary values). When you need what was dropped, use fullseye.ledger.cosmic_ray_reject_stack.raw(...) or call astrostack.cosmic_ray_reject_stack directly.
• Underlying return: (cleaned, masks) -> images
Cosmic-ray removal by comparing frames — it never hits the same place twice.
Cosmic rays are hard to detect in a single frame because "stars are peaked too", but with many aligned frames it becomes simple: a star sits on the same pixel every time, and a cosmic ray comes only once. So per pixel the median and the MAD along the frame axis are taken, and only frames with value > median + kappa * sigma are dropped and filled with the median (the low side is not dropped —— cosmic rays are always outliers in the adding direction, and dropping the low side would erase even missing pixels).
Below min_frames frames neither median nor MAD means anything, so it is refused (3 frames is the minimum: with 2 it is undecidable which one is the outlier). Frames must be aligned —— passed misaligned, stars "come only once" and it is the stars that vanish.
Measured: with few frames the MAD itself cannot be trusted. On a background of 8 frames the MAD estimate was 7.89 against a true σ of 9.22 (-14.5 %), and it scatters strongly from pixel to pixel, so an intended kappa=5 was effectively 4.3, and false positives came out 2.4 times the true positives (546 vs 227 pixels). Moreover the false positives were on the background, not on stars (the median ground truth of the false-positive pixels was 60 = the sky itself), so "because stars are peaked" does not explain it —— the cause is that the MAD of a small sample scatters strongly per pixel. 2 countermeasures are stacked:
1. apply the small-sample correction of _mad_correction (Croux & Rousseeuw 1992).
2. passing *read_sigma* makes the one noise model this module has, sigma = sqrt(median/gain + read_sigma^2), the floor of the scale. This is the same "Poisson(signal) + Gauss(readout)" that synth_starfield uses, not a second theory. With the floor, pixels whose MAD happened to come out small through sampling fluctuation no longer turn into cosmic rays.
Measured effect (128x128 / 25 stars / 10 cosmic rays per frame / 8 frames, true cosmic rays 227 pixels, kappa=5)::
read_sigma=None -> 594 px detected, false positives 367, precision 0.382
read_sigma=5.0 -> 228 px detected, false positives 1, precision 0.996
Recall is 1.000 in both cases (not a single pixel missed). The floor does not act in the direction of missing cosmic rays —— cosmic rays are tens of times the noise, so even if the floor rises by a few e- they still clear it.
Returns (cleaned, masks):
• cleaned —— list of length N, each (H, W) float64 (images vocabulary).
• masks —— (N, H, W) bool, True = pixel judged to be a cosmic ray.
Raises ValueError: if *frames* is not a list / tuple / the number of frames is less than *min_frames* / shapes do not match / *kappa* is non-positive / *gain* is non-positive / *read_sigma* is negative.
• Sample-data catalog (download URLs / licences) — 2-D uses skimage.data (BSD/public domain) plus synthetic images; 3-D lists download URLs for real data sources (Stanford, PDS, …).
• Operator provenance and references — the sources of the research/methods this op family came from.
• The canonical algorithm (author, year) and its uses are named in the family usage guide above.
• astro_stacking — py -3.11 examples/astro_stacking.py
images as input)lucky_select · sigma_clip_stack · drizzle_resample · align_frames
cosmic)*Provenance: astrostack.py — ASTROSTACK operator registry. This per-op note is generated by tools/opdocs.py md (do not hand-edit).*
© 2026 Kazufumi Furuse — Fullseye operator documentation. Licensed under Apache-2.0.