Function detect_hot_pixels 

pub fn detect_hot_pixels(
    data: &Array3<f64>,
    k_mad: f64,
) -> Result<Array2<bool>, IoError>

Detect hot (railed / runaway) pixels via a two-stage criterion: a robust one-sided log-space median + k·MAD screen on per-pixel total counts (stage 1, global), confirmed by a local-neighborhood isolation test (stage 2).

Pipeline-integrity screen only — see the module-level “Pixel masks — pipeline integrity only” section. The algorithm:

1. totals[y, x] = Σ_t data[t, y, x].
2. The statistics sample is ln(totals) over totals > 0 pixels only — dead pixels are excluded before the median/MAD so ln(0) never enters and a large dead population cannot drag the median down. If no live pixels exist, every pixel is unflagged (all-false mask).
3. med = median(ln totals), mad = median(|ln totals − med|).
4. sigma = max(MAD_TO_SIGMA·mad, exp(−med/2)). The second term is the delta-method Poisson floor of ln N: Var[ln N] ≈ 1/N evaluated at N = exp(med) (medians commute with monotone maps, so exp(med) is the median total), giving σ_floor = 1/√exp(med) = exp(−med/2). The robust scale can never legitimately sit below counting noise; this guards mad == 0 on quantized low-count images without becoming a low-count screen — it only ever raises the threshold. Worked check: an image where most totals are 1.0 and some are 2.0 has mad = 0 and floor = 1, so the threshold is e^6 ≈ 403× the median — the 2-count pixels are NOT flagged, while a railed pixel still is.
5. Stage 1 (global): a pixel is a candidate iff totals > 0 && ln(total) > med + k_mad·sigma — upper tail only. A stuck-low pixel is indistinguishable from a low-count-alive pixel and is deliberately kept (masking it would be the banned low-count screen). Railed/always-max pixels are subsumed by the upper tail: no fixed saturation value exists after efficiency correction, so a saturation-constant test would be wrong anyway.
6. Stage 2 (local confirmation), iterated to a fixpoint: a candidate is flagged iff its total also exceeds HOT_LOCAL_FACTOR × the median of its 8-neighborhood reference sample, where each neighbor contributes its total if live (total > 0) and not yet flagged, contributes 0.0 if already flagged (a known defect cannot vouch for its neighbors — see below), and is omitted if dead (a dead pixel carries no scene information). Edge pixels use whatever neighbors exist. A candidate whose reference sample is empty (every neighbor dead — isolated live pixel in a dead field) keeps the global verdict; a candidate whose neighbors are mostly flagged defects likewise stays flagged (its reference median is 0). After each full pass over the fixed stage-1 candidate list, newly confirmed flags are applied and the pass repeats until a pass adds no new flag.

Fixpoint erosion of railed clusters

A single stage-2 pass misses the INTERIOR of a railed cluster ≥2 px wide: an interior pixel’s 8-neighbors are railed too, so its neighbor median is railed and the ratio test refutes the flag. Iterating to a fixpoint erodes such clusters from the boundary inward — once the cluster’s outermost pixels are flagged they stop vouching (each contributes 0.0 instead of its railed total), the reference median of the next ring drops back to the background level, and the next pass flags that ring. Contributing 0.0 (rather than omitting the flagged neighbor) is load-bearing: with omission, a 3-background + 3-railed reference sample has an even-count nereids_core::stats::median midpoint mid-gap (≈ railed/2), the ratio test reads ~2× and the erosion stalls; with the zero contribution the median stays on the background side and erosion completes. Erosion fully consumes clusters up to 3 px wide in their narrower dimension (point defects, 1-px lines, 2–3-px-wide blobs/segments — the physical shapes of railed detector defects) PROVIDED the cluster exposes at least one end cap or convex corner to normal-scene neighbors — erosion must seed somewhere. An EDGE-TO-EDGE railed band ≥2 px wide (both ends off-detector — spanning the full detector width or height) exposes neither: every interior band pixel keeps ≥5 railed of 8 neighbors, and even the on-detector-border band ends keep 3 railed of 5, so every neighbor median stays railed, no pixel ever seeds, and the band is NOT caught (test_detect_hot_pixels_edge_to_edge_2px_band_not_flagged_by_design). This is deliberate, not an oversight: a slit-aperture open beam produces a genuine full-width bright SCENE band that is pixel-for-pixel indistinguishable from such a defect, so a full-span row/column screen would mask it — re-introducing the exact bimodal failure stage 2 exists to prevent. Full-span detector pathologies of width ≥2 belong in a declared/file mask. (A full-span width-1 railed line IS caught — each of its pixels keeps ≥4 normal neighbors, test_detect_hot_pixels_full_railed_column_caught; and a ≥2-px band with even one end cap inside the detector is fully consumed from that cap, test_detect_hot_pixels_2px_band_one_end_on_detector_fully_caught.) A hard-edged railed rectangle ≥4 px wide keeps its interior (only its convex corners flag): it is pixel-for-pixel indistinguishable from a hard-edged bright scene region, which must survive (below).

