309.2. Detecting LSB structures

309.2. Detecting LSB structures¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: large
LSST Science Pipelines version: r30.0.11
Last verified to run: 2026-08-27
Repository: github.com/lsst/tutorial-notebooks

Learning objective: Apply pixel binning and adaptive smoothing to reveal low surface brightness structures in DP2 coadd images, and compare how each technique affects the detectable extent of galaxy outskirts through radial surface brightness profiles.

LSST data products: deep_coadd, skyMap

Packages: lsst.daf.butler, lsst.geom, lsst.afw.display, numpy, matplotlib, scipy, astropy

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

Low surface brightness structures span a wide range of angular scales and physical origins: nearby dwarf satellite galaxies can subtend tens to hundreds of arcseconds; tidal streams and shells can extend over arcminutes; and intracluster light can fill degrees around nearby rich clusters. All share the property that their surface brightness is comparable to, or fainter than, the per-pixel noise. Detecting them requires suppressing pixel-to-pixel noise while preserving diffuse flux.

This notebook applies two complementary noise-suppression techniques to the DP2 Extended Chandra Deep Field South (ECDFS) coadd. Pixel binning block-averages $N\times N$ pixels to reduce noise by $\sqrt{N}$ while exactly preserving surface brightness. Adaptive smoothing adjusts the kernel size pixel-by-pixel so that bright regions retain resolution while faint diffuse regions are smoothed only as much as needed to achieve a target signal-to-noise ratio. The closing section compares radial SNR profiles to show at what radius the azimuthal-median galaxy flux falls to the residual noise level of each technique — the point at which the signal becomes undetectable with that approach.

Related tutorials: Tutorial 309.1 introduces the display stretches used throughout this series; 309.3 covers interactive surface brightness measurement of tidal features.

1.1. Import packages¶

lsst.daf.butler and lsst.geom access the data; lsst.afw.display renders images with an astronomically appropriate stretch. astropy.visualization provides ImageNormalize and AsinhStretch for display; scipy.ndimage.uniform_filter powers the normalized convolution inside adaptive_smooth. Standard scientific Python packages (numpy, matplotlib) handle array operations and figures. lsst.images.get_legacy_deep_coadd_mask_planes maps mask plane names back to their legacy equivalents, and lsst.images.Box builds pixel bounding boxes for cutout retrieval.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import uniform_filter
from astropy.visualization import AsinhStretch, ImageNormalize

import lsst.geom as geom
from lsst.daf.butler import Butler
from lsst.images import Box, get_legacy_deep_coadd_mask_planes

1.2. Define parameters and functions¶

Set plotting defaults and the color constants used for figures, then define TARGET, the sky coordinate of the target galaxy, followed by the helper functions used throughout: build_sky_arr and sky_stats for sky-noise-anchored statistics, nJy_per_pix_to_mu for unit conversion, and sky_norm, bin_image, and adaptive_smooth for display and diffuse-structure enhancement.

In [2]:

plt.rcParams.update(
    {"figure.dpi": 150, "image.origin": "lower", "axes.titlesize": 10})

BLUE = "#648FFF"
PURPLE = "#785EF0"
MAGENTA = "#DC267F"
ORANGE = "#FE6100"
YELLOW = "#FFB000"

In [3]:

TARGET = (53.010, -28.350)  # (RA, Dec)

build_sky_arr masks detected sources and detector artifacts by setting flagged pixels to NaN, preventing bright compact objects from biasing sky statistics. sky_stats returns display limits and the per-pixel sky noise $\hat{\sigma}$ (the same quantity used throughout Tutorial 309.1).

Sky noise is the pixel-to-pixel flux variation on blank patches of sky, measured after masking sources and artifacts. It sets the fundamental detection floor: any structure fainter than a few $\hat{\sigma}$ per pixel is invisible without additional noise suppression. The noise is estimated robustly as $\sigma_\mathrm{sky} = 1.4826 \times \mathrm{MAD}$, where the factor 1.4826 makes the MAD-based estimate consistent with the standard deviation of a Gaussian distribution.

In [4]:

def build_sky_arr(
    coadd,
    arr,
    bad_planes=("SAT", "NO_DATA", "BAD", "CR", "BRIGHT_OBJECT"),
    source_planes=("DETECTED",),
):
    """Return arr with artifact and source pixels set to NaN.

    Parameters
    ----------
    coadd : lsst.images.CellCoadd
        Coadd (or cutout stamp) whose mask attribute is used to identify bad pixels.
    arr : ndarray, shape (ny, nx)
        Pixel array to mask, typically the image or variance plane.
    bad_planes : tuple of str, optional
        Legacy mask plane names treated as artifacts.
    source_planes : tuple of str, optional
        Legacy mask plane names treated as detected sources.

    Returns
    -------
    masked_arr : ndarray, shape (ny, nx)
        Copy of arr with artifact, source, and non-finite pixels set to NaN.
    """
    legacy_to_new = {
        v.name: k for k, v in get_legacy_deep_coadd_mask_planes().items()
    }
    mask = coadd.mask
    schema_names = {plane.name for plane in mask.schema}
    exclude = np.zeros(arr.shape, dtype=bool)
    for plane in list(bad_planes) + list(source_planes):
        new_plane = legacy_to_new.get(plane, plane)
        if new_plane in schema_names:
            exclude |= mask.get(new_plane).astype(bool)
    return np.where(exclude | ~np.isfinite(arr), np.nan, arr)

In [5]:

