Metrics

ADI.Metrics

This module provides code for analyzing the results from ADI in a way that is interpretable statistically. Some of the key functionalities are signal-to-noise, significance, the receiver operating characteristic, and the contrast curve.

detectionmap([method=snr], data, fwhm; fill=0)

Parallel implementation of arbitrary detection mapping applied to each pixel in the input image. Any invalid values will be set to fill.

The following methods are provided in the Metrics module:

• snr - signal-to-noise ratio (S/N) using student-t statistics to account for small sample penalty.
• significance - Gaussian signifance using student-t statistics to account for small samples penalty.
• noise - Standard deviation of apertures in each annulus.

snr(data, position, fwhm)

Calculate the signal to noise ratio (SNR, S/N) for a test point at position using apertures of diameter fwhm in a residual frame.

Uses the method of Mawet et al. (2014) which includes penalties for small sample statistics. These are encoded by using a student's t-test for calculating the SNR.

significance(data, position, fwhm)

Calculates the Gaussian significance from the signal-to-noise ratio (SNR, S/N) for a test point at position using apertures of diameter fwhm in a residual frame.

The Gaussian signifiance is calculated from converting the SNR confidence limit from a student t distribution to a Gaussian via

$\text{sig}(\text{SNR}) = \Phi^{-1}\left[\int_0^\text{SNR}{t_\nu(x)dx}\right]$

where the degrees of freedom $\nu$ is given as $2\pi r / \Gamma - 2$ where r is the radial distance of each pixel from the center of the frame.

See Also

snr

noise(data, position, fwhm)

Calculate the statistical noise for a test point at position using apertures of diameter fwhm in a residual frame.

Uses the biweight middeviance (square root of biweight midvariance, a robust estimator of the standard deviation) of the apertures in the entire annulus. This is distinct from the snr noise calculation, which defines a confidence interval using student-t statistics. This means you cannot simply create a noise map and divide it from the signal to create an equivalent S/N map.

stimmap(residuals, angles)

Calculate the standardized trajectory intensity mean (STIM) map. The inputs are a cube of residuals and the corresponding parallactic angles.

This method seeks to improve upon the typical student-t S/N tests (snr, significance) by calculating statistics in the temporal domain instead of the spatial domain. This is why the full residual cube is required rather than a reduced frame.

In particular, the STIM map is robust to detections with multiple objects or extended sources within the same annuli, which results in very high noise estimates using spatial methods. The STIM map also performs better at small angular separations, since the temporal domain has no limitations from limited resolution elements.

Pairet et al. 2019 derives a detection threshold of τ = 0.5 for the STIM map. The detection threshold can be calculated for a specific dataset using stim_threshold.

Examples

julia> cube, angles = # load data

julia> S = subtract(PCA(10), cube, angles);

julia> sm = stimmap(S, angles);

References

1. Pairet et al. 2019 "STIM map: detection map for exoplanets imaging beyond asymptotic Gaussian residual speckle noise"

See Also

stim_threshold

stim_threshold([stimmap, ] residuals, angles)

Calculate the detection threshold for the standardized trajectory intensity mean (STIM) map. This method uses the same residual cube as stimmap but adds an additional step of estimating the residual noise by derotating the residuals with the opposite parallactic angles.

If the STIM map has already been calculated, it can be passed in, otherwise it will be calculated in addition to the noise map. Note this will not return the STIM map, only the threshold.

The threshold is derived in section 5.1 of Pairet et al. 2019 as the ratio of the number of pixels above the approximated noise map. They found a value of τ = 0.5 to be typical.

Examples

julia> cube, angles = # load data

julia> S = subtract(PCA(10), cube, angles);

julia> sm = stimmap(S, angles);

julia> τ = stim_threshold(sm, S, angles);

References

1. Pairet et al. 2019 "STIM map: detection map for exoplanets imaging beyond asymptotic Gaussian residual speckle noise"

See Also

stimmap

Metrics.stim(cube; dims)

Calculates the STIM map of a derotated cube along the given dims. dims should correspond to the temporal axis of the cube. The STIM statistic is the slice mean divided by the slice standard deviation. Invalid values will become 0.

slimmap(residuals, angles; N)
slimmap(stim_maps; N)

Calculate the STIM largest intensity mask (SLIMask) map. This is an ensemble method which averages the STIM maps for multiple residual cubes. In addition to computing this average, for each STIM map all but the brightest N pixels will be masked out and eventaully averaged to create the SLIMask.

N should represent a cutoff for the number of expected companions. For example, if the FWHM of the companion signal is 5 pixels, then the area under the fwhm is ~20 pixels. If I want to probe the brightest 4 potential companions, I would mask all but the N = 20 * 4 = 80 brightest pixels.

Both the average STIM map and the SLIMask will be returned, and the two can be multiplied together to produce the SLIMask-STIM map. This achieves two things: first, the masking makes most of the pixels 0, providing better visual contrast in the map, and second, by averaging the mask, pixels which are not consistently in the brightest N for each STIM map will have lower probabilities in the corresponding SLIMask-STIM map.

Examples

This example recreates the analysis shown in Pairet, B. (2020) where the SLIM map is computed with the ensemble of residual cubes produced by increasing ranks of PCA subtraction.

julia> cube, angles = # load data

julia> algs = PCA.(5:25);

julia> residual_cubes = subtract.(algs, Ref(cube));

julia> stim_av, slimmask = slimmap(residual_cubes, angles; N=100);