Termination bound: flags are only ever added, and every pass except the last adds at least one, so at most height·width passes can do work; the loop is additionally capped at height·width passes to make the bound structural rather than reasoned. In practice the pass count is on the order of the defect-cluster radius (one pass for point defects and 1-px lines).

Why bright scene regions never erode

Erosion must seed at a bright-region boundary pixel. A boundary pixel of a contiguous bright scene region ≥2 px wide keeps ≥4 of its 8 neighbors on its own (bright) side for any straight or diagonal edge, so its reference median stays bright and its ratio is the scene gradient (≤2–3× across real edges) — far below the ≥10× HOT_LOCAL_FACTOR. Stage 1’s global cut additionally gates which pixels can seed at all. With no seed, the fixpoint is reached with zero flags in the region and it survives intact (test_detect_hot_pixels_large_psf_bright_region_not_eroded). Two documented, test-pinned exceptions — both physically rare on VENUS, where the detector PSF blurs real scene features over ≥2 px so ≥10× single-pixel contrast steps do not occur in scene:

• a width-1 bright line at ≥10× local contrast is spatially indistinguishable from a railed line and IS masked — the accepted trade-off for catching railed rows/columns (test_detect_hot_pixels_1px_bright_line_flagged_by_design);
• the single pixel at a sharp convex (90°) corner of a hard-edged ≥10× region sees only 3 same-side neighbors, its reference median falls on the dark side, and it is flagged (this predates the fixpoint); erosion does NOT propagate past it — the adjacent edge pixels keep ≥4 bright unflagged neighbors (test_detect_hot_pixels_hard_edged_bright_rectangle_corners_only).

The two stages encode complementary definitions of “hot”: stage 1 says statistically implausible for this image, stage 2 says spatially isolated, and a railed/hot pixel is a detector point defect — it must be both. Stage 2 exists because the global cut alone fails catastrophically on bimodal images: with a dark majority holding the median (a sample covering >50 % of the FOV, or an aperture-limited open beam), the MAD reflects only the dark population’s internal spread and the ENTIRE bright minority exceeds med + k·MAD. A contiguous bright REGION is scene, not a defect — masking it would reject statistically plausible pixels, the exact failure the module rules ban. A true point defect beats its neighbor median by ≥100× and is caught by both stages even inside a bright region and even as part of a railed row/column or a small railed cluster (see HOT_LOCAL_FACTOR).

Raw counts required

data must be raw detected counts (unscaled). The Poisson floor in step 4 assumes Var[N] = N; any prior scaling silently breaks that identity — down-scaling (proton-charge-normalized rates ≪ 1, per-pixel gain division) inflates the floor and can suppress real flags, while up-scaling (event weights > 1) deflates it below true counting noise. Run the detectors on unscaled counts and normalize afterwards; the GUI does exactly this (all three of its raw-counts paths — TIFF pair, HDF5 with open beam, HDF5 without open beam — pass the raw sample/open-beam stacks, before any normalization).

Arguments

• data — 3D raw-counts array with shape (n_tof, height, width).
• k_mad — Robust-σ multiplier for the stage-1 upper-tail cut; use HOT_PIXEL_K_MAD unless you have a reason not to.

Returns

2D boolean mask, shape (height, width). true = hot pixel.

Errors

Returns IoError::InvalidParameter if data contains a non-finite or negative value or has an empty TOF axis (shape[0] == 0), or if k_mad is not finite and positive.