cosmicray_median now honors an input mask (#932, PR #979) in a deliberately narrow way: masked pixels are never flagged, are returned unchanged, and are excluded from the noise estimate. They are not excluded from the median filter itself, because scipy.ndimage.median_filter knows nothing about masks. So very bright masked pixels (saturated stars, bad columns) still pull up the local median of their unmasked neighbors, which can hide a cosmic ray next to them or, on a steep background, cause false positives.
What was tried and rejected
An earlier revision of #979 replaced masked pixels with a fill value (global unmasked mean), ran the median filter, refilled with that local median, and filtered again. Measured against an exact masked median (sliding_window_view + np.nanmedian):
- Fine on a flat background (≤ 0.4σ from exact for bad columns, bright blocks, 60% scattered masks; no mis-flags).
- Biased toward the global mean wherever a masked region is wider than ~
mbox/2 and the local background differs from it: 10 false positives at ~15σ next to a 21×11 masked block straddling a 30 ADU amplifier step, a missed 6σ cosmic ray beside a 21×21 block, and a correct pixel inside a masked region flagged at 15.7σ.
- Iterating the refill converges slowly (false flags 10→6→6→4→2→0 at 1/2/3/5/10/20 passes; each pass ~5.6 s on 2048² with
mbox=11).
Options
- Exact chunked
sliding_window_view + np.nanmedian with NaN at masked pixels: ~21 s on 2048²/mbox=11 in 64-row chunks (~0.13 GB temp) vs 5.6 s for one median_filter.
- Exact only near the mask: run (1) on the rows/cols within
binary_dilation(mask, mbox) and ndimage.median_filter everywhere else. Near current cost for typical masks, exact everywhere. Probably the right one.
- Leave as is and document (current state).
Whatever is chosen runs in numpy, so this also belongs with the CPU-only-operations policy in #935. The test fixture should include a sky level and a masked region wider than mbox/2 on a non-flat background; the flat, zero-mean fixtures from the earlier revision could not distinguish any of the approaches.
Follow-up from the review of #979.
cosmicray_mediannow honors an input mask (#932, PR #979) in a deliberately narrow way: masked pixels are never flagged, are returned unchanged, and are excluded from the noise estimate. They are not excluded from the median filter itself, becausescipy.ndimage.median_filterknows nothing about masks. So very bright masked pixels (saturated stars, bad columns) still pull up the local median of their unmasked neighbors, which can hide a cosmic ray next to them or, on a steep background, cause false positives.What was tried and rejected
An earlier revision of #979 replaced masked pixels with a fill value (global unmasked mean), ran the median filter, refilled with that local median, and filtered again. Measured against an exact masked median (
sliding_window_view+np.nanmedian):mbox/2and the local background differs from it: 10 false positives at ~15σ next to a 21×11 masked block straddling a 30 ADU amplifier step, a missed 6σ cosmic ray beside a 21×21 block, and a correct pixel inside a masked region flagged at 15.7σ.mbox=11).Options
sliding_window_view+np.nanmedianwith NaN at masked pixels: ~21 s on 2048²/mbox=11in 64-row chunks (~0.13 GB temp) vs 5.6 s for onemedian_filter.binary_dilation(mask, mbox)andndimage.median_filtereverywhere else. Near current cost for typical masks, exact everywhere. Probably the right one.Whatever is chosen runs in numpy, so this also belongs with the CPU-only-operations policy in #935. The test fixture should include a sky level and a masked region wider than
mbox/2on a non-flat background; the flat, zero-mean fixtures from the earlier revision could not distinguish any of the approaches.Follow-up from the review of #979.