Background Estimation

The module provides tools and algorithms for estimating the background of astronomical data.

Usage

Estimating backgrounds is an important step in performing photometry. Ideally, we could perfectly describe the background with a scalar value or with some distribution. Unfortunately, it's impossible for us to precisely separate the background and foreground signals. Here, we use mixture of robust statistical estimators and meshing to let us get the spatially varying background from an astronomical photo.

Let's show an example

using Photometry
using FITSIO
using Plots

# Download our image, courtesy of astropy
hdu = FITS(download("https://rawcdn.githack.com/astropy/photutils-datasets/8c97b4fa3a6c9e6ea072faeed2d49a20585658ba/data/M6707HH.fits"))
image = read(hdu[1])

# Plot
function imshow(image; kwargs...)
    xs, ys = axes(image)
    data = transpose(image)
    heatmap(xs, ys, data; aspect_ratio=1, xlim=extrema(xs), ylim=extrema(ys), kwargs...)
end

imshow(image)

Now let's try and estimate the background using estimate_background. First, we'll si gma-clip to try and remove the signals from the stars. Then, the background is broken down into boxes, in this case of size (50, 50). Within each box, the given statistical estimators get the background value and RMS. By default, we use SourceExtractorBackground and StdRMS. This creates a low-resolution image, which we then need to resize. We can accomplish this using an interpolator, by default a cubic-spline interpolator via ZoomInterpolator. The end result is a smooth estimate of the spatially varying background and background RMS.

# sigma-clip
clipped = sigma_clip(image, 1, fill=NaN)

# get background and background rms with box-size (50, 50)
bkg, bkg_rms = estimate_background(clipped, 50)

# plot
plot(
    imshow(image, title="Original"),
    imshow(clipped, title="Sigma-Clipped"),
    imshow(bkg, title="Background"),
    imshow(bkg_rms, title="Background RMS"),
    layout=(2, 2), ticks=false
)

We could apply a median filter, too, by specifying filter_size

# get background and background rms with box-size (50, 50) and filter_size (5, 5)
bkg_f, bkg_rms_f = estimate_background(clipped, 50, filter_size=5)

# plot
plot(
    imshow(bkg, title="Unfiltered", ylabel="Background"),
    imshow(bkg_f, title="Filtered"),
    imshow(bkg_rms, ylabel="RMS"),
    imshow(bkg_rms_f);
    layout=(2, 2), ticks=false
)

Now we can see our image after subtracting the filtered background and ready for Aperture Photometry!

subt = image .- bkg_f[axes(image)...]
plot(
    imshow(image, title="Original", colorbar=false),
    imshow(subt, title="Subtracted");
    layout=(1, 2), size=(600, 260),
    xlims=(400, 800), ylims=(400, 800),
    clims=(minimum(subt), maximum(image)),
    ticks=false, aspect_ratio=1
)

IDW Interpolator

Here is a quick example using the IDWInterpolator

b1, r1 = estimate_background(clipped, 50, filter_size=5)
b2, r2 = estimate_background(clipped, 50, itp=IDWInterpolator(50), filter_size=5)

plot(
    imshow(b1, title="ZoomInterpolator", ylabel="Background"),
    imshow(b2, title="IDWInterpolator"),
    imshow(r1, ylabel="RMS"),
    imshow(r2);
    layout=(2, 2), ticks=false
)

API/Reference

Photometry.Background.estimate_background — Function

estimate_background(data;
    location=SourceExtractorBackground(),
    rms=StdRMS(),
    dims=:)

Perform scalar background estimation using the given estimators.

The value returned will be two values corresponding to the estimated background and the estimated background RMS. The dimensionality will depend on the dims keyword.

location and rms can be anything that is callable, for example median, or one of the estimators we provide in Background Estimators.

Examples

julia> data = ones(3, 5);

julia> bkg, bkg_rms = estimate_background(data)
(1.0, 0.0)

julia> using Statistics: median

julia> bkg, bkg_rms = estimate_background(data; location=median, rms=MADStdRMS())
(1.0, 0.0)

See Also

Location Estimators, RMS Estimators

source

estimate_background(data, box_size;
    location=SourceExtractorBackground(),
    rms=StdRMS(),
    itp=ZoomInterpolator(box_size),
    edge_method=:pad,
    [filter_size])

Perform 2D background estimation using the given estimators mapped over windows of the data..

This function will estimate backgrounds in boxes of size box_size. When size(data) is not an integer multiple of the box size, there are two edge methods: :pad and :crop. The default is to pad (and is recommend to avoid losing image data). If box_size is an integer, the implicit shape will be square (eg. box_size=4 is equivalent to box_size=(4,4)).

For evaluating the meshes, each box will be passed into location to estimate the background and then into rms to estimate the background root-mean-square value. These can be anything that is callable, like median or one of our Background Estimators.

Once the meshes are created they will be median filtered if filter_size is given. filter_size can be either an integer or a tuple, with the integer being converted to a tuple the same way box_size is. Filtering is done via ImageFiltering.MapWindow.mapwindow. filter_size must be odd.

After filtering (if applicable), the meshes are passed to the itp to recreate a low-order estimate of the background at the same resolution as the input.

Note

If your box_size is not an integer multiple of the input size, the output background and rms arrays will not have the same size.

See Also

Location Estimators, RMS Estimators, Interpolators

source

Photometry.Background.sigma_clip — Function

sigma_clip(x, sigma; fill=:clamp, center=median(x), std=std(x, corrected=false))
sigma_clip(x, sigma_low, sigma_high; fill=:clamp, center=median(x), std=std(x, corrected=false))

This function returns sigma-clipped values of the input x.

Specify the upper and lower bounds with sigma_low and sigma_high, otherwise assume they are equal. center and std are optional keyword arguments which are functions for finding central element and standard deviation.

If fill === :clamp, this will clamp values in x lower than center - sigma_low * std and values higher than center + sigma_high * std. Otherwise, they will be replaced with fill.

Examples

julia> x = randn(100_000);

julia> extrema(x)
(-4.496308951466683, 4.080724496910187)

julia> x_clip = sigma_clip(x, 1);

julia> extrema(x_clip) # should be close to (-1, 1)
(-1.0042721545326967, 0.9957463910682249)

source

Photometry.Background.sigma_clip! — Function

sigma_clip!(x, sigma; fill=:clamp, center=median(x), std=std(x))
sigma_clip!(x, sigma_low, sigma_high; fill=:clamp, center=median(x), std=std(x))

In-place version of sigma_clip

Warning

sigma_clip! mutates the element in place and mutation cannot lead to change in type. Please be considerate of your input type, because if you are using Int64 and we try to clip it to 0.5 an InexactError will be thrown.

To avoid this, we recommend converting to float before clipping, or using sigma_clip which does this internally.

source