ngmix masked pixels: neighbours and defects as separate questions

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.<tile>.r.seg.fits.fz, on the tile's
           pixel grid), so both tile_detection options mask the same way.
@@ -1376,22 +1377,26 @@ 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
           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". In the committed default (`BLEND_HANDLING = noisefill`),
-          masked pixels are replaced with independent noise at the per-pixel
-          background RMS when supplied, or the stamp's robust noise scale
-          otherwise. With `BLEND_HANDLING = uberseg`, the image is left
-          untouched and weights are zeroed on defects and neighbour-side
-          pixels. [LINT] the prepare_ngmix_weights docstring says noisefill
-          keeps the weight of filled pixels (the code zeroes it), and the
-          ngmix_runner comment says noisefill fills neighbour pixels (it fills
-          flagged pixels and leaves neighbours untouched). The committed fill
+          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. 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
@@ -1413,7 +1418,11 @@ analyses:
           interpolate:
             label: Interpolate short bounded runs; noise-fill the rest
             description: >-
-              Not implemented. PR #916 measures c1 = -1.3e-3 for a column
+              DEFECT_FILL = interpolate: row or column runs of at most
+              MAX_INTERPOLATED_RUN = 3 defect pixels that stop short of the
+              stamp border are Clough-Tocher interpolated from the clean pixels
+              within SUPPORT_RADIUS = 4 px; other defects are noise-filled.
+              PR #916 measures c1 = -1.3e-3 for a column
               8 px from a 0.5 arcsec galaxy and +2e-6 when the weight is also
               zeroed on its quarter-turn orbit. For a 3-px bleed 6 px from that
               galaxy, PR #916 measures m11 = +0.89% when the fill is
@@ -1431,11 +1440,11 @@ analyses:
               with symmetrization against 0.7e-4 without it.
             insights: [mask_bad_column_symmetrize, mask_des_defect_practice, mask_fixed_orientation]
           raw:
-            label: "No fill: raw defect values (BLEND_HANDLING = uberseg)"
+            label: "No fill: raw defect values"
             description: >-
-              With BLEND_HANDLING = uberseg, defect pixels keep weight 0 but
-              their image values stay untouched in the image that metacal
-              deconvolves, shears and reconvolves.
+              Defect pixels keep weight 0 but their image values stay
+              untouched in the image that metacal deconvolves, shears and
+              reconvolves.
             excluded: true
             excluded_reason: >-
               Metacal acts on every pixel regardless of weight, so raw defects
@@ -1445,10 +1454,18 @@ analyses:
       blend_handling:
         label: Neighbour treatment before metacal
         rationale: >-
-          Covers only pixels shared with a neighbour; defect_fill is coupled to
-          it through 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
+          Covers only pixels shared with a neighbour; defect_fill does not
+          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 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
+          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
@@ -1469,15 +1486,16 @@ 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 = noisefill)
+          noisefill:
+            label: Noise-fill the 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:
@@ -1508,11 +1526,13 @@ analyses:
       central_defect_veto:
         label: Per-epoch veto on a defect near the stamp centre
         rationale: >-
-          Not on develop; implemented on feat/defect-fill-veto (7777181b).
-          There 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. The veto reads only the
-          defect mask, so it selects on nothing shear-responsive; for the same
+          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 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
           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,
@@ -1523,16 +1543,18 @@ analyses:
           wide defects through the elliptical PSF sit at the 1% bound (m11 =
           -0.98%; -0.24% at 11 px), and noise fill needs 14 px for 0.7 and 0.9
           arcsec galaxies. Measured on feat/defect-fill-veto and
-          feat/defect-interpolation.
-        default: disabled
+          feat/defect-interpolation. The 7 px radius applies only under
+          DEFECT_FILL = interpolate; with the committed noise fill every defect
+          takes the 10 px radius.
+          Values:
+          EPOCH_CENTRAL_DEFECT_RADIUS = 10;
+          EPOCH_INTERPOLATED_DEFECT_RADIUS = 7.
+        default: fixed_radii
         options:
           disabled:
-            label: No central veto (committed code)
+            label: No central veto (EPOCH_CENTRAL_DEFECT_RADIUS = 0)
           fixed_radii:
             label: Fixed radii, 10 px for noise-filled and 7 px for interpolated defects
-            description: >-
-              Implemented on feat/defect-fill-veto and feat/defect-interpolation,
-              not on develop.
           size_scaled_radius:
             label: Veto radius scaled by galaxy size
             excluded: true
@@ -1542,14 +1564,16 @@ analyses:
       epoch_masked_fraction_cut:
         label: Per-epoch masked-fraction cut
         rationale: >-
-          [HARDCODED] on develop, an epoch whose stamp has more than 1/3 of its
-          pixels flagged (any nonzero flag bit, including the tile-coverage bit
-          2**10 set where the tile vignet is off-image) is dropped from the
-          multi-epoch fit; an object with no surviving epoch has no shape.
-          Zero-weight and invalid-RMS pixels are not counted. On
-          feat/defect-fill-veto the cut counts the raw, unsymmetrized defect
-          set (flagged, zero-weight and invalid-RMS pixels) against
-          EPOCH_MASKED_FRACTION_CUT, default 1/3. Before the cut, an epoch is
+          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, invalid-RMS pixels and
+          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
+          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:
           Y1 rejected any epoch with a masked or zero-weight pixel, or with
@@ -1565,12 +1589,12 @@ analyses:
             label: 1/3 of the stamp in the defect set
             description: >-
               Drop an epoch only if more than 1/3 of the stamp pixels are
-              defects (on develop, flagged pixels).
+              defects.
           ten_percent:
             label: 10% (DES Y3 / Y6)
             description: >-
-              Not a default; EPOCH_MASKED_FRACTION_CUT = 0.1 on
-              feat/defect-fill-veto. Matches DES Y3 max_zero_weight_frac and
+              Not a default; EPOCH_MASKED_FRACTION_CUT = 0.1. Matches DES Y3
+              max_zero_weight_frac and
               Y6 max_masked_fraction.
             insights: [mask_multi_epoch_drop]
           any_masked:

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)

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, 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,
+    ) & ~excluded
+    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, excluded, target):
+    """Replace the ``target`` pixels of every plane by a Clough-Tocher
+    interpolant of the kept pixels around them.
+
+    @sc [decision:shape_measurement.defect_fill] shared-rotation-averaged-interpolant
+    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
+    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)``.
+    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`).
+
+    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)
+    excluded = np.asarray(excluded, 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(excluded, k),
+                np.rot90(target, k),
+            ),
+            -k, axes=(1, 2),
+        )[:, target]
+        for k in range(4)
+    ]
+    out[:, target] = np.mean(turns, axis=0)
+    return out