def sky_stats(arr, lo_sigma=2.0, hi_sigma=8.0):
    """Return (vmin, vmax) anchored to the sky background.

    Uses the robust median as sky level and 1.4826 x MAD as pixel noise.
    Pass a source-masked array (e.g. from build_sky_arr) so that detected
    sources do not bias the estimates.

    Parameters
    ----------
    arr : ndarray
        Pixel array, typically source-masked (NaN at bad/source pixels).
    lo_sigma : float, optional
        Lower stretch limit in units of sky noise below the sky median.
    hi_sigma : float, optional
        Upper stretch limit in units of sky noise above the sky median.

    Returns
    -------
    vmin : float
        Lower display limit (sky - lo_sigma * noise).
    vmax : float
        Upper display limit (sky + hi_sigma * noise).
    sigma : float
        Per-pixel noise estimate (1.4826 * MAD).
    """
    flat = arr[np.isfinite(arr)]
    sky = np.median(flat)
    sigma = 1.4826 * np.median(np.abs(flat - sky))
    return sky - lo_sigma * sigma, sky + hi_sigma * sigma, sigma

nJy_per_pix_to_mu converts pixel flux in nJy to AB surface brightness in mag arcsec$^{-2}$, using PIXEL_AREA.

PIXEL_SCALE and PIXEL_AREA are not passed to the function as arguments: they are assigned as global variables in Section 2, immediately after the coadd is loaded. Python looks up a name used inside a function body in the global namespace at the time the function is called, not when it is defined, so nJy_per_pix_to_mu can be defined here, before PIXEL_AREA exists, and still work correctly as long as the notebook is run top to bottom before the function is used.

In [6]:

def nJy_per_pix_to_mu(flux_nJy):
    """Convert nJy pixel**-1 to AB mag arcsec**-2 for the LSST pixel scale.

    Parameters
    ----------
    flux_nJy : float or array_like
        Pixel flux in nJy pixel**-1.

    Returns
    -------
    mu : float or array_like
        Surface brightness in AB mag arcsec**-2.
    """
    return -2.5 * np.log10(np.abs(flux_nJy) * 1e-9 / PIXEL_AREA / 3631.0)

sky_norm wraps sky_stats into an ImageNormalize object for use with imshow.

In [7]:

def sky_norm(arr, lo_sigma=2.0, hi_sigma=8.0):
    """Return an ImageNormalize anchored to the sky background.

    Wraps sky_stats. Pass a source-masked array so that detected objects do
    not bias the sky median or noise estimate.

    Parameters
    ----------
    arr : ndarray
        Pixel array, typically source-masked (NaN at bad/source pixels).
    lo_sigma : float, optional
        Lower stretch limit in units of sky noise below the sky median.
    hi_sigma : float, optional
        Upper stretch limit in units of sky noise above the sky median.

    Returns
    -------
    norm : astropy.visualization.ImageNormalize
        Normalization with asinh stretch and sky-anchored limits.
    """
    vmin, vmax, _ = sky_stats(arr, lo_sigma=lo_sigma, hi_sigma=hi_sigma)
    return ImageNormalize(vmin=vmin, vmax=vmax, stretch=AsinhStretch())

bin_image performs $N\times N$ block averaging with NaN propagation, reducing pixel noise while preserving surface brightness.

In [8]:

