From 28b0b929e76a29533e00457b8cc813230482744b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:07:02 +0200 Subject: [PATCH 01/14] test(ngmix): specify defect fill, epoch cuts and shear recovery near defects Defects (nonzero flag, zero weight, invalid RMS) must be zero-weighted and noise-filled under every BLEND_HANDLING, with neighbour pixels left raw under uberseg. The masked-fraction cut must count that same raw set, and an epoch with a defect closer than EPOCH_CENTRAL_DEFECT_RADIUS to the stamp centre must be dropped, at a strict boundary and at the configured radius. The runner must pass both cut options on, the HSM centroid must ignore raw defect values, each epoch must keep its own OFFSET, and every tile must log its epoch-cut tally. The science test recovers the full 2x2 response for columns, 3-px bleeds, single pixels and edge bands at the veto radius, on 0.3" and 0.5" galaxies through round and elliptical 0.7" PSFs, and requires |c1|, |c2| < 5e-4 and |m11|, |m22| < 1%. Flagged pixels in the simulation hold a hot value, so an unfilled defect shows up. A column two pixels inside the radius is the positive control. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01QAc1ywyd9Sbbra6po2iuaA --- tests/helpers/defect_response.py | 89 +++++ tests/module/test_ngmix_defect_fill.py | 532 +++++++++++++++++++++++++ tests/module/test_ngmix_uberseg.py | 69 ++-- tests/science/test_defect_veto.py | 104 +++++ 4 files changed, 755 insertions(+), 39 deletions(-) create mode 100644 tests/helpers/defect_response.py create mode 100644 tests/module/test_ngmix_defect_fill.py create mode 100644 tests/science/test_defect_veto.py diff --git a/tests/helpers/defect_response.py b/tests/helpers/defect_response.py new file mode 100644 index 000000000..c3e3ac157 --- /dev/null +++ b/tests/helpers/defect_response.py @@ -0,0 +1,89 @@ +"""Full-matrix metacal recovery for fixed detector defects.""" + +import numpy as np + +from shapepipe.modules.ngmix_package import ngmix as ngm +from shapepipe.testing.simulate import make_data +from tests.helpers.metacal_sim import build_stamp + + +def defect_response(bad, hlr=0.5, psf=0.7, seeds=range(4), + psf_shear=(0.0, 0.0), options=None, known_rms=True, + defect_value=1e3): + """Recover c and M = inverse(mean R) A - I with paired seed errors. + + Null pairs rotate the pixels by 90 degrees while pre-rotating the PSF + ellipticity oppositely, so the final PSF stays fixed in detector + coordinates. + This does NOT average away elliptical-PSF leakage. + + The flagged pixels ``bad`` hold ``defect_value`` (far above the galaxy's + peak) rather than sky, as a hot column or bleed would, so any defect value + that reaches metacal shows up as a bias. + """ + gamma = 0.02 + options = {} if options is None else options + samples = [] + for seed in seeds: + def arm(axis, sign, rotation=0): + shear = [0.0, 0.0] + if axis >= 0: + shear[axis] = sign * gamma + ps = tuple(v * (-1 if rotation else 1) for v in psf_shear) + data = list(make_data( + rng=np.random.RandomState(seed + 100), shear=shear, + psf_shear=ps, noise=1e-4, n_epochs=1, img_size=51, + gal_hlr=hlr, psf_fwhm=psf, return_centers=True, + )) + centre = data.pop()[0] + offset = np.array([centre.y - 26, centre.x - 26]) + if rotation: + data[0] = [np.rot90(a).copy() for a in data[0]] + data[1] = [np.rot90(a).copy() for a in data[1]] + offset = np.array([-offset[1], offset[0]]) + data[0] = [np.where(bad, defect_value, a) for a in data[0]] + data[4] = [bad.astype(np.int32)] + stamp = build_stamp(data) + stamp.offsets = [offset] + if known_rms: + stamp.bkg_rms = [np.full((51, 51), 1e-4)] + rng = np.random.RandomState(seed) + result, _, _ = ngm.do_ngmix_metacal( + stamp, ngm.get_prior(0.1857, rng), 1.0, rng, + centroid_source="wcs", **options, + ) + assert all(result[t]["flags"] == 0 for t in ngm.METACAL_TYPES) + e = np.asarray(result["noshear"]["g"]) + response = np.column_stack([ + (np.asarray(result[p]["g"]) - result[m]["g"]) / 0.02 + for p, m in (("1p", "1m"), ("2p", "2m")) + ]) + assert np.all(np.isfinite(e)) and np.all(np.isfinite(response)) + return e, response + null = [arm(-1, 0, k) for k in (0, 1)] + e0 = np.mean([a[0] for a in null], axis=0) + r0 = np.mean([a[1] for a in null], axis=0) + derivatives, responses = [], [] + for axis in (0, 1): + plus, minus = arm(axis, 1), arm(axis, -1) + derivatives.append((plus[0] - minus[0]) / (2 * gamma)) + responses.extend([plus[1], minus[1]]) + samples.append((e0, r0, np.column_stack(derivatives), + np.mean(responses, axis=0))) + e, r, a, rm = [np.array([s[i] for s in samples]) for i in range(4)] + assert np.linalg.svd(rm.mean(axis=0), compute_uv=False).min() > 0.1 + c = np.linalg.solve(r.mean(axis=0), e.mean(axis=0)) + matrix = np.linalg.solve(rm.mean(axis=0), a.mean(axis=0)) - np.eye(2) + draw = np.random.RandomState(91).randint(len(e), size=(1000, len(e))) + cb = np.linalg.solve(r[draw].mean(axis=1), + e[draw].mean(axis=1)[..., None])[..., 0] + mb = np.linalg.solve(rm[draw].mean(axis=1), a[draw].mean(axis=1)) + mb -= np.eye(2) + return dict(c=c.tolist(), m=np.diag(matrix).tolist(), + matrix=matrix.tolist(), + c_err=cb.std(axis=0).tolist(), + m_err=np.diag(mb.std(axis=0)).tolist(), + response=r.mean(axis=0).tolist(), + seed_groups=[dict(e=s[0].tolist(), R=s[1].tolist(), + A=s[2].tolist(), Rm=s[3].tolist()) + for s in samples]) diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py new file mode 100644 index 000000000..0480c7ef7 --- /dev/null +++ b/tests/module/test_ngmix_defect_fill.py @@ -0,0 +1,532 @@ +"""Defect fill and the epoch cuts (ngmix module). + +A defect is a stamp pixel with a nonzero flag, zero exposure weight or an +invalid background RMS. Before metacal, :func:`prepare_ngmix_weights` gives +every defect weight 0 and fills it with noise at its background RMS, whatever +``BLEND_HANDLING`` is. Under uberseg, pixels on the neighbour side only lose +their weight, and their image values stay raw. The per-epoch cuts in +:func:`prepare_postage_stamps` act on the same defect set: the masked-fraction +cut counts it, and the central-defect veto drops an epoch with a defect near +the stamp centre. +""" + +import re +from types import SimpleNamespace + +import numpy as np +import numpy.testing as npt +from astropy.io import fits +from astropy.wcs import WCS +from hypothesis import given +from hypothesis import strategies as st +from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package import ngmix as ngmix_module + +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + Ngmix, + make_ngmix_observation, + prepare_ngmix_weights, + prepare_postage_stamps, + uberseg_weight, +) + + +# --- prepare_ngmix_weights: the filled set is the defect set --------------- + +@st.composite +def defect_stamps(draw): + """Square stamp with random flagged, zero-weight and bad-RMS pixels.""" + n = draw(st.integers(min_value=5, max_value=21)) + pixels = st.lists( + st.tuples(st.integers(0, n - 1), st.integers(0, n - 1)), max_size=n + ) + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + bkg_rms = np.ones((n, n)) + for i, j in draw(pixels): + weight[i, j] = 0.0 + for i, j in draw(pixels): + flag[i, j] = draw(st.sampled_from([1, 2, 2**10])) + for i, j in draw(pixels): + bkg_rms[i, j] = draw(st.sampled_from([0.0, -1.0, np.nan, np.inf])) + # Raw values far outside the unit-RMS noise, so a filled pixel is + # unambiguous. + gal = 1.0e3 + np.arange(n * n, dtype=float).reshape(n, n) + return gal, weight, flag, bkg_rms + + +def _uberseg_seg(n): + """Central object on the centre pixel, neighbour footprint in a corner.""" + seg = np.zeros((n, n), dtype=np.int32) + seg[n // 2, n // 2] = 1 + seg[:2, :2] = 2 + return seg + + +@given( + stamp=defect_stamps(), + blend_handling=st.sampled_from(["none", "uberseg"]), + seed=st.integers(0, 2**31 - 1), +) +def test_filled_set_is_the_defect_set(stamp, blend_handling, seed): + """Filled pixels are exactly the defects (flag, zero weight, bad RMS). + They carry zero weight and look like noise. Every other pixel keeps its + raw value, under either BLEND_HANDLING. + + Failure modes: + * a defect source (flag, zero weight, bad RMS) is left out of the fill; + * the filled set grows beyond the defects (for example symmetrized); + * the filled set and the zero-weight defect set differ; + * the fill is skipped under uberseg; + * neighbour-side pixels are filled. + """ + gal, weight, flag, bkg_rms = stamp + n = gal.shape[0] + kwargs = ( + dict(seg=_uberseg_seg(n), object_number=1, dilate_neighbour=1) + if blend_handling == "uberseg" + else {} + ) + + gal_out, w_out, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(seed), bkg_rms=bkg_rms, + blend_handling=blend_handling, **kwargs, + ) + + defect = ( + (weight == 0) | (flag != 0) | ~(np.isfinite(bkg_rms) & (bkg_rms > 0)) + ) + if defect.all(): + return # fully masked: no clean pixel left to set the noise level + filled = gal_out != gal + npt.assert_array_equal(filled, defect, "filled set is not the defect set") + assert np.all(np.abs(gal_out[filled]) < 10.0), "fill is not unit noise" + neighbour = ( + uberseg_weight(np.ones((n, n)), kwargs["seg"], 1, dilate_neighbour=1) + == 0.0 + if blend_handling == "uberseg" + else np.zeros((n, n), dtype=bool) + ) + # Zero weight on defects and neighbour side; neighbour pixels that are + # not defects keep their raw values (filled-set equality above). + npt.assert_array_equal(w_out == 0.0, defect | neighbour) + npt.assert_array_equal(w_out[~(defect | neighbour)], 1.0) + + +def test_uberseg_defect_in_neighbour_region_is_filled(): + """A defect pixel that also lies on the neighbour side is filled. + + Failure mode: the fill is restricted to pixels uberseg keeps, so a raw bad + pixel on the neighbour side still reaches metacal. + """ + n = 21 + gal = 1.0e3 + np.arange(n * n, dtype=float).reshape(n, n) + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + flag[1, 1] = 1 # inside the neighbour footprint + gal_out, w_out, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), bkg_rms=np.ones((n, n)), + blend_handling="uberseg", seg=_uberseg_seg(n), object_number=1, + ) + assert w_out[1, 1] == 0.0 and gal_out[1, 1] != gal[1, 1] + # A neighbour-side pixel that is not a defect: zero weight, raw value. + assert w_out[0, 3] == 0.0 and gal_out[0, 3] == gal[0, 3] + + +# --- prepare_postage_stamps: the fraction cut counts the defect set -------- + +N_STAMP = 51 +RA, DEC = 150.0, 2.0 + + +def _fake_inputs(epochs): + """Minimal vignet / tile-catalogue stand-ins for prepare_postage_stamps. + + ``epochs`` maps ``"-"`` to ``(flag, weight)`` or + ``(flag, weight, bkg_rms)``. Returns ``(vignet, tile_cat, psf_obj, + gal_obj)``. + """ + rng = np.random.default_rng(1) + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [RA, DEC] + wcs.wcs.crpix = [N_STAMP / 2, N_STAMP / 2] + wcs.wcs.cdelt = [-0.187 / 3600, 0.187 / 3600] + header = fits.Header({"FSCALE": 1.0}).tostring() + + def per_epoch(make): + return {k: {"VIGNET": make(k)} for k in epochs} + + with_rms = any(len(v) == 3 for v in epochs.values()) + psf_obj = per_epoch(lambda k: np.ones((N_STAMP, N_STAMP))) + gal_obj = per_epoch(lambda k: rng.normal(0.0, 1.0, (N_STAMP, N_STAMP))) + vignet = SimpleNamespace( + gal_vign_cat={"1": gal_obj}, + bkg_vign_cat=None, + bkg_rms_vign_cat=( + {"1": per_epoch( + lambda k: epochs[k][2] if len(epochs[k]) == 3 + else np.ones((N_STAMP, N_STAMP)) + )} + if with_rms + else None + ), + flag_vign_cat={"1": per_epoch(lambda k: epochs[k][0])}, + weight_vign_cat={"1": per_epoch(lambda k: epochs[k][1])}, + f_wcs_file={ + k.split("-")[0]: { + int(k.split("-")[1]): {"WCS": wcs, "header": header} + } + for k in epochs + }, + ) + tile_cat = SimpleNamespace( + vign=None, seg=None, ra=np.array([RA]), dec=np.array([DEC]) + ) + return vignet, tile_cat, psf_obj, gal_obj + + +def _two_sided_band(width): + """Mask of ``width`` columns on each side of the stamp, far from the + centre (at least 17 px for width 9).""" + band = np.zeros((N_STAMP, N_STAMP), dtype=bool) + band[:, :width] = True + band[:, -width:] = True + return band + + +def _epochs(): + """A clean epoch and four masked ones, all defects far from the centre. + + * band10: 10 flagged columns on one side, 19.6% (39% if symmetrized); + * flag18: 2 x 9 flagged columns, 35.3%; + * dead18: the same columns at zero weight, no flags; + * rms18: the same columns with an invalid background RMS. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + band10 = clean.copy() + band10[:, -10:] = 1 + two = _two_sided_band(9) + dead = ones.copy() + dead[two] = 0.0 + rms = ones.copy() + rms[two] = np.nan + return { + "2100001-10": (clean, ones), + "2100002-11": (band10, ones), + "2100003-12": (two.astype(np.int32), ones), + "2100004-13": (clean.copy(), dead), + "2100005-14": (clean.copy(), ones, rms), + } + + +def _surviving(epochs, **kwargs): + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, **kwargs, + ) + names = {id(v[0]): name for name, v in epochs.items()} + return sorted(names[id(flag)] for flag in stamp.flags) + + +def test_epoch_cut_counts_the_defect_set(): + """At the default 1/3 cut, the three 35% epochs are dropped, whichever + defect source masks them, and the 20% band survives. + + Failure modes: the cut counts flags only (keeps dead18 and rms18), omits + one defect source, or counts a symmetrized set (drops band10); in each + case the cut disagrees with the set that is zero-weighted and filled + (epoch-cut-on-defect-mask). + """ + assert _surviving(_epochs()) == ["2100001-10", "2100002-11"] + + +def test_epoch_cut_threshold_is_the_configured_fraction(): + """At a 10% cut (the DES Y3 and Y6 value), the 19.6% band is dropped as + well, and only the clean epoch survives. + + Failure mode: the configured threshold is ignored. + """ + assert _surviving(_epochs(), epoch_masked_fraction_cut=0.1) == [ + "2100001-10" + ] + + +# --- prepare_postage_stamps: the central-defect veto ----------------------- + +def _veto_epochs(radius): + """A clean epoch, one with a single flagged pixel just inside + ``radius`` of the stamp centre, and one with a flagged column exactly + ``radius`` away. Both masked epochs are far below the fraction cut. + """ + centre = N_STAMP // 2 + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + near, far = clean.copy(), clean.copy() + near[centre, centre + radius - 1] = 1 + far[:, centre + radius] = 1 + return { + "2100001-10": (clean, ones), + "2100002-11": (near, ones), + "2100003-12": (far, ones), + } + + +def test_central_defect_vetoes_the_epoch(): + """At the default radius, a single defect pixel inside it drops the + epoch, and a column at the radius does not. + + Failure mode: the veto is skipped, so an epoch whose filled hole overlaps + the object's light enters the fit; or the boundary is inclusive, dropping + the column at the radius (epoch-central-defect-veto). + """ + epochs = _veto_epochs(EPOCH_CENTRAL_DEFECT_RADIUS) + assert _surviving(epochs) == ["2100001-10", "2100003-12"] + + +def test_central_defect_radius_is_the_configured_value(): + """Radius 0 disables the veto; a radius one pixel larger than the far + column's distance drops that epoch too. + + Failure mode: the configured radius is ignored. + """ + epochs = _veto_epochs(EPOCH_CENTRAL_DEFECT_RADIUS) + assert _surviving(epochs, epoch_central_defect_radius=0) == [ + "2100001-10", "2100002-11", "2100003-12" + ] + assert _surviving( + epochs, epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS + 1 + ) == ["2100001-10"] + + +# --- prepare_postage_stamps: per-epoch OFFSET ------------------------------ + +def test_each_surviving_epoch_carries_its_own_offset(): + """``stamp.offsets`` holds each surviving epoch's vignette OFFSET, in the + order of ``stamp.flags``. + + Failure mode: the offset is dropped or read from another epoch, so the + default "wcs" centroid raises or puts the Jacobian origin off the object. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + epochs = {f"210000{i}-1{i}": (clean.copy(), ones) for i in range(3)} + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + del vignet.gal_vign_cat # OFFSET must come from the supplied gal_obj. + for i, name in enumerate(epochs): + gal_obj[name]["OFFSET"] = np.array([0.1 * i, -0.1 * i]) + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, + ) + names = {id(flag): name for name, (flag, _) in epochs.items()} + assert len(stamp.offsets) == len(stamp.flags) == 3 + for flag, offset in zip(stamp.flags, stamp.offsets): + npt.assert_array_equal(offset, gal_obj[names[id(flag)]]["OFFSET"]) + + +# --- Ngmix.process: the per-tile epoch-cut tally --------------------------- + +class _RecordingLogger: + def __init__(self): + self.messages = [] + + def info(self, msg, *_args, **_kwargs): + self.messages.append(msg) + + +def test_process_logs_the_epoch_cut_tally(tmp_path, monkeypatch): + """One tile, four objects; the end-of-tile line counts each cut's drops. + + * object 1: clean, 18-column edge band, defect inside the radius -> one + epoch each for considered, masked_fraction, central_veto; survives. + * object 2: edge band and central defect -> both epochs dropped; emptied. + * object 3: one all-zero stamp, skipped before the cuts -> not considered, + and not emptied by the cuts. + * object 4: no PSF ('empty') -> never reaches the cuts. + + Failure modes: a cut's drops are not counted or land in the wrong + counter; epochs skipped before the cuts are counted as considered; an + object with no epoch at all is reported as emptied by the cuts; counts + from one object overwrite another's. + """ + radius = EPOCH_CENTRAL_DEFECT_RADIUS + centre = N_STAMP // 2 + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + wide18, near = clean.copy(), clean.copy() + wide18[:, -18:] = 1 + near[centre, centre + radius - 1] = 1 + objects = { + 1: {"2100001-10": (clean, ones), "2100002-11": (wide18, ones), + "2100003-12": (near, ones)}, + 2: {"2100004-13": (wide18.copy(), ones), + "2100005-14": (near.copy(), ones)}, + 3: {"2100006-15": (clean.copy(), ones)}, + } + stores = {} + for obj_id, epochs in objects.items(): + vignet, _, psf_obj, gal_obj = _fake_inputs(epochs) + if obj_id == 3: + gal_obj["2100006-15"]["VIGNET"] = np.zeros((N_STAMP, N_STAMP)) + stores[obj_id] = (vignet, psf_obj, gal_obj) + vignet = SimpleNamespace( + bkg_vign_cat=None, + bkg_rms_vign_cat=None, + psf_vign_cat={ + "4": "empty", **{str(i): s[1] for i, s in stores.items()} + }, + gal_vign_cat={ + "4": "empty", **{str(i): s[2] for i, s in stores.items()} + }, + flag_vign_cat={ + str(i): s[0].flag_vign_cat["1"] for i, s in stores.items() + }, + weight_vign_cat={ + str(i): s[0].weight_vign_cat["1"] for i, s in stores.items() + }, + f_wcs_file={ + k: v for s in stores.values() for k, v in s[0].f_wcs_file.items() + }, + close=lambda: None, + ) + tile_cat = SimpleNamespace( + obj_id=np.array([1, 2, 3, 4]), ra=np.full(4, RA), dec=np.full(4, DEC), + vign=None, seg=None, flux=None, + ) + + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + log = _RecordingLogger() + ngmix = Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, 0.186, str(paths[4]), log, + bkg_sub=False, + ) + ngmix._vignet_cat.close() + ngmix._vignet_cat = vignet + monkeypatch.setattr(ngmix_module, "Tile_cat", lambda *a, **k: tile_cat) + + def no_fit(*_args, **_kwargs): + raise RuntimeError("metacal is not under test") + + monkeypatch.setattr(ngmix_module, "do_ngmix_metacal", no_fit) + for method in ("compile_results", "save_results", "log_mean_ellipticity"): + monkeypatch.setattr(Ngmix, method, lambda *_a, **_k: None) + + ngmix.process() + + lines = [m for m in log.messages if m.startswith("epoch cuts:")] + assert len(lines) == 1, log.messages + tally = dict( + (k, int(v)) for k, v in re.findall(r"(\w+)=(\d+)", lines[0]) + ) + assert tally == dict( + considered=5, masked_fraction=2, central_veto=2, objects_emptied=1 + ), lines[0] + + +# --- ngmix_runner: the epoch-cut options reach Ngmix ------------------------ + +class _OptionConfig: + """Config stub answering from a dict; absent options take the caller's + fallback.""" + + def __init__(self, options): + self._options = {"MAG_ZP": "30.0", "ID_OBJ_MIN": "-1", + "ID_OBJ_MAX": "-1", **options} + + def has_option(self, _sec, key): + return key in self._options + + def get(self, _sec, key): + return self._options[key] + + def getexpanded(self, _sec, key): + return self._options[key] + + def getfloat(self, _sec, key): + return float(self._options[key]) + + def getint(self, _sec, key): + return int(self._options[key]) + + def getboolean(self, _sec, key, fallback=False): + return fallback + + +def test_runner_threads_the_epoch_cut_options(tmp_path, monkeypatch): + """EPOCH_CENTRAL_DEFECT_RADIUS and EPOCH_MASKED_FRACTION_CUT reach Ngmix + as configured, and default to the module constants when absent. + + Failure mode: the runner drops or ignores an option, so a configured A/B + arm silently runs the default cuts. + """ + from shapepipe.modules import ngmix_runner as runner_module + + captured = [] + + class _Capture: + def __init__(self, *_args, **kwargs): + captured.append(kwargs) + + def process(self): + pass + + monkeypatch.setattr(runner_module, "Ngmix", _Capture) + inputs = [str(tmp_path / f"in{i}.sqlite") for i in range(7)] + for path in inputs: + SqliteDict(path).close() + + for options, radius, fraction in ( + ({"EPOCH_CENTRAL_DEFECT_RADIUS": "7.5", + "EPOCH_MASKED_FRACTION_CUT": "0.1"}, 7.5, 0.1), + ({}, EPOCH_CENTRAL_DEFECT_RADIUS, + ngmix_module.EPOCH_MASKED_FRACTION_CUT), + ): + runner_module.ngmix_runner( + inputs, {"output": str(tmp_path)}, "-001-001", + _OptionConfig(options), "NGMIX_RUNNER", _RecordingLogger(), + ) + assert captured[-1]["epoch_central_defect_radius"] == radius + assert captured[-1]["epoch_masked_fraction_cut"] == fraction + + +# --- make_ngmix_observation: the HSM centroid reads the filled image ------- + +def test_hsm_centroid_ignores_raw_defect_values(): + """With ``centroid_source="hsm"``, the Jacobian origin lands on the + object even when a flagged column a few pixels away holds a raw value + ten times the object's peak. + + Failure mode: HSM measures the raw stamp, so the defect drags the centroid + off the object (or HSM fails and falls back to the stamp centre). + """ + import galsim + + n = 51 + centre = (n - 1) / 2 + d_row, d_col = 1.3, -0.7 + rows, cols = np.mgrid[:n, :n] + gal = np.exp( + -((rows - centre - d_row) ** 2 + (cols - centre - d_col) ** 2) + / (2 * 2.0 ** 2) + ) + flag = np.zeros((n, n), dtype=np.int32) + flag[:, n // 2 + 6] = 1 + gal[flag != 0] = 10.0 + psf = np.exp(-((rows - centre) ** 2 + (cols - centre) ** 2) / 2.0) + obs = make_ngmix_observation( + gal, np.ones((n, n)), flag, psf / psf.sum(), + galsim.PixelScale(0.1857).jacobian(), np.random.RandomState(0), + bkg_rms=np.full((n, n), 1e-3), centroid_source="hsm", + ) + row, col = obs.jacobian.get_cen() + npt.assert_allclose( + [row - centre, col - centre], [d_row, d_col], atol=0.05 + ) diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 0c80486b2..86dad33bc 100644 --- a/tests/module/test_ngmix_uberseg.py +++ b/tests/module/test_ngmix_uberseg.py @@ -6,10 +6,10 @@ Geometry assertions on a synthetic two-object stamp: neighbour-side pixels are zeroed, the surviving central core is a *single connected* region (the emergent "circularisation"), and the neighbour footprint is fully removed. -* :func:`prepare_ngmix_weights` — the ``noisefill`` default is byte-for-byte - unchanged (asserted against an independent recomputation of the legacy - three-line noise-fill on a shared RNG), while ``uberseg`` hard-masks the - weight (weight -> 0) and leaves the image untouched. +* :func:`prepare_ngmix_weights` under ``uberseg`` — neighbour-side pixels + lose their weight and keep their raw image values, while defect pixels are + noise-filled as under any blend handling (the defect fill itself is covered + in ``test_ngmix_defect_fill.py``). * The error contract when ``uberseg`` is selected without a segmentation map (the seg-map source is plumbing-gated upstream). """ @@ -226,7 +226,7 @@ def test_uberseg_matches_bruteforce_nearest_segment(): npt.assert_array_equal(out, brute) -# --- prepare_ngmix_weights: default unchanged, uberseg hard-masks ---------- +# --- prepare_ngmix_weights: uberseg zeroes neighbour weights only --------- def _gal_flag_weight(npix=41, seed=7): rng = np.random.default_rng(seed) @@ -239,32 +239,6 @@ def _gal_flag_weight(npix=41, seed=7): return gal, flag, weight -def test_noisefill_default_is_byte_identical_to_legacy(): - """The default path reproduces the legacy three-line noise-fill exactly - (same RNG stream): masked pixels replaced by noise, weight 1/sigma^2.""" - gal, flag, weight = _gal_flag_weight() - - gal_out, w_out, noise_out = prepare_ngmix_weights( - gal, weight, flag, np.random.RandomState(123), - ) - - # Independent recomputation of the legacy algorithm on the same stream. - from modopt.math.stats import sigma_mad - rng = np.random.RandomState(123) - mask = np.copy(weight) != 0 - mask[flag != 0] = False - sig = sigma_mad(gal) - w_exp = mask.astype(float) / sig ** 2 - noise_exp = rng.standard_normal(gal.shape) * sig - noise_gal = rng.standard_normal(gal.shape) * sig - gal_exp = np.copy(gal) - gal_exp[~mask] = noise_gal[~mask] - - npt.assert_array_equal(gal_out, gal_exp) - npt.assert_array_equal(w_out, w_exp) - npt.assert_array_equal(noise_out, noise_exp) - - def test_noisefill_ignores_seg_and_dilate_kwargs(): """Under noisefill, passing seg / dilate_neighbour changes nothing: the result matches the plain default call on the same RNG stream.""" @@ -282,9 +256,15 @@ def test_noisefill_ignores_seg_and_dilate_kwargs(): npt.assert_array_equal(a, b) -def test_uberseg_hard_masks_weight_and_leaves_image_untouched(): - """uberseg: image returned untouched, weight zeroed on neighbour-side and - flagged pixels, positive on the central core.""" +def test_uberseg_fills_defects_and_leaves_neighbour_pixels_raw(): + """uberseg: flagged pixels are noise-filled at weight 0; neighbour-side + pixels get weight 0 and keep their raw image values; the central core + keeps weight and image. + + Failure modes: the defect fill is skipped under uberseg, so raw bad + pixels reach metacal; or the neighbour side is noise-filled + (defects-filled-neighbours-raw). + """ npix = 41 gal, flag, weight = _gal_flag_weight(npix=npix) seg, centre, neigh = two_object_seg(npix=npix, sep=12) @@ -294,13 +274,16 @@ def test_uberseg_hard_masks_weight_and_leaves_image_untouched(): blend_handling="uberseg", seg=seg, object_number=1, ) - # Image untouched under uberseg (no noise fill). - npt.assert_array_equal(gal_out, gal) - # Neighbour footprint hard-masked; central centre kept. + # Flagged pixels: zero weight, image replaced by noise. + for pix in [(5, 5), (30, 12)]: + assert w_out[pix] == 0.0 + assert gal_out[pix] != gal[pix] + # Neighbour footprint: zero weight, raw image (never noise-filled). assert np.all(w_out[seg == 2] == 0.0) + npt.assert_array_equal(gal_out[seg == 2], gal[seg == 2]) + # Central core: weight and image untouched. assert w_out[centre] > 0.0 - # Flagged bad pixels remain at weight 0 (folded into the base mask). - assert w_out[5, 5] == 0.0 + assert gal_out[centre] == gal[centre] def test_uberseg_requires_seg_and_object_number(): @@ -422,3 +405,11 @@ def test_runner_missing_seg_vignet_file_raises(tmp_path): "NGMIX_RUNNER", _RecordingLogger(), ) + + +def test_retired_noisefill_name_is_rejected(): + """Do not silently interpret the retired neighbour option.""" + with pytest.raises(ValueError, match="noisefill"): + prepare_ngmix_weights(np.ones((5, 5)), np.ones((5, 5)), + np.zeros((5, 5)), np.random.RandomState(3), + blend_handling="noisefill") diff --git a/tests/science/test_defect_veto.py b/tests/science/test_defect_veto.py new file mode 100644 index 000000000..3384263e4 --- /dev/null +++ b/tests/science/test_defect_veto.py @@ -0,0 +1,104 @@ +"""Shear recovery for the defects the central-defect veto keeps (noise fill). + +Physics invariant: a defect the veto keeps -- a column, a 3-px bleed or a +single pixel at the veto radius or beyond, or an edge band -- leaves both +additive terms |c1|, |c2| < 5e-4 and both diagonal multiplicative terms +|m11|, |m22| < 1%, from the full 2x2 response matrix. The grid spans galaxies +with half-light radius 0.3" and 0.5" through a 0.7" PSF, round and with +ellipticity (0.05, 0.02). A noise-filled defect biases m anisotropically +(m11 and m22 differ by up to a factor of ten), so a scalar m would hide it. +The cases sit at ``EPOCH_CENTRAL_DEFECT_RADIUS`` itself, so lowering the +radius below the calibrated value turns this red. + +Positive control: the same column two pixels inside the radius breaks the +bound, so the recovery check can see the bias the veto removes. +""" + +import json + +import numpy as np +import pytest + +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + defect_mask, + has_central_defect, +) +from tests.helpers.defect_response import defect_response + +N = 51 +CENTRE = N // 2 +R = int(np.ceil(EPOCH_CENTRAL_DEFECT_RADIUS)) +ELLIPTICAL = (0.05, 0.02) +SEEDS = range(6) + + +def geometry(kind, distance): + """A detector defect whose nearest pixel is ``distance`` px from the + stamp centre; an edge band reaches in to that distance.""" + bad = np.zeros((N, N), dtype=bool) + if kind == "pixel": + bad[CENTRE, CENTRE + distance] = True + elif kind == "column": + bad[:, CENTRE + distance] = True + elif kind == "bleed": + bad[:, CENTRE + distance:CENTRE + distance + 3] = True + elif kind == "edge": + bad[:, CENTRE + distance:] = True + return bad + + +# Wide defects (a 3-px bleed, the widest edge band the veto keeps) at the +# radius itself on the 0.5" galaxy through the elliptical PSF sit at the +# bound (m11 = -0.98% +/- 0.04% and -0.94% +/- 0.06%; see +# :func:`has_central_defect`), so those two are checked one pixel out. +CASES = ( + [(k, d, h, (0.0, 0.0)) for k in ("column", "bleed", "pixel") + for d in (R, R + 2) for h in (0.3, 0.5)] + + [(k, R, 0.5, ELLIPTICAL) for k in ("column", "pixel")] + + [(k, R + 1, 0.5, ELLIPTICAL) for k in ("bleed", "edge")] + + [("edge", R, 0.5, (0.0, 0.0))] + + [("edge", N - CENTRE - 5, 0.5, psf) for psf in ((0.0, 0.0), ELLIPTICAL)] +) + + +def recover(bad, hlr, psf_shear, tmp_path): + result = defect_response(bad, hlr=hlr, psf=0.7, seeds=SEEDS, + psf_shear=psf_shear) + (tmp_path / "recovery.json").write_text(json.dumps(result, indent=2)) + return np.abs(result["m"]).max(), np.abs(result["c"]).max(), result + + +def case_id(case): + kind, distance, hlr, psf_shear = case + psf = "elliptical" if any(psf_shear) else "round" + return f"{kind}-{distance}px-hlr{hlr}-{psf}" + + +@pytest.mark.parametrize("kind,distance,hlr,psf_shear", CASES, + ids=[case_id(c) for c in CASES]) +def test_kept_defects_recover_shear_on_both_axes(kind, distance, hlr, + psf_shear, tmp_path): + """Failure modes: the veto radius is below the calibrated value; the + fill leaks raw defect values or leaves more of the object's light + missing; the check reads m11 alone and misses m22, or c1 alone. + """ + bad = geometry(kind, distance) + masked = defect_mask(np.ones((N, N)), bad.astype(np.int32)) + np.testing.assert_array_equal(masked, bad) + assert masked.mean() <= EPOCH_MASKED_FRACTION_CUT + assert not has_central_defect(masked, EPOCH_CENTRAL_DEFECT_RADIUS) + m, c, result = recover(bad, hlr, psf_shear, tmp_path) + assert m < 0.01, result + assert c < 5e-4, result + + +def test_vetoed_column_breaks_the_bound(tmp_path): + """Positive control: a column two pixels inside the radius, on the + 0.5" galaxy, gives |m| > 2% and |c| > 1e-3. The veto drops it.""" + bad = geometry("column", R - 2) + assert has_central_defect(bad, EPOCH_CENTRAL_DEFECT_RADIUS) + m, c, result = recover(bad, 0.5, (0.0, 0.0), tmp_path) + assert m > 0.02, result + assert c > 1e-3, result From cd8b321b2bc81e5920e9a2a9b7d501c846fff9a0 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:07:40 +0200 Subject: [PATCH 02/14] feat(ngmix): noise-fill defects under every BLEND_HANDLING, veto central ones Metacal shears the whole stamp without reading the weights, so a defect left raw under uberseg reached the fit. Every defect (nonzero flag, zero weight, invalid RMS) is now zero-weighted and noise-filled whatever BLEND_HANDLING is; uberseg only zeroes the weight of neighbour-side pixels and keeps their light. The option that meant "no neighbour treatment" is renamed from noisefill to none, and the retired name is rejected. The defect set is not symmetrized: at the veto radius the one-sided fill leaves |c| <= 3e-4, while a four-fold mask would quadruple m (-2.7% against -0.64% for a 3-px bleed at 10 px on a 0.5" galaxy) and, through an elliptical PSF, still leave c1 = 7.3e-4. An epoch is dropped when a defect lies closer than EPOCH_CENTRAL_DEFECT_RADIUS (default 10 px) to the stamp centre, the smallest radius at which columns, 3-px bleeds, single pixels and edge bands give |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 on 0.3" and 0.5" galaxies through a 0.7" PSF. The masked-fraction cut counts the same raw defect set; EPOCH_MASKED_FRACTION_CUT (default 1/3) makes it configurable. Each tile logs "epoch cuts: considered=N masked_fraction=A central_veto=B objects_emptied=C". The HSM centroid is measured on the filled image. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01QAc1ywyd9Sbbra6po2iuaA --- src/shapepipe/modules/ngmix_package/ngmix.py | 353 +++++++++++++++---- src/shapepipe/modules/ngmix_runner.py | 38 +- 2 files changed, 317 insertions(+), 74 deletions(-) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index b2e092240..807b99ff0 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -8,6 +8,7 @@ import os import re +from collections import Counter from typing import NamedTuple import ngmix @@ -25,7 +26,17 @@ from shapepipe.pipeline import file_io # Neighbour treatments selectable with the BLEND_HANDLING option. -BLEND_HANDLINGS = ("noisefill", "uberseg") +BLEND_HANDLINGS = ("none", "uberseg") + +# Default of the EPOCH_MASKED_FRACTION_CUT option: an epoch is dropped when +# more than this fraction of its stamp is in :func:`defect_mask` (see +# :func:`prepare_postage_stamps`). +EPOCH_MASKED_FRACTION_CUT = 1 / 3 + +# Default of the EPOCH_CENTRAL_DEFECT_RADIUS option (pixels): an epoch is +# dropped when a pixel of :func:`defect_mask` lies closer than this to the +# stamp centre (see :func:`has_central_defect`). +EPOCH_CENTRAL_DEFECT_RADIUS = 10 METACAL_TYPES = ('noshear', '1p', '1m', '2p', '2m') @@ -191,7 +202,7 @@ class Tile_cat(): vignetmaker output cut from the tile ``SEGMENTATION`` check image), row-aligned to ``cat_path``. When given, ``self.seg`` holds one integer seg stamp per object for the ``"uberseg"`` blend handling; ``None`` - leaves ``self.seg`` unset (the noise-fill path is unaffected). + leaves ``self.seg`` unset, which only ``"uberseg"`` needs. """ def __init__( @@ -226,7 +237,7 @@ def get_data(self, cat_path): # Coadd-frame SExtractor segmentation stamp (integer labels), one per # object and row-aligned to the tile catalogue, overlaid unchanged on # every epoch for uberseg neighbour masking (shapepipe#776). None -> - # uberseg unavailable; the noise-fill path is unaffected. + # uberseg unavailable; ``"none"`` does not read it. self.seg = None if self.seg_cat_path: seg_cat = file_io.FITSCatalogue( @@ -281,7 +292,7 @@ def __init__( self.flags = [] self.bkg_rms = [] # Segmentation stamps, one per epoch, used only by the "uberseg" blend - # handling; empty for the default noise-fill path. All epochs carry the + # handling; empty under the default "none". All epochs carry the # SAME coadd-frame seg stamp (shapepipe#776: one coadd seg per object, # no per-epoch reprojection), each MegaCam-flipped to match its galaxy # stamp so the overlay stays registered. @@ -298,6 +309,9 @@ def __init__( # CCD number of the first epoch, used only to build the per-object # position seed (see :func:`position_seed`). self.ccd = None + self.epoch_cuts = Counter( + considered=0, masked_fraction=0, central_veto=0 + ) self.bkg_sub = bkg_sub self.megacam_flip = megacam_flip @@ -426,16 +440,26 @@ class Ngmix(object): propagated on the vignette. ``"hsm"`` re-centers on the adaptive-moment centroid measured from the stamp pixels. See :func:`make_ngmix_observation`. - blend_handling : {"noisefill", "uberseg"}, optional - Neighbour treatment; ``"noisefill"`` (default) is the historical - noise-fill, ``"uberseg"`` hard-masks neighbour-side pixels from the - coadd segmentation map and requires ``seg_cat_path``. + blend_handling : {"none", "uberseg"}, optional + Neighbour treatment. ``"none"`` (default) leaves neighbour pixels + untouched; ``"uberseg"`` zeroes the weight of neighbour-side pixels + from the coadd segmentation map and requires ``seg_cat_path``. Defect + pixels are noise-filled under both (see + :func:`prepare_ngmix_weights`). seg_cat_path : str, optional Path to the coadd-frame segmentation VIGNET catalogue (see :class:`Tile_cat`). Required when ``blend_handling="uberseg"``. dilate_neighbour : int, optional Neighbour-mask dilation iterations for ``"uberseg"`` (see :func:`uberseg_weight`); the default is ``1``. + epoch_central_defect_radius : float, optional + Drop an epoch when a masked pixel lies closer than this many pixels + to the stamp centre (see :func:`has_central_defect`); the default is + ``EPOCH_CENTRAL_DEFECT_RADIUS``. + epoch_masked_fraction_cut : float, optional + Drop an epoch when more than this fraction of its stamp is masked + (see :func:`prepare_postage_stamps`); the default is + ``EPOCH_MASKED_FRACTION_CUT``. Notes ----- @@ -466,10 +490,12 @@ def __init__( id_obj_max=-1, bkg_sub=True, centroid_source="wcs", - blend_handling="noisefill", + blend_handling="none", seg_cat_path=None, dilate_neighbour=1, metacal_psf="fitgauss", + epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, + epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, ): # Base count = catalogue + vignets, excluding the f_wcs headers (passed @@ -539,6 +565,8 @@ def __init__( self._seg_cat_path = seg_cat_path self._dilate_neighbour = dilate_neighbour self._metacal_psf = metacal_psf + self._epoch_central_defect_radius = epoch_central_defect_radius + self._epoch_masked_fraction_cut = epoch_masked_fraction_cut self._w_log = w_log @@ -985,6 +1013,8 @@ def process(self): id_first = -1 id_last = -1 count_batch = 0 + epoch_cuts = Counter(considered=0, masked_fraction=0, central_veto=0) + n_emptied = 0 saved_batch_cumul = 0 for i_tile, obj_id in enumerate(tile_cat.obj_id): @@ -1014,10 +1044,14 @@ def process(self): self._bkg_sub, psf_obj, gal_obj, + epoch_central_defect_radius=self._epoch_central_defect_radius, + epoch_masked_fraction_cut=self._epoch_masked_fraction_cut, ) + epoch_cuts.update(stamp.epoch_cuts) if len(stamp.gals) == 0: n_no_epoch += 1 + n_emptied += stamp.epoch_cuts["considered"] > 0 continue # Per-object RNG, seeded from (ra, dec, ccd) — see @@ -1131,6 +1165,13 @@ def process(self): + f" {n_fitted} fitted" ) + self._w_log.info( + "epoch cuts:" + + f" considered={epoch_cuts['considered']}" + + f" masked_fraction={epoch_cuts['masked_fraction']}" + + f" central_veto={epoch_cuts['central_veto']}" + + f" objects_emptied={n_emptied}" + ) vignet_cat.close() # Put all results together @@ -1150,7 +1191,46 @@ def prepare_postage_stamps( bkg_sub=True, psf_obj=None, gal_obj=None, + epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, + epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, ): + """Gather one object's epoch stamps, dropping epochs its defects spoil. + + @sc [decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.defect_fill] epoch-cut-on-defect-mask + An epoch is dropped when more than ``epoch_masked_fraction_cut`` of its + stamp lies in :func:`defect_mask`: flagged, zero-weight and invalid-RMS + pixels, the set that :func:`prepare_ngmix_weights` zero-weights and + noise-fills. Counting flags alone would keep epochs whose filled area + exceeds the cut. The default is 1/3; 10%, the DES Y3 and Y6 value, is + the alternative to test. + + Parameters + ---------- + vignet : Vignet + Per-object vignet stores. + obj_id : int + Object ID (SExtractor ``NUMBER``). + i_tile : int + Row of the object in ``tile_cat``. + tile_cat : Tile_cat + Tile catalogue. + bkg_sub : bool, optional + Subtract the background vignet; the default is ``True``. + psf_obj, gal_obj : dict, optional + The object's PSF and galaxy vignet dicts, if already read. + epoch_central_defect_radius : float, optional + Drop an epoch with a defect pixel closer than this many pixels to the + stamp centre (:func:`has_central_defect`); 0 disables the veto. The + default is ``EPOCH_CENTRAL_DEFECT_RADIUS``. + epoch_masked_fraction_cut : float, optional + Drop an epoch with more than this fraction of its stamp in + :func:`defect_mask`. The default is ``EPOCH_MASKED_FRACTION_CUT``. + + Returns + ------- + Postage_stamp + The surviving epochs' stamps. + """ # define per-object lists of individual exposures to go into ngmix stamp = Postage_stamp(bkg_sub=bkg_sub) # Read each store's per-object dict ONCE: every sqlitedict access @@ -1220,17 +1300,23 @@ def prepare_postage_stamps( flag_vign = flag_obj[expccd_name]['VIGNET'] if tile_vign is not None: flag_vign[np.where(tile_vign == -1e30)] = 2**10 - v_flag_tmp = flag_vign.ravel() - # remove objects that are more than 1/3 masked - if len(np.where(v_flag_tmp != 0)[0]) / v_flag_tmp.size > 1 / 3.0: - continue - weight_vign = weight_obj[expccd_name]['VIGNET'] bkg_rms_vign = ( bkg_rms_obj[expccd_name]['VIGNET'] if bkg_rms_obj is not None else None ) + # Drop the epoch when too much of it would be zero-weighted and + # noise-filled (epoch-cut-on-defect-mask), or when a filled pixel + # would sit near the object (epoch-central-defect-veto). + masked = defect_mask(weight_vign, flag_vign, bkg_rms_vign) + stamp.epoch_cuts["considered"] += 1 + if masked.mean() > epoch_masked_fraction_cut: + stamp.epoch_cuts["masked_fraction"] += 1 + continue + if has_central_defect(masked, epoch_central_defect_radius): + stamp.epoch_cuts["central_veto"] += 1 + continue # One unpickle per exposure (all CCDs), reused across this object's # epochs; the cache is per-call, so bounded by the object's exposures @@ -1503,9 +1589,9 @@ def uberseg_weight(weight, seg, object_number, dilate_neighbour=0): Because the partition is by distance to the nearest footprint, the pixels surviving around a compact central object form a single connected, roughly circular core; the "circularisation" is emergent geometry, not a - separate aperture. Unlike the noise-fill treatment, masked pixels are - handed to ngmix as a hard mask (weight = 0), never replaced by a noise - realisation. + separate aperture. Only the weight changes: neighbour-side pixels are + handed to ngmix as a hard mask (weight = 0) and keep their image values, + so metacal shears the neighbour's light along with the target's. Parameters ---------- @@ -1566,33 +1652,163 @@ def uberseg_weight(weight, seg, object_number, dilate_neighbour=0): return weight + +def defect_mask(weight, flag, bkg_rms=None): + """Defect pixels of one epoch stamp. + + @sc [decision:shape_measurement.defect_fill,label:physics] defect-set-unsymmetrized + A defect is a pixel with zero exposure weight, a nonzero flag, or (when a + background RMS map is given) a non-finite or non-positive RMS. This one + set is zero-weighted and noise-filled by :func:`prepare_ngmix_weights` + and counted by the epoch cuts in :func:`prepare_postage_stamps`. It is + not ORed with its rotations. For defects the central-defect veto keeps + (:func:`has_central_defect`), the unsymmetrized fill leaves + |c| <= 3e-4 per affected epoch on galaxies with half-light radius 0.3" + and 0.5" through a 0.7" PSF, round or elliptical. Symmetrizing would + quadruple the filled area near the object, and with it m: a 3-px bleed at + 10 px on the 0.5" galaxy gives m11 = -2.7% four-fold against -0.64% + unsymmetrized. Through a PSF with ellipticity (0.05, 0.02) it would not + cancel c either: c1 = 7.3e-4 four-fold against 0.7e-4. It would also + turn a 5-px edge band (9.8% of the stamp), which biases nothing, into a + 35.4% frame that fails the masked-fraction cut. + + Parameters + ---------- + weight : numpy.ndarray + Exposure weight stamp. + flag : numpy.ndarray + Flag stamp. + bkg_rms : numpy.ndarray, optional + Background RMS stamp. + + Returns + ------- + numpy.ndarray of bool + ``True`` on defect pixels. + """ + defect = (weight == 0) | (flag != 0) + if bkg_rms is not None: + defect |= ~(np.isfinite(bkg_rms) & (bkg_rms > 0)) + return defect + + +def has_central_defect(defect, radius): + """Whether a defect pixel lies closer than ``radius`` to the stamp centre. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.defect_fill] epoch-central-defect-veto + An epoch is dropped when any pixel of :func:`defect_mask` lies closer + than ``radius`` pixels to the stamp centre. A noise-filled hole in the + object's light is sheared by metacal but not by the sky, so the metacal + response is wrong for that epoch, and a one-sided hole adds an additive + term. The default radius (``EPOCH_CENTRAL_DEFECT_RADIUS``, 10 px or + 1.9") is the smallest at which columns, 3-px bleeds, single pixels and + edge bands give |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 (full response + matrix) on galaxies with half-light radius 0.3" and 0.5" through a 0.7" + PSF. The bias is anisotropic: a column at 10 px on the 0.5" galaxy gives + m11 = -0.60%, m22 = +0.01% and c1 = -1.4e-4; at 9 px, m11 = -1.75% and + c1 = -1.1e-3. Through a PSF with ellipticity (0.05, 0.02), wide defects + at 10 px on the 0.5" galaxy sit at the bound: a 3-px bleed gives + m11 = -0.98% +/- 0.04% and the widest edge band the veto keeps (16 px) + -0.94% +/- 0.06%; at 11 px the bleed gives -0.24%. Larger galaxies need + more: at hlr 0.7" through a 0.9" PSF a 3-px bleed at 10 px gives + m11 = -6.8% and c1 = -4.9e-3, and passes from 14 px. The distance is + measured from the stamp centre, where the extractor places the object to + within half a pixel. The veto reads only the defect mask, never the + object's pixels, so it selects on nothing that responds to shear. + Radius 0 disables it. Guarded by ``tests/science/test_defect_veto.py``. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp (:func:`defect_mask`). + radius : float + Veto radius in pixels. + + Returns + ------- + bool + ``True`` if the epoch should be dropped. + """ + rows, cols = np.nonzero(defect) + centre_row = (defect.shape[0] - 1) / 2 + centre_col = (defect.shape[1] - 1) / 2 + distance = np.hypot(rows - centre_row, cols - centre_col) + return bool(np.any(distance < radius)) + + +def fill_defects(image, defect, noise): + """Replace an image's defect pixels by a noise realisation. + + @sc [decision:shape_measurement.defect_fill,label:physics] defect-noise-fill + Defect pixels are replaced by ``noise``, an independent realisation at + the per-pixel background RMS. Metacal deconvolves, shears and reconvolves + the whole image without reading the weights, so a raw defect value (bad + column, bleed, cosmic ray) would leak into the fit. The noise fill leaves + a hole in the object's light that metacal shears and the sky does not, + which biases m for an epoch with a defect near the object; the + central-defect veto drops those epochs (epoch-central-defect-veto). + + Parameters + ---------- + image : numpy.ndarray + Stamp image. + defect : numpy.ndarray of bool + Pixels to replace (:func:`defect_mask`). + noise : numpy.ndarray + Noise realisation on the stamp grid. + + Returns + ------- + numpy.ndarray + ``image`` with ``defect`` pixels taken from ``noise``. + """ + return np.where(defect, noise, image) + + def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, - blend_handling="noisefill", seg=None, object_number=None, + blend_handling="none", seg=None, object_number=None, dilate_neighbour=0, ): - """bookkeeping for ngmix weights. runs on a single galaxy and epoch - pixel scale and galaxy guess - TO DO: decide if we want galaxy guess stuff + """Build one epoch's image, weight map and noise image for ngmix. + + Defect pixels (:func:`defect_mask`: flagged, zero-weight or invalid-RMS + pixels) get weight 0 and are replaced by an independent noise + realisation at their background RMS (:func:`fill_defects`), under either + ``blend_handling``. ``blend_handling`` decides only the neighbour side. + + @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-neighbours-raw + Every pixel of ``defect_mask`` is zero-weighted and noise-filled whatever + ``blend_handling`` is, and the filled set equals the zero-weight defect + set. ``blend_handling`` acts on the neighbour side only: uberseg zeroes + the weight of pixels nearer a neighbour's footprint and leaves their + image values raw, because that light is real sky that metacal shears + along with the target; filling it with noise would cut the target along + an unsheared edge. Neighbour pixels are never noise-filled. The noise + image covers the whole stamp and is independent of both. Parameters ---------- gal : numpy.ndarray + Background-subtracted galaxy stamp. weight : numpy.ndarray + Exposure weight stamp; zero marks a defect. flag : numpy.ndarray + Flag stamp; nonzero marks a defect. rng : numpy.random.RandomState Random state for the noise realisations (seeded per object; see :func:`position_seed`). bkg_rms : numpy.ndarray, optional - Per-pixel background RMS map. If supplied, unmasked pixels use - ``1 / bkg_rms**2`` as the ngmix inverse variance. - blend_handling : {"noisefill", "uberseg"}, optional - How to treat pixels shared with a neighbour. ``"noisefill"`` (default) - replaces flagged pixels with a noise realisation and keeps their - inverse-variance weight — the historical behaviour. ``"uberseg"`` - instead hard-masks (weight = 0) every pixel closer to a neighbour's - segmentation footprint than to the central object's, leaving the - image untouched (see :func:`uberseg_weight`). + Per-pixel background RMS map. If supplied, clean pixels use + ``1 / bkg_rms**2`` as the ngmix inverse variance, and non-finite or + non-positive values mark defects. Otherwise every clean pixel gets + ``1 / sigma_mad(gal)**2``. + blend_handling : {"none", "uberseg"}, optional + Neighbour treatment. ``"none"`` (default) leaves neighbour + pixels weighted and untouched. ``"uberseg"`` zeroes the weight of + every pixel closer to a neighbour's segmentation footprint than to the + central object's and keeps its raw image value (see + :func:`uberseg_weight`). seg : numpy.ndarray, optional Segmentation map on the stamp grid (object NUMBERs). Required for ``blend_handling="uberseg"``; ignored otherwise. @@ -1606,44 +1822,57 @@ def prepare_ngmix_weights( Returns ------- numpy.ndarray - Galaxy image. For ``"noisefill"`` masked pixels are replaced by noise; - for ``"uberseg"`` the image is returned untouched. + Galaxy image with defect pixels noise-filled. numpy.ndarray - Variance map for NGMIX. + Inverse-variance weight map for ngmix. numpy.ndarray - Noise image. + Noise image: an independent realisation over the whole stamp, for + metacal's ``fixnoise``. + + Raises + ------ + ValueError + If ``blend_handling`` is unknown, or ``"uberseg"`` lacks ``seg`` or + ``object_number``. """ if blend_handling not in BLEND_HANDLINGS: raise ValueError( f"Unknown blend_handling '{blend_handling}'; expected one of" + f" {BLEND_HANDLINGS}" ) + if blend_handling == "uberseg" and (seg is None or object_number is None): + raise ValueError( + "blend_handling='uberseg' requires a segmentation map and the" + + " central object_number; none reached prepare_ngmix_weights." + + " Set SEG_VIGNET_PATH on the ngmix run (see" + + " CosmoStat/shapepipe#776)." + ) - mask = np.copy(weight) != 0 - mask[flag != 0] = False + defect = defect_mask(weight, flag, bkg_rms) + clean = ~defect if bkg_rms is None: sig_noise = sigma_mad(gal) # Guard the degenerate constant stamp (sigma_mad == 0): 0 * inf # would otherwise put NaN in a fully-masked weight map. weight_map = ( - mask.astype(float) / sig_noise ** 2 + clean.astype(float) / sig_noise ** 2 if sig_noise > 0 else np.zeros_like(gal, dtype=float) ) else: - valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) - mask &= valid_rms weight_map = np.zeros_like(gal, dtype=float) - weight_map[mask] = 1.0 / bkg_rms[mask] ** 2 + weight_map[clean] = 1.0 / bkg_rms[clean] ** 2 # Per-pixel noise sigma for the realisations below: metacal's # fixnoise bookkeeping (1/w + 1/w_noise) assumes the noise image # is a faithful realisation of the per-pixel variance the weights # claim; a scalar sigma there mis-reports errors and erodes the # inverse-variance advantage whenever the RMS map actually varies. + # Pixels without a valid RMS take the median over clean pixels. + valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) sig_noise = ( - np.where(valid_rms, bkg_rms, np.median(bkg_rms[mask])) - if mask.any() + np.where(valid_rms, bkg_rms, np.median(bkg_rms[clean])) + if clean.any() else sigma_mad(gal) ) @@ -1656,33 +1885,20 @@ def prepare_ngmix_weights( noise_img = rng.standard_normal(gal.shape) * sig_noise noise_img_gal = rng.standard_normal(gal.shape) * sig_noise + gal_filled = fill_defects(gal, defect, noise_img_gal) - gal_masked = np.copy(gal) if blend_handling == "uberseg": - # Hard-mask neighbour-side pixels (weight -> 0) from the segmentation - # geometry; the image is left untouched (the masked pixels carry no - # weight, so ngmix ignores them in the likelihood). Bad/flagged - # pixels already sit at weight 0 from the mask above. - if seg is None or object_number is None: - raise ValueError( - "blend_handling='uberseg' requires a segmentation map and the" - + " central object_number; none reached prepare_ngmix_weights." - + " Set SEG_VIGNET_PATH on the ngmix run (see" - + " CosmoStat/shapepipe#776)." - ) weight_map = uberseg_weight( weight_map, seg, object_number, dilate_neighbour=dilate_neighbour ) - elif (~mask).any(): - # noisefill (default): replace masked pixels with a noise realisation. - gal_masked[~mask] = noise_img_gal[~mask] - return gal_masked, weight_map, noise_img + return gal_filled, weight_map, noise_img + def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, - blend_handling="noisefill", seg=None, object_number=None, + blend_handling="none", seg=None, object_number=None, dilate_neighbour=0, ): """Build an ngmix Observation for a single galaxy epoch. @@ -1725,9 +1941,9 @@ def make_ngmix_observation( Sub-pixel ``[row, col]`` coadd-centroid offset propagated from the stamp extractor. Required for ``centroid_source="wcs"`` (ignored for ``"hsm"``). - blend_handling : {"noisefill", "uberseg"}, optional + blend_handling : {"none", "uberseg"}, optional Neighbour treatment passed through to :func:`prepare_ngmix_weights`; - the default ``"noisefill"`` is the historical behaviour. + the default ``"none"`` leaves neighbour pixels untouched. seg : numpy.ndarray, optional Segmentation map on the stamp grid. Required for ``blend_handling="uberseg"`` (ignored otherwise). @@ -1769,11 +1985,12 @@ def make_ngmix_observation( if centroid_source == "hsm": # Re-center the Jacobian on the HSM adaptive-moment centroid (pixel - # offset from the stamp center); fall back to the stamp center if - # HSM fails. + # offset from the stamp center), measured on the filled image so no + # raw defect value pulls it; fall back to the stamp center if HSM + # fails. try: _hsm = galsim.hsm.FindAdaptiveMom( - galsim.Image(gal, scale=1.0), strict=False + galsim.Image(gal_masked, scale=1.0), strict=False ) if _hsm.error_message != "": raise galsim.hsm.GalSimHSMError(_hsm.error_message) @@ -1988,7 +2205,7 @@ def make_runners(prior, flux_guess, rng): def do_ngmix_metacal( stamp, prior, flux_guess, rng, centroid_source="wcs", - blend_handling="noisefill", object_number=None, dilate_neighbour=0, + blend_handling="none", object_number=None, dilate_neighbour=0, metacal_psf="fitgauss", ): """Do Ngmix Metacal. @@ -2012,10 +2229,10 @@ def do_ngmix_metacal( coadd-centroid offset in ``stamp.offsets``, propagated from the stamp extractor); ``"hsm"`` uses the adaptive-moment centroid from the stamp pixels — see that function. - blend_handling : {"noisefill", "uberseg"}, optional + blend_handling : {"none", "uberseg"}, optional Neighbour treatment passed through to - :func:`make_ngmix_observation`; the default ``"noisefill"`` is the - historical behaviour. ``"uberseg"`` consumes ``stamp.segs`` and + :func:`make_ngmix_observation`; the default ``"none"`` leaves + neighbour pixels untouched. ``"uberseg"`` consumes ``stamp.segs`` and ``object_number``. object_number : int, optional Central object's segmentation label — its SExtractor ``NUMBER`` diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 27648947f..ffc4de083 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -11,7 +11,11 @@ from sqlitedict import SqliteDict from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.ngmix_package.ngmix import Ngmix +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + Ngmix, +) @module_runner( @@ -125,14 +129,16 @@ def ngmix_runner( else: centroid_source = "wcs" - # Neighbour treatment: "noisefill" (default, historical) replaces a - # neighbour's pixels with a noise realisation; "uberseg" hard-masks - # (weight -> 0) every pixel closer to a neighbour than to the central - # object, from the segmentation map. See the ngmix module docstrings. + # Neighbour treatment: "none" (default) leaves neighbour pixels + # weighted and untouched; "uberseg" zeroes the weight of every pixel + # closer to a neighbour than to the central object, from the segmentation + # map, and leaves its image raw. Defect pixels (flagged, zero-weight or + # invalid-RMS) are zero-weighted and noise-filled under both; see + # prepare_ngmix_weights. if config.has_option(module_config_sec, "BLEND_HANDLING"): blend_handling = config.get(module_config_sec, "BLEND_HANDLING") else: - blend_handling = "noisefill" + blend_handling = "none" # DILATE_NEIGHBOUR (optional): binary-dilation iterations enlarging the # uberseg neighbour mask, to absorb the few-pixel coadd-vs-epoch seg-overlay @@ -142,6 +148,24 @@ def ngmix_runner( else: dilate_neighbour = 1 + # EPOCH_CENTRAL_DEFECT_RADIUS (optional, pixels): drop an epoch when a + # masked pixel lies closer than this to the stamp centre; 0 disables. + if config.has_option(module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS"): + epoch_central_defect_radius = config.getfloat( + module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS" + ) + else: + epoch_central_defect_radius = EPOCH_CENTRAL_DEFECT_RADIUS + + # EPOCH_MASKED_FRACTION_CUT (optional): drop an epoch when more than this + # fraction of its stamp is masked (flagged, zero-weight or invalid-RMS). + if config.has_option(module_config_sec, "EPOCH_MASKED_FRACTION_CUT"): + epoch_masked_fraction_cut = config.getfloat( + module_config_sec, "EPOCH_MASKED_FRACTION_CUT" + ) + else: + epoch_masked_fraction_cut = EPOCH_MASKED_FRACTION_CUT + # Check PSF vignets first: if all are empty dicts {}, the exposures for this # tile are absent from the PSF dictionary and no shape measurement is possible. # This check must come before reading image vignets to avoid a C-level malloc @@ -197,6 +221,8 @@ def ngmix_runner( seg_cat_path=seg_vignet_path, dilate_neighbour=dilate_neighbour, metacal_psf=metacal_psf, + epoch_central_defect_radius=epoch_central_defect_radius, + epoch_masked_fraction_cut=epoch_masked_fraction_cut, ) # Process ngmix shape measurement and metacalibration From 7777181bc6cd76f3c0b7a41d4b98bccb218442b1 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:07:40 +0200 Subject: [PATCH 03/14] ngmix: read each epoch's OFFSET from the object dict in hand prepare_postage_stamps already holds the object's galaxy vignet dict, so the per-epoch OFFSET is read from it instead of a second sqlitedict lookup. The pixel scale sets only the centroid-prior width; the comments that also credited it with a noise window are corrected. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01QAc1ywyd9Sbbra6po2iuaA --- .../modules/ngmix_package/__init__.py | 4 ++-- src/shapepipe/modules/ngmix_package/ngmix.py | 18 ++++++++---------- src/shapepipe/modules/ngmix_runner.py | 2 +- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/shapepipe/modules/ngmix_package/__init__.py b/src/shapepipe/modules/ngmix_package/__init__.py index 66323822c..012b3264b 100644 --- a/src/shapepipe/modules/ngmix_package/__init__.py +++ b/src/shapepipe/modules/ngmix_package/__init__.py @@ -36,8 +36,8 @@ PIXEL_SCALE : float, optional Pixel scale in arcsec. Optional override; when omitted (or non-positive) it is read from the image WCS so it cannot drift from the pixels. Only - sets the centroid-prior width and noise window -- the fit Jacobian is - built per object from the full WCS. + sets the centroid-prior width -- the fit Jacobian is built per object + from the full WCS. SAVE_BATCH : int, optional Save the output catalogue in batches of this size; default is ``-1`` (no batch saving) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 807b99ff0..a19d2a83e 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -367,11 +367,11 @@ def pixel_scale_from_wcs(f_wcs_file): """Representative pixel scale (arcsec) from the tile's image WCS. The ngmix fit builds each object's Jacobian from the full per-epoch WCS, - so this scalar only sets the centroid-prior width and the noise-window - scale (see :func:`get_prior`, :func:`get_noise`). A single value read from - the astrometry is therefore sufficient -- and, unlike a hard-coded config - constant, it cannot silently drift from the pixels it describes (this - mirrors SExtractor's ``PIXEL_SCALE 0`` convention). + so this scalar only sets the centroid-prior width (see :func:`get_prior`). + A single value read from the astrometry is therefore sufficient -- and, + unlike a hard-coded config constant, it cannot silently drift from the + pixels it describes (this mirrors SExtractor's ``PIXEL_SCALE 0`` + convention). Parameters ---------- @@ -578,8 +578,8 @@ def __init__( # Pixel scale: an explicit positive PIXEL_SCALE overrides; otherwise # derive it from the image WCS so it can never drift from the pixels # (mirrors SExtractor's ``PIXEL_SCALE 0`` convention). Only the - # centroid-prior width and noise window use it -- the fit Jacobian is - # built per object from the full WCS. + # centroid-prior width uses it -- the fit Jacobian is built per + # object from the full WCS. if pixel_scale is None or pixel_scale <= 0: self._pixel_scale = pixel_scale_from_wcs( self._vignet_cat.f_wcs_file @@ -1359,9 +1359,7 @@ def prepare_postage_stamps( # make_ngmix_observation), which raises if it is missing; the "hsm" # path ignores it, so read it leniently rather than coupling hsm to a # field it never uses. - stamp.offsets.append( - vignet.gal_vign_cat[str(obj_id)][expccd_name].get('OFFSET') - ) + stamp.offsets.append(gal_obj[expccd_name].get('OFFSET')) stamp.ra.append(tile_cat.ra[i_tile]) stamp.dec.append(tile_cat.dec[i_tile]) # CCD of the first surviving epoch — Fabian's coord_list[0] convention diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index ffc4de083..605262395 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -54,7 +54,7 @@ def ngmix_runner( # Pixel scale -- optional override. When absent (or non-positive) it is # derived from the image WCS inside Ngmix, so it cannot drift from the - # pixels. Only the centroid-prior width and noise window use it. + # pixels. Only the centroid-prior width uses it. if config.has_option(module_config_sec, "PIXEL_SCALE"): pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") else: From 9f0ea4430f32f998a3a9e83fd1466cbeab69546d Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:44:32 +0200 Subject: [PATCH 04/14] tests: name the none-mode uberseg test for none, not noisefill Co-Authored-By: Claude Opus 5.5 --- tests/module/test_ngmix_uberseg.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 86dad33bc..4d6ebb32a 100644 --- a/tests/module/test_ngmix_uberseg.py +++ b/tests/module/test_ngmix_uberseg.py @@ -239,9 +239,9 @@ def _gal_flag_weight(npix=41, seed=7): return gal, flag, weight -def test_noisefill_ignores_seg_and_dilate_kwargs(): - """Under noisefill, passing seg / dilate_neighbour changes nothing: the - result matches the plain default call on the same RNG stream.""" +def test_none_ignores_seg_and_dilate_kwargs(): + """Under BLEND_HANDLING = none, passing seg / dilate_neighbour changes + nothing: the result matches the plain default call on the same RNG stream.""" gal, flag, weight = _gal_flag_weight() seg, _, _ = two_object_seg(npix=gal.shape[0], sep=12) From a92a7471fbb09a067512ffb8f13a7c775426cac1 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:22:49 +0200 Subject: [PATCH 05/14] test(ngmix): specify DEFECT_FILL = interpolate and its shear recovery The interpolation fills only short defect runs (at most 3 px along a row or column, with clean pixels at both ends); anything wider, and runs reaching the stamp border, stays noise-filled. The interpolant must reproduce linear planes without reading defect values, commute with quarter turns of the stamp, and apply one operator to the science and metacal noise images. The weight is zero on the defects and on the quarter-turn orbit of the interpolated pixels, whose light stays. Each defect is vetoed at the radius of its own fill: EPOCH_INTERPOLATED_DEFECT_RADIUS for interpolated pixels, EPOCH_CENTRAL_DEFECT_RADIUS for noise-filled ones. The runner and Ngmix pass DEFECT_FILL to the epoch cuts and to metacal; noise stays the default. The science test recovers the full response matrix for interpolated columns, full and finite 3-px bleeds and single pixels at the interpolated radius, on 0.3" and 0.5" galaxies through round and elliptical 0.7" PSFs and on a 0.7" galaxy through a 0.9" PSF, and for noise-filled wide defects at the noise radius. A bleed three pixels inside the radius is the positive control. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01QAc1ywyd9Sbbra6po2iuaA --- tests/module/test_defect_interpolation.py | 237 +++++++++++++++++++++ tests/module/test_ngmix_defect_fill.py | 223 ++++++++++++++++++- tests/science/test_defect_interpolation.py | 119 +++++++++++ 3 files changed, 568 insertions(+), 11 deletions(-) create mode 100644 tests/module/test_defect_interpolation.py create mode 100644 tests/science/test_defect_interpolation.py diff --git a/tests/module/test_defect_interpolation.py b/tests/module/test_defect_interpolation.py new file mode 100644 index 000000000..1368e78cc --- /dev/null +++ b/tests/module/test_defect_interpolation.py @@ -0,0 +1,237 @@ +"""Defect interpolation (``DEFECT_FILL = interpolate``): which pixels are +interpolated, and the properties of the interpolant. + +:func:`interpolable_defects` picks the defect pixels that lie in a short row +or column run with clean pixels at both ends; :func:`interpolate_defects` +fills them from the nearby clean pixels with a Clough-Tocher interpolant, +averaged over the four quarter turns of the stamp and shared by every plane +(science image and metacal noise image). +""" + +import numpy as np +import numpy.testing as npt +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from shapepipe.modules.ngmix_package.defect_interpolation import ( + MAX_INTERPOLATED_RUN, + fourfold, + interpolable_defects, + interpolate_defects, +) + + +# --- interpolable_defects --------------------------------------------------- + +def _oracle(defect, max_run=MAX_INTERPOLATED_RUN): + """Brute force: walk each defect pixel's row and column run.""" + n0, n1 = defect.shape + out = np.zeros_like(defect) + for i, j in zip(*np.nonzero(defect)): + for di, dj in ((0, 1), (1, 0)): + lo, hi = (i, j), (i, j) + while (0 <= lo[0] - di and 0 <= lo[1] - dj + and defect[lo[0] - di, lo[1] - dj]): + lo = (lo[0] - di, lo[1] - dj) + while (hi[0] + di < n0 and hi[1] + dj < n1 + and defect[hi[0] + di, hi[1] + dj]): + hi = (hi[0] + di, hi[1] + dj) + length = hi[0] - lo[0] + hi[1] - lo[1] + 1 + bounded = (lo[0] - di >= 0 and lo[1] - dj >= 0 + and hi[0] + di < n0 and hi[1] + dj < n1) + if bounded and length <= max_run: + out[i, j] = True + return out + + +@given( + n=st.integers(5, 21), + density=st.floats(0.0, 0.6), + seed=st.integers(0, 2**31 - 1), +) +@settings(max_examples=60, deadline=None) +def test_interpolable_defects_are_the_short_bounded_runs(n, density, seed): + """A defect pixel is interpolated exactly when its row or column run is + at most MAX_INTERPOLATED_RUN long and has clean pixels at both ends; the + rule commutes with quarter turns. + + Failure modes: a run touching the stamp border (edge band, corner) or a + wide hole is interpolated from one side; a 3-px bleed is left to noise; + a clean pixel is selected; one axis is ignored, so the rule has a + preferred direction. + """ + defect = np.random.RandomState(seed).uniform(size=(n, n)) < density + out = interpolable_defects(defect) + npt.assert_array_equal(out, _oracle(defect)) + assert not out[~defect].any() + for k in range(1, 4): + npt.assert_array_equal( + interpolable_defects(np.rot90(defect, k)), np.rot90(out, k) + ) + + +@pytest.mark.parametrize("kind,expected", [ + ("column", True), ("bleed3", True), ("finite_bleed3", True), + ("pixel", True), ("bleed4", False), ("blob5", False), + ("edge_band", False), ("corner", False), +]) +def test_calibrated_widths_are_interpolated_and_wider_holes_are_not( + kind, expected, +): + """Columns, 3-px bleeds and single pixels are interpolated; 4-px bleeds, + blobs, edge bands and corners are not.""" + n, c = 51, 25 + defect = np.zeros((n, n), dtype=bool) + region = { + "column": np.s_[:, c + 8], + "bleed3": np.s_[:, c + 8:c + 11], + "finite_bleed3": np.s_[c - 5:c + 6, c + 8:c + 11], + "pixel": np.s_[c, c + 8], + "bleed4": np.s_[:, c + 8:c + 12], + "blob5": np.s_[c - 2:c + 3, c + 8:c + 13], + "edge_band": np.s_[:, -3:], + "corner": np.s_[:2, :2], + }[kind] + defect[region] = True + out = interpolable_defects(defect) + assert out[defect].all() if expected else not out.any() + + +# --- fourfold --------------------------------------------------------------- + +@given(n=st.integers(2, 20), seed=st.integers(0, 2**31 - 1)) +@settings(max_examples=30, deadline=None) +def test_fourfold_is_the_quarter_turn_orbit(n, seed): + """The union contains the mask, is invariant under quarter turns, and is + the smallest such set (idempotent).""" + mask = np.random.RandomState(seed).uniform(size=(n, n)) < 0.1 + out = fourfold(mask) + assert out[mask].all() + for k in range(1, 4): + npt.assert_array_equal(np.rot90(out, k), out) + npt.assert_array_equal(fourfold(out), out) + npt.assert_array_equal( + out, mask | np.rot90(mask) | np.rot90(mask, 2) | np.rot90(mask, 3) + ) + + +def test_fourfold_rejects_rectangular_stamps(): + with pytest.raises(ValueError, match="square"): + fourfold(np.zeros((11, 12), dtype=bool)) + + +# --- interpolate_defects ---------------------------------------------------- + +def _mask(n=31): + """A column, a finite 3-px bleed and a single pixel, all bounded.""" + defect = np.zeros((n, n), dtype=bool) + defect[:, 19] = True + defect[4:12, 7:10] = True + defect[22, 11] = True + return defect + + +@given(seed=st.integers(0, 2**31 - 1)) +@settings(max_examples=20, deadline=None) +def test_interpolation_reproduces_planes_without_reading_defects(seed): + """Linear planes are reproduced at the interpolated pixels; clean pixels, + and defect pixels not selected, are returned untouched; the inputs are + not modified; and no defect value (NaN or a sentinel) is ever read. + + Failure modes: defect pixels enter the support; coordinates are + transposed or mis-rotated; the fill smooths clean pixels. + """ + n = 31 + rng = np.random.RandomState(seed) + rows, cols = np.indices((n, n)) + coeff = rng.uniform(-2, 2, (2, 3)) + planes = np.array([a + b * rows + c * cols for a, b, c in coeff]) + defect = _mask(n) + defect[:, -2:] = True # an edge band and a blob next to the column: + defect[13:18, 21:26] = True # defects that are not targets + target = interpolable_defects(defect) + assert target[15, 19] + assert not target[13:18, 21:26].any() and not target[:, -2:].any() + contaminated = planes.copy() + contaminated[:, defect] = np.nan + saved = contaminated.copy() + + out = interpolate_defects(contaminated, defect, target) + + npt.assert_allclose(out[:, target], planes[:, target], atol=1e-5) + npt.assert_array_equal(out[:, ~target], contaminated[:, ~target]) + npt.assert_array_equal(contaminated, saved) + contaminated[:, defect] = 1e30 + npt.assert_array_equal( + interpolate_defects(contaminated, defect, target)[:, target], + out[:, target], + ) + + +def test_interpolation_commutes_with_quarter_turns(): + """Rotating the stamp and its mask rotates the fill. A regular grid's + Delaunay triangulation has degenerate diagonals, so a single-orientation + interpolant does not commute; the four-orientation average does. + + Failure mode: the average over orientations is skipped, so the fill has a + preferred direction. + """ + n = 31 + rows, cols = np.indices((n, n)) + image = np.exp(-((rows - 15.3) ** 2 + (cols - 14.6) ** 2) / 10.0) + image += 0.05 * np.random.RandomState(3).normal(size=(n, n)) + planes = image[None] + defect = _mask(n) + target = interpolable_defects(defect) + out = interpolate_defects(planes, defect, target) + for k in range(1, 4): + rotated = interpolate_defects( + np.rot90(planes, k, axes=(1, 2)), np.rot90(defect, k), + np.rot90(target, k), + ) + npt.assert_allclose( + rotated, np.rot90(out, k, axes=(1, 2)), atol=1e-12, rtol=0 + ) + + +def test_every_plane_sees_the_same_operator(): + """The fill is one linear operator applied to every plane: filling + a * image + b * noise gives a * fill(image) + b * fill(noise), up to the + Clough-Tocher gradient solver's tolerance. + + Failure mode: the noise image is filled differently from the science + image (another support, triangulation or orientation set), so metacal's + fixnoise no longer mirrors the science image's correlated noise. + """ + n = 31 + rng = np.random.RandomState(11) + image, noise = rng.normal(size=(2, n, n)) + defect = _mask(n) + target = interpolable_defects(defect) + both = interpolate_defects(np.array([image, noise]), defect, target) + mixed = interpolate_defects( + np.array([2.0 * image - 3.0 * noise]), defect, target + ) + npt.assert_allclose(mixed[0], 2.0 * both[0] - 3.0 * both[1], atol=1e-5) + alone = interpolate_defects(noise[None], defect, target) + npt.assert_allclose(alone[0], both[1], atol=1e-5) + + +def test_no_target_is_a_no_op(): + planes = np.random.RandomState(5).normal(size=(2, 15, 15)) + defect = np.zeros((15, 15), dtype=bool) + defect[:, -3:] = True + out = interpolate_defects(planes, defect, np.zeros_like(defect)) + npt.assert_array_equal(out, planes) + + +def test_degenerate_support_gives_nan_not_an_error(): + """With fewer than three non-collinear clean pixels the target is NaN, + which the caller replaces by noise.""" + defect = np.ones((5, 5), dtype=bool) + defect[2, 1] = defect[2, 3] = False + target = np.zeros_like(defect) + target[2, 2] = True + out = interpolate_defects(np.ones((1, 5, 5)), defect, target) + assert np.isnan(out[0, 2, 2]) diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 0480c7ef7..cbbbaa3b0 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -2,12 +2,14 @@ A defect is a stamp pixel with a nonzero flag, zero exposure weight or an invalid background RMS. Before metacal, :func:`prepare_ngmix_weights` gives -every defect weight 0 and fills it with noise at its background RMS, whatever -``BLEND_HANDLING`` is. Under uberseg, pixels on the neighbour side only lose -their weight, and their image values stay raw. The per-epoch cuts in -:func:`prepare_postage_stamps` act on the same defect set: the masked-fraction -cut counts it, and the central-defect veto drops an epoch with a defect near -the stamp centre. +every defect weight 0 and fills it, whatever ``BLEND_HANDLING`` is: with noise +at its background RMS (``DEFECT_FILL = noise``, the default) or, for short +defect runs, with an interpolant of the clean pixels around them +(``DEFECT_FILL = interpolate``). Under uberseg, pixels on the neighbour side +only lose their weight, and their image values stay raw. The per-epoch cuts +in :func:`prepare_postage_stamps` act on the same defect set: the +masked-fraction cut counts it, and the central-defect veto drops an epoch +with a defect near the stamp centre, at a radius set by the defect's fill. """ import re @@ -15,6 +17,7 @@ import numpy as np import numpy.testing as npt +import pytest from astropy.io import fits from astropy.wcs import WCS from hypothesis import given @@ -22,8 +25,14 @@ from sqlitedict import SqliteDict from shapepipe.modules.ngmix_package import ngmix as ngmix_module +from shapepipe.modules.ngmix_package.defect_interpolation import ( + fourfold, + interpolable_defects, + interpolate_defects, +) from shapepipe.modules.ngmix_package.ngmix import ( EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, Ngmix, make_ngmix_observation, prepare_ngmix_weights, @@ -461,8 +470,9 @@ def getboolean(self, _sec, key, fallback=False): def test_runner_threads_the_epoch_cut_options(tmp_path, monkeypatch): - """EPOCH_CENTRAL_DEFECT_RADIUS and EPOCH_MASKED_FRACTION_CUT reach Ngmix - as configured, and default to the module constants when absent. + """EPOCH_CENTRAL_DEFECT_RADIUS, EPOCH_MASKED_FRACTION_CUT, DEFECT_FILL and + EPOCH_INTERPOLATED_DEFECT_RADIUS reach Ngmix as configured, and default + to the module constants (and the noise fill) when absent. Failure mode: the runner drops or ignores an option, so a configured A/B arm silently runs the default cuts. @@ -483,11 +493,15 @@ def process(self): for path in inputs: SqliteDict(path).close() - for options, radius, fraction in ( + for options, radius, fraction, fill, interpolated_radius in ( ({"EPOCH_CENTRAL_DEFECT_RADIUS": "7.5", - "EPOCH_MASKED_FRACTION_CUT": "0.1"}, 7.5, 0.1), + "EPOCH_MASKED_FRACTION_CUT": "0.1", + "DEFECT_FILL": "interpolate", + "EPOCH_INTERPOLATED_DEFECT_RADIUS": "5.5"}, + 7.5, 0.1, "interpolate", 5.5), ({}, EPOCH_CENTRAL_DEFECT_RADIUS, - ngmix_module.EPOCH_MASKED_FRACTION_CUT), + ngmix_module.EPOCH_MASKED_FRACTION_CUT, "noise", + EPOCH_INTERPOLATED_DEFECT_RADIUS), ): runner_module.ngmix_runner( inputs, {"output": str(tmp_path)}, "-001-001", @@ -495,6 +509,9 @@ def process(self): ) assert captured[-1]["epoch_central_defect_radius"] == radius assert captured[-1]["epoch_masked_fraction_cut"] == fraction + assert captured[-1]["defect_fill"] == fill + assert (captured[-1]["epoch_interpolated_defect_radius"] + == interpolated_radius) # --- make_ngmix_observation: the HSM centroid reads the filled image ------- @@ -530,3 +547,187 @@ def test_hsm_centroid_ignores_raw_defect_values(): npt.assert_allclose( [row - centre, col - centre], [d_row, d_col], atol=0.05 ) + + +# --- DEFECT_FILL = interpolate ---------------------------------------------- + +def _hot_stamp(): + """A bright object with hot defects: a column and a finite 3-px bleed + (interpolated), an edge band and a 5x5 blob (noise-filled).""" + centre = N_STAMP // 2 + rows, cols = np.mgrid[:N_STAMP, :N_STAMP] + gal = 1e3 * np.exp( + -((rows - centre) ** 2 + (cols - centre) ** 2) / (2 * 3.0 ** 2) + ) + flag = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + flag[:, centre + 8] = 1 + flag[centre - 5:centre + 6, centre - 11:centre - 8] = 1 + flag[:, :4] = 1 + flag[centre + 9:centre + 14, centre + 12:centre + 17] = 1 + gal[flag != 0] = 5e4 + return gal, np.ones((N_STAMP, N_STAMP)), flag + + +@pytest.mark.parametrize("blend_handling", ["none", "uberseg"]) +def test_interpolated_fill_and_its_weights(blend_handling): + """Short defect runs take the interpolant of the clean image; other + defects take noise; the weight is zero on the defects and on the + quarter-turn orbit of the interpolated pixels, whose light stays; the + metacal noise image is interpolated with the same operator. + + Failure modes: the orbit is not zero-weighted (a one-sided hole in the + likelihood biases c), or its light is replaced; the fill mask is + symmetrized; wide defects or edge bands are extrapolated; raw defect + values leak; the noise image keeps independent noise where the science + image is smooth. + """ + gal, weight, flag = _hot_stamp() + defect = flag != 0 + target = interpolable_defects(defect) + assert target.any() and (defect & ~target).any() + kwargs = ( + dict(seg=_uberseg_seg(N_STAMP), object_number=1, dilate_neighbour=1) + if blend_handling == "uberseg" + else {} + ) + gal_out, w_out, noise_out = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(4), + bkg_rms=np.ones((N_STAMP, N_STAMP)), blend_handling=blend_handling, + defect_fill="interpolate", **kwargs, + ) + + neighbour = ( + uberseg_weight(np.ones_like(gal), kwargs["seg"], 1, dilate_neighbour=1) + == 0.0 + if kwargs else np.zeros_like(defect) + ) + npt.assert_array_equal(w_out == 0.0, defect | fourfold(target) | neighbour) + npt.assert_array_equal(gal_out[~defect], gal[~defect]) + expected = interpolate_defects(gal[None], defect, target)[0] + npt.assert_allclose(gal_out[target], expected[target], rtol=1e-5) + assert np.all(np.abs(gal_out[defect & ~target]) < 10.0) + refilled = interpolate_defects(noise_out[None], defect, target)[0] + npt.assert_allclose(noise_out[target], refilled[target], atol=1e-5) + + +def test_noise_is_the_default_fill(): + """Leaving DEFECT_FILL unset is bit-identical to the noise fill.""" + gal, weight, flag = _hot_stamp() + rms = np.ones_like(gal) + default = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(9), bkg_rms=rms + ) + noise = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(9), bkg_rms=rms, + defect_fill="noise", + ) + for a, b in zip(default, noise): + npt.assert_array_equal(a, b) + + +def test_unknown_defect_fill_is_rejected(): + gal, weight, flag = _hot_stamp() + with pytest.raises(ValueError, match="DEFECT_FILL"): + prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + defect_fill="interp", + ) + + +def _fill_veto_epochs(): + """A clean epoch; a column at the interpolated-defect radius; a column + one pixel inside it; a 5-px bleed (noise-filled) at the same radius.""" + centre = N_STAMP // 2 + radius = int(EPOCH_INTERPOLATED_DEFECT_RADIUS) + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + at, inside, wide = clean.copy(), clean.copy(), clean.copy() + at[:, centre + radius] = 1 + inside[:, centre + radius - 1] = 1 + wide[:, centre + radius:centre + radius + 5] = 1 + return { + "2100001-10": (clean, ones), + "2100002-11": (at, ones), + "2100003-12": (inside, ones), + "2100004-13": (wide, ones), + } + + +def test_the_veto_radius_follows_the_fill(): + """Under interpolation, an interpolated defect is vetoed inside + EPOCH_INTERPOLATED_DEFECT_RADIUS and a noise-filled one inside + EPOCH_CENTRAL_DEFECT_RADIUS. Under the noise fill every defect keeps the + noise radius. The masked-fraction cut counts the raw defect set in + both modes. + + Failure modes: the interpolated radius is applied to noise-filled pixels + (a wide hole next to the object survives), or not applied at all; the + configured radius is ignored; the boundary is inclusive; the fraction + cut counts the zero-weight orbit. + """ + epochs = _fill_veto_epochs() + assert _surviving(epochs, defect_fill="interpolate") == [ + "2100001-10", "2100002-11" + ] + assert _surviving(epochs) == ["2100001-10"] + assert _surviving( + epochs, defect_fill="interpolate", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS + 1, + ) == ["2100001-10"] + assert _surviving(_epochs(), defect_fill="interpolate") == [ + "2100001-10", "2100002-11" + ] + + +def test_process_threads_the_defect_fill(tmp_path, monkeypatch): + """Ngmix hands DEFECT_FILL to both the epoch cuts and metacal: under + interpolation the column at the interpolated-defect radius survives and + metacal is asked to interpolate; under noise it is vetoed. + + Failure mode: the option is dropped on the way to either, so the tile + silently runs the noise fill or the noise cuts. + """ + epochs = {k: v for k, v in _fill_veto_epochs().items() + if k in ("2100001-10", "2100002-11")} + vignet, tile_cat, psf_obj, _ = _fake_inputs(epochs) + vignet.psf_vign_cat = {"1": psf_obj} + vignet.close = lambda: None + tile_cat.obj_id = np.array([1]) + tile_cat.flux = None + monkeypatch.setattr(ngmix_module, "Tile_cat", lambda *a, **k: tile_cat) + for method in ("compile_results", "save_results", "log_mean_ellipticity"): + monkeypatch.setattr(Ngmix, method, lambda *_a, **_k: None) + calls = [] + + def record(stamp, *_args, **kwargs): + calls.append((len(stamp.gals), kwargs.get("defect_fill"))) + raise RuntimeError("metacal is not under test") + + monkeypatch.setattr(ngmix_module, "do_ngmix_metacal", record) + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + for fill in ("interpolate", "noise"): + ngmix = Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, 0.186, str(paths[4]), + _RecordingLogger(), bkg_sub=False, defect_fill=fill, + ) + ngmix._vignet_cat.close() + ngmix._vignet_cat = vignet + ngmix.process() + assert calls == [(2, "interpolate"), (1, "noise")] + + +def test_ngmix_rejects_an_unknown_defect_fill(tmp_path): + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + with pytest.raises(ValueError, match="DEFECT_FILL"): + Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, 0.186, str(paths[4]), + _RecordingLogger(), bkg_sub=False, defect_fill="interp", + ) diff --git a/tests/science/test_defect_interpolation.py b/tests/science/test_defect_interpolation.py new file mode 100644 index 000000000..d36be39b0 --- /dev/null +++ b/tests/science/test_defect_interpolation.py @@ -0,0 +1,119 @@ +"""Shear recovery for the defects kept under ``DEFECT_FILL = interpolate``. + +Physics invariant: every defect the veto keeps leaves both additive terms +|c1|, |c2| < 5e-4 and both diagonal multiplicative terms |m11|, |m22| < 1%, +from the full 2x2 response matrix. Interpolated defects (columns, full and +finite 3-px bleeds, single pixels) are checked at +``EPOCH_INTERPOLATED_DEFECT_RADIUS`` and two pixels beyond it, on galaxies +with half-light radius 0.3" and 0.5" through a 0.7" PSF, round and with +ellipticity (0.05, 0.02), and on a 0.7" galaxy through a 0.9" PSF. Defects +too wide to interpolate are noise-filled, and are checked at +``EPOCH_CENTRAL_DEFECT_RADIUS``. The cases sit at the radii themselves, so +lowering either below its calibrated value turns this red. + +Positive control: a 3-px bleed three pixels inside the interpolated-defect +radius breaks the bound. +""" + +import json + +import numpy as np +import pytest + +from shapepipe.modules.ngmix_package.defect_interpolation import ( + interpolable_defects, +) +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + central_defect_vetoes, + defect_mask, +) +from tests.helpers.defect_response import defect_response + +N = 51 +CENTRE = N // 2 +RI = int(np.ceil(EPOCH_INTERPOLATED_DEFECT_RADIUS)) +RN = int(np.ceil(EPOCH_CENTRAL_DEFECT_RADIUS)) +ROUND = (0.0, 0.0) +ELLIPTICAL = (0.05, 0.02) +SEEDS = range(6) +INTERPOLATE = {"defect_fill": "interpolate"} + + +def geometry(kind, distance): + """A detector defect whose nearest pixel is ``distance`` px from the + stamp centre.""" + bad = np.zeros((N, N), dtype=bool) + near = CENTRE + distance + if kind == "pixel": + bad[CENTRE, near] = True + elif kind == "column": + bad[:, near] = True + elif kind == "bleed": + bad[:, near:near + 3] = True + elif kind == "finite_bleed": + bad[CENTRE - 5:CENTRE + 6, near:near + 3] = True + elif kind == "wide_bleed": + bad[:, near:near + 5] = True + elif kind == "edge": + bad[:, near:] = True + return bad + + +INTERPOLATED = ("column", "bleed", "finite_bleed", "pixel") +CASES = ( + [(k, d, h, 0.7, ROUND) for k in INTERPOLATED + for d in (RI, RI + 2) for h in (0.3, 0.5)] + + [(k, RI, 0.5, 0.7, ELLIPTICAL) for k in INTERPOLATED] + + [(k, RI, 0.7, 0.9, psf) for k in ("bleed", "finite_bleed") + for psf in (ROUND, ELLIPTICAL)] + + [(k, RN, h, 0.7, ROUND) for k in ("wide_bleed", "edge") + for h in (0.3, 0.5)] +) + + +def case_id(case): + kind, distance, hlr, fwhm, psf_shear = case + psf = "elliptical" if any(psf_shear) else "round" + return f"{kind}-{distance}px-hlr{hlr}-psf{fwhm}-{psf}" + + +def recover(bad, hlr, fwhm, psf_shear, tmp_path): + result = defect_response(bad, hlr=hlr, psf=fwhm, seeds=SEEDS, + psf_shear=psf_shear, options=INTERPOLATE) + (tmp_path / "recovery.json").write_text(json.dumps(result, indent=2)) + return np.abs(result["m"]).max(), np.abs(result["c"]).max(), result + + +@pytest.mark.parametrize("kind,distance,hlr,fwhm,psf_shear", CASES, + ids=[case_id(c) for c in CASES]) +def test_kept_defects_recover_shear_on_both_axes(kind, distance, hlr, fwhm, + psf_shear, tmp_path): + """Failure modes: a veto radius is below its calibrated value; the + interpolated pixels' quarter-turn orbit keeps its weight (a one-sided + hole in the likelihood); the fill mask is symmetrized; a wide hole is + interpolated; raw defect values leak into metacal.""" + bad = geometry(kind, distance) + masked = defect_mask(np.ones((N, N)), bad.astype(np.int32)) + interpolated = interpolable_defects(masked) + assert interpolated.any() == (kind in INTERPOLATED) + assert masked.mean() <= EPOCH_MASKED_FRACTION_CUT + assert not central_defect_vetoes(masked, EPOCH_CENTRAL_DEFECT_RADIUS, + "interpolate") + m, c, result = recover(bad, hlr, fwhm, psf_shear, tmp_path) + assert m < 0.01, result + assert c < 5e-4, result + + +def test_vetoed_bleed_breaks_the_bound(tmp_path): + """Positive control: a 3-px bleed three pixels inside the + interpolated-defect radius, on the 0.3" galaxy, gives |c| > 1e-3 and + |m| > 1.5%. The veto drops it.""" + bad = geometry("bleed", RI - 3) + assert central_defect_vetoes(bad, EPOCH_CENTRAL_DEFECT_RADIUS, + "interpolate") + m, c, result = recover(bad, 0.3, 0.7, ROUND, tmp_path) + assert c > 1e-3, result + assert m > 0.015, result From 1bb6437afb0286ec6a084758e0549603f53b5b0a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 26 Sep 2026 07:22:49 +0200 Subject: [PATCH 06/14] feat(ngmix): DEFECT_FILL = interpolate, with a fill-dependent central veto DEFECT_FILL = interpolate fills short defect runs (at most 3 px along a row or column, with clean pixels at both ends: columns, 3-px bleeds, single pixels) with a Clough-Tocher interpolant of the clean pixels within 4 px. The interpolant is averaged over the four quarter turns of the stamp, and each turn's triangulation serves the science image and the metacal noise image alike. Wider holes and runs reaching the stamp border keep the noise fill. noise stays the default and is unchanged. Interpolation restores the object's light, so the remaining bias is the hole in the likelihood: a one-sided zero-weight column 8 px from a 0.5" galaxy still gives c1 = -1.3e-3. The weight is therefore also zero on the quarter-turn orbit of the interpolated pixels (their light stays), which brings that column to c1 = +2e-6. Symmetrizing the fill instead would trade true light for interpolated light over four times the area. Each defect pixel is vetoed at the radius of its own fill: EPOCH_INTERPOLATED_DEFECT_RADIUS (default 7 px) for interpolated pixels and EPOCH_CENTRAL_DEFECT_RADIUS (10 px) for noise-filled ones. 7 px is the smallest radius at which interpolated columns, 3-px bleeds and single pixels give |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 on 0.3", 0.5" and 0.7" galaxies, round and elliptical PSFs; it does not scale with galaxy size. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01QAc1ywyd9Sbbra6po2iuaA --- .../ngmix_package/defect_interpolation.py | 168 ++++++++++++++ src/shapepipe/modules/ngmix_package/ngmix.py | 208 +++++++++++++++--- src/shapepipe/modules/ngmix_runner.py | 28 ++- 3 files changed, 375 insertions(+), 29 deletions(-) create mode 100644 src/shapepipe/modules/ngmix_package/defect_interpolation.py diff --git a/src/shapepipe/modules/ngmix_package/defect_interpolation.py b/src/shapepipe/modules/ngmix_package/defect_interpolation.py new file mode 100644 index 000000000..713e326f7 --- /dev/null +++ b/src/shapepipe/modules/ngmix_package/defect_interpolation.py @@ -0,0 +1,168 @@ +"""DEFECT INTERPOLATION. + +Clough-Tocher interpolation of short defect runs before metacal, selected +with ``DEFECT_FILL = interpolate`` (see +:func:`shapepipe.modules.ngmix_package.ngmix.prepare_ngmix_weights`). + +""" + +import numpy as np +from scipy.interpolate import CloughTocher2DInterpolator +from scipy.ndimage import binary_dilation, label +from scipy.spatial import QhullError + +# Longest row or column run of defect pixels that is interpolated. +MAX_INTERPOLATED_RUN = 3 + +# Clean pixels within this Chebyshev distance (pixels) of the interpolated +# pixels support the interpolant. +SUPPORT_RADIUS = 4 + +_ROW_RUNS = np.array([[0, 0, 0], [1, 1, 1], [0, 0, 0]]) + + +def _short_row_runs(defect, max_run): + """Pixels in row runs of at most ``max_run`` defects that stop short of + both stamp borders.""" + labels, n_runs = label(defect, structure=_ROW_RUNS) + if n_runs == 0: + return np.zeros_like(defect) + short = np.bincount(labels.ravel(), minlength=n_runs + 1) <= max_run + short[0] = False + short[labels[:, 0]] = False + short[labels[:, -1]] = False + return short[labels] + + +def interpolable_defects(defect, max_run=MAX_INTERPOLATED_RUN): + """Defect pixels that ``DEFECT_FILL = interpolate`` interpolates. + + @sc [decision:shape_measurement.defect_fill] interpolable-defects + A defect pixel is interpolated when its row or its column run of defect + pixels is at most ``max_run`` (3) long and has clean pixels at both + ends. That covers columns, 3-px bleeds and isolated pixels, the widths + whose shear recovery is calibrated. Wider holes, and runs that reach the + stamp border (edge bands, corners), have clean light on one side only; + they are noise-filled and vetoed at the noise-fill radius (see + :func:`~shapepipe.modules.ngmix_package.ngmix.central_defect_vetoes`). + The rule reads only the mask and commutes with quarter turns of the + stamp. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp. + max_run : int, optional + Longest interpolated run; the default is ``MAX_INTERPOLATED_RUN``. + + Returns + ------- + numpy.ndarray of bool + ``True`` on the defect pixels to interpolate. + """ + defect = np.asarray(defect, dtype=bool) + return ( + _short_row_runs(defect, max_run) + | _short_row_runs(defect.T, max_run).T + ) + + +def fourfold(mask): + """Union of a square stamp mask with its quarter turns. + + Parameters + ---------- + mask : numpy.ndarray of bool + Square mask. + + Returns + ------- + numpy.ndarray of bool + ``mask`` ORed with its rotations by 90, 180 and 270 degrees about + the stamp centre. + + Raises + ------ + ValueError + If ``mask`` is not square. + """ + mask = np.asarray(mask, dtype=bool) + if mask.ndim != 2 or mask.shape[0] != mask.shape[1]: + raise ValueError( + f"A quarter-turn orbit needs a square stamp, not {mask.shape}" + ) + return mask | np.rot90(mask) | np.rot90(mask, 2) | np.rot90(mask, 3) + + +def _interpolate_once(planes, defect, target): + """One Clough-Tocher interpolant of every plane at ``target``; NaN + elsewhere and where the support cannot reach.""" + out = np.full(planes.shape, np.nan) + support = binary_dilation( + target, structure=np.ones((3, 3), dtype=bool), + iterations=SUPPORT_RADIUS, + ) & ~defect + points = np.argwhere(support).astype(float) + if len(points) < 3: + return out + query = np.argwhere(target) + try: + interpolant = CloughTocher2DInterpolator( + points, planes[:, support].T, fill_value=np.nan, + ) + except QhullError: + return out + out[:, query[:, 0], query[:, 1]] = interpolant(query.astype(float)).T + return out + + +def interpolate_defects(planes, defect, target): + """Replace the ``target`` pixels of every plane by a Clough-Tocher + interpolant of the clean pixels around them. + + @sc [decision:shape_measurement.defect_fill] shared-rotation-averaged-interpolant + The support is the clean pixels within ``SUPPORT_RADIUS`` (4 px) of the + target; no defect pixel enters it, so defect values are never read. For + each quarter turn of the stamp, one Delaunay triangulation of the support + serves every plane, so the science image and the metacal noise image + see the same linear operator and fixnoise mirrors the science image's + interpolated noise. A regular grid's triangulation has degenerate + diagonals, so one orientation has a preferred direction; averaging the + four quarter-turned operators makes the fill commute with rotations of + the stamp. That removes up to 6e-5 of c2 for a column or 3-px bleed + 6 px from the object. + + Parameters + ---------- + planes : array_like + Stamp planes, shape ``(n, ny, nx)``. + defect : numpy.ndarray of bool + Every defect pixel, shape ``(ny, nx)``; none of them supports the + interpolant. + target : numpy.ndarray of bool + Defect pixels to interpolate (:func:`interpolable_defects`). + + Returns + ------- + numpy.ndarray + A copy of ``planes`` with ``target`` pixels interpolated; NaN at a + target pixel whose support is degenerate in some orientation. + """ + planes = np.asarray(planes, dtype=float) + defect = np.asarray(defect, dtype=bool) + target = np.asarray(target, dtype=bool) + out = planes.copy() + if not target.any(): + return out + turns = [ + np.rot90( + _interpolate_once( + np.rot90(planes, k, axes=(1, 2)), np.rot90(defect, k), + np.rot90(target, k), + ), + -k, axes=(1, 2), + )[:, target] + for k in range(4) + ] + out[:, target] = np.mean(turns, axis=0) + return out diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index a19d2a83e..acb425835 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -23,11 +23,20 @@ from scipy.spatial import cKDTree from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package.defect_interpolation import ( + fourfold, + interpolable_defects, + interpolate_defects, +) from shapepipe.pipeline import file_io # Neighbour treatments selectable with the BLEND_HANDLING option. BLEND_HANDLINGS = ("none", "uberseg") +# Defect fills selectable with the DEFECT_FILL option (see +# :func:`prepare_ngmix_weights`). +DEFECT_FILLS = ("noise", "interpolate") + # Default of the EPOCH_MASKED_FRACTION_CUT option: an epoch is dropped when # more than this fraction of its stamp is in :func:`defect_mask` (see # :func:`prepare_postage_stamps`). @@ -38,6 +47,11 @@ # stamp centre (see :func:`has_central_defect`). EPOCH_CENTRAL_DEFECT_RADIUS = 10 +# Default of the EPOCH_INTERPOLATED_DEFECT_RADIUS option (pixels): under +# DEFECT_FILL = interpolate, the veto radius for interpolated defect pixels +# (see :func:`central_defect_vetoes`). +EPOCH_INTERPOLATED_DEFECT_RADIUS = 7 + METACAL_TYPES = ('noshear', '1p', '1m', '2p', '2m') # Noise budget for the PSF observation's flat weight map (psf_wt = @@ -444,8 +458,7 @@ class Ngmix(object): Neighbour treatment. ``"none"`` (default) leaves neighbour pixels untouched; ``"uberseg"`` zeroes the weight of neighbour-side pixels from the coadd segmentation map and requires ``seg_cat_path``. Defect - pixels are noise-filled under both (see - :func:`prepare_ngmix_weights`). + pixels are filled under both (see :func:`prepare_ngmix_weights`). seg_cat_path : str, optional Path to the coadd-frame segmentation VIGNET catalogue (see :class:`Tile_cat`). Required when ``blend_handling="uberseg"``. @@ -460,6 +473,14 @@ class Ngmix(object): Drop an epoch when more than this fraction of its stamp is masked (see :func:`prepare_postage_stamps`); the default is ``EPOCH_MASKED_FRACTION_CUT``. + defect_fill : {"noise", "interpolate"}, optional + How defect pixels are filled before metacal (see + :func:`prepare_ngmix_weights`); the default is ``"noise"``. + epoch_interpolated_defect_radius : float, optional + Under ``defect_fill="interpolate"``, drop an epoch when an + interpolated defect pixel lies closer than this many pixels to the + stamp centre (see :func:`central_defect_vetoes`); the default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. Notes ----- @@ -471,8 +492,8 @@ class Ngmix(object): IndexError If the length of the input file list is incorrect ValueError - If ``blend_handling`` is unknown, or ``"uberseg"`` is selected without - ``seg_cat_path``. + If ``blend_handling`` or ``defect_fill`` is unknown, or ``"uberseg"`` + is selected without ``seg_cat_path``. """ @@ -496,6 +517,8 @@ def __init__( metacal_psf="fitgauss", epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, + defect_fill="noise", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, ): # Base count = catalogue + vignets, excluding the f_wcs headers (passed @@ -514,6 +537,11 @@ def __init__( f"Unknown BLEND_HANDLING '{blend_handling}'; expected one of" + f" {BLEND_HANDLINGS}" ) + if defect_fill not in DEFECT_FILLS: + raise ValueError( + f"Unknown DEFECT_FILL '{defect_fill}'; expected one of" + + f" {DEFECT_FILLS}" + ) # Fail fast at construction (not deep in the per-epoch loop) when # uberseg is requested without its segmentation input (shapepipe#776). @@ -567,6 +595,10 @@ def __init__( self._metacal_psf = metacal_psf self._epoch_central_defect_radius = epoch_central_defect_radius self._epoch_masked_fraction_cut = epoch_masked_fraction_cut + self._defect_fill = defect_fill + self._epoch_interpolated_defect_radius = ( + epoch_interpolated_defect_radius + ) self._w_log = w_log @@ -1046,6 +1078,10 @@ def process(self): gal_obj, epoch_central_defect_radius=self._epoch_central_defect_radius, epoch_masked_fraction_cut=self._epoch_masked_fraction_cut, + defect_fill=self._defect_fill, + epoch_interpolated_defect_radius=( + self._epoch_interpolated_defect_radius + ), ) epoch_cuts.update(stamp.epoch_cuts) @@ -1092,6 +1128,7 @@ def process(self): object_number=obj_id, dilate_neighbour=self._dilate_neighbour, metacal_psf=self._metacal_psf, + defect_fill=self._defect_fill, ) except Exception as ee: self._w_log.info( @@ -1193,6 +1230,8 @@ def prepare_postage_stamps( gal_obj=None, epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, + defect_fill="noise", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, ): """Gather one object's epoch stamps, dropping epochs its defects spoil. @@ -1200,8 +1239,9 @@ def prepare_postage_stamps( An epoch is dropped when more than ``epoch_masked_fraction_cut`` of its stamp lies in :func:`defect_mask`: flagged, zero-weight and invalid-RMS pixels, the set that :func:`prepare_ngmix_weights` zero-weights and - noise-fills. Counting flags alone would keep epochs whose filled area - exceeds the cut. The default is 1/3; 10%, the DES Y3 and Y6 value, is + fills. Counting flags alone would keep epochs whose filled area exceeds + the cut. The zero-weight orbit of interpolated pixels keeps its light and + is not counted. The default is 1/3; 10%, the DES Y3 and Y6 value, is the alternative to test. Parameters @@ -1225,6 +1265,14 @@ def prepare_postage_stamps( epoch_masked_fraction_cut : float, optional Drop an epoch with more than this fraction of its stamp in :func:`defect_mask`. The default is ``EPOCH_MASKED_FRACTION_CUT``. + defect_fill : {"noise", "interpolate"}, optional + The fill :func:`prepare_ngmix_weights` will apply; it sets the veto + radius of each defect pixel (:func:`central_defect_vetoes`). The + default is ``"noise"``. + epoch_interpolated_defect_radius : float, optional + Veto radius for interpolated defect pixels under + ``defect_fill="interpolate"``. The default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. Returns ------- @@ -1307,14 +1355,17 @@ def prepare_postage_stamps( else None ) # Drop the epoch when too much of it would be zero-weighted and - # noise-filled (epoch-cut-on-defect-mask), or when a filled pixel - # would sit near the object (epoch-central-defect-veto). + # filled (epoch-cut-on-defect-mask), or when a filled pixel would sit + # near the object (veto-radius-follows-the-fill). masked = defect_mask(weight_vign, flag_vign, bkg_rms_vign) stamp.epoch_cuts["considered"] += 1 if masked.mean() > epoch_masked_fraction_cut: stamp.epoch_cuts["masked_fraction"] += 1 continue - if has_central_defect(masked, epoch_central_defect_radius): + if central_defect_vetoes( + masked, epoch_central_defect_radius, defect_fill, + epoch_interpolated_defect_radius, + ): stamp.epoch_cuts["central_veto"] += 1 continue @@ -1734,6 +1785,59 @@ def has_central_defect(defect, radius): return bool(np.any(distance < radius)) +def central_defect_vetoes( + defect, radius, defect_fill="noise", + interpolated_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, +): + """Whether the central-defect veto drops an epoch under ``defect_fill``. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.defect_fill] veto-radius-follows-the-fill + Each defect pixel is vetoed (:func:`has_central_defect`) at the radius + calibrated for the operator that fills it. Under ``"noise"`` every + defect keeps ``radius`` (``EPOCH_CENTRAL_DEFECT_RADIUS``). Under + ``"interpolate"`` the pixels :func:`interpolable_defects` selects are + interpolated and vetoed at ``interpolated_radius`` + (``EPOCH_INTERPOLATED_DEFECT_RADIUS``); the others are noise-filled and + keep ``radius``. 7 px is the smallest interpolated radius at which + columns, full and finite 3-px bleeds and single pixels give + |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 (full response matrix) on + galaxies with half-light radius 0.3" and 0.5" through a 0.7" PSF, round + or with ellipticity (0.05, 0.02), and on a 0.7" galaxy through a 0.9" + PSF. There the worst case is c1 = 3.6e-4 (3-px bleed, 0.7" galaxy, + elliptical PSF); on the 0.3" and 0.5" galaxies |c| <= 5e-5 and + |m| <= 0.21%. At 6 px those two still pass (|c| <= 3.1e-4), but a 3-px + bleed on the 0.7" galaxy gives c1 = 7.5e-4; at 5 px a 3-px bleed on the + 0.5" galaxy gives 9.5e-4. The radius does not scale with galaxy size: + interpolation needs one pixel more for the 0.7" galaxy (noise fill needs + four, 10 to 14 px), and a size-dependent veto would select on measured + size, which responds to shear. Guarded by + ``tests/science/test_defect_interpolation.py``. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp (:func:`defect_mask`). + radius : float + Veto radius in pixels for noise-filled defect pixels. + defect_fill : {"noise", "interpolate"}, optional + The fill the epoch will get; the default is ``"noise"``. + interpolated_radius : float, optional + Veto radius in pixels for interpolated defect pixels; the default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. + + Returns + ------- + bool + ``True`` if the epoch should be dropped. + """ + if defect_fill == "interpolate": + interpolated = interpolable_defects(defect) + return has_central_defect( + interpolated, interpolated_radius + ) or has_central_defect(defect & ~interpolated, radius) + return has_central_defect(defect, radius) + + def fill_defects(image, defect, noise): """Replace an image's defect pixels by a noise realisation. @@ -1766,24 +1870,42 @@ def fill_defects(image, defect, noise): def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, blend_handling="none", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", ): """Build one epoch's image, weight map and noise image for ngmix. Defect pixels (:func:`defect_mask`: flagged, zero-weight or invalid-RMS - pixels) get weight 0 and are replaced by an independent noise - realisation at their background RMS (:func:`fill_defects`), under either - ``blend_handling``. ``blend_handling`` decides only the neighbour side. + pixels) get weight 0 and are filled under either ``blend_handling``: + by an independent noise realisation at their background RMS + (:func:`fill_defects`), or under ``defect_fill="interpolate"``, for + short defect runs, by an interpolant of the clean pixels around them. + ``blend_handling`` decides only the neighbour side. @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-neighbours-raw - Every pixel of ``defect_mask`` is zero-weighted and noise-filled whatever - ``blend_handling`` is, and the filled set equals the zero-weight defect - set. ``blend_handling`` acts on the neighbour side only: uberseg zeroes + Every pixel of ``defect_mask`` is zero-weighted and filled whatever + ``blend_handling`` is, and the filled set equals the defect set. + ``blend_handling`` acts on the neighbour side only: uberseg zeroes the weight of pixels nearer a neighbour's footprint and leaves their image values raw, because that light is real sky that metacal shears along with the target; filling it with noise would cut the target along - an unsheared edge. Neighbour pixels are never noise-filled. The noise - image covers the whole stamp and is independent of both. + an unsheared edge. Neighbour pixels are never filled. The noise image + covers the whole stamp and is independent of both. + + @sc [decision:shape_measurement.defect_fill,label:physics] interpolated-defect-weight-orbit + Under ``defect_fill="interpolate"`` the short defect runs + (:func:`interpolable_defects`) take the interpolant of the clean pixels + around them (:func:`interpolate_defects`) in the image and in the noise + image alike, and the other defects are noise-filled. Interpolation + restores the object's light, which leaves the hole in the likelihood: + ngmix fits a Gaussian to a non-Gaussian profile, and a one-sided + zero-weight hole pulls that fit. Once interpolated, a column 8 px from a + galaxy with half-light radius 0.5" through a 0.7" PSF still gives + c1 = -1.3e-3. The weight is therefore also zero on the quarter-turn + orbit of the interpolated pixels, whose light stays, which cancels that + spin-2 term: c1 = +2e-6 for the same column. Only the weights are + symmetrized, not the fill: interpolating the orbit as well would trade + true light for interpolated light over four times the area, raising m11 + for a 3-px bleed 6 px from the same galaxy from +0.19% to +0.89%. Parameters ---------- @@ -1816,23 +1938,31 @@ def prepare_ngmix_weights( dilate_neighbour : int, optional Neighbour-mask dilation iterations, passed to :func:`uberseg_weight` under ``blend_handling="uberseg"``; ignored otherwise. + defect_fill : {"noise", "interpolate"}, optional + ``"noise"`` (default) noise-fills every defect; ``"interpolate"`` + interpolates the short defect runs and noise-fills the rest. Returns ------- numpy.ndarray - Galaxy image with defect pixels noise-filled. + Galaxy image with defect pixels filled. numpy.ndarray Inverse-variance weight map for ngmix. numpy.ndarray Noise image: an independent realisation over the whole stamp, for - metacal's ``fixnoise``. + metacal's ``fixnoise``, interpolated where the galaxy image is. Raises ------ ValueError - If ``blend_handling`` is unknown, or ``"uberseg"`` lacks ``seg`` or - ``object_number``. + If ``blend_handling`` or ``defect_fill`` is unknown, or + ``"uberseg"`` lacks ``seg`` or ``object_number``. """ + if defect_fill not in DEFECT_FILLS: + raise ValueError( + f"Unknown DEFECT_FILL '{defect_fill}'; expected one of" + + f" {DEFECT_FILLS}" + ) if blend_handling not in BLEND_HANDLINGS: raise ValueError( f"Unknown blend_handling '{blend_handling}'; expected one of" @@ -1848,19 +1978,29 @@ def prepare_ngmix_weights( defect = defect_mask(weight, flag, bkg_rms) clean = ~defect + interpolated = ( + interpolable_defects(defect) + if defect_fill == "interpolate" + else np.zeros_like(defect) + ) + # The quarter-turn orbit of the interpolated pixels carries no weight + # (interpolated-defect-weight-orbit). + weighted = ( + ~(defect | fourfold(interpolated)) if interpolated.any() else clean + ) if bkg_rms is None: sig_noise = sigma_mad(gal) # Guard the degenerate constant stamp (sigma_mad == 0): 0 * inf # would otherwise put NaN in a fully-masked weight map. weight_map = ( - clean.astype(float) / sig_noise ** 2 + weighted.astype(float) / sig_noise ** 2 if sig_noise > 0 else np.zeros_like(gal, dtype=float) ) else: weight_map = np.zeros_like(gal, dtype=float) - weight_map[clean] = 1.0 / bkg_rms[clean] ** 2 + weight_map[weighted] = 1.0 / bkg_rms[weighted] ** 2 # Per-pixel noise sigma for the realisations below: metacal's # fixnoise bookkeeping (1/w + 1/w_noise) assumes the noise image # is a faithful realisation of the per-pixel variance the weights @@ -1884,6 +2024,13 @@ def prepare_ngmix_weights( noise_img = rng.standard_normal(gal.shape) * sig_noise noise_img_gal = rng.standard_normal(gal.shape) * sig_noise gal_filled = fill_defects(gal, defect, noise_img_gal) + if interpolated.any(): + # One operator for the image and the noise image; a pixel whose + # support is degenerate keeps its noise fill. + filled = interpolate_defects([gal, noise_img], defect, interpolated) + done = interpolated & np.all(np.isfinite(filled), axis=0) + gal_filled[done] = filled[0][done] + noise_img = np.where(done, filled[1], noise_img) if blend_handling == "uberseg": weight_map = uberseg_weight( @@ -1897,7 +2044,7 @@ def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, blend_handling="none", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", ): """Build an ngmix Observation for a single galaxy epoch. @@ -1951,6 +2098,9 @@ def make_ngmix_observation( dilate_neighbour : int, optional Neighbour-mask dilation iterations passed through to :func:`prepare_ngmix_weights` under ``blend_handling="uberseg"``. + defect_fill : {"noise", "interpolate"}, optional + Defect fill passed through to :func:`prepare_ngmix_weights`; the + default is ``"noise"``. Returns ------- @@ -1978,7 +2128,7 @@ def make_ngmix_observation( gal_masked, weight_map, noise_img = prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=bkg_rms, blend_handling=blend_handling, seg=seg, object_number=object_number, - dilate_neighbour=dilate_neighbour, + dilate_neighbour=dilate_neighbour, defect_fill=defect_fill, ) if centroid_source == "hsm": @@ -2204,7 +2354,7 @@ def make_runners(prior, flux_guess, rng): def do_ngmix_metacal( stamp, prior, flux_guess, rng, centroid_source="wcs", blend_handling="none", object_number=None, dilate_neighbour=0, - metacal_psf="fitgauss", + metacal_psf="fitgauss", defect_fill="noise", ): """Do Ngmix Metacal. @@ -2251,6 +2401,9 @@ def do_ngmix_metacal( round PSF that metacal reconvolves with after shearing, so it moves the metacal *response* (and therefore the recovered shear) but never reaches the deconvolution, which is by the PSF image. + defect_fill : {"noise", "interpolate"}, optional + Defect fill passed through to :func:`make_ngmix_observation`; the + default is ``"noise"``. Returns ------- @@ -2283,6 +2436,7 @@ def do_ngmix_metacal( seg=stamp.segs[n_e] if n_e < len(stamp.segs) else None, object_number=object_number, dilate_neighbour=dilate_neighbour, + defect_fill=defect_fill, ) gal_obs_list.append(gal_obs) diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 605262395..158e14ef8 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -13,6 +13,7 @@ from shapepipe.modules.module_decorator import module_runner from shapepipe.modules.ngmix_package.ngmix import ( EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, EPOCH_MASKED_FRACTION_CUT, Ngmix, ) @@ -133,8 +134,8 @@ def ngmix_runner( # weighted and untouched; "uberseg" zeroes the weight of every pixel # closer to a neighbour than to the central object, from the segmentation # map, and leaves its image raw. Defect pixels (flagged, zero-weight or - # invalid-RMS) are zero-weighted and noise-filled under both; see - # prepare_ngmix_weights. + # invalid-RMS) are zero-weighted and filled under both (DEFECT_FILL + # below); see prepare_ngmix_weights. if config.has_option(module_config_sec, "BLEND_HANDLING"): blend_handling = config.get(module_config_sec, "BLEND_HANDLING") else: @@ -166,6 +167,27 @@ def ngmix_runner( else: epoch_masked_fraction_cut = EPOCH_MASKED_FRACTION_CUT + # DEFECT_FILL (optional): "noise" (default) noise-fills every defect; + # "interpolate" interpolates short defect runs (at most 3 px along a row + # or column) from the clean pixels around them and noise-fills the rest. + if config.has_option(module_config_sec, "DEFECT_FILL"): + defect_fill = config.get(module_config_sec, "DEFECT_FILL") + else: + defect_fill = "noise" + + # EPOCH_INTERPOLATED_DEFECT_RADIUS (optional, pixels): under + # DEFECT_FILL = interpolate, drop an epoch when an interpolated defect + # pixel lies closer than this to the stamp centre; noise-filled defect + # pixels keep EPOCH_CENTRAL_DEFECT_RADIUS. + if config.has_option( + module_config_sec, "EPOCH_INTERPOLATED_DEFECT_RADIUS" + ): + epoch_interpolated_defect_radius = config.getfloat( + module_config_sec, "EPOCH_INTERPOLATED_DEFECT_RADIUS" + ) + else: + epoch_interpolated_defect_radius = EPOCH_INTERPOLATED_DEFECT_RADIUS + # Check PSF vignets first: if all are empty dicts {}, the exposures for this # tile are absent from the PSF dictionary and no shape measurement is possible. # This check must come before reading image vignets to avoid a C-level malloc @@ -223,6 +245,8 @@ def ngmix_runner( metacal_psf=metacal_psf, epoch_central_defect_radius=epoch_central_defect_radius, epoch_masked_fraction_cut=epoch_masked_fraction_cut, + defect_fill=defect_fill, + epoch_interpolated_defect_radius=epoch_interpolated_defect_radius, ) # Process ngmix shape measurement and metacalibration From b9d92a9643538346289e16145834d0ef8d6d35f0 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:14:13 +0200 Subject: [PATCH 07/14] ngmix: restore the BLEND_HANDLING name noisefill The neighbour option keeps develop's name, noisefill, in the code, the runner config key, tests, astra.yaml and the committed universe; the name none is dropped. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 8 +++--- src/shapepipe/modules/ngmix_package/ngmix.py | 30 ++++++++++---------- src/shapepipe/modules/ngmix_runner.py | 4 +-- tests/module/test_ngmix.py | 2 +- tests/module/test_ngmix_defect_fill.py | 4 +-- tests/module/test_ngmix_uberseg.py | 11 ++----- universes/committed.yaml | 2 +- 7 files changed, 27 insertions(+), 34 deletions(-) diff --git a/astra.yaml b/astra.yaml index 85ee5f48c..978fe94a7 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1392,7 +1392,7 @@ analyses: rationale: >- Covers only pixels shared with a neighbour; defect_fill does not depend on it: defects are noise-filled under every BLEND_HANDLING. - none (the default; the committed config + noisefill (the default; the committed config sets no key) leaves neighbours fully weighted and untouched. uberseg zeroes the weight of pixels nearer a neighbour's coadd segmentation footprint than the target's (DILATE_NEIGHBOUR, default 1) @@ -1415,10 +1415,10 @@ analyses: instead cut the target's own light along an unsheared boundary, a sharp edge that rings in the FFTs. The recommended comparison arm is uberseg (weight-only), with defect_fill held equal across arms. - default: none + default: noisefill options: - none: - label: No neighbour treatment (BLEND_HANDLING = none) + noisefill: + label: No neighbour treatment (BLEND_HANDLING = noisefill) description: >- Neighbour pixels keep their full weight and image values, so the fit sees all neighbour light, which biases shapes toward diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 9ebb018d3..5bd942913 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -31,7 +31,7 @@ from shapepipe.pipeline import file_io # Neighbour treatments selectable with the BLEND_HANDLING option. -BLEND_HANDLINGS = ("none", "uberseg") +BLEND_HANDLINGS = ("noisefill", "uberseg") # Defect fills selectable with the DEFECT_FILL option (see # :func:`prepare_ngmix_weights`). @@ -550,7 +550,7 @@ def get_data(self, cat_path): # Coadd-frame SExtractor segmentation stamp (integer labels), one per # object and row-aligned to the tile catalogue, overlaid unchanged on # every epoch for uberseg neighbour masking (shapepipe#776). None -> - # uberseg unavailable; ``"none"`` does not read it. + # uberseg unavailable; ``"noisefill"`` does not read it. self.seg = None if self.seg_cat_path: seg_cat = file_io.FITSCatalogue( @@ -605,7 +605,7 @@ def __init__( self.flags = [] self.bkg_rms = [] # Segmentation stamps, one per epoch, used only by the "uberseg" blend - # handling; empty under the default "none". All epochs carry the + # handling; empty under the default "noisefill". All epochs carry the # SAME coadd-frame seg stamp (shapepipe#776: one coadd seg per object, # no per-epoch reprojection), each MegaCam-flipped to match its galaxy # stamp so the overlay stays registered. @@ -753,8 +753,8 @@ class Ngmix(object): propagated on the vignette. ``"hsm"`` re-centers on the adaptive-moment centroid measured from the stamp pixels. See :func:`make_ngmix_observation`. - blend_handling : {"none", "uberseg"}, optional - Neighbour treatment. ``"none"`` (default) leaves neighbour pixels + blend_handling : {"noisefill", "uberseg"}, optional + Neighbour treatment. ``"noisefill"`` (default) leaves neighbour pixels untouched; ``"uberseg"`` zeroes the weight of neighbour-side pixels from the coadd segmentation map and requires ``seg_cat_path``. Defect pixels are filled under both (see :func:`prepare_ngmix_weights`). @@ -810,7 +810,7 @@ def __init__( id_obj_max=-1, bkg_sub=True, centroid_source="wcs", - blend_handling="none", + blend_handling="noisefill", seg_cat_path=None, dilate_neighbour=1, metacal_psf="fitgauss", @@ -2144,7 +2144,7 @@ def fill_defects(image, defect, noise): def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, - blend_handling="none", seg=None, object_number=None, + blend_handling="noisefill", seg=None, object_number=None, dilate_neighbour=0, defect_fill="noise", ): """Build one epoch's image, weight map and noise image for ngmix. @@ -2198,8 +2198,8 @@ def prepare_ngmix_weights( ``1 / bkg_rms**2`` as the ngmix inverse variance, and non-finite or non-positive values mark defects. Otherwise every clean pixel gets ``1 / sigma_mad(gal)**2``. - blend_handling : {"none", "uberseg"}, optional - Neighbour treatment. ``"none"`` (default) leaves neighbour + blend_handling : {"noisefill", "uberseg"}, optional + Neighbour treatment. ``"noisefill"`` (default) leaves neighbour pixels weighted and untouched. ``"uberseg"`` zeroes the weight of every pixel closer to a neighbour's segmentation footprint than to the central object's and keeps its raw image value (see @@ -2319,7 +2319,7 @@ def prepare_ngmix_weights( def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, - blend_handling="none", seg=None, object_number=None, + blend_handling="noisefill", seg=None, object_number=None, dilate_neighbour=0, defect_fill="noise", ): """Build an ngmix Observation for a single galaxy epoch. @@ -2362,9 +2362,9 @@ def make_ngmix_observation( Sub-pixel ``[row, col]`` coadd-centroid offset propagated from the stamp extractor. Required for ``centroid_source="wcs"`` (ignored for ``"hsm"``). - blend_handling : {"none", "uberseg"}, optional + blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to :func:`prepare_ngmix_weights`; - the default ``"none"`` leaves neighbour pixels untouched. + the default ``"noisefill"`` leaves neighbour pixels untouched. seg : numpy.ndarray, optional Segmentation map on the stamp grid. Required for ``blend_handling="uberseg"`` (ignored otherwise). @@ -2635,7 +2635,7 @@ def make_runners(prior, flux_guess, rng): def do_ngmix_metacal( stamp, prior, flux_guess, rng, centroid_source="wcs", - blend_handling="none", object_number=None, dilate_neighbour=0, + blend_handling="noisefill", object_number=None, dilate_neighbour=0, metacal_psf="fitgauss", defect_fill="noise", ): """Do Ngmix Metacal. @@ -2659,9 +2659,9 @@ def do_ngmix_metacal( coadd-centroid offset in ``stamp.offsets``, propagated from the stamp extractor); ``"hsm"`` uses the adaptive-moment centroid from the stamp pixels — see that function. - blend_handling : {"none", "uberseg"}, optional + blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to - :func:`make_ngmix_observation`; the default ``"none"`` leaves + :func:`make_ngmix_observation`; the default ``"noisefill"`` leaves neighbour pixels untouched. ``"uberseg"`` consumes ``stamp.segs`` and ``object_number``. object_number : int, optional diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index d350ba6fb..166e8cae5 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -144,7 +144,7 @@ def ngmix_runner( else: centroid_source = "wcs" - # Neighbour treatment: "none" (default) leaves neighbour pixels + # Neighbour treatment: "noisefill" (default) leaves neighbour pixels # weighted and untouched; "uberseg" zeroes the weight of every pixel # closer to a neighbour than to the central object, from the segmentation # map, and leaves its image raw. Defect pixels (flagged, zero-weight or @@ -153,7 +153,7 @@ def ngmix_runner( if config.has_option(module_config_sec, "BLEND_HANDLING"): blend_handling = config.get(module_config_sec, "BLEND_HANDLING") else: - blend_handling = "none" + blend_handling = "noisefill" # DILATE_NEIGHBOUR (optional): binary-dilation iterations enlarging the # uberseg neighbour mask, to absorb the few-pixel coadd-vs-epoch seg-overlay diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index edd54c118..a69aa0e69 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -660,7 +660,7 @@ def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags inst._id_obj_min = inst._id_obj_max = -1 inst._bkg_sub = True inst._pixel_scale = .186 - inst._blend_handling = "none" + inst._blend_handling = "noisefill" inst._dilate_neighbour = 1 inst._metacal_psf = "fitgauss" inst._epoch_central_defect_radius = module.EPOCH_CENTRAL_DEFECT_RADIUS diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index ac10057b7..74cb3045b 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -75,7 +75,7 @@ def _uberseg_seg(n): @given( stamp=defect_stamps(), - blend_handling=st.sampled_from(["none", "uberseg"]), + blend_handling=st.sampled_from(["noisefill", "uberseg"]), seed=st.integers(0, 2**31 - 1), ) def test_filled_set_is_the_defect_set(stamp, blend_handling, seed): @@ -572,7 +572,7 @@ def _hot_stamp(): return gal, np.ones((N_STAMP, N_STAMP)), flag -@pytest.mark.parametrize("blend_handling", ["none", "uberseg"]) +@pytest.mark.parametrize("blend_handling", ["noisefill", "uberseg"]) def test_interpolated_fill_and_its_weights(blend_handling): """Short defect runs take the interpolant of the clean image; other defects take noise; the weight is zero on the defects and on the diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 4d6ebb32a..92616fde0 100644 --- a/tests/module/test_ngmix_uberseg.py +++ b/tests/module/test_ngmix_uberseg.py @@ -239,8 +239,8 @@ def _gal_flag_weight(npix=41, seed=7): return gal, flag, weight -def test_none_ignores_seg_and_dilate_kwargs(): - """Under BLEND_HANDLING = none, passing seg / dilate_neighbour changes +def test_noisefill_ignores_seg_and_dilate_kwargs(): + """Under BLEND_HANDLING = noisefill, passing seg / dilate_neighbour changes nothing: the result matches the plain default call on the same RNG stream.""" gal, flag, weight = _gal_flag_weight() seg, _, _ = two_object_seg(npix=gal.shape[0], sep=12) @@ -406,10 +406,3 @@ def test_runner_missing_seg_vignet_file_raises(tmp_path): _RecordingLogger(), ) - -def test_retired_noisefill_name_is_rejected(): - """Do not silently interpret the retired neighbour option.""" - with pytest.raises(ValueError, match="noisefill"): - prepare_ngmix_weights(np.ones((5, 5)), np.ones((5, 5)), - np.zeros((5, 5)), np.random.RandomState(3), - blend_handling="noisefill") diff --git a/universes/committed.yaml b/universes/committed.yaml index 4e3e02bec..6f2cf8da7 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -52,7 +52,7 @@ analyses: psf_likelihood_noise: psf_noise_1em5 megacam_ccd_flip: megapipe_flip defect_fill: noise - blend_handling: none + blend_handling: noisefill epoch_masked_fraction_cut: one_third central_defect_veto: fixed_radii catalogue_assembly: From 23cc5cb29e1035ac766d4cd1a047216a0548ad6c Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:16:34 +0200 Subject: [PATCH 08/14] test(ngmix): specify neighbour markers apart from defects SExtractor's -1e30 markers in the tile VIGNET form their own per-epoch mask (stamp.neighbours): a footprint 3 px from the centre, or covering 45% of the stamp, drops no epoch; noisefill zero-weights and noise-fills exactly the marked pixels, bit for bit as develop does; uberseg leaves them raw and weighted; a flagged column near the centre is still vetoed; and do_ngmix_metacal hands each epoch's mask to make_ngmix_observation. Co-Authored-By: Claude Opus 5.5 --- tests/module/test_ngmix_defect_fill.py | 267 +++++++++++++++++++++++++ 1 file changed, 267 insertions(+) diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 74cb3045b..b9e3096c5 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -10,6 +10,9 @@ in :func:`prepare_postage_stamps` act on the same defect set: the masked-fraction cut counts it, and the central-defect veto drops an epoch with a defect near the stamp centre, at a radius set by the defect's fill. +SExtractor's -1e30 neighbour markers in the tile VIGNET are not defects: +noisefill zero-weights and noise-fills them, uberseg ignores them, and the +epoch cuts never count them. """ import re @@ -22,6 +25,7 @@ from astropy.wcs import WCS from hypothesis import given from hypothesis import strategies as st +from modopt.math.stats import sigma_mad from sqlitedict import SqliteDict from shapepipe.modules.ngmix_package import ngmix as ngmix_module @@ -735,3 +739,266 @@ def test_ngmix_rejects_an_unknown_defect_fill(tmp_path): str(tmp_path), "-001-001", 30.0, 0.186, str(paths[4]), _RecordingLogger(), bkg_sub=False, defect_fill="interp", ) + + +# --- SExtractor's -1e30 neighbour markers are not defects ------------------- +# +# The tile VIGNET carries -1e30 on the footprints of other detections. Every +# epoch shares that tile stamp, so a marker counted as a defect would drop +# every epoch of an object with a neighbour inside the veto radius. The +# markers are their own per-epoch mask (``stamp.neighbours``): noisefill +# zero-weights and noise-fills them, uberseg ignores them, and the epoch cuts +# never read them. + +_MARKER = -1.0e30 +_CENTRE = N_STAMP // 2 +# Flipped (ccd < 18) and unflipped (ccd >= 18) MegaCam CCDs. +_MARKER_EPOCH_NAMES = ["2100001-10", "2100002-20", "2100003-11"] + + +def _tile_with_neighbour(columns_from=_CENTRE + 3, rows=(_CENTRE - 2, + _CENTRE + 3)): + """Tile VIGNET with a -1e30 neighbour footprint whose nearest pixel is + 3 px from the stamp centre. The footprint is off-centre, so the MegaCam + flip moves it.""" + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[rows[0]:rows[1], columns_from:columns_from + 6] = _MARKER + return tile + + +def _marker_stamp(tile, epochs=None, **kwargs): + """Run prepare_postage_stamps on defect-free epochs under ``tile``.""" + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + if epochs is None: + epochs = {name: (clean.copy(), ones) for name in _MARKER_EPOCH_NAMES} + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + tile_cat.vign = tile[np.newaxis] + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, **kwargs, + ) + return stamp, epochs, gal_obj + + +def _expected_neighbours(tile, name): + return Ngmix.MegaCamFlip(tile, int(name.split("-")[1])) == _MARKER + + +def test_neighbour_markers_near_the_centre_keep_every_epoch(): + """A neighbour footprint 3 px from the centre, and one covering 45% of + the stamp, drop no epoch: the masked-fraction cut and the central veto + count no marker. + + Failure mode: the markers are written into the flag stamp and counted as + defects, so every epoch (they all share the tile VIGNET) is dropped by + the veto or the fraction cut and the object loses its shape + (neighbour-markers-are-not-defects). + """ + small = _tile_with_neighbour() + large = _tile_with_neighbour(rows=(0, N_STAMP)) + large[:, _CENTRE + 3:] = _MARKER + assert (large == _MARKER).mean() > 1 / 3 + for tile in (small, large): + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + assert stamp.epoch_cuts["considered"] == len(_MARKER_EPOCH_NAMES) + assert stamp.epoch_cuts["masked_fraction"] == 0 + assert stamp.epoch_cuts["central_veto"] == 0 + + +def test_neighbour_markers_are_their_own_per_epoch_mask(): + """``stamp.neighbours`` holds the MegaCam-flipped marker mask of each + surviving epoch, and the flag stamps stay the exposure's own. + + Failure mode: the markers are merged into the flags, or the neighbour + mask is not flipped with its epoch and lands on the wrong pixels. + """ + tile = _tile_with_neighbour() + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.neighbours) == len(stamp.flags) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + name = names[id(flag)] + npt.assert_array_equal(flag, 0) + npt.assert_array_equal(neighbour, _expected_neighbours(tile, name)) + assert not np.array_equal(stamp.neighbours[0], stamp.neighbours[1]) + + +def test_a_flagged_column_near_the_centre_is_still_vetoed(): + """With a neighbour footprint present, an epoch with a genuinely flagged + column 3 px from the centre is still dropped, and only that epoch. + + Failure mode: handling the markers apart also exempts real defects from + the central veto. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + column = clean.copy() + column[:, _CENTRE - 3] = 1 + epochs = { + "2100001-10": (clean, ones), + "2100002-20": (column, ones), + "2100003-11": (clean.copy(), ones), + } + stamp, _, _ = _marker_stamp(_tile_with_neighbour(), epochs) + names = {id(v[0]): name for name, v in epochs.items()} + assert sorted(names[id(f)] for f in stamp.flags) == [ + "2100001-10", "2100003-11" + ] + assert stamp.epoch_cuts["central_veto"] == 1 + assert stamp.epoch_cuts["masked_fraction"] == 0 + + +def _stamp_epoch_weights(stamp, i, seed, **kwargs): + return prepare_ngmix_weights( + 1.0e3 + stamp.gals[i], stamp.weights[i], stamp.flags[i], + np.random.RandomState(seed), bkg_rms=stamp.bkg_rms[i], + neighbour=stamp.neighbours[i], **kwargs, + ) + + +def test_noisefill_fills_exactly_the_marked_pixels(): + """Under noisefill, a defect-free epoch has zero weight and noise + exactly on the marked neighbour pixels; every other pixel keeps its raw + value and its weight. + + Failure mode: the neighbour markers are dropped with the flags, so + noisefill no longer removes neighbour light (noisefill-fills-markers). + """ + tile = _tile_with_neighbour() + stamp, epochs, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling="noisefill", + ) + neighbour = stamp.neighbours[i] + assert neighbour.any() + npt.assert_array_equal(gal_out != gal, neighbour) + npt.assert_array_equal(w_out == 0.0, neighbour) + assert np.all(np.abs(gal_out[neighbour]) < 10.0) + + +def test_uberseg_leaves_the_marked_pixels_raw(): + """Under uberseg, the markers mask nothing: with a seg map holding only + the central object, a defect-free epoch keeps every pixel raw and + weighted, marked or not. + + Failure mode: the markers reach the defect set or the fill, so uberseg + noise-fills the neighbour's light instead of leaving it to the seg-based + weight (uberseg-ignores-markers). + """ + tile = _tile_with_neighbour() + stamp, _, _ = _marker_stamp(tile) + seg = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + seg[_CENTRE - 1:_CENTRE + 2, _CENTRE - 1:_CENTRE + 2] = 1 + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling="uberseg", seg=seg, + object_number=1, + ) + npt.assert_array_equal(gal_out, gal) + assert np.all(w_out > 0.0) + + +def _develop_noisefill(gal, weight, flag, rng, bkg_rms=None): + """prepare_ngmix_weights under BLEND_HANDLING = noisefill on develop + (240b37e4), where the markers arrive as flag 2**10: the reference the + noisefill output is pinned to.""" + mask = np.copy(weight) != 0 + mask[flag != 0] = False + if bkg_rms is None: + sig_noise = sigma_mad(gal) + weight_map = mask.astype(float) / sig_noise ** 2 + else: + valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) + mask &= valid_rms + weight_map = np.zeros_like(gal, dtype=float) + weight_map[mask] = 1.0 / bkg_rms[mask] ** 2 + sig_noise = np.where(valid_rms, bkg_rms, np.median(bkg_rms[mask])) + noise_img = rng.standard_normal(gal.shape) * sig_noise + noise_img_gal = rng.standard_normal(gal.shape) * sig_noise + gal_masked = np.copy(gal) + gal_masked[~mask] = noise_img_gal[~mask] + return gal_masked, weight_map, noise_img + + +@given( + seed=st.integers(0, 2**31 - 1), + with_rms=st.booleans(), + rms_scale=st.floats(0.5, 2.0), + with_defects=st.booleans(), +) +def test_noisefill_matches_develop_on_marked_neighbours( + seed, with_rms, rms_scale, with_defects, +): + """For an epoch with a marked neighbour, noisefill returns the image, + weight and noise image develop returned, bit for bit: with no defect, + and with a flagged column and a dead pixel under the default noise fill. + + Failure mode: carrying the markers apart from the flags changes what + noisefill does to neighbour pixels, their weights, the noise level or + the RNG stream (noisefill-matches-develop). + """ + rng = np.random.default_rng(seed) + n = 31 + gal = rng.normal(0.0, 1.0, (n, n)) + gal[n // 2 - 2:n // 2 + 3, n // 2 - 2:n // 2 + 3] += 50.0 + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + neighbour = np.zeros((n, n), dtype=bool) + neighbour[n // 2 - 3:n // 2 + 4, n // 2 + 3:n // 2 + 9] = True + gal[neighbour] += 30.0 + if with_defects: + flag[:, 2] = 1 + weight[n - 3, n // 2] = 0.0 + bkg_rms = ( + rms_scale * (1.0 + 0.1 * rng.random((n, n))) if with_rms else None + ) + + new = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(seed), bkg_rms=bkg_rms, + blend_handling="noisefill", neighbour=neighbour, + ) + old = _develop_noisefill( + gal, weight, np.where(neighbour, 2**10, flag), + np.random.RandomState(seed), bkg_rms=bkg_rms, + ) + for a, b in zip(new, old): + npt.assert_array_equal(a, b) + + +def test_do_ngmix_metacal_threads_each_epochs_neighbour_mask(monkeypatch): + """Each epoch's neighbour mask reaches make_ngmix_observation. + + Failure mode: the mask is built but never used, so noisefill silently + stops filling neighbours. + """ + tile = _tile_with_neighbour() + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + seen = [] + + class _Stop(Exception): + pass + + def fake_observation(*args, **kwargs): + seen.append(kwargs["neighbour"]) + if len(seen) == len(stamp.gals): + raise _Stop + return None + + monkeypatch.setattr( + ngmix_module, "make_ngmix_observation", fake_observation, + ) + monkeypatch.setattr(ngmix_module, "ObsList", list) + with pytest.raises(_Stop): + ngmix_module.do_ngmix_metacal( + stamp, None, 1.0, np.random.RandomState(0), + ) + for got, want in zip(seen, stamp.neighbours): + assert got is want From 2517887cd49194969034df6948ba5bde44ba0a17 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:22:58 +0200 Subject: [PATCH 09/14] fix(ngmix): keep SExtractor's neighbour markers out of the defect set prepare_postage_stamps wrote the tile VIGNET's -1e30 markers (other detections' footprints) into each epoch's flag stamp as 2**10, so defect_mask counted them: they fed the masked-fraction cut and the central veto. Every epoch shares the tile VIGNET, so a neighbour within EPOCH_CENTRAL_DEFECT_RADIUS dropped every epoch of the object. The markers are now their own per-epoch mask, stamp.neighbours (MegaCam-flipped with the epoch), threaded through do_ngmix_metacal and make_ngmix_observation to prepare_ngmix_weights. The epoch cuts count genuine defects only. Under noisefill the marked pixels get weight 0 and the defects' noise realisation, the output noisefill gave when they arrived as flags; under uberseg they are ignored and the seg-based weight handles neighbours. Defects are treated identically under both. Docstrings, @sc contracts, runner comments and the astra.yaml decisions say so. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 48 +++--- src/shapepipe/modules/ngmix_package/ngmix.py | 150 +++++++++++++------ src/shapepipe/modules/ngmix_runner.py | 15 +- tests/module/test_ngmix_uberseg.py | 3 +- 4 files changed, 143 insertions(+), 73 deletions(-) diff --git a/astra.yaml b/astra.yaml index 978fe94a7..493922a10 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1329,10 +1329,11 @@ analyses: weighted pixels within about a PSF width. DES's ngmixer fills for that reason: "it may be important for codes that take moments or use FFTs". Under every BLEND_HANDLING and the default DEFECT_FILL = - noise, defect pixels are replaced with independent noise at the per-pixel background RMS when supplied, or - the stamp's robust noise scale otherwise; `BLEND_HANDLING = uberseg` - additionally zeroes the weight of neighbour-side pixels and leaves - their image values untouched. The committed fill + noise, defect pixels are replaced with independent noise at the + per-pixel background RMS when supplied, or the stamp's robust noise + scale otherwise. SExtractor's -1e30 neighbour markers in the tile + VIGNET are not defects; blend_handling decides their treatment. The + committed fill uses the unsymmetrized defect set: DES symmetrized its masks, but four-fold symmetrization quadruples m and still leaves an additive c1 (symmetrized_4fold_noise). The cost of not symmetrizing, a hole in @@ -1391,10 +1392,14 @@ analyses: label: Neighbour treatment before metacal rationale: >- Covers only pixels shared with a neighbour; defect_fill does not - depend on it: defects are noise-filled under every BLEND_HANDLING. - noisefill (the default; the committed config - sets no key) leaves neighbours fully weighted and untouched. - uberseg zeroes the weight of pixels nearer a neighbour's coadd + depend on it: defects are filled the same way under every + BLEND_HANDLING, and the epoch cuts never count neighbour pixels. + noisefill (the default; the committed config sets no key) gives + weight 0 to the pixels SExtractor marks -1e30 in the tile VIGNET + (other detections' footprints; about 93% of neighbour-footprint + pixels on a SExtractor-mode sims tile) and replaces them with noise; unmarked + neighbour pixels keep their weight and light. uberseg ignores the + markers and zeroes the weight of pixels nearer a neighbour's coadd segmentation footprint than the target's (DILATE_NEIGHBOUR, default 1) and leaves the image untouched, as official uberseg does: esheldon/meds get_uberseg returns a weight map (a @@ -1418,12 +1423,13 @@ analyses: default: noisefill options: noisefill: - label: No neighbour treatment (BLEND_HANDLING = noisefill) + label: Noise-fill SExtractor-marked neighbour pixels (BLEND_HANDLING = noisefill) description: >- - Neighbour pixels keep their full weight and image values, so the - fit sees all neighbour light, which biases shapes toward - neighbours (Jarvis et al. 2016); this masks less than even their - plain segmentation map. A candidate cause of the FLAGS=2 B-modes + Pixels marked -1e30 in the tile VIGNET get weight 0 and noise. + The fill stops at the marked footprint, so unmarked neighbour + pixels and the neighbour's wings keep their weight and light, + which biases shapes toward neighbours (Jarvis et al. 2016). A + candidate cause of the FLAGS=2 B-modes investigated in #814. insights: [mask_uberseg_neighbour_bias] uberseg: @@ -1456,9 +1462,10 @@ analyses: rationale: >- An epoch is dropped when a defect pixel lies strictly closer to the stamp centre than its fill's radius, beside the masked-fraction cut - in the epoch loop; EPOCH_CENTRAL_DEFECT_RADIUS = 0 disables it. The - veto reads only the - defect mask, so it selects on nothing shear-responsive; for the same + in the epoch loop; EPOCH_CENTRAL_DEFECT_RADIUS = 0 disables it. + SExtractor's -1e30 neighbour markers are not defects: all epochs + share the tile VIGNET, so a neighbour inside the radius would drop + every epoch. The veto reads only the defect mask, so it selects on nothing shear-responsive; for the same reason the radii are fixed rather than scaled by galaxy size: 10 px for noise-filled pixels (EPOCH_CENTRAL_DEFECT_RADIUS) and 7 px for interpolated ones (EPOCH_INTERPOLATED_DEFECT_RADIUS, @@ -1492,10 +1499,11 @@ analyses: rationale: >- An epoch whose stamp has more than EPOCH_MASKED_FRACTION_CUT (default 1/3) of its pixels in the raw, unsymmetrized defect set is - dropped from the multi-epoch fit: flagged pixels (any nonzero flag - bit, including the tile-coverage bit 2**10 set where the tile vignet - is off-image), zero-weight pixels and invalid-RMS pixels, the set - defect_fill noise-fills. An object with no surviving epoch has no + dropped from the multi-epoch fit: flagged pixels (any nonzero + exposure flag bit), zero-weight pixels and invalid-RMS pixels, the + set defect_fill fills. SExtractor's -1e30 neighbour markers in the + tile VIGNET are not counted: every epoch shares the tile VIGNET, so + a large neighbour would drop them all. An object with no surviving epoch has no shape. Before the cut, an epoch is dropped silently if its galaxy stamp is all zeros or its background-subtracted sigma_mad is not positive. DES was stricter: diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 5bd942913..1e9d10797 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -603,6 +603,11 @@ def __init__( self.psfs = [] self.weights = [] self.flags = [] + # Neighbour masks, one per epoch: the pixels SExtractor marks -1e30 in + # the tile VIGNET (other detections' footprints), MegaCam-flipped to + # the epoch. noisefill zero-weights and noise-fills them; uberseg and + # the epoch cuts do not read them (see prepare_ngmix_weights). + self.neighbours = [] self.bkg_rms = [] # Segmentation stamps, one per epoch, used only by the "uberseg" blend # handling; empty under the default "noisefill". All epochs carry the @@ -754,10 +759,12 @@ class Ngmix(object): adaptive-moment centroid measured from the stamp pixels. See :func:`make_ngmix_observation`. blend_handling : {"noisefill", "uberseg"}, optional - Neighbour treatment. ``"noisefill"`` (default) leaves neighbour pixels - untouched; ``"uberseg"`` zeroes the weight of neighbour-side pixels - from the coadd segmentation map and requires ``seg_cat_path``. Defect - pixels are filled under both (see :func:`prepare_ngmix_weights`). + Neighbour treatment. ``"noisefill"`` (default) zero-weights and + noise-fills the pixels SExtractor marks -1e30 in the tile VIGNET + (other detections' footprints); ``"uberseg"`` ignores those markers, + zeroes the weight of neighbour-side pixels from the coadd + segmentation map and requires ``seg_cat_path``. Defect pixels are + filled under both (see :func:`prepare_ngmix_weights`). seg_cat_path : str, optional Path to the coadd-frame segmentation VIGNET catalogue (see :class:`Tile_cat`). Required when ``blend_handling="uberseg"``. @@ -765,11 +772,11 @@ class Ngmix(object): Neighbour-mask dilation iterations for ``"uberseg"`` (see :func:`uberseg_weight`); the default is ``1``. epoch_central_defect_radius : float, optional - Drop an epoch when a masked pixel lies closer than this many pixels + Drop an epoch when a defect pixel lies closer than this many pixels to the stamp centre (see :func:`has_central_defect`); the default is ``EPOCH_CENTRAL_DEFECT_RADIUS``. epoch_masked_fraction_cut : float, optional - Drop an epoch when more than this fraction of its stamp is masked + Drop an epoch when more than this fraction of its stamp is defects (see :func:`prepare_postage_stamps`); the default is ``EPOCH_MASKED_FRACTION_CUT``. defect_fill : {"noise", "interpolate"}, optional @@ -1511,10 +1518,19 @@ def prepare_postage_stamps( An epoch is dropped when more than ``epoch_masked_fraction_cut`` of its stamp lies in :func:`defect_mask`: flagged, zero-weight and invalid-RMS pixels, the set that :func:`prepare_ngmix_weights` zero-weights and - fills. Counting flags alone would keep epochs whose filled area exceeds - the cut. The zero-weight orbit of interpolated pixels keeps its light and - is not counted. The default is 1/3; 10%, the DES Y3 and Y6 value, is - the alternative to test. + fills under every ``blend_handling``. Counting flags alone would keep + epochs whose filled area exceeds the cut. The zero-weight orbit of + interpolated pixels keeps its light and is not counted. The default is + 1/3; 10%, the DES Y3 and Y6 value, is the alternative to test. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.blend_handling] neighbour-markers-are-not-defects + SExtractor writes -1e30 into the tile VIGNET on the footprints of other + detections (and beyond the tile's edge). These markers form the epoch's neighbour mask + (``stamp.neighbours``, MegaCam-flipped like the epoch), kept apart from + its flag stamp, so neither the masked-fraction cut nor the central + veto counts them. Every epoch shares the tile VIGNET: counting a + neighbour within the veto radius as a defect would drop every epoch of + the object. Parameters ---------- @@ -1618,8 +1634,12 @@ def prepare_postage_stamps( tile_seg = Ngmix.MegaCamFlip(tile_seg, int(ccd_n)) flag_vign = flag_obj[expccd_name]['VIGNET'] - if tile_vign is not None: - flag_vign[np.where(tile_vign == -1e30)] = 2**10 + # Neighbour markers (neighbour-markers-are-not-defects). + neighbour = ( + tile_vign == -1e30 + if tile_vign is not None + else np.zeros(np.shape(gal_vign), dtype=bool) + ) weight_vign = weight_obj[expccd_name]['VIGNET'] bkg_rms_vign = ( bkg_rms_obj[expccd_name]['VIGNET'] @@ -1673,6 +1693,7 @@ def prepare_postage_stamps( stamp.psfs.append(psf_obj[expccd_name]['VIGNET']) stamp.weights.append(weight_vign_scaled) stamp.flags.append(flag_vign) + stamp.neighbours.append(neighbour) stamp.bkg_rms.append(bkg_rms_vign_scaled) if tile_seg is not None: stamp.segs.append(tile_seg) @@ -1981,11 +2002,13 @@ def defect_mask(weight, flag, bkg_rms=None): """Defect pixels of one epoch stamp. @sc [decision:shape_measurement.defect_fill,label:physics] defect-set-unsymmetrized - A defect is a pixel with zero exposure weight, a nonzero flag, or (when a - background RMS map is given) a non-finite or non-positive RMS. This one - set is zero-weighted and noise-filled by :func:`prepare_ngmix_weights` - and counted by the epoch cuts in :func:`prepare_postage_stamps`. It is - not ORed with its rotations. For defects the central-defect veto keeps + A defect is a pixel with zero exposure weight, a nonzero exposure flag, + or (when a background RMS map is given) a non-finite or non-positive RMS. + This one set is zero-weighted and filled by :func:`prepare_ngmix_weights` + under every ``blend_handling`` and counted by the epoch cuts in + :func:`prepare_postage_stamps`. SExtractor's neighbour markers are not + in it (neighbour-markers-are-not-defects). It is not ORed with its + rotations. For defects the central-defect veto keeps (:func:`has_central_defect`), the unsymmetrized fill leaves |c| <= 3e-4 per affected epoch on galaxies with half-light radius 0.3" and 0.5" through a 0.7" PSF, round or elliptical. Symmetrizing would @@ -2001,7 +2024,7 @@ def defect_mask(weight, flag, bkg_rms=None): weight : numpy.ndarray Exposure weight stamp. flag : numpy.ndarray - Flag stamp. + Exposure flag stamp. bkg_rms : numpy.ndarray, optional Background RMS stamp. @@ -2130,7 +2153,8 @@ def fill_defects(image, defect, noise): image : numpy.ndarray Stamp image. defect : numpy.ndarray of bool - Pixels to replace (:func:`defect_mask`). + Pixels to replace: the defects (:func:`defect_mask`), and under + noisefill the marked neighbour pixels. noise : numpy.ndarray Noise realisation on the stamp grid. @@ -2145,7 +2169,7 @@ def fill_defects(image, defect, noise): def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, blend_handling="noisefill", seg=None, object_number=None, - dilate_neighbour=0, defect_fill="noise", + dilate_neighbour=0, defect_fill="noise", neighbour=None, ): """Build one epoch's image, weight map and noise image for ngmix. @@ -2154,17 +2178,29 @@ def prepare_ngmix_weights( by an independent noise realisation at their background RMS (:func:`fill_defects`), or under ``defect_fill="interpolate"``, for short defect runs, by an interpolant of the clean pixels around them. - ``blend_handling`` decides only the neighbour side. + ``blend_handling`` decides only the neighbour pixels. - @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-neighbours-raw + @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-whatever-the-blend-handling Every pixel of ``defect_mask`` is zero-weighted and filled whatever - ``blend_handling`` is, and the filled set equals the defect set. - ``blend_handling`` acts on the neighbour side only: uberseg zeroes - the weight of pixels nearer a neighbour's footprint and leaves their - image values raw, because that light is real sky that metacal shears - along with the target; filling it with noise would cut the target along - an unsheared edge. Neighbour pixels are never filled. The noise image - covers the whole stamp and is independent of both. + ``blend_handling`` is, by the same operator. ``blend_handling`` acts on + the neighbour pixels only, and the noise image covers the whole stamp. + + @sc [decision:shape_measurement.blend_handling] noisefill-fills-markers + Under ``"noisefill"`` the pixels of ``neighbour``, SExtractor's -1e30 + neighbour markers, get weight 0 and are replaced by the same noise + realisation as the noise-filled defects, so no marked neighbour light + reaches metacal. Under the default noise fill, the image, weight map and + noise image are those the marked pixels would get as flagged defects; + only the epoch cuts treat them differently, by not counting them + (neighbour-markers-are-not-defects). Unmarked neighbour light stays + raw and weighted. + + @sc [decision:shape_measurement.blend_handling] uberseg-ignores-markers + Under ``"uberseg"`` the markers are ignored: :func:`uberseg_weight` + zeroes the weight of pixels nearer a neighbour's segmentation footprint + than the target's and leaves their image values raw, because that light + is real sky that metacal shears along with the target; filling it with + noise would cut the target along an unsheared edge. @sc [decision:shape_measurement.defect_fill,label:physics] interpolated-defect-weight-orbit Under ``defect_fill="interpolate"`` the short defect runs @@ -2189,7 +2225,7 @@ def prepare_ngmix_weights( weight : numpy.ndarray Exposure weight stamp; zero marks a defect. flag : numpy.ndarray - Flag stamp; nonzero marks a defect. + Exposure flag stamp; nonzero marks a defect. rng : numpy.random.RandomState Random state for the noise realisations (seeded per object; see :func:`position_seed`). @@ -2199,11 +2235,11 @@ def prepare_ngmix_weights( non-positive values mark defects. Otherwise every clean pixel gets ``1 / sigma_mad(gal)**2``. blend_handling : {"noisefill", "uberseg"}, optional - Neighbour treatment. ``"noisefill"`` (default) leaves neighbour - pixels weighted and untouched. ``"uberseg"`` zeroes the weight of - every pixel closer to a neighbour's segmentation footprint than to the - central object's and keeps its raw image value (see - :func:`uberseg_weight`). + Neighbour treatment. ``"noisefill"`` (default) zero-weights and + noise-fills the ``neighbour`` pixels. ``"uberseg"`` ignores + ``neighbour``, zeroes the weight of every pixel closer to a + neighbour's segmentation footprint than to the central object's and + keeps its raw image value (see :func:`uberseg_weight`). seg : numpy.ndarray, optional Segmentation map on the stamp grid (object NUMBERs). Required for ``blend_handling="uberseg"``; ignored otherwise. @@ -2216,11 +2252,15 @@ def prepare_ngmix_weights( defect_fill : {"noise", "interpolate"}, optional ``"noise"`` (default) noise-fills every defect; ``"interpolate"`` interpolates the short defect runs and noise-fills the rest. + neighbour : numpy.ndarray of bool, optional + The epoch's neighbour mask (``Postage_stamp.neighbours``); read only + under ``blend_handling="noisefill"``. ``None`` marks no pixel. Returns ------- numpy.ndarray - Galaxy image with defect pixels filled. + Galaxy image with defect pixels, and under noisefill marked + neighbour pixels, filled. numpy.ndarray Inverse-variance weight map for ngmix. numpy.ndarray @@ -2253,7 +2293,13 @@ def prepare_ngmix_weights( ) defect = defect_mask(weight, flag, bkg_rms) - clean = ~defect + # Marked neighbour pixels that noisefill removes (noisefill-fills-markers). + removed_neighbour = ( + np.asarray(neighbour, dtype=bool) + if blend_handling == "noisefill" and neighbour is not None + else np.zeros_like(defect) + ) + clean = ~(defect | removed_neighbour) interpolated = ( interpolable_defects(defect) if defect_fill == "interpolate" @@ -2262,7 +2308,7 @@ def prepare_ngmix_weights( # The quarter-turn orbit of the interpolated pixels carries no weight # (interpolated-defect-weight-orbit). weighted = ( - ~(defect | fourfold(interpolated)) if interpolated.any() else clean + clean & ~fourfold(interpolated) if interpolated.any() else clean ) if bkg_rms is None: @@ -2299,12 +2345,17 @@ def prepare_ngmix_weights( noise_img = rng.standard_normal(gal.shape) * sig_noise noise_img_gal = rng.standard_normal(gal.shape) * sig_noise - gal_filled = fill_defects(gal, defect, noise_img_gal) + gal_filled = fill_defects(gal, ~clean, noise_img_gal) if interpolated.any(): # One operator for the image and the noise image; a pixel whose - # support is degenerate keeps its noise fill. + # support is degenerate, or that noisefill removes as a neighbour, + # keeps its noise fill. filled = interpolate_defects([gal, noise_img], defect, interpolated) - done = interpolated & np.all(np.isfinite(filled), axis=0) + done = ( + interpolated + & ~removed_neighbour + & np.all(np.isfinite(filled), axis=0) + ) gal_filled[done] = filled[0][done] noise_img = np.where(done, filled[1], noise_img) @@ -2320,7 +2371,7 @@ def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, blend_handling="noisefill", seg=None, object_number=None, - dilate_neighbour=0, defect_fill="noise", + dilate_neighbour=0, defect_fill="noise", neighbour=None, ): """Build an ngmix Observation for a single galaxy epoch. @@ -2364,7 +2415,8 @@ def make_ngmix_observation( ``"hsm"``). blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to :func:`prepare_ngmix_weights`; - the default ``"noisefill"`` leaves neighbour pixels untouched. + the default ``"noisefill"`` zero-weights and noise-fills the + ``neighbour`` pixels. seg : numpy.ndarray, optional Segmentation map on the stamp grid. Required for ``blend_handling="uberseg"`` (ignored otherwise). @@ -2377,6 +2429,8 @@ def make_ngmix_observation( defect_fill : {"noise", "interpolate"}, optional Defect fill passed through to :func:`prepare_ngmix_weights`; the default is ``"noise"``. + neighbour : numpy.ndarray of bool, optional + Neighbour mask passed through to :func:`prepare_ngmix_weights`. Returns ------- @@ -2406,6 +2460,7 @@ def make_ngmix_observation( gal, weight, flag, rng, bkg_rms=bkg_rms, blend_handling=blend_handling, seg=seg, object_number=object_number, dilate_neighbour=dilate_neighbour, defect_fill=defect_fill, + neighbour=neighbour, ) if centroid_source == "hsm": @@ -2661,9 +2716,9 @@ def do_ngmix_metacal( stamp pixels — see that function. blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to - :func:`make_ngmix_observation`; the default ``"noisefill"`` leaves - neighbour pixels untouched. ``"uberseg"`` consumes ``stamp.segs`` and - ``object_number``. + :func:`make_ngmix_observation`; the default ``"noisefill"`` + zero-weights and noise-fills the pixels of ``stamp.neighbours``. + ``"uberseg"`` consumes ``stamp.segs`` and ``object_number``. object_number : int, optional Central object's segmentation label — its SExtractor ``NUMBER`` (``obj_id``), authoritative because seg labels are the NUMBERs of the @@ -2720,6 +2775,9 @@ def do_ngmix_metacal( object_number=object_number, dilate_neighbour=dilate_neighbour, defect_fill=defect_fill, + neighbour=( + stamp.neighbours[n_e] if n_e < len(stamp.neighbours) else None + ), ) gal_obs_list.append(gal_obs) diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 166e8cae5..dc27b8afc 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -144,10 +144,11 @@ def ngmix_runner( else: centroid_source = "wcs" - # Neighbour treatment: "noisefill" (default) leaves neighbour pixels - # weighted and untouched; "uberseg" zeroes the weight of every pixel - # closer to a neighbour than to the central object, from the segmentation - # map, and leaves its image raw. Defect pixels (flagged, zero-weight or + # Neighbour treatment: "noisefill" (default) zero-weights and noise-fills + # the pixels SExtractor marks -1e30 in the tile VIGNET (other detections' + # footprints); "uberseg" ignores those markers, zeroes the weight of every + # pixel closer to a neighbour than to the central object, from the + # segmentation map, and leaves its image raw. Defect pixels (flagged, zero-weight or # invalid-RMS) are zero-weighted and filled under both (DEFECT_FILL # below); see prepare_ngmix_weights. if config.has_option(module_config_sec, "BLEND_HANDLING"): @@ -164,7 +165,8 @@ def ngmix_runner( dilate_neighbour = 1 # EPOCH_CENTRAL_DEFECT_RADIUS (optional, pixels): drop an epoch when a - # masked pixel lies closer than this to the stamp centre; 0 disables. + # defect pixel (flagged, zero-weight or invalid-RMS; not a neighbour + # marker) lies closer than this to the stamp centre; 0 disables. if config.has_option(module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS"): epoch_central_defect_radius = config.getfloat( module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS" @@ -173,7 +175,8 @@ def ngmix_runner( epoch_central_defect_radius = EPOCH_CENTRAL_DEFECT_RADIUS # EPOCH_MASKED_FRACTION_CUT (optional): drop an epoch when more than this - # fraction of its stamp is masked (flagged, zero-weight or invalid-RMS). + # fraction of its stamp is defects (flagged, zero-weight or invalid-RMS; + # neighbour markers are not counted). if config.has_option(module_config_sec, "EPOCH_MASKED_FRACTION_CUT"): epoch_masked_fraction_cut = config.getfloat( module_config_sec, "EPOCH_MASKED_FRACTION_CUT" diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 92616fde0..5c20789f6 100644 --- a/tests/module/test_ngmix_uberseg.py +++ b/tests/module/test_ngmix_uberseg.py @@ -263,7 +263,8 @@ def test_uberseg_fills_defects_and_leaves_neighbour_pixels_raw(): Failure modes: the defect fill is skipped under uberseg, so raw bad pixels reach metacal; or the neighbour side is noise-filled - (defects-filled-neighbours-raw). + (defects-filled-whatever-the-blend-handling, + uberseg-ignores-markers). """ npix = 41 gal, flag, weight = _gal_flag_weight(npix=npix) From 7faa9e78dcd7118a1308563a7c4c322741b08fc5 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:45:43 +0200 Subject: [PATCH 10/14] fix(ngmix): off-tile pixels are defects, not neighbours SExtractor writes -1e30 into the tile VIGNET beyond the tile's edge as well as on other detections' footprints. Beyond the edge the epoch holds the object's own light, cut off, so treating those pixels as neighbours exempted edge objects from both epoch cuts and measured them with a noise-filled band near their centre. split_tile_markers classifies the markers geometrically on the flipped tile VIGNET: the union of stamp rows and columns that are entirely -1e30 (the off-image part of a rectangle clip) is off-tile, the rest is the neighbour mask. Off-tile pixels are flagged OFF_TILE_FLAG (2**10) in the epoch's flag stamp, so under both blend handlings they are zero-weighted, filled, and counted by the masked-fraction cut and the central veto. Docstrings, @sc contracts, runner comments and astra.yaml say so. Tests: an object 3 px from the tile edge loses every epoch; a band 12 px away passes the default cuts and is vetoed at radius 13; a corner's L-shaped off-tile region is flagged exactly and a border-touching neighbour stays a neighbour; off-tile pixels are filled under both blend handlings. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 33 +++-- src/shapepipe/modules/ngmix_package/ngmix.py | 77 +++++++++--- src/shapepipe/modules/ngmix_runner.py | 14 +-- tests/module/test_ngmix_defect_fill.py | 121 ++++++++++++++++++- 4 files changed, 206 insertions(+), 39 deletions(-) diff --git a/astra.yaml b/astra.yaml index 493922a10..fde5b8a70 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1321,8 +1321,8 @@ analyses: defect_fill: label: Image content of defect pixels before metacal rationale: >- - A defect pixel (nonzero instrument flag, zero exposure weight or - invalid background RMS) gets weight 0 in prepare_ngmix_weights. Its + A defect pixel (nonzero instrument flag, zero exposure weight, + invalid background RMS, or off-tile) gets weight 0 in prepare_ngmix_weights. Its image value still matters: ngmix's metacal deconvolves, shears and reconvolves an InterpolatedImage of the whole image and copies the weights through, so a zero-weight pixel's content spreads into the @@ -1331,9 +1331,13 @@ analyses: FFTs". Under every BLEND_HANDLING and the default DEFECT_FILL = noise, defect pixels are replaced with independent noise at the per-pixel background RMS when supplied, or the stamp's robust noise - scale otherwise. SExtractor's -1e30 neighbour markers in the tile - VIGNET are not defects; blend_handling decides their treatment. The - committed fill + scale otherwise. SExtractor writes -1e30 into the tile VIGNET on + other detections' footprints and beyond the tile's edge. The stamp + rows and columns that are entirely -1e30 (the off-image part of a + rectangle clip) are off-tile: the epoch holds the object's own light + there, so they are defects, flagged 2**10. The other markers are + neighbour pixels, not defects; blend_handling decides their + treatment. The committed fill uses the unsymmetrized defect set: DES symmetrized its masks, but four-fold symmetrization quadruples m and still leaves an additive c1 (symmetrized_4fold_noise). The cost of not symmetrizing, a hole in @@ -1396,8 +1400,10 @@ analyses: BLEND_HANDLING, and the epoch cuts never count neighbour pixels. noisefill (the default; the committed config sets no key) gives weight 0 to the pixels SExtractor marks -1e30 in the tile VIGNET - (other detections' footprints; about 93% of neighbour-footprint - pixels on a SExtractor-mode sims tile) and replaces them with noise; unmarked + on other detections' footprints (about 93% of neighbour-footprint + pixels on a SExtractor-mode sims tile), excluding the off-tile + rows and columns, which are defects, and replaces them with noise; + unmarked neighbour pixels keep their weight and light. uberseg ignores the markers and zeroes the weight of pixels nearer a neighbour's coadd segmentation footprint than the target's (DILATE_NEIGHBOUR, default 1) @@ -1465,7 +1471,8 @@ analyses: in the epoch loop; EPOCH_CENTRAL_DEFECT_RADIUS = 0 disables it. SExtractor's -1e30 neighbour markers are not defects: all epochs share the tile VIGNET, so a neighbour inside the radius would drop - every epoch. The veto reads only the defect mask, so it selects on nothing shear-responsive; for the same + every epoch. Off-tile pixels are defects, so an object near the + tile edge is vetoed. The veto reads only the defect mask, so it selects on nothing shear-responsive; for the same reason the radii are fixed rather than scaled by galaxy size: 10 px for noise-filled pixels (EPOCH_CENTRAL_DEFECT_RADIUS) and 7 px for interpolated ones (EPOCH_INTERPOLATED_DEFECT_RADIUS, @@ -1500,10 +1507,12 @@ analyses: An epoch whose stamp has more than EPOCH_MASKED_FRACTION_CUT (default 1/3) of its pixels in the raw, unsymmetrized defect set is dropped from the multi-epoch fit: flagged pixels (any nonzero - exposure flag bit), zero-weight pixels and invalid-RMS pixels, the - set defect_fill fills. SExtractor's -1e30 neighbour markers in the - tile VIGNET are not counted: every epoch shares the tile VIGNET, so - a large neighbour would drop them all. An object with no surviving epoch has no + exposure flag bit), zero-weight pixels, invalid-RMS pixels and + off-tile pixels (whole -1e30 rows and columns of the tile VIGNET, + flag 2**10), the set defect_fill fills. On a 51-px stamp an object + within about 8.5 px of the tile edge fails the 1/3 cut. The other + -1e30 markers, neighbour footprints, are not counted: every epoch + shares the tile VIGNET, so a large neighbour would drop them all. An object with no surviving epoch has no shape. Before the cut, an epoch is dropped silently if its galaxy stamp is all zeros or its background-subtracted sigma_mad is not positive. DES was stricter: diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 1e9d10797..9cd196743 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -30,6 +30,9 @@ ) from shapepipe.pipeline import file_io +# Flag bit set on an epoch's off-tile pixels (see :func:`split_tile_markers`). +OFF_TILE_FLAG = 2**10 + # Neighbour treatments selectable with the BLEND_HANDLING option. BLEND_HANDLINGS = ("noisefill", "uberseg") @@ -604,8 +607,9 @@ def __init__( self.weights = [] self.flags = [] # Neighbour masks, one per epoch: the pixels SExtractor marks -1e30 in - # the tile VIGNET (other detections' footprints), MegaCam-flipped to - # the epoch. noisefill zero-weights and noise-fills them; uberseg and + # the tile VIGNET on other detections' footprints (off-tile markers + # are flagged as defects instead; see split_tile_markers), + # MegaCam-flipped to the epoch. noisefill zero-weights and noise-fills them; uberseg and # the epoch cuts do not read them (see prepare_ngmix_weights). self.neighbours = [] self.bkg_rms = [] @@ -761,7 +765,7 @@ class Ngmix(object): blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment. ``"noisefill"`` (default) zero-weights and noise-fills the pixels SExtractor marks -1e30 in the tile VIGNET - (other detections' footprints); ``"uberseg"`` ignores those markers, + on other detections' footprints; ``"uberseg"`` ignores those markers, zeroes the weight of neighbour-side pixels from the coadd segmentation map and requires ``seg_cat_path``. Defect pixels are filled under both (see :func:`prepare_ngmix_weights`). @@ -1525,12 +1529,20 @@ def prepare_postage_stamps( @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.blend_handling] neighbour-markers-are-not-defects SExtractor writes -1e30 into the tile VIGNET on the footprints of other - detections (and beyond the tile's edge). These markers form the epoch's neighbour mask - (``stamp.neighbours``, MegaCam-flipped like the epoch), kept apart from - its flag stamp, so neither the masked-fraction cut nor the central - veto counts them. Every epoch shares the tile VIGNET: counting a - neighbour within the veto radius as a defect would drop every epoch of - the object. + detections and beyond the tile's edge (:func:`split_tile_markers`). The + footprint markers form the epoch's neighbour mask (``stamp.neighbours``, + MegaCam-flipped like the epoch), kept apart from its flag stamp, so + neither the masked-fraction cut nor the central veto counts them. Every + epoch shares the tile VIGNET: counting a neighbour within the veto + radius as a defect would drop every epoch of the object. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.defect_fill] off-tile-pixels-are-defects + Beyond the tile's edge the epoch holds real data, the object's own light + cut off by the tile. Those pixels are flagged ``OFF_TILE_FLAG`` (2**10) + in the epoch's flag stamp and so join its defect set under every + ``blend_handling``: zero weight, the defect fill, and both epoch cuts. + On a 51-px stamp an object within about 8.5 px of the tile edge fails + the 1/3 cut. Parameters ---------- @@ -1634,12 +1646,10 @@ def prepare_postage_stamps( tile_seg = Ngmix.MegaCamFlip(tile_seg, int(ccd_n)) flag_vign = flag_obj[expccd_name]['VIGNET'] - # Neighbour markers (neighbour-markers-are-not-defects). - neighbour = ( - tile_vign == -1e30 - if tile_vign is not None - else np.zeros(np.shape(gal_vign), dtype=bool) - ) + # Off-tile pixels are defects (off-tile-pixels-are-defects); the + # other -1e30 markers are neighbours (neighbour-markers-are-not-defects). + neighbour, off_tile = split_tile_markers(tile_vign, np.shape(gal_vign)) + flag_vign[off_tile] = OFF_TILE_FLAG weight_vign = weight_obj[expccd_name]['VIGNET'] bkg_rms_vign = ( bkg_rms_obj[expccd_name]['VIGNET'] @@ -1714,6 +1724,38 @@ def prepare_postage_stamps( return stamp +def split_tile_markers(tile_vign, shape): + """Split the tile VIGNET's -1e30 markers into neighbour and off-tile. + + @sc [decision:shape_measurement.blend_handling,decision:shape_measurement.defect_fill] off-tile-is-whole-marked-rows-and-columns + SExtractor writes -1e30 on the footprints of other detections and on + stamp pixels beyond the tile's edge. A stamp clipped by the tile's + rectangle loses whole rows and whole columns, so the off-tile pixels are + the union of the stamp rows and columns that are entirely -1e30. The + remaining markers are neighbour pixels; a footprint touching the stamp + border stays a neighbour unless it fills a whole row or column. + + Parameters + ---------- + tile_vign : numpy.ndarray or None + Tile VIGNET stamp, oriented like the epoch; ``None`` marks nothing. + shape : tuple of int + Stamp shape, used when ``tile_vign`` is ``None``. + + Returns + ------- + numpy.ndarray of bool + Neighbour pixels. + numpy.ndarray of bool + Off-tile pixels. + """ + if tile_vign is None: + return np.zeros(shape, dtype=bool), np.zeros(shape, dtype=bool) + marker = tile_vign == -1e30 + off_tile = marker.all(axis=1)[:, None] | marker.all(axis=0)[None, :] + return marker & ~off_tile, off_tile + + def background_subtract(gal,bkg): """background subtraction. @@ -2007,7 +2049,8 @@ def defect_mask(weight, flag, bkg_rms=None): This one set is zero-weighted and filled by :func:`prepare_ngmix_weights` under every ``blend_handling`` and counted by the epoch cuts in :func:`prepare_postage_stamps`. SExtractor's neighbour markers are not - in it (neighbour-markers-are-not-defects). It is not ORed with its + in it (neighbour-markers-are-not-defects); off-tile pixels are, as flag + ``OFF_TILE_FLAG`` (off-tile-pixels-are-defects). It is not ORed with its rotations. For defects the central-defect veto keeps (:func:`has_central_defect`), the unsymmetrized fill leaves |c| <= 3e-4 per affected epoch on galaxies with half-light radius 0.3" @@ -2196,7 +2239,7 @@ def prepare_ngmix_weights( raw and weighted. @sc [decision:shape_measurement.blend_handling] uberseg-ignores-markers - Under ``"uberseg"`` the markers are ignored: :func:`uberseg_weight` + Under ``"uberseg"`` the neighbour markers are ignored: :func:`uberseg_weight` zeroes the weight of pixels nearer a neighbour's segmentation footprint than the target's and leaves their image values raw, because that light is real sky that metacal shears along with the target; filling it with diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index dc27b8afc..e217995c2 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -148,9 +148,9 @@ def ngmix_runner( # the pixels SExtractor marks -1e30 in the tile VIGNET (other detections' # footprints); "uberseg" ignores those markers, zeroes the weight of every # pixel closer to a neighbour than to the central object, from the - # segmentation map, and leaves its image raw. Defect pixels (flagged, zero-weight or - # invalid-RMS) are zero-weighted and filled under both (DEFECT_FILL - # below); see prepare_ngmix_weights. + # segmentation map, and leaves its image raw. Defect pixels (flagged, + # zero-weight, invalid-RMS or off-tile) are zero-weighted and filled + # under both (DEFECT_FILL below); see prepare_ngmix_weights. if config.has_option(module_config_sec, "BLEND_HANDLING"): blend_handling = config.get(module_config_sec, "BLEND_HANDLING") else: @@ -165,8 +165,8 @@ def ngmix_runner( dilate_neighbour = 1 # EPOCH_CENTRAL_DEFECT_RADIUS (optional, pixels): drop an epoch when a - # defect pixel (flagged, zero-weight or invalid-RMS; not a neighbour - # marker) lies closer than this to the stamp centre; 0 disables. + # defect pixel (flagged, zero-weight, invalid-RMS or off-tile; not a + # neighbour marker) lies closer than this to the stamp centre; 0 disables. if config.has_option(module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS"): epoch_central_defect_radius = config.getfloat( module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS" @@ -175,8 +175,8 @@ def ngmix_runner( epoch_central_defect_radius = EPOCH_CENTRAL_DEFECT_RADIUS # EPOCH_MASKED_FRACTION_CUT (optional): drop an epoch when more than this - # fraction of its stamp is defects (flagged, zero-weight or invalid-RMS; - # neighbour markers are not counted). + # fraction of its stamp is defects (flagged, zero-weight, invalid-RMS or + # off-tile; neighbour markers are not counted). if config.has_option(module_config_sec, "EPOCH_MASKED_FRACTION_CUT"): epoch_masked_fraction_cut = config.getfloat( module_config_sec, "EPOCH_MASKED_FRACTION_CUT" diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index b9e3096c5..9c15d1807 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -786,7 +786,7 @@ def _expected_neighbours(tile, name): def test_neighbour_markers_near_the_centre_keep_every_epoch(): - """A neighbour footprint 3 px from the centre, and one covering 45% of + """A neighbour footprint 3 px from the centre, and one covering 41% of the stamp, drop no epoch: the masked-fraction cut and the central veto count no marker. @@ -796,8 +796,10 @@ def test_neighbour_markers_near_the_centre_keep_every_epoch(): (neighbour-markers-are-not-defects). """ small = _tile_with_neighbour() - large = _tile_with_neighbour(rows=(0, N_STAMP)) - large[:, _CENTRE + 3:] = _MARKER + # Large, but short of the stamp border: no row or column is entirely + # marked, so it is a neighbour, not off-tile. + large = _tile_with_neighbour() + large[1:-1, _CENTRE + 3:-1] = _MARKER assert (large == _MARKER).mean() > 1 / 3 for tile in (small, large): stamp, _, _ = _marker_stamp(tile) @@ -939,6 +941,11 @@ def test_noisefill_matches_develop_on_marked_neighbours( """For an epoch with a marked neighbour, noisefill returns the image, weight and noise image develop returned, bit for bit: with no defect, and with a flagged column and a dead pixel under the default noise fill. + The stamps carry no off-tile pixels, and the equivalence is claimed for + such stamps only. Off-tile pixels reach this function as flag 2**10 + defects (off-tile-pixels-are-defects), as every marker did on develop; + which epochs survive the cuts differs from develop wherever there are + neighbour markers. Failure mode: carrying the markers apart from the flags changes what noisefill does to neighbour pixels, their weights, the noise level or @@ -1002,3 +1009,111 @@ def fake_observation(*args, **kwargs): ) for got, want in zip(seen, stamp.neighbours): assert got is want + + +# --- Off-tile pixels are defects ------------------------------------------- +# +# SExtractor also writes -1e30 beyond the tile's edge, where the epoch holds +# the object's own light, cut off. Those pixels are the stamp rows and +# columns that are entirely -1e30 (the off-image part of a rectangle clip); +# they join the epoch's defect set as flag 2**10. The other markers are the +# neighbour mask. + +_OFF_TILE = 2**10 + + +def _off_tile_expected(tile, name): + flipped = Ngmix.MegaCamFlip(tile, int(name.split("-")[1])) == _MARKER + return flipped.all(axis=1)[:, None] | flipped.all(axis=0)[None, :] + + +def test_object_three_px_from_the_tile_edge_is_dropped(): + """With the tile edge 3 px from the object, every epoch is dropped: the + off-tile band counts toward the epoch cuts. + + Failure mode: off-tile pixels are treated as neighbour markers, so an + edge object is measured with a noise-filled band through its own light + (off-tile-pixels-are-defects). + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:, :_CENTRE - 2] = _MARKER + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == 0 + assert stamp.epoch_cuts["considered"] == len(_MARKER_EPOCH_NAMES) + assert ( + stamp.epoch_cuts["masked_fraction"] + stamp.epoch_cuts["central_veto"] + == len(_MARKER_EPOCH_NAMES) + ) + + +def test_the_central_veto_sees_off_tile_pixels(): + """An off-tile band 12 px from the object (27% of the stamp) passes the + default cuts and is vetoed at radius 13. + + Failure mode: the central veto does not read the off-tile set. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:, :_CENTRE - 11] = _MARKER + kept, _, _ = _marker_stamp(tile) + assert len(kept.gals) == len(_MARKER_EPOCH_NAMES) + vetoed, _, _ = _marker_stamp(tile, epoch_central_defect_radius=13) + assert len(vetoed.gals) == 0 + assert vetoed.epoch_cuts["central_veto"] == len(_MARKER_EPOCH_NAMES) + + +def test_corner_off_tile_region_and_border_neighbour_are_classified(): + """At a tile corner, the L-shaped off-tile region is flagged 2**10 + exactly, and a neighbour footprint touching the stamp border without + filling a row or column stays in the neighbour mask. + + Failure modes: off-tile pixels are classified by something other than + whole marked rows and columns (the L is missed or a border-touching + neighbour is swallowed); the classification ignores the MegaCam flip. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:5, :] = _MARKER + tile[:, -5:] = _MARKER + tile[40:, :6] = _MARKER # neighbour on the bottom-left border + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + name = names[id(flag)] + off_tile = _off_tile_expected(tile, name) + assert off_tile.sum() == 5 * N_STAMP * 2 - 25 + npt.assert_array_equal(flag == _OFF_TILE, off_tile) + npt.assert_array_equal( + neighbour, _expected_neighbours(tile, name) & ~off_tile + ) + assert neighbour.sum() == 11 * 6 + + +@pytest.mark.parametrize("blend_handling", ["noisefill", "uberseg"]) +def test_off_tile_pixels_are_zero_weighted_and_filled(blend_handling): + """Off-tile pixels are zero-weighted and noise-filled under either + BLEND_HANDLING; under uberseg the neighbour markers stay raw. + + Failure mode: under uberseg the off-tile band keeps its weight, or the + fill differs between blend handlings. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:5, :] = _MARKER + tile[20:24, 35:40] = _MARKER + stamp, _, _ = _marker_stamp(tile) + seg = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + seg[_CENTRE - 1:_CENTRE + 2, _CENTRE - 1:_CENTRE + 2] = 1 + kwargs = ( + dict(seg=seg, object_number=1) if blend_handling == "uberseg" else {} + ) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling=blend_handling, **kwargs, + ) + off_tile = stamp.flags[i] == _OFF_TILE + assert off_tile.sum() == 5 * N_STAMP + removed = off_tile | ( + stamp.neighbours[i] if blend_handling == "noisefill" else False + ) + npt.assert_array_equal(gal_out != gal, removed) + npt.assert_array_equal(w_out == 0.0, removed) From f0654053f5e300d9ea6293afb2b5022c0895b144 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:47:26 +0200 Subject: [PATCH 11/14] fix(ngmix): interpolate defects from the pixels the image keeps Under DEFECT_FILL = interpolate and BLEND_HANDLING = noisefill, the Clough-Tocher support included marked neighbour pixels that noisefill replaces with noise, so a defect beside a bright neighbour was filled with light the image no longer contains. The support now excludes ~clean (the defects and the removed neighbour pixels); the NaN fallback to the noise fill covers degenerate supports. Under uberseg the neighbour light stays in the image and the support is unchanged. interpolate_defects names its second argument `excluded`. Test: a flagged pixel next to a 1e4 marked neighbour interpolates to the kept-pixel interpolant (1.4, against 3102 with the neighbour in the support) under noisefill, and from the raw neighbour light under uberseg. Co-Authored-By: Claude Opus 5.5 --- .../ngmix_package/defect_interpolation.py | 24 ++++----- src/shapepipe/modules/ngmix_package/ngmix.py | 18 +++++-- tests/module/test_ngmix_defect_fill.py | 52 +++++++++++++++++++ 3 files changed, 78 insertions(+), 16 deletions(-) diff --git a/src/shapepipe/modules/ngmix_package/defect_interpolation.py b/src/shapepipe/modules/ngmix_package/defect_interpolation.py index 713e326f7..ae0d53c41 100644 --- a/src/shapepipe/modules/ngmix_package/defect_interpolation.py +++ b/src/shapepipe/modules/ngmix_package/defect_interpolation.py @@ -94,14 +94,14 @@ def fourfold(mask): return mask | np.rot90(mask) | np.rot90(mask, 2) | np.rot90(mask, 3) -def _interpolate_once(planes, defect, target): +def _interpolate_once(planes, excluded, target): """One Clough-Tocher interpolant of every plane at ``target``; NaN elsewhere and where the support cannot reach.""" out = np.full(planes.shape, np.nan) support = binary_dilation( target, structure=np.ones((3, 3), dtype=bool), iterations=SUPPORT_RADIUS, - ) & ~defect + ) & ~excluded points = np.argwhere(support).astype(float) if len(points) < 3: return out @@ -116,14 +116,14 @@ def _interpolate_once(planes, defect, target): return out -def interpolate_defects(planes, defect, target): +def interpolate_defects(planes, excluded, target): """Replace the ``target`` pixels of every plane by a Clough-Tocher - interpolant of the clean pixels around them. + interpolant of the kept pixels around them. @sc [decision:shape_measurement.defect_fill] shared-rotation-averaged-interpolant - The support is the clean pixels within ``SUPPORT_RADIUS`` (4 px) of the - target; no defect pixel enters it, so defect values are never read. For - each quarter turn of the stamp, one Delaunay triangulation of the support + The support is the pixels within ``SUPPORT_RADIUS`` (4 px) of the target + outside ``excluded``; no excluded pixel enters it, so their values are + never read. For each quarter turn of the stamp, one Delaunay triangulation of the support serves every plane, so the science image and the metacal noise image see the same linear operator and fixnoise mirrors the science image's interpolated noise. A regular grid's triangulation has degenerate @@ -136,9 +136,9 @@ def interpolate_defects(planes, defect, target): ---------- planes : array_like Stamp planes, shape ``(n, ny, nx)``. - defect : numpy.ndarray of bool - Every defect pixel, shape ``(ny, nx)``; none of them supports the - interpolant. + excluded : numpy.ndarray of bool + Pixels that never support the interpolant, shape ``(ny, nx)``: every + defect, and any pixel whose light the image does not keep. target : numpy.ndarray of bool Defect pixels to interpolate (:func:`interpolable_defects`). @@ -149,7 +149,7 @@ def interpolate_defects(planes, defect, target): target pixel whose support is degenerate in some orientation. """ planes = np.asarray(planes, dtype=float) - defect = np.asarray(defect, dtype=bool) + excluded = np.asarray(excluded, dtype=bool) target = np.asarray(target, dtype=bool) out = planes.copy() if not target.any(): @@ -157,7 +157,7 @@ def interpolate_defects(planes, defect, target): turns = [ np.rot90( _interpolate_once( - np.rot90(planes, k, axes=(1, 2)), np.rot90(defect, k), + np.rot90(planes, k, axes=(1, 2)), np.rot90(excluded, k), np.rot90(target, k), ), -k, axes=(1, 2), diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 9cd196743..bf9bf0b12 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -2228,6 +2228,15 @@ def prepare_ngmix_weights( ``blend_handling`` is, by the same operator. ``blend_handling`` acts on the neighbour pixels only, and the noise image covers the whole stamp. + @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] interpolation-support-is-the-kept-image + Under ``defect_fill="interpolate"`` the interpolant is supported on the + pixels whose light the image keeps: never a defect, and under + ``"noisefill"`` never a marked neighbour pixel, whose light is replaced + by noise; a Clough-Tocher interpolant next to a bright neighbour would + otherwise build the defect's value from light the image no longer + contains. Under ``"uberseg"`` the neighbour light stays in the image and + supports the interpolant. + @sc [decision:shape_measurement.blend_handling] noisefill-fills-markers Under ``"noisefill"`` the pixels of ``neighbour``, SExtractor's -1e30 neighbour markers, get weight 0 and are replaced by the same noise @@ -2390,10 +2399,11 @@ def prepare_ngmix_weights( noise_img_gal = rng.standard_normal(gal.shape) * sig_noise gal_filled = fill_defects(gal, ~clean, noise_img_gal) if interpolated.any(): - # One operator for the image and the noise image; a pixel whose - # support is degenerate, or that noisefill removes as a neighbour, - # keeps its noise fill. - filled = interpolate_defects([gal, noise_img], defect, interpolated) + # One operator for the image and the noise image, supported on the + # pixels the image keeps (interpolation-support-is-the-kept-image); a + # pixel whose support is degenerate, or that noisefill removes as a + # neighbour, keeps its noise fill. + filled = interpolate_defects([gal, noise_img], ~clean, interpolated) done = ( interpolated & ~removed_neighbour diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 9c15d1807..67aff0d20 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -1117,3 +1117,55 @@ def test_off_tile_pixels_are_zero_weighted_and_filled(blend_handling): ) npt.assert_array_equal(gal_out != gal, removed) npt.assert_array_equal(w_out == 0.0, removed) + + +# --- DEFECT_FILL = interpolate beside a removed neighbour ------------------- + +def _defect_beside_neighbour(): + """A single flagged pixel 6 px right of the centre, with a bright marked + neighbour footprint starting on the next column.""" + n = 31 + c = n // 2 + gal = np.random.default_rng(3).normal(0.0, 1.0, (n, n)) + flag = np.zeros((n, n), dtype=np.int32) + flag[c, c + 6] = 1 + neighbour = np.zeros((n, n), dtype=bool) + neighbour[c - 2:c + 3, c + 7:c + 10] = True + gal[neighbour] += 1.0e4 + seg = np.zeros((n, n), dtype=np.int32) + seg[c - 1:c + 2, c - 1:c + 2] = 1 + return gal, flag, neighbour, seg, (c, c + 6) + + +def test_noisefill_interpolation_does_not_read_removed_neighbour_light(): + """Under noisefill, the interpolant of a defect beside a marked neighbour + is built from the pixels the image keeps: the neighbour's light, which + noisefill removes, is not in its support. Under uberseg the neighbour + light is raw in the image and supports the interpolant. + + Failure mode: the support includes the removed neighbour pixels, so the + defect is filled with light the image no longer contains + (interpolation-support-is-the-kept-image). + """ + gal, flag, neighbour, seg, pix = _defect_beside_neighbour() + weight = np.ones_like(gal) + target = np.zeros_like(neighbour) + target[pix] = True + + out, _, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + blend_handling="noisefill", neighbour=neighbour, + defect_fill="interpolate", + ) + expected = interpolate_defects(gal[None], (flag != 0) | neighbour, target) + assert out[pix] == pytest.approx(expected[0][pix]) + assert abs(out[pix]) < 100.0 + + out, _, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + blend_handling="uberseg", seg=seg, object_number=1, + neighbour=neighbour, defect_fill="interpolate", + ) + expected = interpolate_defects(gal[None], flag != 0, target) + assert out[pix] == pytest.approx(expected[0][pix]) + assert out[pix] > 1000.0 From b27eca849426564bee89ce814bd788db728807e0 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:48:52 +0200 Subject: [PATCH 12/14] fix(ngmix): the defect fill keeps the stamp's dtype fill_defects used np.where, which promotes a float32 stamp to float64; develop's noisefill kept the galaxy stamp's dtype. The filled image is cast back to the input dtype, so noisefill stays bit-identical to develop on float32 stamps. The develop-equivalence test gains a float32 case and checks dtypes. Co-Authored-By: Claude Opus 5.5 --- src/shapepipe/modules/ngmix_package/ngmix.py | 5 +++-- tests/module/test_ngmix_defect_fill.py | 9 +++++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index bf9bf0b12..8f0ef66e1 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -2204,9 +2204,10 @@ def fill_defects(image, defect, noise): Returns ------- numpy.ndarray - ``image`` with ``defect`` pixels taken from ``noise``. + ``image`` with ``defect`` pixels taken from ``noise``, in the dtype + of ``image``. """ - return np.where(defect, noise, image) + return np.where(defect, noise, image).astype(image.dtype, copy=False) def prepare_ngmix_weights( diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 67aff0d20..17f04a2a8 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -934,13 +934,16 @@ def _develop_noisefill(gal, weight, flag, rng, bkg_rms=None): with_rms=st.booleans(), rms_scale=st.floats(0.5, 2.0), with_defects=st.booleans(), + dtype=st.sampled_from([np.float64, np.float32]), ) def test_noisefill_matches_develop_on_marked_neighbours( - seed, with_rms, rms_scale, with_defects, + seed, with_rms, rms_scale, with_defects, dtype, ): """For an epoch with a marked neighbour, noisefill returns the image, weight and noise image develop returned, bit for bit: with no defect, - and with a flagged column and a dead pixel under the default noise fill. + and with a flagged column and a dead pixel under the default noise fill, + for float64 and float32 stamps (the filled image keeps the stamp's + dtype). The stamps carry no off-tile pixels, and the equivalence is claimed for such stamps only. Off-tile pixels reach this function as flag 2**10 defects (off-tile-pixels-are-defects), as every marker did on develop; @@ -963,6 +966,7 @@ def test_noisefill_matches_develop_on_marked_neighbours( if with_defects: flag[:, 2] = 1 weight[n - 3, n // 2] = 0.0 + gal = gal.astype(dtype) bkg_rms = ( rms_scale * (1.0 + 0.1 * rng.random((n, n))) if with_rms else None ) @@ -976,6 +980,7 @@ def test_noisefill_matches_develop_on_marked_neighbours( np.random.RandomState(seed), bkg_rms=bkg_rms, ) for a, b in zip(new, old): + assert a.dtype == b.dtype npt.assert_array_equal(a, b) From a5bfef10ca0fb905d1aa5658295208d57f573ff9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 02:01:47 +0200 Subject: [PATCH 13/14] fix(ngmix): off-tile pixels are the marked rows and columns at the stamp border Every entirely -1e30 stamp row or column counted as off-tile. In an edge stamp, a neighbour footprint that covers the on-image part of a row (a wide DR6 footprint running from the tile edge) completes that row, so the neighbour's pixels became OFF_TILE_FLAG defects: counted by the epoch cuts, interpolated under DEFECT_FILL = interpolate, filled under uberseg. Only the runs of marked rows and columns that start at a stamp border are off-tile now. Tests build tile VIGNETs with the catalogue-mode converter: on the real 202.301 DR6 patch every stamp splits exactly into off-image pixels and seg-map neighbour pixels, and the edge case above keeps its footprint in the neighbour mask through prepare_postage_stamps. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BbT81hykcvtuefzB3FhYvY --- astra.yaml | 17 ++-- src/shapepipe/modules/ngmix_package/ngmix.py | 30 +++++-- tests/module/test_ngmix_defect_fill.py | 91 ++++++++++++++++++++ 3 files changed, 123 insertions(+), 15 deletions(-) diff --git a/astra.yaml b/astra.yaml index 8eea6bc4c..f83377c80 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1386,13 +1386,16 @@ analyses: FFTs". Under every BLEND_HANDLING and the default DEFECT_FILL = noise, defect pixels are replaced with independent noise at the per-pixel background RMS when supplied, or the stamp's robust noise - scale otherwise. SExtractor writes -1e30 into the tile VIGNET on - other detections' footprints and beyond the tile's edge. The stamp - rows and columns that are entirely -1e30 (the off-image part of a - rectangle clip) are off-tile: the epoch holds the object's own light - there, so they are defects, flagged 2**10. The other markers are - neighbour pixels, not defects; blend_handling decides their - treatment. The committed fill + scale otherwise. The tile VIGNET holds -1e30 on other detections' + footprints and beyond the tile's edge (written by SExtractor, or by + the unions_catalogue converter from the catalogue's segmentation + map; catalogue_neighbour_marking). The runs of entirely -1e30 stamp + rows and columns that start at a stamp border (the off-image part of + a rectangle clip) are off-tile: the epoch holds the object's own + light there, so they are defects, flagged 2**10. The other markers, + including a neighbour footprint that completes an interior row + beside an off-tile band, are neighbour pixels, not defects; + blend_handling decides their treatment. The committed fill uses the unsymmetrized defect set: DES symmetrized its masks, but four-fold symmetrization quadruples m and still leaves an additive c1 (symmetrized_4fold_noise). The cost of not symmetrizing, a hole in diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 2f172b4a9..4908fcd20 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -1754,13 +1754,17 @@ def prepare_postage_stamps( def split_tile_markers(tile_vign, shape): """Split the tile VIGNET's -1e30 markers into neighbour and off-tile. - @sc [decision:shape_measurement.blend_handling,decision:shape_measurement.defect_fill] off-tile-is-whole-marked-rows-and-columns - SExtractor writes -1e30 on the footprints of other detections and on - stamp pixels beyond the tile's edge. A stamp clipped by the tile's - rectangle loses whole rows and whole columns, so the off-tile pixels are - the union of the stamp rows and columns that are entirely -1e30. The - remaining markers are neighbour pixels; a footprint touching the stamp - border stays a neighbour unless it fills a whole row or column. + @sc [decision:shape_measurement.blend_handling,decision:shape_measurement.defect_fill] off-tile-is-marked-border-rows-and-columns + The tile VIGNET holds -1e30 on the footprints of other detections and on + stamp pixels beyond the tile's edge; SExtractor writes it, and in + catalogue mode the converter paints it from the catalogue's segmentation + map (``read_ext_sexcat._extract_vignets``). A stamp clipped by the tile's + rectangle loses whole rows and whole columns from its border, so the + off-tile pixels are the union of the runs of entirely -1e30 rows and + columns that start at a stamp border. The remaining markers are + neighbour pixels: a footprint touching the stamp border, and a footprint + that completes an interior row or column beside an off-tile band, stay + neighbours. Parameters ---------- @@ -1779,7 +1783,17 @@ def split_tile_markers(tile_vign, shape): if tile_vign is None: return np.zeros(shape, dtype=bool), np.zeros(shape, dtype=bool) marker = tile_vign == -1e30 - off_tile = marker.all(axis=1)[:, None] | marker.all(axis=0)[None, :] + + def border_runs(full): + # Lines in the unbroken run of marked lines from either border. + lead = np.logical_and.accumulate(full) + trail = np.logical_and.accumulate(full[::-1])[::-1] + return lead | trail + + off_tile = ( + border_runs(marker.all(axis=1))[:, None] + | border_runs(marker.all(axis=0))[None, :] + ) return marker & ~off_tile, off_tile diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 17f04a2a8..8919453c2 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -16,6 +16,7 @@ """ import re +from pathlib import Path from types import SimpleNamespace import numpy as np @@ -41,8 +42,10 @@ make_ngmix_observation, prepare_ngmix_weights, prepare_postage_stamps, + split_tile_markers, uberseg_weight, ) +from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat # --- prepare_ngmix_weights: the filled set is the defect set --------------- @@ -1124,6 +1127,94 @@ def test_off_tile_pixels_are_zero_weighted_and_filled(blend_handling): npt.assert_array_equal(w_out == 0.0, removed) +# --- Catalogue mode: converter-built tile VIGNETs --------------------------- +# +# In catalogue mode the tile VIGNET comes from read_ext_sexcat, which paints +# -1e30 off the image and on other footprints of the DR6 segmentation map. +# The off-image pixels are the off-tile defects and the footprint pixels are +# the neighbour mask, exactly. + +DR6_PATCH = Path(__file__).parent / "data" / "dr6_202.301_seg_patch.fits" + + +def _converter_stamps(seg, number, x, y): + """Converter tile VIGNETs on a unit image, and each stamp's off-image + mask.""" + relabelled, _ = read_ext_sexcat.relabel_seg(seg, number, x, y) + vignets = read_ext_sexcat._extract_vignets( + np.ones(seg.shape, np.float32), x, y, N_STAMP, seg=relabelled, + number=number, + ) + half = N_STAMP // 2 + ny, nx = seg.shape + off_image = [] + for xi, yi in zip(x, y): + rows = int(np.rint(yi)) - 1 - half + np.arange(N_STAMP) + cols = int(np.rint(xi)) - 1 - half + np.arange(N_STAMP) + off_image.append( + ((rows < 0) | (rows >= ny))[:, None] + | ((cols < 0) | (cols >= nx))[None, :] + ) + return vignets, off_image + + +def test_dr6_converter_stamps_split_into_off_image_and_neighbours(): + """On the real 202.301 patch, every converter stamp splits into its + off-image pixels (off-tile) and its other -1e30 pixels (neighbours). + + Failure modes: the converter's float32 -1e30 is not recognised as a + marker; off-image pixels of an edge stamp land in the neighbour mask; + neighbour-footprint pixels become off-tile defects. + """ + with fits.open(DR6_PATCH) as hdul: + seg = hdul["SEG"].data + objects = hdul["OBJECTS"].data + number = np.array(objects["NUMBER"]) + x, y = np.array(objects["X_IMAGE"]), np.array(objects["Y_IMAGE"]) + vignets, off_image = _converter_stamps(seg, number, x, y) + assert vignets.dtype == np.float32 + n_edge = 0 + for vign, off in zip(vignets, off_image): + neighbour, off_tile = split_tile_markers(vign, vign.shape) + npt.assert_array_equal(off_tile, off) + npt.assert_array_equal(neighbour, (vign == _MARKER) & ~off) + n_edge += off.any() + assert n_edge >= 5 + assert sum( + split_tile_markers(v, v.shape)[0].sum() for v in vignets + ) > 1000 + + +def test_a_neighbour_completing_rows_beside_the_tile_edge_stays_a_neighbour(): + """An object 14 px from the tile's left edge, with a wide neighbour + footprint that runs from the tile edge across the stamp: in the rows of + that footprint every stamp pixel is -1e30, off the image or on the + neighbour. Only the off-image columns are off-tile; the footprint is the + neighbour mask, through prepare_postage_stamps. + + Failure mode: every entirely -1e30 row counts as off-tile, so the + neighbour's rows become defects that the epoch cuts count and the defect + fill interpolates (off-tile-is-marked-border-rows-and-columns). + """ + seg = np.zeros((80, 80), np.int32) + seg[20:24, 0:45] = 5 + seg[27:32, 13:18] = 1 + number, x, y = np.array([1, 5]), np.array([15.0, 30.0]), np.array([30.0, 22.0]) + vignets, off_image = _converter_stamps(seg, number, x, y) + tile, off = vignets[0], off_image[0] + footprint = (tile == _MARKER) & ~off + assert footprint.sum() == 4 * (N_STAMP - 11) + assert ((tile == _MARKER).all(axis=1) & ~off.all(axis=1)).sum() == 4 + + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.flags) == len(_MARKER_EPOCH_NAMES) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + ccd = int(names[id(flag)].split("-")[1]) + npt.assert_array_equal(flag == _OFF_TILE, Ngmix.MegaCamFlip(off, ccd)) + npt.assert_array_equal(neighbour, Ngmix.MegaCamFlip(footprint, ccd)) + + # --- DEFECT_FILL = interpolate beside a removed neighbour ------------------- def _defect_beside_neighbour(): From 633425d94973521b1f1f7627cc69116b2e23fcc2 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 02:02:35 +0200 Subject: [PATCH 14/14] docs: -1e30 tile-VIGNET markers come from SExtractor or the catalogue converter The neighbour-marking rationale and the converter's BIG comment describe the markers as ngmix neighbours under flag 2**10; ngmix splits them into off-tile defects and a neighbour mask. The masking text names the tile VIGNET, not SExtractor, as the markers' carrier, since catalogue mode paints them from the DR6 segmentation map. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BbT81hykcvtuefzB3FhYvY --- astra.yaml | 22 ++++++++++--------- src/shapepipe/modules/ngmix_package/ngmix.py | 18 +++++++-------- src/shapepipe/modules/ngmix_runner.py | 2 +- .../read_ext_sexcat.py | 5 +++-- tests/module/test_ngmix_defect_fill.py | 11 +++++----- 5 files changed, 31 insertions(+), 27 deletions(-) diff --git a/astra.yaml b/astra.yaml index f83377c80..bf019aa51 100644 --- a/astra.yaml +++ b/astra.yaml @@ -637,9 +637,10 @@ analyses: catalogue_neighbour_marking: label: Neighbour pixels in the catalogue path's VIGNET rationale: >- - ngmix masks a neighbour only where the tile VIGNET is -1e30 (flag - 2**10: zero weight, noise-filled under noisefill), and SExtractor - writes -1e30 on neighbours' footprints and off the image. The + ngmix's noisefill masks a neighbour only where the tile VIGNET is + -1e30 (zero weight, noise-filled; the off-tile rows and columns are + defects instead), and SExtractor writes -1e30 on neighbours' + footprints and off the image. The unions_catalogue converter reproduces that from the catalogue's own r-band segmentation map (CFIS..r.seg.fits.fz, on the tile's pixel grid), so both tile_detection options mask the same way. @@ -1457,9 +1458,10 @@ analyses: depend on it: defects are filled the same way under every BLEND_HANDLING, and the epoch cuts never count neighbour pixels. noisefill (the default; the committed config sets no key) gives - weight 0 to the pixels SExtractor marks -1e30 in the tile VIGNET - on other detections' footprints (about 93% of neighbour-footprint - pixels on a SExtractor-mode sims tile), excluding the off-tile + weight 0 to the pixels marked -1e30 in the tile VIGNET on other + detections' footprints (about 93% of neighbour-footprint pixels on + a SExtractor-mode sims tile; every one in catalogue mode, from the + segmentation map), excluding the off-tile rows and columns, which are defects, and replaces them with noise; unmarked neighbour pixels keep their weight and light. uberseg ignores the @@ -1487,7 +1489,7 @@ analyses: default: noisefill options: noisefill: - label: Noise-fill SExtractor-marked neighbour pixels (BLEND_HANDLING = noisefill) + label: Noise-fill the marked neighbour pixels (BLEND_HANDLING = noisefill) description: >- Pixels marked -1e30 in the tile VIGNET get weight 0 and noise. The fill stops at the marked footprint, so unmarked neighbour @@ -1527,7 +1529,7 @@ analyses: An epoch is dropped when a defect pixel lies strictly closer to the stamp centre than its fill's radius, beside the masked-fraction cut in the epoch loop; EPOCH_CENTRAL_DEFECT_RADIUS = 0 disables it. - SExtractor's -1e30 neighbour markers are not defects: all epochs + The tile VIGNET's -1e30 neighbour markers are not defects: all epochs share the tile VIGNET, so a neighbour inside the radius would drop every epoch. Off-tile pixels are defects, so an object near the tile edge is vetoed. The veto reads only the defect mask, so it selects on nothing shear-responsive; for the same @@ -1566,8 +1568,8 @@ analyses: (default 1/3) of its pixels in the raw, unsymmetrized defect set is dropped from the multi-epoch fit: flagged pixels (any nonzero exposure flag bit), zero-weight pixels, invalid-RMS pixels and - off-tile pixels (whole -1e30 rows and columns of the tile VIGNET, - flag 2**10), the set defect_fill fills. On a 51-px stamp an object + off-tile pixels (whole -1e30 rows and columns at the tile VIGNET's + border, flag 2**10), the set defect_fill fills. On a 51-px stamp an object within about 8.5 px of the tile edge fails the 1/3 cut. The other -1e30 markers, neighbour footprints, are not counted: every epoch shares the tile VIGNET, so a large neighbour would drop them all. An object with no surviving epoch has no diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 4908fcd20..471e62b2c 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -633,8 +633,8 @@ def __init__( self.psfs = [] self.weights = [] self.flags = [] - # Neighbour masks, one per epoch: the pixels SExtractor marks -1e30 in - # the tile VIGNET on other detections' footprints (off-tile markers + # Neighbour masks, one per epoch: the pixels marked -1e30 in the tile + # VIGNET on other detections' footprints (off-tile markers # are flagged as defects instead; see split_tile_markers), # MegaCam-flipped to the epoch. noisefill zero-weights and noise-fills them; uberseg and # the epoch cuts do not read them (see prepare_ngmix_weights). @@ -791,8 +791,8 @@ class Ngmix(object): :func:`make_ngmix_observation`. blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment. ``"noisefill"`` (default) zero-weights and - noise-fills the pixels SExtractor marks -1e30 in the tile VIGNET - on other detections' footprints; ``"uberseg"`` ignores those markers, + noise-fills the pixels marked -1e30 in the tile VIGNET on other + detections' footprints; ``"uberseg"`` ignores those markers, zeroes the weight of neighbour-side pixels from the coadd segmentation map and requires ``seg_cat_path``. Defect pixels are filled under both (see :func:`prepare_ngmix_weights`). @@ -1555,8 +1555,8 @@ def prepare_postage_stamps( 1/3; 10%, the DES Y3 and Y6 value, is the alternative to test. @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.blend_handling] neighbour-markers-are-not-defects - SExtractor writes -1e30 into the tile VIGNET on the footprints of other - detections and beyond the tile's edge (:func:`split_tile_markers`). The + The tile VIGNET holds -1e30 on the footprints of other detections and + beyond the tile's edge (:func:`split_tile_markers`). The footprint markers form the epoch's neighbour mask (``stamp.neighbours``, MegaCam-flipped like the epoch), kept apart from its flag stamp, so neither the masked-fraction cut nor the central veto counts them. Every @@ -2089,7 +2089,7 @@ def defect_mask(weight, flag, bkg_rms=None): or (when a background RMS map is given) a non-finite or non-positive RMS. This one set is zero-weighted and filled by :func:`prepare_ngmix_weights` under every ``blend_handling`` and counted by the epoch cuts in - :func:`prepare_postage_stamps`. SExtractor's neighbour markers are not + :func:`prepare_postage_stamps`. The tile VIGNET's neighbour markers are not in it (neighbour-markers-are-not-defects); off-tile pixels are, as flag ``OFF_TILE_FLAG`` (off-tile-pixels-are-defects). It is not ORed with its rotations. For defects the central-defect veto keeps @@ -2280,8 +2280,8 @@ def prepare_ngmix_weights( supports the interpolant. @sc [decision:shape_measurement.blend_handling] noisefill-fills-markers - Under ``"noisefill"`` the pixels of ``neighbour``, SExtractor's -1e30 - neighbour markers, get weight 0 and are replaced by the same noise + Under ``"noisefill"`` the pixels of ``neighbour``, the tile VIGNET's + -1e30 neighbour markers, get weight 0 and are replaced by the same noise realisation as the noise-filled defects, so no marked neighbour light reaches metacal. Under the default noise fill, the image, weight map and noise image are those the marked pixels would get as flagged defects; diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 205a38cab..08113aa84 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -145,7 +145,7 @@ def ngmix_runner( centroid_source = "wcs" # Neighbour treatment: "noisefill" (default) zero-weights and noise-fills - # the pixels SExtractor marks -1e30 in the tile VIGNET (other detections' + # the pixels marked -1e30 in the tile VIGNET (other detections' # footprints); "uberseg" ignores those markers, zeroes the weight of every # pixel closer to a neighbour than to the central object, from the # segmentation map, and leaves its image raw. Defect pixels (flagged, diff --git a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py index ff4f58cc8..4b1d66012 100644 --- a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py +++ b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py @@ -48,8 +48,9 @@ def _build_ldac_imhead(img_header): # SExtractor's VIGNET value for pixels that are not the object's: off the -# image and on a neighbour's segmentation footprint. ngmix flags these pixels -# (tile VIGNET == -1e30) and noise-fills them at zero weight. +# image and on a neighbour's segmentation footprint. ngmix splits these pixels +# (tile VIGNET == -1e30, ``ngmix.split_tile_markers``): the off-image rows and +# columns are defects, the footprint pixels its neighbour mask. BIG = -1e30 # The label of a footprint no catalogue object claims in the relabelled diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 8919453c2..f068c1290 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -10,7 +10,7 @@ in :func:`prepare_postage_stamps` act on the same defect set: the masked-fraction cut counts it, and the central-defect veto drops an epoch with a defect near the stamp centre, at a radius set by the defect's fill. -SExtractor's -1e30 neighbour markers in the tile VIGNET are not defects: +The tile VIGNET's -1e30 neighbour markers are not defects: noisefill zero-weights and noise-fills them, uberseg ignores them, and the epoch cuts never count them. """ @@ -744,7 +744,7 @@ def test_ngmix_rejects_an_unknown_defect_fill(tmp_path): ) -# --- SExtractor's -1e30 neighbour markers are not defects ------------------- +# --- The tile VIGNET's -1e30 neighbour markers are not defects ------------- # # The tile VIGNET carries -1e30 on the footprints of other detections. Every # epoch shares that tile stamp, so a marker counted as a defect would drop @@ -1021,9 +1021,10 @@ def fake_observation(*args, **kwargs): # --- Off-tile pixels are defects ------------------------------------------- # -# SExtractor also writes -1e30 beyond the tile's edge, where the epoch holds -# the object's own light, cut off. Those pixels are the stamp rows and -# columns that are entirely -1e30 (the off-image part of a rectangle clip); +# The tile VIGNET also holds -1e30 beyond the tile's edge, where the epoch +# holds the object's own light, cut off. Those pixels are the runs of +# entirely -1e30 stamp rows and columns that start at a stamp border (the +# off-image part of a rectangle clip); # they join the epoch's defect set as flag 2**10. The other markers are the # neighbour mask.