julia> slim_prob_map = stim_av .* slimmask;

References

1. Pairet, B. 2020 "Signal processing methods for high-contrast observations of planetary systems"

See also

stimmap, stim

throughput(alg, cube, angles, psf;
           fwhm, nbranch=1, theta=0, inner_rad=1,
           fc_rad_sep=3, snr=100, kwargs...)

Calculate the throughput of alg by injecting fake companions into cube and measuring the relative photometry of each companion in the reduced frame. The photometry is measured using a circular aperture with a diameter matching the fwhm. Any additional kwargs will be passed to alg when it is called.

Keyword Arguments

• nbranch - number of azimuthal branches to use
• theta - position angle of initial branch
• inner_rad - position of innermost planet in FWHM
• fc_rad_sep - the separation between planets in FWHM for a single reduction
• snr - the target signal to noise ratio of the injected planets
• reduced_empty - the collapsed residual frame for estimating the noise. Will process using alg if not provided.

throughput(alg, cube, angles, psf, position;
           fwhm, snr=100,
           reduced_empty=nothing,
           verbose=true, kwargs...)

Calculate the throughput of alg by injecting psf into cube at the given position and measuring the relative photometry of the companion in the reduced frame. The photometry is measured using a circular aperture with a diameter matching the fwhm. Any additional kwargs will be passed to alg when it is called.

If position is a tuple or a vector, it will be parsed as the cartesian coordinate (x, y). If position is a CoordinateTransformations.Polar it will be parsed as polar coordinates from the center of the cube. Note the Polar type expects angles in radians.

Keyword Arguments

• fwhm - the full width at half-maximum
• snr - the target signal to noise ratio of the injected planets
• reduced_empty - the collapsed residual frame for estimating the noise. Will process using alg if not provided.
• verbose - show informative messages during process

contrast_curve(alg, cube, angles, psf;
               fwhm, sigma=5, nbranch=1, theta=0, inner_rad=1,
               starphot=Metrics.estimate_starphot(cube, fwhm),
               fc_rad_sep=3, snr=100, k=2, smooth=true,
               subsample=true, kwargs...)

Calculate the throughput-calibrated contrast. This first processes the algorithmic throughput by injecting instances of psf into cube. These are processed through alg and the ratio of the recovered flux to the injected flux is calculated. These companions are injected in resolution elements across the frame, which can be changed via the various keyword arguments.

The throughput can only be calculated for discrete resolution elements, but we typically want a much smoother curve. To accomplish this, we measure the noise (the standard deviation of all resolution elements in an annulus at a given radius) for every pixel in increasing radii. We then interpolate the throughput to this grid and return the subsampled curves.

Returned Fields

• distance - The radial distance (in pixels) for each measurement
• contrast - The Gaussian sensitivity
• contrast_corr - The Student-t sensitivity
• noise - The noise measured for each distance
• throughput - The throughput measured for each distance.

Keyword Arguments

• sigma - The confidence limit in terms of Gaussian standard deviations
• starphot - The flux of the star. By default calculates the flux in the central core.

Injection Options (See also throughput)

• nbranch - number of azimuthal branches to use
• theta - position angle of initial branch
• inner_rad - position of innermost planet in FWHM
• fc_rad_sep - the separation between planets in FWHM for a single reduction
• snr - the target signal to noise ratio of the injected planets

Subsampling Options (See also Metrics.subsample_contrast)

• subsample - If true, subsamples the throughput measurements to increase density of curve
• k - The order of the BSpline used for subsampling the throughput
• smooth - If true, will smooth the subsampled noise measurements with a 2nd order Savitzky-Golay filter

Any additional kwargs will be passed to alg when it is called.

contrast_curve(alg, cube, angles, psf; fwhm=fwhm) |> DataFrame

Metrics.subsample_contrast(empty_frame, distance, throughput;
                           fwhm, starphot, sigma=5, inner_rad=1,
                           theta=0, smooth=true, k=2)

Helper function to subsample and smooth a contrast curve.

Contrast curves, by definition, are calculated with discrete resolution elements. This can cause contrast curves to have very few points instead of appearing as a continuously measured statistic across the pixels. We alleviate this by sub-sampling the throughput (via BSpline interpolation) across each pixel (instead of each resolution element).

The noise can be found efficiently enough, so rather than interpolate we measure the noise in annuli of width fwhm increasing in distance by 1 pixel. We measure this noise in empty_frame, which should be a 2D reduced ADI frame.

The noise measurements can be noisy, so a 2nd order Savitzky-Golay filter can be applied via smooth. This fits a quadratic polynomial over a window of fwhm/2 points together to reduce high-frequency jitter.

Examples

Here is an example which calculates the exact contrast curve in addition to a subsampled curve without re-calculating the throughput.

cube, angles, psf = # load data

alg = PCA(10)
cc = contrast_curve(alg, cube, angles, psf; fwhm=8.4, subsample=false)
reduced_empty = alg(cube, angles)
sp = Metrics.estimate_starphot(cube, 8.4)
cc_sub = Metrics.subsample_contrast(reduced_empty, cc.distance, cc.throughput;
                                    fwhm=8.4, starphot=sp)

Metrics.estimate_starphot(cube, fwhm)
Metrics.estimate_starphot(frame, fwhm)

Simple utility to estimate the stellar photometry by placing a circular aperture with fwhm diameter in the center of the frame. If a cube is provided, first the median frame will be found.