"""
Background estimation routines for astronomical images.
This module provides various methods for estimating and subtracting background
from astronomical images, including grid-based methods (SEP, photutils) and
non-grid methods (percentile filtering, morphological opening).
"""
import numpy as np
from astropy.stats import mad_std
[docs]
def get_background_percentile(
image, mask=None, size=51, percentile=25.0, maxiters=3, sigma=3.0, use_mode=False
):
"""
Background estimation via percentile filtering (non-grid method).
Uses a sliding window to compute low percentile value at each pixel.
Stars (positive outliers) are naturally above the percentile and thus rejected.
This is a truly non-grid method: operates on all pixels with sliding window,
not limited by grid cell resolution.
Parameters
----------
image : ndarray
Input image
mask : ndarray, optional
Boolean mask (True = masked pixels to exclude)
size : int
Filter kernel size (should be larger than typical star size, e.g., 51-101 pixels).
Larger values are more robust but may miss fine background structure.
percentile : float
Percentile value to use (0-100). Typical values:
- 25: 25th percentile (default, good for moderate star density)
- 10: More aggressive star rejection (for crowded fields)
- 50: Median (less star rejection but more robust to noise)
maxiters : int
Number of sigma-clipping iterations to refine percentile estimate
sigma : float
Sigma threshold for clipping
use_mode : bool
If True, estimate mode instead of percentile (experimental)
Returns
-------
background : ndarray
Estimated background map
Notes
-----
- Computationally intensive for large kernels (scales as O(N * kernel^2))
- Consider using smaller kernels (31-51) for large images
- For very crowded fields, use lower percentiles (10-15)
- Handles arbitrary background shapes (not limited to polynomial forms)
"""
from scipy.ndimage import percentile_filter, median_filter
# Work with a copy to avoid modifying input
img = image.copy()
# Apply mask by replacing with NaN (will be ignored in percentile computation)
has_nan = False
if mask is not None:
img[mask] = np.nan
has_nan = True
# Fix #3: Define masked_percentile closure once before the loop
def masked_percentile(values):
valid = values[~np.isnan(values)]
if len(valid) < 3:
return np.nan
return np.percentile(valid, percentile)
# Iterative sigma-clipping approach
for iteration in range(maxiters):
if use_mode:
# Mode estimation via histogram (in sliding window)
# This is expensive, so use median as approximation
back = median_filter(img, size=size)
else:
# Fix #2 & #4: Use fast percentile_filter with NaN-filling approach
# instead of slow generic_filter with Python callback
if has_nan:
# Two-pass approach: fill NaN with local median, then run fast filter
# Pass 1: Fill NaN pixels with local median
img_filled = img.copy()
nan_mask = np.isnan(img_filled)
if np.any(nan_mask):
# Use median filter to get fill values
# Temporarily replace NaN with image median for the filter
fill_val = np.nanmedian(img_filled)
img_filled[nan_mask] = fill_val
local_med = median_filter(img_filled, size=min(size, 15))
img_filled[nan_mask] = local_med[nan_mask]
# Pass 2: Run fast C-based percentile_filter on the filled image
back = percentile_filter(img_filled, percentile=percentile, size=size)
else:
back = percentile_filter(img, percentile=percentile, size=size)
# Sigma-clip for next iteration
if iteration < maxiters - 1:
residual = img - back
# Estimate RMS from negative residuals (not contaminated by stars)
neg_residuals = residual[residual < 0]
if len(neg_residuals) > 100:
rms = np.std(neg_residuals) * np.sqrt(2) # Correct for one-sided
else:
rms = mad_std(residual[np.isfinite(residual)])
# Clip positive outliers (stars)
outliers = residual > sigma * rms
img[outliers] = np.nan
if np.any(outliers):
has_nan = True
# Fill any remaining NaN values
if np.any(np.isnan(back)):
valid = np.isfinite(back)
if np.sum(valid) > 0:
# Use median of valid pixels
fill_value = np.median(back[valid])
back[~valid] = fill_value
return back
[docs]
def get_background_morphology(
image, mask=None, size=25, iterations=1, smooth=True, smooth_size=None
):
"""
Background estimation via morphological opening (non-grid method).
Performs morphological opening (erosion followed by dilation) which removes
structures smaller than the structuring element (stars). Fast and simple.
This is a truly non-grid method: operates on all pixels with sliding window.
Parameters
----------
image : ndarray
Input image
mask : ndarray, optional
Boolean mask (True = masked pixels to exclude)
size : int
Size of structuring element (should be larger than typical star size).
Typical values: 15-35 pixels depending on seeing.
iterations : int
Number of opening iterations (default: 1, more iterations = more aggressive)
smooth : bool
If True, apply additional smoothing to reduce kernel-sized artifacts.
This significantly improves the smoothness of the result.
smooth_size : int, optional
Size of smoothing kernel in pixels. If None, automatically set based on
the structuring element size. Use larger values for smoother results.
Returns
-------
background : ndarray
Estimated background map
Notes
-----
- Very fast (much faster than percentile filtering)
- Simple and robust
- Can create kernel-sized artifacts without smoothing
- May suppress legitimate background features similar in size to stars
- Best for images with clear point sources on smooth backgrounds
- Use smooth=True (default) to eliminate spotty kernel artifacts
"""
from scipy.ndimage import grey_opening, median_filter, gaussian_filter
# Work with a copy
img = image.copy()
# Handle mask by replacing with local median
if mask is not None:
# Fill masked regions with local median to avoid artifacts
if np.any(mask):
# Simple approach: use median filter to fill masked regions
img_filled = median_filter(img, size=5)
img[mask] = img_filled[mask]
# Morphological opening with circular structuring element
# Create circular structuring element
y, x = np.ogrid[-size // 2 : size // 2 + 1, -size // 2 : size // 2 + 1]
structure = (x * x + y * y) <= (size / 2) ** 2
# Apply opening
back = grey_opening(img, structure=structure)
# Additional iterations if requested
for i in range(iterations - 1):
back = grey_opening(back, structure=structure)
# Apply smoothing to reduce kernel-sized artifacts
# Fix #6: smooth_size consistently means filter size in pixels for both paths
if smooth:
if mask is not None and np.any(mask):
# Use median filter for masked images (more robust to artifacts)
if smooth_size is None:
filter_size = max(5, size // 5)
else:
filter_size = int(smooth_size)
back = median_filter(back, size=filter_size)
else:
# Use Gaussian smoothing for clean images (smoother results)
if smooth_size is None:
# Default: adaptive smoothing based on kernel size
sigma_val = max(2.0, size / 5.0)
else:
# Convert filter size to sigma (size ~ 3*sigma is a reasonable heuristic)
sigma_val = smooth_size / 3.0
back = gaussian_filter(back, sigma=sigma_val, mode='reflect')
return back
[docs]
def estimate_background_rms_percentile(image, mask=None, size=51):
"""
Estimate background RMS using percentile range.
Uses interquartile range (IQR) as robust RMS estimator.
Parameters
----------
image : ndarray
Input image
mask : ndarray, optional
Boolean mask (True = masked pixels)
size : int
Kernel size for local percentile estimation
Returns
-------
rms : ndarray
Background RMS map
"""
from scipy.ndimage import percentile_filter, median_filter
# Fix #4: Replace slow generic_filter with fast two-pass approach
if mask is not None:
img = image.copy()
img[mask] = np.nan
# Two-pass: fill NaN with local median, then run fast percentile_filter
nan_mask = np.isnan(img)
if np.any(nan_mask):
fill_val = np.nanmedian(img)
img_filled = img.copy()
img_filled[nan_mask] = fill_val
local_med = median_filter(img_filled, size=min(size, 15))
img_filled[nan_mask] = local_med[nan_mask]
else:
img_filled = img
# Fix #9: Always use p25/p75 for IQR regardless of caller's percentile
p25 = percentile_filter(img_filled, percentile=25, size=size)
p75 = percentile_filter(img_filled, percentile=75, size=size)
else:
p25 = percentile_filter(image, percentile=25, size=size)
p75 = percentile_filter(image, percentile=75, size=size)
# IQR / 1.349 ≈ sigma for Gaussian
rms = (p75 - p25) / 1.349
return rms
[docs]
def estimate_background_rms_local(image, background, mask=None, size=51):
"""
Estimate background RMS from residuals.
Computes local standard deviation of (image - background).
Parameters
----------
image : ndarray
Input image
background : ndarray
Estimated background
mask : ndarray, optional
Boolean mask (True = masked pixels)
size : int
Kernel size for local RMS estimation
Returns
-------
rms : ndarray
Background RMS map
"""
from scipy.ndimage import uniform_filter
residual = image - background
if mask is not None:
residual = residual.copy()
residual[mask] = np.nan
# Fix #4: Replace slow generic_filter with fast uniform_filter-based local std
# For masked data, fill NaN first then use fast filters
has_nan = np.any(np.isnan(residual))
if has_nan:
nan_mask = np.isnan(residual)
res_filled = residual.copy()
fill_val = np.nanmedian(res_filled)
res_filled[nan_mask] = fill_val
# Local mean and mean of squares via uniform_filter
local_mean = uniform_filter(res_filled, size=size, mode='reflect')
local_mean_sq = uniform_filter(res_filled**2, size=size, mode='reflect')
else:
local_mean = uniform_filter(residual, size=size, mode='reflect')
local_mean_sq = uniform_filter(residual**2, size=size, mode='reflect')
# Var = E[X^2] - E[X]^2
local_var = local_mean_sq - local_mean**2
# Clamp negative values from numerical precision
local_var = np.maximum(local_var, 0.0)
rms = np.sqrt(local_var)
# Fill NaN values
if np.any(np.isnan(rms)):
valid = np.isfinite(rms)
if np.sum(valid) > 0:
fill_value = np.median(rms[valid])
rms[~valid] = fill_value
return rms
[docs]
def get_background_polynomial(
image,
mask=None,
order=1,
sigma=3.0,
maxiters=5,
directional=False,
direction=None,
):
"""Estimate background as a low-order 2D polynomial.
Fits a smooth polynomial surface to the image using iterative
sigma-clipping to reject stars and other outliers. By construction,
low-order polynomials cannot represent localised features such as
nebulae, so only the smooth gradient component is captured.
Parameters
----------
image : 2D ndarray
Input image.
mask : 2D bool ndarray, optional
Pixels to exclude (True = masked).
order : int
Polynomial order: 1 = plane, 2 = quadratic, 3 = cubic, etc.
Default 1.
sigma : float
Sigma threshold for iterative outlier clipping. Default 3.0.
maxiters : int
Maximum number of sigma-clipping iterations. Default 5.
directional : bool or float
Controls directional constraint for higher-order (≥ 2) terms.
- ``False`` (default) — use all 2D polynomial terms.
- ``True`` — auto-detect the gradient direction from a linear
fit, then restrict order ≥ 2 terms to powers of the
directional coordinate *t* = *x* cos θ + *y* sin θ.
This prevents high-order polynomials from fitting localised
features while still capturing curvature along the gradient.
- A ``float`` value — same as ``True`` but uses the given angle
(degrees, 0 = x-axis, 90 = y-axis) instead of auto-detecting.
For order ≤ 1 this parameter has no effect (a plane already
has only 3 coefficients).
direction : float or None
*Deprecated* synonym for ``directional=<angle>``. If both are
given, *directional* takes precedence.
Returns
-------
background : 2D ndarray
Fitted polynomial background, same shape as *image*.
rms : float
Scalar RMS of the residuals (clipped pixels excluded).
Notes
-----
**Full 2D mode** (``directional=False``):
Number of polynomial coefficients = (order+1)(order+2)/2:
order 1 → 3, order 2 → 6, order 3 → 10, order 4 → 15.
**Directional mode** (``directional=True`` or angle):
Number of coefficients = order + 2 (full linear + higher-order
directional terms): order 2 → 4, order 3 → 5, order 4 → 6.
Coordinates are normalised to [-1, 1] internally for numerical
stability. The fit uses random sub-sampling (100 000 points) for
images larger than ~300×300 to keep memory and CPU bounded.
Examples
--------
Remove a linear gradient:
>>> bg, rms = get_background_polynomial(image, order=1)
Fit a cubic curve along the gradient direction (5 parameters):
>>> bg, rms = get_background_polynomial(image, order=3, directional=True)
Specify the gradient direction explicitly (e.g. 45°):
>>> bg, rms = get_background_polynomial(image, order=3, directional=45.0)
"""
ny, nx = image.shape
# --- resolve directional / direction arguments -----------------------
if isinstance(directional, (int, float)) and not isinstance(directional, bool):
# directional=45.0 → direction angle in degrees
direction_deg = float(directional)
directional = True
elif directional and direction is not None:
direction_deg = float(direction)
elif direction is not None:
direction_deg = float(direction)
directional = True
else:
direction_deg = None # auto-detect later if directional=True
# Directional constraint only matters for order >= 2
if order < 2:
directional = False
# --- build a fixed-size working sample -------------------------------
# All fitting and sigma-clipping operates on this subsample only;
# the full-grid evaluation is a single pass at the end.
MAX_SAMPLE = 100_000
data_flat = image.ravel()
valid_idx = np.where(np.isfinite(data_flat))[0]
if mask is not None:
valid_idx = valid_idx[~mask.ravel()[valid_idx]]
if len(valid_idx) == 0:
return np.full_like(image, np.nan, dtype=np.float64), np.nan
# Subsample once — reused for every clipping iteration
if len(valid_idx) > MAX_SAMPLE:
rng = np.random.RandomState(0)
sample_idx = rng.choice(valid_idx, MAX_SAMPLE, replace=False)
else:
sample_idx = valid_idx
# Coordinates and data for the working sample (normalised to [-1, 1])
x_scale = max(nx / 2, 1)
y_scale = max(ny / 2, 1)
sample_y, sample_x = np.divmod(sample_idx, nx)
xs = (sample_x.astype(np.float64) - nx / 2) / x_scale
ys = (sample_y.astype(np.float64) - ny / 2) / y_scale
ds = data_flat[sample_idx].astype(np.float64)
# --- helpers ---------------------------------------------------------
def _poly_terms(ord_):
"""Return list of (p, q) exponent pairs for full 2D polynomial."""
return [(p, q) for p in range(ord_ + 1) for q in range(ord_ + 1 - p)]
terms = _poly_terms(order)
def _build_basis(x, y, theta_rad):
"""Build design matrix for the current model."""
if directional:
t = x * np.cos(theta_rad) + y * np.sin(theta_rad)
cols = [np.ones_like(x), x, y]
tk = t.copy()
for _ in range(2, order + 1):
tk = tk * t
cols.append(tk)
return np.column_stack(cols)
cols = []
for p, q in terms:
col = np.ones_like(x)
if p:
col = col * x**p
if q:
col = col * y**q
cols.append(col)
return np.column_stack(cols)
def _eval_2d(coeffs, theta_rad):
"""Evaluate polynomial on the full pixel grid via broadcasting."""
x_1d = (np.arange(nx, dtype=np.float64) - nx / 2) / x_scale
y_1d = (np.arange(ny, dtype=np.float64) - ny / 2) / y_scale
bg = np.full((ny, nx), coeffs[0], dtype=np.float64)
if directional:
bg += coeffs[1] * x_1d[np.newaxis, :]
bg += coeffs[2] * y_1d[:, np.newaxis]
t = np.cos(theta_rad) * x_1d[np.newaxis, :] + np.sin(theta_rad) * y_1d[:, np.newaxis]
tk = t.copy()
for k in range(2, order + 1):
tk = tk * t
bg += coeffs[k + 1] * tk
else:
idx = 0
for p, q in terms:
if idx == 0:
idx += 1
continue # constant already handled
term = np.ones((ny, nx), dtype=np.float64)
if p:
term *= x_1d[np.newaxis, :] ** p
if q:
term *= y_1d[:, np.newaxis] ** q
bg += coeffs[idx] * term
idx += 1
return bg
# --- auto-detect gradient direction ----------------------------------
theta_rad = 0.0
if directional and direction_deg is None:
A = np.column_stack([np.ones_like(xs), xs, ys])
c, _, _, _ = np.linalg.lstsq(A, ds, rcond=None)
theta_rad = np.arctan2(c[2], c[1])
elif direction_deg is not None:
theta_rad = np.radians(direction_deg)
# --- iterative sigma-clipped polynomial fit --------------------------
# All iterations work on the fixed subsample (xs, ys, ds).
# Note: some BLAS implementations (e.g. Apple Accelerate) raise
# spurious divide-by-zero warnings during matmul/lstsq even when
# all inputs and outputs are finite.
good = np.ones(len(xs), dtype=bool)
coeffs = None
rms_val = np.nan
n_params = (order + 2) if directional else len(terms)
with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
for iteration in range(maxiters + 1):
n_good = int(np.sum(good))
if n_good < n_params + 1:
break
A = _build_basis(xs[good], ys[good], theta_rad)
coeffs, _, _, _ = np.linalg.lstsq(A, ds[good], rcond=None)
if iteration == maxiters:
break
# Evaluate model on the subsample for clipping
model = _build_basis(xs, ys, theta_rad) @ coeffs
res = ds - model
# Robust RMS via MAD
med_res = np.median(res)
rms_val = 1.4826 * np.median(np.abs(res - med_res))
if rms_val <= 0:
break
good = np.abs(res - med_res) < sigma * rms_val
# --- final evaluation on the full grid (single pass) -----------------
if coeffs is None:
return np.full_like(image, np.nan, dtype=np.float64), np.nan
background = _eval_2d(coeffs, theta_rad)
# Final RMS from clipped subsample residuals
if np.any(good) and np.all(np.isfinite(coeffs)):
with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
model_good = _build_basis(xs[good], ys[good], theta_rad) @ coeffs
rms_val = float(mad_std(ds[good] - model_good))
else:
rms_val = np.nan
return background, rms_val
[docs]
def get_background_gp(
image,
mask=None,
max_points=3000,
grid_step=None,
clip_sigma=3.5,
n_clip_iter=2,
length_scale=64.0,
length_scale_bounds=(8.0, 512.0),
matern_nu=1.5,
white_noise=1.0,
white_noise_bounds=(1e-3, 1e3),
normalize_y=True,
n_restarts_optimizer=1,
random_state=0,
get_uncertainty=False,
chunk_rows=256,
):
"""
Background estimation via Gaussian Process regression (non-grid method).
Uses a Gaussian Process with a Matern kernel to model the background as a
smooth function of spatial coordinates (x, y). This method is very flexible
and can handle complex background variations including gradients, vignetting,
and other large-scale structures.
Workflow:
1. Select candidate background pixels (masked + sigma clipping to remove stars)
2. Subsample them randomly or on a grid (to keep computation tractable)
3. Fit a GP to the (x, y) -> background mapping
4. Predict background map everywhere (with optional uncertainty)
Parameters
----------
image : ndarray
Input image
mask : ndarray, optional
Boolean mask (True = masked pixels to exclude)
max_points : int
Maximum number of training points for the GP. More points = better fit
but much slower (O(N^3) cost). Typical values: 1000-5000.
grid_step : int, optional
If set, sample every N pixels on a grid (after masking). If None,
random sampling is used.
clip_sigma : float
Sigma threshold for iterative clipping of outliers (stars) in intensity space
n_clip_iter : int
Number of sigma-clipping iterations to remove stars
length_scale : float
Initial guess for GP length scale in pixels. Should be larger than
typical star size. Typical values: 32-128 pixels.
length_scale_bounds : tuple
Bounds for length scale optimization (min, max)
matern_nu : float
Smoothness parameter for Matern kernel. Typical values:
- 0.5: Exponential kernel (rough)
- 1.5: Once differentiable (default, good balance)
- 2.5: Twice differentiable (very smooth)
white_noise : float
Initial guess for white noise level (in image units, e.g., ADU)
white_noise_bounds : tuple
Bounds for white noise optimization (min, max)
normalize_y : bool
If True, normalize target values during GP training (recommended)
n_restarts_optimizer : int
Number of random restarts for hyperparameter optimization (more = better
but slower)
random_state : int
Random seed for reproducibility
get_uncertainty : bool
If True, return uncertainty map. If False, return None for RMS.
chunk_rows : int
Process prediction in chunks of this many rows to reduce memory usage
Returns
-------
background : ndarray
Estimated background map
uncertainty : ndarray or None
If get_uncertainty=True: 2D array of GP posterior standard deviation
If get_uncertainty=False: None
Notes
-----
- Requires scikit-learn to be installed
- Computationally expensive for large images (use smaller max_points)
- Very flexible: can handle arbitrary background shapes
- Length scale should be >> PSF size (otherwise may fit stars)
- For very large images, consider using grid_step to reduce training points
Examples
--------
Basic usage with default parameters:
>>> bg = get_background_gp(image, mask=mask)
Custom length scale for wide-field imaging:
>>> bg = get_background_gp(image, mask=mask, length_scale=128.0)
Get uncertainty map:
>>> bg, bg_std = get_background_gp(image, mask=mask, get_uncertainty=True)
Faster computation with grid sampling:
>>> bg = get_background_gp(image, mask=mask, max_points=1000, grid_step=10)
"""
try:
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, Matern, WhiteKernel
except ImportError:
raise ImportError(
"Gaussian Process background estimation requires scikit-learn. "
"Install with: pip install scikit-learn"
)
# Prepare image and mask
img = np.asarray(image, dtype=float)
if img.ndim != 2:
raise ValueError("image must be 2D")
ny, nx = img.shape
if mask is None:
mask = np.zeros(img.shape, dtype=bool)
else:
mask = np.asarray(mask, dtype=bool)
if mask.shape != img.shape:
raise ValueError(f"mask shape {mask.shape} != image shape {img.shape}")
# Step 1: Select candidate background pixels via sigma clipping
valid = ~mask & np.isfinite(img)
y_vals = img[valid]
if y_vals.size < 50:
raise RuntimeError("Too few valid pixels to estimate background")
# Iterative sigma clipping to remove stars
# Fix #7: Use mad_std directly instead of inner robust_sigma function
keep = np.ones(y_vals.shape, dtype=bool)
for _ in range(max(0, int(n_clip_iter))):
yy = y_vals[keep]
if yy.size < 50:
break
med = np.median(yy)
sig = mad_std(yy, ignore_nan=True)
if not np.isfinite(sig) or sig <= 0:
break
lo = med - clip_sigma * sig
hi = med + clip_sigma * sig
keep = keep & (y_vals >= lo) & (y_vals <= hi)
# Get coordinates for kept pixels
ys_idx, xs_idx = np.nonzero(valid)
xs_idx = xs_idx[keep]
ys_idx = ys_idx[keep]
y_vals = y_vals[keep]
# Step 2: Subsample training points
rng = np.random.default_rng(random_state)
if grid_step is not None and grid_step > 1:
# Grid-based subsampling
gx = xs_idx // grid_step
gy = ys_idx // grid_step
# Unique grid cells
key = gx.astype(np.int64) + (gy.astype(np.int64) << 32)
_, idx = np.unique(key, return_index=True)
xs_sub, ys_sub, y_sub = xs_idx[idx], ys_idx[idx], y_vals[idx]
else:
xs_sub, ys_sub, y_sub = xs_idx, ys_idx, y_vals
# Further random subsampling if needed
n_candidates = y_sub.size
if n_candidates > max_points:
idx = rng.choice(n_candidates, size=max_points, replace=False)
xs_sub, ys_sub, y_sub = xs_sub[idx], ys_sub[idx], y_sub[idx]
# Prepare training data
X_train = np.column_stack([xs_sub, ys_sub]).astype(float)
y_train = y_sub.astype(float)
# Step 3: Build and fit GP
kernel = ConstantKernel(1.0, (1e-3, 1e3)) * Matern(
length_scale=length_scale,
length_scale_bounds=length_scale_bounds,
nu=matern_nu,
) + WhiteKernel(
noise_level=white_noise,
noise_level_bounds=white_noise_bounds,
)
gp = GaussianProcessRegressor(
kernel=kernel,
normalize_y=normalize_y,
n_restarts_optimizer=n_restarts_optimizer,
random_state=random_state,
)
gp.fit(X_train, y_train)
# Step 4: Predict full background map in chunks
background = np.empty((ny, nx), dtype=float)
if get_uncertainty:
uncertainty = np.empty((ny, nx), dtype=float)
else:
uncertainty = None
x_coords = np.arange(nx, dtype=float)
for y0 in range(0, ny, chunk_rows):
y1 = min(ny, y0 + chunk_rows)
rows = np.arange(y0, y1, dtype=float)
# Create (x, y) coordinate pairs for this chunk
X_pred = np.column_stack(
[
np.tile(x_coords, y1 - y0),
np.repeat(rows, nx),
]
)
if get_uncertainty:
pred, pred_std = gp.predict(X_pred, return_std=True)
background[y0:y1, :] = pred.reshape((y1 - y0, nx))
uncertainty[y0:y1, :] = pred_std.reshape((y1 - y0, nx))
else:
pred = gp.predict(X_pred, return_std=False)
background[y0:y1, :] = pred.reshape((y1 - y0, nx))
# Fix #8: Always return (background, uncertainty_or_None)
# When get_uncertainty=False, return None so caller can compute RMS if needed
return background, uncertainty
[docs]
def get_background(image, mask=None, method='sep', size=128, get_rms=False, **kwargs):
"""
Estimate background using various methods.
This is the main interface for background estimation, supporting multiple
algorithms optimized for different scenarios.
Parameters
----------
image : ndarray
Input image
mask : ndarray, optional
Boolean mask (True = masked pixels)
method : str
Background estimation method:
- 'sep': SEP grid-based (default) - Fast, accurate for flat/linear backgrounds
- 'photutils': photutils grid-based - Similar to SEP
- 'polynomial': Low-order polynomial fit - Captures smooth gradients, preserves nebulae
- 'percentile': Percentile filtering (non-grid, robust to stars) - Slow but very robust
- 'morphology': Morphological opening (non-grid, fast) - Good for moderate gradients
- 'gp': Gaussian Process regression (non-grid, very flexible) - Best for complex backgrounds
size : int
Grid size for sep/photutils, or kernel size for percentile/morphology.
For 'gp' method, this parameter is used as the initial length_scale.
Typical values:
- SEP/photutils: 64-128 pixels
- Percentile: 31-51 pixels (smaller for gradients)
- Morphology: 15-35 pixels (should be > star size)
- GP: 32-128 pixels (length scale, should be >> PSF size)
get_rms : bool
If True, return (background, background_rms). For GP method with get_rms=True,
returns the full uncertainty map instead of a scalar RMS.
**kwargs : dict
Additional parameters passed to the method:
- For polynomial: order, sigma, maxiters, directional, direction
(see get_background_polynomial for full list)
- For percentile: percentile, maxiters, sigma
- For morphology: iterations, smooth, smooth_size
- For gp: max_points, grid_step, clip_sigma, n_clip_iter, matern_nu, etc.
(see get_background_gp for full list)
Returns
-------
background : ndarray
Estimated background
background_rms : ndarray or float, optional
Background RMS (if get_rms=True). For GP method, this is a 2D uncertainty map.
Examples
--------
Default grid-based background estimation:
>>> bg = get_background(image, method='sep', size=128)
Non-grid morphological method for moderate gradients:
>>> bg = get_background(image, method='morphology', size=25)
Morphology with extra smoothing to eliminate kernel artifacts:
>>> bg = get_background(image, method='morphology', size=25, smooth_size=15)
Percentile method for very crowded fields:
>>> bg = get_background(image, method='percentile', size=31, percentile=15.0)
Gaussian Process for complex backgrounds with vignetting:
>>> bg = get_background(image, method='gp', size=64, max_points=2000)
Get both background and RMS:
>>> bg, bg_rms = get_background(image, method='sep', get_rms=True)
Get GP background with uncertainty map:
>>> bg, bg_std = get_background(image, method='gp', get_rms=True)
Polynomial fit for gradient removal (preserves nebulae):
>>> bg = get_background(image, method='polynomial', order=2)
Directional gradient with cubic curvature:
>>> bg = get_background(image, method='polynomial', order=3, directional=True)
See Also
--------
get_background_polynomial : Polynomial fit details
get_background_percentile : Percentile filtering details
get_background_morphology : Morphological opening details
get_background_gp : Gaussian Process regression details
Notes
-----
**Method Selection Guide**:
- **Flat backgrounds**: Use 'sep' (default) - fastest and most accurate
- **Linear gradients**: Use 'sep' with size=64-128
- **Smooth gradients with nebulae**: Use 'polynomial' — only captures the
gradient, preserves localised background features by construction
- **Directional gradients**: Use 'polynomial' with ``directional=True`` —
higher-order terms confined to the gradient direction
- **Quadratic gradients**: Use 'morphology' or local gradient fitting (bkg_order)
- **Complex backgrounds** (vignetting, strong curvature): Use 'gp'
- **Very crowded fields**: Use 'percentile' with percentile=10-15
- **Large images**: Avoid 'percentile' (very slow), use 'sep' or 'morphology'
**Performance**:
- SEP: ~1-3 ms (512x512 image)
- Morphology: ~300 ms (512x512 image)
- Percentile: ~2-5 seconds (512x512 image)
- GP: ~1-10 seconds (512x512 image, depends on max_points)
For complex backgrounds (strong curvature, vignetting), consider using
local gradient fitting instead (bkg_order parameter in measure_objects),
or the GP method for maximum flexibility.
"""
if method == 'sep':
# Fix #1: Lazy import
import sep
# Fix #5: Avoid unnecessary copy if already float64 and contiguous
image = np.ascontiguousarray(image, dtype=np.double)
bg = sep.Background(image, mask=mask, bw=size, bh=size, **kwargs)
back, backrms = bg.back(), bg.rms()
elif method == 'photutils':
# Fix #1: Lazy import
import photutils.background
bg = photutils.background.Background2D(image, size, mask=mask, **kwargs)
back, backrms = bg.background, bg.background_rms
elif method == 'percentile':
back = get_background_percentile(image, mask=mask, size=size, **kwargs)
# Estimate RMS from percentile range if needed
# Fix #9: Don't pass **kwargs to RMS estimator (percentile param would
# incorrectly override the fixed p25/p75 used for IQR)
if get_rms:
backrms = estimate_background_rms_percentile(image, mask=mask, size=size)
else:
backrms = None
elif method == 'morphology':
back = get_background_morphology(image, mask=mask, size=size, **kwargs)
# Estimate RMS from local statistics if needed
if get_rms:
backrms = estimate_background_rms_local(image, back, mask=mask, size=size)
else:
backrms = None
elif method == 'gp':
# Use size as length_scale if not explicitly provided
if 'length_scale' not in kwargs:
kwargs['length_scale'] = float(size)
# Set get_uncertainty based on get_rms
kwargs['get_uncertainty'] = get_rms
back, gp_uncertainty = get_background_gp(image, mask=mask, **kwargs)
# Fix #8: When get_uncertainty=False, gp_uncertainty is None.
# Compute RMS via estimate_background_rms_local if get_rms is requested.
if get_rms:
if gp_uncertainty is not None:
backrms = gp_uncertainty
else:
backrms = estimate_background_rms_local(image, back, mask=mask, size=size)
else:
backrms = None
elif method == 'polynomial':
back, poly_rms = get_background_polynomial(image, mask=mask, **kwargs)
backrms = np.full_like(back, poly_rms) if get_rms else None
else:
raise ValueError(f"Unknown background method: {method}")
if get_rms:
return back, backrms
else:
return back