ngmix: opt-in DEFECT_FILL = interpolate for narrow defects, with a 7 px veto

astra.yaml

@@ -1329,8 +1329,8 @@ analyses:
           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". Under every BLEND_HANDLING, defect pixels are replaced with
-          independent noise at the per-pixel background RMS when supplied, or
+          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
@@ -1355,7 +1355,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
@@ -1466,11 +1470,12 @@ 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. The 7 px interpolated radius arrives with
-          defect_fill's interpolate option (feat/defect-interpolation); until
-          then every defect is noise-filled and the 10 px radius applies.
+          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_CENTRAL_DEFECT_RADIUS = 10;
+          EPOCH_INTERPOLATED_DEFECT_RADIUS = 7.
         default: fixed_radii
         options:
           disabled:

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