def bin_image(arr, bin_factor):
    """Block-average arr by bin_factor in each dimension, ignoring NaN (masked) pixels.

    Parameters
    ----------
    arr : ndarray, shape (ny, nx)
        Input pixel array; NaN at masked pixels.
    bin_factor : int
        Downsampling factor; output shape is (ny // bin_factor, nx // bin_factor).

    Returns
    -------
    binned : ndarray, shape (ny // bin_factor, nx // bin_factor)
        Block-averaged array.
    """
    ny, nx = arr.shape
    ny_trim = (ny // bin_factor) * bin_factor
    nx_trim = (nx // bin_factor) * bin_factor
    trimmed = arr[:ny_trim, :nx_trim]
    with np.errstate(invalid="ignore"):
        return np.nanmean(
            trimmed.reshape(
                ny_trim
                // bin_factor,
                bin_factor,
                nx_trim
                // bin_factor,
                bin_factor),
            axis=(
                1,
                3),
        )

adaptive_smooth applies a spatially varying kernel via normalized convolution, using a larger kernel where the local signal-to-noise ratio is low.

In [9]:

def adaptive_smooth(arr, noise_sigma, target_snr=3.0, kernel_sizes=None):
    """Adaptive kernel smoothing using normalized convolution to handle masked pixels.

    At each pixel, selects the smallest kernel for which |smoothed / noise_eff| >= target_snr.
    Masked (NaN) pixels are excluded from each kernel window via normalized convolution, so
    they do not bleed into adjacent valid pixels.

    Parameters
    ----------
    arr : ndarray, shape (ny, nx)
        Input pixel array; NaN at bad or source-masked pixels.
    noise_sigma : float
        Per-pixel sky noise from clean (source-masked) pixels, in nJy pixel**-1.
    target_snr : float, optional
        Target |signal/noise| ratio used to select the smoothing kernel.
    kernel_sizes : list of int, optional
        Uniform-filter widths in pixels; defaults to [3, 7, 15, 31, 61].

    Returns
    -------
    result : ndarray, shape (ny, nx)
        Adaptively smoothed array; NaN at bad pixels.
    kernel_map : ndarray, shape (ny, nx)
        Kernel width applied at each pixel; NaN at bad pixels.
    """
    if kernel_sizes is None:
        kernel_sizes = [3, 7, 15, 31, 61]

    valid = np.isfinite(arr)
    arr_z = np.where(valid, arr, 0.0)
    valid_f = valid.astype(np.float64)

    def _nconv(k):
        num = uniform_filter(arr_z, size=k)
        den = uniform_filter(valid_f, size=k)
        with np.errstate(invalid="ignore", divide="ignore"):
            return np.where(den > 0.1, num / den, np.nan)

    result = _nconv(kernel_sizes[-1])
    kernel_map = np.full(arr.shape, float(kernel_sizes[-1]))

    for k in reversed(kernel_sizes[:-1]):
        sm = _nconv(k)
        noise_eff = noise_sigma / k
        snr = np.abs(sm) / noise_eff
        meets = (snr >= target_snr) & np.isfinite(sm)
        result = np.where(meets, sm, result)
        kernel_map = np.where(meets, float(k), kernel_map)

    return np.where(valid, result, np.nan), np.where(valid, kernel_map, np.nan)

2. Load a coadd patch¶

Instantiate the Butler with the DP2 repository, load the $i$-band deep_coadd for the ECDFS field centred on the target galaxy, and apply the 'pretty' background (deep_coadd pixel values are not background-subtracted by default). Display the full patch with a sky-noise anchored stretch.

In [10]:

butler = Butler.from_config('dp2', collections='dp2')
assert butler is not None

skymap_name = "lsst_cells_v2"

target_ra, target_dec = TARGET
band = "i"

skymap = butler.get("skyMap", skymap=skymap_name, collections="skymaps")
sky_point = geom.SpherePoint(
    target_ra * geom.degrees, target_dec * geom.degrees)
tract_info = skymap.findTract(sky_point)
patch_info = tract_info.findPatch(sky_point)
tract_id = tract_info.getId()
patch_id = patch_info.getSequentialIndex()
patch_wcs = tract_info.getWcs()
patch_bbox = patch_info.getOuterBBox()

data_id = {"band": band, "tract": tract_id, "patch": patch_id}
coadd = butler.get("deep_coadd", dataId=data_id)
coadd.apply_background("pretty")
x0, y0 = coadd.bbox.x.start, coadd.bbox.y.start
sky_00 = coadd.sky_projection.pixel_to_sky(x=x0, y=y0)
sky_10 = coadd.sky_projection.pixel_to_sky(x=x0 + 1, y=y0)
PIXEL_SCALE = sky_00.separation(sky_10).arcsec
PIXEL_AREA = PIXEL_SCALE**2
print(f"Pixel scale: {PIXEL_SCALE: .4f} arcsec/pixel")
coadd_arr = coadd.image.array.copy()
print(f"Loaded {band}-band coadd: {coadd_arr.shape}")

Pixel scale:  0.2000 arcsec/pixel
Loaded i-band coadd: (3300, 3300)

Display the full coadd patch with a sky-noise anchored asinh stretch.

In [11]:

coadd_sky_full = build_sky_arr(coadd, coadd_arr)
norm_patch = sky_norm(coadd_sky_full)

fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(coadd_arr, norm=norm_patch, cmap="Greys_r", origin="lower")
ax.set_title(f"{band}-band deep_coadd: tract {tract_id}, patch {patch_id}")
ax.set_xlabel("x (pixels)")
ax.set_ylabel("y (pixels)")

Out[11]:

Text(0, 0.5, 'y (pixels)')

Figure 1: Full $i$-band deep_coadd patch containing the target galaxy (RA = 53.010°, Dec = −28.350°), displayed with an asinh stretch anchored to the sky noise estimated from source-masked pixels.

Build display and analysis masks from the pipeline mask planes. coadd_display retains detected-source pixels for a natural-looking image; coadd_sky masks detected sources and artifacts, leaving only blank sky pixels. sky_stats estimates the per-pixel noise sky_sigma from coadd_sky.

In [12]:

coadd_display = build_sky_arr(
    coadd, coadd_arr, source_planes=())

coadd_sky = build_sky_arr(coadd, coadd_arr)
_, _, sky_sigma = sky_stats(coadd_sky)

sky_val = np.nanmedian(coadd_sky)
mu_sigma = nJy_per_pix_to_mu(sky_sigma)
mu_val = nJy_per_pix_to_mu(sky_val)
print(f"Sky level: {sky_val: .1f} nJy/pix = {mu_val: .1f} mag/arcsec^2")
print(f"Per-pixel 1-sigma noise: {sky_sigma: .1f} nJy/pix = {mu_sigma: .1f} mag/arcsec^2")

Sky level:  0.0 nJy/pix =  31.8 mag/arcsec^2
Per-pixel 1-sigma noise:  3.4 nJy/pix =  26.6 mag/arcsec^2

3. Pixel binning to reveal diffuse structure¶

Block-averaging $N\times N$ pixels reduces pixel-to-pixel noise by $\sqrt{N}$ while exactly preserving surface brightness — the mean of $N$ sky pixels estimates the sky level equally well regardless of bin size, and the surface brightness units (nJy per original pixel area) remain unchanged.

Masked (NaN) pixels are excluded from each block average automatically.

3.1. Noise scaling¶

The table below converts the $\sqrt{N}$ noise scaling to surface brightness depth, assuming uncorrelated pixel noise. Surface brightness is an intensive quantity: it measures flux per unit area and does not depend on how many pixels are averaged, so the sky level in $\mathrm{mag\ arcsec}^{-2}$ is unchanged by binning. Only the noise floor — the minimum surface brightness detectable at $1\hat{\sigma}$ per binned pixel — drops by $\sqrt{N}$ with each factor-of-$N$ binning step.

Warning: This cell may produce a RuntimeWarning: Mean of empty slice from numpy. This is expected: bin blocks that fall entirely within masked regions (detector edges, saturated columns) contain only NaN pixels, so np.nanmean has nothing to average and returns NaN for those blocks. The warning can be safely ignored — the output array is correct.

In [13]:

mu_1x = nJy_per_pix_to_mu(sky_sigma)
print(f"Per-pixel 1-sigma noise: {sky_sigma: .1f} nJy/pix -> {mu_1x: .1f} mag/arcsec^2\n")
cols = [
    f"{'Bin factor': >10}",
    f"{'Resolution': >12}",
    f"{'sigma_eff (nJy/pix)': >20}",
    f"{'1-sigma depth (mag/arcsec^2)': >28}",
    f"{'Gain (mag)': >10}",
]
print("  ".join(cols))
print("-" * 85)
for bf in [1, 2, 4, 8, 16, 32]:
    sigma_eff = sky_sigma / bf
    mu = nJy_per_pix_to_mu(sigma_eff)
    gain = mu - mu_1x
    res_str = f'{bf * PIXEL_SCALE: .1f}"'
    row = [
        f"{bf: >10}",
        f"{res_str: >12}",
        f"{sigma_eff: >20.1f}",
        f"{mu: >28.1f}",
        f"{gain: >+10.2f}",
    ]
    print("  ".join(row))

bin_factors = [1, 4, 8, 16]
binned_arrays = [bin_image(coadd_display, bf) for bf in bin_factors]
shared_norm = sky_norm(bin_image(coadd_sky, 1), lo_sigma=2.0, hi_sigma=8.0)

Per-pixel 1-sigma noise:  3.4 nJy/pix ->  26.6 mag/arcsec^2

Bin factor    Resolution   sigma_eff (nJy/pix)  1-sigma depth (mag/arcsec^2)  Gain (mag)
-------------------------------------------------------------------------------------
         1          0.2"                   3.4                          26.6       +0.00
         2          0.4"                   1.7                          27.3       +0.75
         4          0.8"                   0.8                          28.1       +1.51
         8          1.6"                   0.4                          28.8       +2.26
        16          3.2"                   0.2                          29.6       +3.01
        32          6.4"                   0.1                          30.3       +3.76

/tmp/ipykernel_1667/2531585028.py:21: RuntimeWarning: Mean of empty slice
  return np.nanmean(

Quantify the noise reduction as a function of bin factor. Since all real astrophysical emission contributes positive flux, pixels with negative values can only arise from noise. Taking the MAD of these negative-valued pixels — scaled by 1.4826 to match the standard deviation of a Gaussian — gives a noise estimate uncontaminated by faint galaxy outskirts, unmasked PSF wings, or diffuse emission. This negative-half MAD estimator is compared to the $\sqrt{N}$ expectation in Figure 2 right.

Warning: This cell may produce a RuntimeWarning: Mean of empty slice from numpy. This is expected: bin blocks that fall entirely within masked regions (detector edges, saturated columns) contain only NaN pixels, so np.nanmean has nothing to average and returns NaN for those blocks. The warning can be safely ignored — the output array is correct.

In [14]:

x_range = 4.0 * sky_sigma
measured_sigmas = []
colors = [YELLOW, ORANGE, MAGENTA, BLUE]

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

for bf, color in zip(bin_factors, colors):
    binned_lsb = bin_image(coadd_sky, bf)
    pixels = binned_lsb[np.isfinite(binned_lsb)].ravel()
    neg = pixels[pixels < 0]
    meas_sigma = 1.4826 * np.median(np.abs(neg))
    measured_sigmas.append(meas_sigma)

    axes[0].hist(
        pixels,
        bins=120,
        range=(-x_range, x_range),
        density=True,
        histtype="step",
        color=color,
        lw=1.5,
        label=f"Bin {bf}x | sigma = {meas_sigma: .2f} nJy",
    )

axes[0].set_xlabel("Pixel value (nJy)")
axes[0].set_ylabel("Probability density")
axes[0].set_title("Sky pixel distributions at each bin factor")
axes[0].legend(fontsize=8)
axes[0].axvline(0, color="k", ls=":", lw=0.8)

expected_sigmas = [sky_sigma / bf for bf in bin_factors]
axes[1].plot(bin_factors, expected_sigmas, "k--o",
             ms=6, label="Expected: sigma / bin_factor")
axes[1].plot(
    bin_factors,
    measured_sigmas,
    "s-",
    ms=6,
    color=ORANGE,
    label="Measured (negative-half MAD)",
)
axes[1].set_xlabel("Bin factor")
axes[1].set_ylabel("sigma (nJy)")
axes[1].set_title("Noise reduction: expected vs measured")
axes[1].set_xscale("log", base=2)
axes[1].set_yscale("log")
axes[1].set_xticks(bin_factors)
axes[1].set_xticklabels([str(b) for b in bin_factors])
axes[1].legend()

/tmp/ipykernel_1667/2531585028.py:21: RuntimeWarning: Mean of empty slice
  return np.nanmean(

Out[14]:

<matplotlib.legend.Legend at 0x7b2bb5fc2210>

Figure 2: Left — sky pixel histograms at each bin factor; the distribution narrows as noise decreases. Right — measured $\sigma$ versus the $\sqrt{N}$ expectation. The ratio exceeds 1 at large bin factors, indicating spatially correlated background residuals at arcminute scales.

At bin factor 1 the measured and expected noise agree by construction. At larger bin factors the measured $\sigma$ exceeds the $\sqrt{N}$ prediction, reflecting spatially correlated background residuals left by the pipeline sky subtraction. The divergence begins at bin factor 4 (0.8 arcsec) and grows steadily, indicating that background correlations extend over several arcseconds.

3.2. Visual comparison¶

Zoom in on the target galaxy at each bin factor to show how resolution and noise trade off at the scale of the galaxy outskirts (field of view $\pm$120 arcsec around the target).

In [15]:

pixel_target = patch_wcs.skyToPixel(sky_point)
target_px = int(pixel_target.getX()) - patch_bbox.getMinX()
target_py = int(pixel_target.getY()) - patch_bbox.getMinY()
print(f"Target galaxy patch-relative position: ({target_px}, {target_py})")

zoom_half_arcsec = 120.0
zoom_half_pix = int(zoom_half_arcsec / PIXEL_SCALE)

fig, axes = plt.subplots(
    1, len(bin_factors), figsize=(5 * len(bin_factors), 5), squeeze=False
)
for col, (bf, barr) in enumerate(zip(bin_factors, binned_arrays)):
    bx = target_px // bf
    by = target_py // bf
    bh = max(1, zoom_half_pix // bf)

    x0 = max(0, bx - bh)
    x1 = min(barr.shape[1], bx + bh)
    y0 = max(0, by - bh)
    y1 = min(barr.shape[0], by + bh)

    extent = [
        (x0 - bx) * bf * PIXEL_SCALE,
        (x1 - bx) * bf * PIXEL_SCALE,
        (y0 - by) * bf * PIXEL_SCALE,
        (y1 - by) * bf * PIXEL_SCALE,
    ]
    axes[0][col].imshow(
        barr[y0:y1, x0:x1],
        norm=shared_norm,
        cmap="Greys_r",
        origin="lower",
        extent=extent,
    )
    axes[0][col].set_title(
        f'Bin {bf}x  ({
            bf
            * PIXEL_SCALE: .1f}" res.)\nsigma = {
            sky_sigma
            / bf: .0f} nJy')
    axes[0][col].set_xlabel("Delta arcsec")
    if col == 0:
        axes[0][col].set_ylabel("Delta arcsec")

Target galaxy patch-relative position: (1168, 2953)

Figure 3: Zoom-in ($\pm120$ arcsec) around the target galaxy at each bin factor on a shared stretch. Cutouts are clipped to the patch boundary, so a panel may not reach $\pm120$ arcsec on all sides if the target lies near the patch edge.

4. Adaptive smoothing¶

Pixel binning treats every location in the image identically. Adaptive smoothing instead applies a spatially varying kernel: small where the signal is bright enough to be detected at fine resolution, large where the signal is faint and noise dominates. The result preserves compact structure while still revealing diffuse emission.

The adaptive_smooth helper (defined in Section 1.2) implements this approach. It accepts a list of candidate kernel sizes as the kernel_sizes parameter; here kernel_sizes = [1, 3, 7, 15, 31, 61] pixels, chosen to span roughly 0.2 arcsec to 20 arcsec in steps of approximately a factor of two. At each pixel it selects the smallest kernel for which $|\text{smoothed value}| / (\sigma_{\rm sky} / k) \geq 3$. Masked pixels are excluded via normalized convolution (zeroing masked pixels in both numerator and denominator), so they do not bleed into adjacent valid pixels.

First, locate the host galaxy within the patch and extract a postage stamp centred on it; the stamp is used for the algorithm demonstration in Section 4.1.

In [16]:

host_sky_point = sky_point
host_pix = patch_wcs.skyToPixel(host_sky_point)

host_x = int(host_pix.getX()) - patch_bbox.getMinX()
host_y = int(host_pix.getY()) - patch_bbox.getMinY()
print(f"Host galaxy pixel position (patch-relative): ({host_x}, {host_y})")

coadd_ny, coadd_nx = coadd_arr.shape
if 0 <= host_x < coadd_nx and 0 <= host_y < coadd_ny:
    print(f"Pixel value at host: {coadd_arr[host_y, host_x]: .1f} nJy/pix")
else:
    print("WARNING: position outside patch")

Host galaxy pixel position (patch-relative): (1168, 2953)
Pixel value at host:  279.2 nJy/pix

Extract the postage stamp centred on the host galaxy by slicing it from the already background-subtracted patch coadd. Two versions of the stamp array are built with build_sky_arr: stamp_display masks only genuinely bad pixels (saturated columns, cosmic rays, missing data) and is used for display; stamp_sky additionally masks pixels flagged as detected sources, leaving only blank sky for noise estimation. stamp_sigma is the per-pixel sky noise derived from stamp_sky and is used as the noise threshold in both adaptive smoothing and the SNR profile.

In [17]:

stamp_half = 600
host_x_full, host_y_full = int(host_pix.getX()), int(host_pix.getY())
patch_box_new = Box.factory[
    patch_bbox.getMinY(): patch_bbox.getMaxY() + 1,
    patch_bbox.getMinX(): patch_bbox.getMaxX() + 1,
]
stamp_bbox = Box.factory[
    host_y_full - stamp_half: host_y_full + stamp_half,
    host_x_full - stamp_half: host_x_full + stamp_half,
].intersection(patch_box_new)

stamp = coadd[stamp_bbox].copy()
stamp_arr = stamp.image.array.copy()

stamp_display = build_sky_arr(
    stamp, stamp_arr, source_planes=())
stamp_bad = ~np.isfinite(stamp_display)

stamp_sky = build_sky_arr(stamp, stamp_arr)
_, _, stamp_sigma = sky_stats(stamp_sky)
mu_stamp = nJy_per_pix_to_mu(stamp_sigma)
print(f"Stamp per-pixel 1-sigma: {stamp_sigma: .1f} nJy/pix -> {mu_stamp: .1f} mag/arcsec^2")

Stamp per-pixel 1-sigma:  3.9 nJy/pix ->  26.4 mag/arcsec^2

4.1. Algorithm demonstration¶

Apply adaptive_smooth (Section 1.2) to the postage stamp centred on the host galaxy. The stamp is used here rather than the full patch so the kernel size map can be inspected at the scale of an individual galaxy. Print the fraction of stamp pixels assigned to each kernel size, then display the kernel size map alongside the original stamp to show spatially where fine or coarse smoothing was applied. A direct comparison with fixed 8x pixel binning follows.

In [18]:

kernel_sizes = [1, 3, 7, 15, 31, 61]
smooth_result, kernel_map = adaptive_smooth(
    stamp_display, stamp_sigma, target_snr=3.0, kernel_sizes=kernel_sizes
)

kvals = kernel_map[np.isfinite(kernel_map)].ravel()
for ks in kernel_sizes:
    frac = 100 * np.mean(kvals == float(ks))
    print(f"  {ks: 2d} px ({ks * PIXEL_SCALE: .1f}\"): {frac: .1f}% of stamp pixels")

   1 px ( 0.2"):  13.8% of stamp pixels
   3 px ( 0.6"):  16.3% of stamp pixels
   7 px ( 1.4"):  21.2% of stamp pixels
   15 px ( 3.0"):  19.0% of stamp pixels
   31 px ( 6.2"):  16.2% of stamp pixels
   61 px ( 12.2"):  13.5% of stamp pixels

Display the original stamp alongside the kernel size map: bright compact sources receive small kernels; noise-dominated regions receive large kernels.

In [19]:

norm = sky_norm(stamp_sky, lo_sigma=2.0, hi_sigma=8.0)

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
axes[0].imshow(stamp_display, norm=norm, cmap="Greys_r", origin="lower")
axes[0].set_title("Original (bad pixels masked)")
axes[0].set_xlabel("x (pixels)")
axes[0].set_ylabel("y (pixels)")

im = axes[1].imshow(kernel_map, cmap="cividis",
                    origin="lower", vmin=0, vmax=61)
axes[1].set_title("Kernel size map (pixels)")
axes[1].set_xlabel("x (pixels)")
axes[1].set_ylabel("y (pixels)")
plt.colorbar(im, ax=axes[1], label="Kernel size (pixels)")

Out[19]:

<matplotlib.colorbar.Colorbar at 0x7b2bccc7f0e0>

Figure 4: Left -- original stamp with bad pixels masked. Right -- kernel size map in pixels; blue indicates small kernels applied to compact, well-detected sources; yellow indicates large kernels applied to noise-dominated and low-surface-brightness regions.

Compare a fixed 8x pixel binning (uniform across the entire stamp) with the adaptive result on the same stretch.

In [20]:

bin_factor_ref = 8
stamp_binned = bin_image(stamp_display, bin_factor_ref)
ny_full, nx_full = stamp_display.shape

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(
    stamp_binned,
    norm=norm,
    cmap="Greys_r",
    origin="lower",
    extent=[0, nx_full, 0, ny_full],
)
axes[0].set_title(
    f'Pixel binning {bin_factor_ref}x ({
        bin_factor_ref
        * PIXEL_SCALE: .1f}" res.)')
axes[0].set_xlabel("x (pixels, native scale)")
axes[0].set_ylabel("y (pixels, native scale)")

axes[1].imshow(smooth_result, norm=norm, cmap="Greys_r", origin="lower")
axes[1].set_title("Adaptive smoothing (variable resolution)")
axes[1].set_xlabel("x (pixels, native scale)")
axes[1].set_ylabel("y (pixels, native scale)")

Out[20]:

Text(0, 0.5, 'y (pixels, native scale)')

Figure 5: Pixel binning at 8× (left) versus adaptive smoothing (right) on the same stretch. Binning achieves uniform noise reduction at the cost of spatial resolution everywhere; adaptive smoothing retains compact source morphology while suppressing noise in the diffuse outskirts. Ring-like artifacts visible around bright sources in the adaptive-smoothed panel are a known side effect: the kernel size changes abruptly at source boundaries, creating a visible discontinuity between the finely-smoothed core and the more heavily smoothed surrounding halo.

4.2. Noise quantification¶

Apply adaptive_smooth to the full coadd patch rather than the stamp, so that the sky pixel sample is large enough for robust noise statistics. Print the before and after sky noise, then produce two diagnostic plots: a bar chart showing the fraction of stamp pixels assigned to each kernel size (indicating the typical angular scale at which the field becomes noise-dominated), and a histogram comparing the sky pixel distributions before and after smoothing (a narrower distribution after smoothing confirms genuine noise reduction).

In [21]:

kernel_sizes = [1, 3, 7, 15, 31, 61]
coadd_smooth, _ = adaptive_smooth(
    coadd_sky, sky_sigma, target_snr=3.0, kernel_sizes=kernel_sizes
)

sky_valid = np.isfinite(coadd_sky)
sky_raw = coadd_sky[sky_valid]
sky_smoothed = coadd_smooth[sky_valid]

_, _, sky_sigma_raw = sky_stats(sky_raw)
_, _, sky_sigma_smooth = sky_stats(sky_smoothed)
noise_reduction = sky_sigma_raw / sky_sigma_smooth
depth_gain_mag = (
    nJy_per_pix_to_mu(sky_sigma_smooth) - nJy_per_pix_to_mu(sky_sigma_raw)
)

mu_raw = nJy_per_pix_to_mu(sky_sigma_raw)
mu_smooth = nJy_per_pix_to_mu(sky_sigma_smooth)
gain_str = format(depth_gain_mag, '+.2f')
noise_equiv = noise_reduction**2
print("Sky noise in inter-source regions (full patch):")
print(f"  Before adaptive smoothing: {sky_sigma_raw: .1f} nJy/pix = {mu_raw: .1f} mag/arcsec^2")
smooth_str = f"{sky_sigma_smooth: .1f}"
mu_smooth_str = f"{mu_smooth: .1f}"
print(f"  After adaptive smoothing: {smooth_str} nJy/pix = {mu_smooth_str} mag/arcsec^2")
nr_str = f"{noise_reduction: .1f}"
ne_str = f"{noise_equiv: .0f}"
print(f"  Noise reduction: {nr_str}x (equivalent to binning ~{ne_str} pixels)")
print(f"  Surface brightness gain: {gain_str} mag/arcsec^2")

Sky noise in inter-source regions (full patch):
  Before adaptive smoothing:  3.4 nJy/pix =  26.6 mag/arcsec^2
  After adaptive smoothing:  0.5 nJy/pix =  28.7 mag/arcsec^2
  Noise reduction:  7.2x (equivalent to binning ~ 52 pixels)
  Surface brightness gain: +2.14 mag/arcsec^2

Plot the kernel size distribution from the stamp and the sky pixel histograms before and after smoothing.

In [22]:

kvals = kernel_map[np.isfinite(kernel_map)].ravel()
counts = [np.sum(kvals == float(k)) for k in kernel_sizes]
fracs = [c / len(kvals) for c in counts]

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].bar(
    range(len(kernel_sizes)),
    [100 * f for f in fracs],
    tick_label=[f'{k} px\n({k * PIXEL_SCALE: .1f}")' for k in kernel_sizes],
    color=["gray", YELLOW, ORANGE, MAGENTA, PURPLE, BLUE],
)
axes[0].set_xlabel("Kernel size selected")
axes[0].set_ylabel("Fraction of pixels (%)")
axes[0].set_title("Kernel size distribution")
for i, f in enumerate(fracs):
    axes[0].text(i, 100 * f + 0.5, f"{100 * f: .1f}%", ha="center", fontsize=8)

x_range = 4.0 * sky_sigma
axes[1].hist(
    sky_raw,
    bins=100,
    range=(-x_range, x_range),
    density=True,
    histtype="step",
    color="gray",
    lw=1.5,
    label="Unsmoothed (source-masked)",
)
axes[1].hist(
    sky_smoothed,
    bins=100,
    range=(-x_range, x_range),
    density=True,
    histtype="step",
    color=ORANGE,
    lw=1.5,
    label="Adaptive smooth (source-masked)",
)
axes[1].axvline(0, color="k", ls=":", lw=0.8)
axes[1].set_xlabel("Pixel value (nJy)")
axes[1].set_ylabel("Probability density")
axes[1].set_title("Sky pixel distributions before/after adaptive smoothing")
axes[1].legend(fontsize=8)

Out[22]:

<matplotlib.legend.Legend at 0x7b2bb43b1810>

Figure 6: Left -- fraction of stamp pixels assigned to each smoothing kernel; the dominant kernel indicates the typical angular scale at which the field becomes noise-dominated. Right -- pixel value distributions from the full patch, restricted to source-masked pixels (i.e., locations with no detected sources), before (grey) and after (orange) adaptive smoothing. These distributions describe the noise: a narrower peak after smoothing confirms genuine noise reduction.

5. Signal-to-noise profile¶

A radial SNR profile shows how the azimuthal-median galaxy flux compares to the residual noise level of each technique as a function of radius. Where SNR approaches or falls below 1, the surface brightness in that annulus is comparable to the noise floor. The signal-to-noise ratio in each annulus is the azimuthal median flux divided by the effective per-pixel noise $\sigma_\mathrm{eff}$.

For pixel binning at factor $N$, averaging $N^2$ independent pixels reduces noise by $N$, so $\sigma_\mathrm{eff} = \hat{\sigma} / N$. For adaptive smoothing, the kernel size $k$ varies pixel by pixel, giving a spatially varying $\sigma_\mathrm{eff}(i,j) = \hat{\sigma} / k_{i,j}$, where $(i,j)$ are pixel coordinates and $k_{i,j}$ is the value of kernel_map at that pixel.

Compute the galaxy center position in the stamp, build circular annuli out to 50 arcsec, and measure the median flux in each. For binned methods, $\sigma_\mathrm{eff} = \hat{\sigma} / N$ is uniform across the annulus. For adaptive smoothing, the per-pixel SNR is computed using the local $\sigma_\mathrm{eff}(i,j) = \hat{\sigma} / k_{i,j}$, and the annulus value is the median of those per-pixel SNRs.

In [23]:

stamp_ny, stamp_nx = stamp_arr.shape
center_x = int(host_pix.getX()) - stamp_bbox.x.start
center_y = int(host_pix.getY()) - stamp_bbox.y.start
print(f"Host position in stamp: ({center_x}, {center_y})")

if not (5 < center_x < stamp_nx - 5 and 5 < center_y < stamp_ny - 5):
    print("WARNING: host near edge -- using stamp center as fallback.")
    center_x, center_y = stamp_nx // 2, stamp_ny // 2

grid_y, grid_x = np.mgrid[:stamp_ny, :stamp_nx]
r_arcsec = np.sqrt((grid_x - center_x) ** 2
                   + (grid_y - center_y) ** 2) * PIXEL_SCALE

r_max = 50.0
r_bins = np.linspace(0, r_max, 40)
r_centers = 0.5 * (r_bins[:-1] + r_bins[1:])

Host position in stamp: (600, 600)

expand_binned tiles each block-averaged pixel back to its original bf x bf footprint, returning an array at native resolution so binned and unbinned images can be compared pixel-for-pixel in the same annuli.

In [24]:

def expand_binned(arr_binned, bf):
    """Repeat block-averaged values back to native pixel resolution.

    Parameters
    ----------
    arr_binned : ndarray, shape (ny // bf, nx // bf)
        Block-averaged array from bin_image.
    bf : int
        Bin factor used to create arr_binned.

    Returns
    -------
    out : ndarray, shape (stamp_ny, stamp_nx)
        Array expanded to native pixel resolution by pixel repetition.
    """
    expanded = np.repeat(np.repeat(arr_binned, bf, axis=0), bf, axis=1)
    out = np.full((stamp_ny, stamp_nx), np.nan)
    h = min(expanded.shape[0], stamp_ny)
    w = min(expanded.shape[1], stamp_nx)
    out[:h, :w] = expanded[:h, :w]
    return out

snr_profile steps through each radial annulus and divides the median flux by a scalar sigma_eff, giving the SNR assuming uniform noise. Used for the unbinned image and each fixed bin factor.

In [25]:

def snr_profile(arr, sigma_eff):
    """Median flux in each annulus divided by scalar sigma_eff.

    Parameters
    ----------
    arr : ndarray, shape (stamp_ny, stamp_nx)
        Pixel array at native or expanded-binned resolution.
    sigma_eff : float
        Effective noise per pixel in the same units as arr.

    Returns
    -------
    prof : ndarray, shape (n_bins,)
        Median SNR in each radial annulus; NaN where fewer than 10 pixels.
    """
    prof = []
    for r_lo, r_hi in zip(r_bins[:-1], r_bins[1:]):
        ann = (r_arcsec >= r_lo) & (
            r_arcsec < r_hi) & ~stamp_bad & np.isfinite(arr)
        prof.append(np.median(arr[ann])
                    / sigma_eff if ann.sum() > 10 else np.nan)
    return np.array(prof)

snr_profile_adaptive uses the spatially varying effective noise $\hat{\sigma} / k_{i,j}$ at each pixel rather than a scalar, then takes the annulus median of those per-pixel SNRs.

In [26]:

def snr_profile_adaptive(smooth, kmap, sigma):
    """Median per-pixel SNR in each annulus using spatially varying sigma = sigma / k.

    Parameters
    ----------
    smooth : ndarray, shape (stamp_ny, stamp_nx)
        Adaptively smoothed image from adaptive_smooth.
    kmap : ndarray, shape (stamp_ny, stamp_nx)
        Kernel-size map from adaptive_smooth.
    sigma : float
        Per-pixel sky noise before smoothing, in nJy pixel**-1.

    Returns
    -------
    prof : ndarray, shape (n_bins,)
        Median SNR in each radial annulus; NaN where fewer than 10 pixels.
    """
    prof = []
    for r_lo, r_hi in zip(r_bins[:-1], r_bins[1:]):
        ann = (
            (r_arcsec >= r_lo)
            & (r_arcsec < r_hi)
            & ~stamp_bad
            & np.isfinite(smooth)
            & np.isfinite(kmap)
        )
        if ann.sum() > 10:
            prof.append(np.median(smooth[ann] * kmap[ann] / sigma))
        else:
            prof.append(np.nan)
    return np.array(prof)

Compute the radial SNR profile for the unbinned image, three bin factors, and adaptive smoothing.

Warning: This cell may produce a RuntimeWarning: Mean of empty slice from numpy. This is expected: bin blocks that fall entirely within masked regions (detector edges, saturated columns) contain only NaN pixels, so np.nanmean has nothing to average and returns NaN for those blocks. The warning can be safely ignored — the output array is correct.

In [27]:

snr_unbinned = snr_profile(stamp_arr, stamp_sigma)
snr_bin4 = snr_profile(expand_binned(
    bin_image(stamp_display, 4), 4), stamp_sigma / 4)
snr_bin8 = snr_profile(expand_binned(
    bin_image(stamp_display, 8), 8), stamp_sigma / 8)
snr_bin16 = snr_profile(
    expand_binned(bin_image(stamp_display, 16), 16), stamp_sigma / 16
)
snr_adaptive = snr_profile_adaptive(smooth_result, kernel_map, stamp_sigma)

Plot the radial SNR profiles.

In [28]:

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(r_centers, snr_unbinned, "-", lw=1.5,
        color="gray", alpha=0.8, label="Unbinned")
ax.plot(r_centers, snr_bin4, "-", lw=1.5, color=YELLOW, label="Bin 4x")
ax.plot(r_centers, snr_bin8, "-", lw=1.5, color=ORANGE, label="Bin 8x")
ax.plot(r_centers, snr_bin16, "-", lw=1.5, color=MAGENTA, label="Bin 16x")
ax.plot(r_centers, snr_adaptive, "-", lw=2.0,
        color=BLUE, label="Adaptive smooth")
ax.set_xlim(0, r_max)
ax.set_yscale("log")
ax.set_xlabel("Radius (arcsec)")
ax.set_ylabel(r"SNR (median flux / $\sigma_{\rm eff}$)")
ax.legend(fontsize=9)

Out[28]:

<matplotlib.legend.Legend at 0x7b2bb4278190>

Figure 7: Radial SNR profile (azimuthal median flux divided by the effective per-pixel noise $\sigma_{\rm eff}$) for five noise-reduction methods, out to 50 arcsec. Fixed binning at $16\times$ (magenta) achieves far higher SNR than necessary in the bright inner regions, sacrificing spatial resolution everywhere for a gain that is redundant where the signal is already strong. Adaptive smoothing (blue) instead applies only the smallest kernel that reaches the target SNR of 3 at each pixel, tracking a noise floor rather than exceeding it: resolution is preserved in bright regions and smoothing is applied only in the faint outskirts where it is genuinely needed.

6. Summary¶

• Pixel binning reduces noise by $\sqrt{N}$ while exactly preserving surface brightness. At large bin factors (beyond 8–16x) correlated background residuals cause the measured $\sigma$ to exceed the $\sqrt{N}$ prediction.
• Adaptive smoothing via normalized convolution selects the finest kernel that meets a target SNR at each pixel, combining resolution preservation in bright regions with aggressive noise suppression in diffuse outskirts.
• A radial SNR profile, dividing the azimuthal median by the effective noise of each method, shows directly at what radius the galaxy signal drops below detection — and how far each technique extends that limit.