Annulus background

This notebook estimates a local background from a circular annulus and subtracts it from exact aperture sums.

import numpy as np
import matplotlib.pyplot as plt
import astroapers as aap

Synthetic image with a sloped background

ny, nx = 60, 60
y, x = np.mgrid[:ny, :nx]

x0, y0 = 34.7, 33.9
background = 20.0 + 0.03 * x - 0.02 * y
source = 10.0 * np.exp(-0.5 * ((x - x0) ** 2 + (y - y0) ** 2) / 3.2**2)
data = background + source

Source aperture and sky annulus

Use exact fractional weights for the source sum, and center-selected annulus data values for robust background samples. The center weights define which pixels are selected; sampled_values(data) returns the corresponding raw image values.

ap = aap.CircAp((x0, y0), r=5.0)
an = aap.CircAn((x0, y0), r_in=9.0, r_out=15.0)

apsum, src_npix = ap.apsum_exact(data)
sky_values = an.sampled_values(data)[0]
# you can do complicated math to get best sky value.
# Here I just use median:
sky_median = np.median(sky_values)

# sky subtraction
srcsum = apsum - src_npix * sky_median

# apsum, src_npix, sky_values.size, sky_median, srcsum
print(f"Aperture sum: {apsum:.3f}")
print(f"Sky median: {sky_median:.3f}")
print(f"Source npix: {src_npix:.3f}")
print(f"Source sum: {srcsum:.3f} = apsum - npix * sky")

Aperture sum: 2050.868
Sky median: 20.385
Source npix: 78.540
Source sum: 449.827 = apsum - npix * sky

Visual check

fig, ax = plt.subplots(1, 2, figsize=(7, 4))
ax[0].imshow(data, origin="lower", cmap="viridis")
ap.plot(ax=ax[0], color="r", lw=1)
an.plot(ax=ax[0], color="w", lw=1)
ax[0].set_title("Source aperture and sky annulus")

ax[1].hist(sky_values, bins=20, color="gray", edgecolor="k")
ax[1].set_title("Sky values")
ax[1].axvline(sky_median, color="r", ls="--", label=f"Median={sky_median:.2f}")
ax[1].legend()

plt.show()

A masked-pixel example

Bad-pixel masks are full image-shaped boolean arrays where True means excluded. For the sky sample, pass the same mask to sampled_values() before taking the median.

bad_pixels = np.zeros(data.shape, dtype=bool)
bad_pixels[40:45, 45:52] = True

sky_values_masked = an.sampled_values(data, mask=bad_pixels)[0]
sky_median_masked = np.median(sky_values_masked)
apsum_masked, src_npix_masked = ap.apsum_exact(data, mask=bad_pixels)
srcsum_masked = apsum_masked - src_npix_masked * sky_median_masked

print(f"Unmasked sky pixels: {sky_values.size}")
print(f"Masked sky pixels: {sky_values_masked.size}")
print(f"Unmasked sky median: {sky_median:.3f}")
print(f"Masked sky median: {sky_median_masked:.3f}")
print(f"Masked source npix: {src_npix_masked:.3f}")
print(f"Masked source sum: {srcsum_masked:.3f}")

Unmasked sky pixels: 451
Masked sky pixels: 438
Unmasked sky median: 20.385
Masked sky median: 20.372
Masked source npix: 78.540
Masked source sum: 450.871

fig, ax = plt.subplots(figsize=(5, 4))
ax.imshow(data, origin="lower", cmap="viridis")
ax.imshow(bad_pixels, origin="lower", cmap="gray_r", alpha=0.4)
ap.plot(ax=ax, color="r", lw=1)
an.plot(ax=ax, color="w", lw=1)
ax.set_title("Bad pixels excluded from aperture and annulus")
plt.show()

Maximum Performance

For fixed unmasked work, import the raw Rust extension as aapr. Raw Rust functions do not accept the bad_pixels mask; if masked background subtraction is required, keep using the object or Python wrapper layer. For more raw-call patterns, inspect astroapers.kernels.

import astroapers._rust as aapr

raw_data = np.ascontiguousarray(data, dtype=np.float64)
xpos = np.ascontiguousarray([x0], dtype=np.float64)
ypos = np.ascontiguousarray([y0], dtype=np.float64)

raw_src_apsum, raw_src_npix = aapr.apsum_circ_exact(
    raw_data,
    xpos,
    ypos,
    5.0,
)
raw_sky_weights, ixmins, ixmaxs, iymins, iymaxs = aapr.weights_circ_ann_center(
    xpos,
    ypos,
    9.0,
    15.0,
)

sky_shape = (iymaxs[0] - iymins[0], ixmaxs[0] - ixmins[0])
sky_weights0 = np.asarray(raw_sky_weights[0], dtype=np.float64).reshape(sky_shape)
sky_cutout = raw_data[iymins[0]:iymaxs[0], ixmins[0]:ixmaxs[0]]
raw_sky_values = sky_cutout[sky_weights0 > 0]
raw_sky_median = np.median(raw_sky_values)
raw_srcsum = raw_src_apsum[0] - raw_src_npix[0] * raw_sky_median

assert np.allclose(raw_src_apsum[0], apsum)
assert np.allclose(raw_src_npix[0], src_npix)
assert np.allclose(raw_sky_values, sky_values)
assert np.allclose(raw_sky_median, sky_median)

raw_sky_median, raw_srcsum

(np.float64(20.385090278577795), np.float64(449.8267084524746))