ADI.ADIAlgorithm

This abstract type is used for defining ADI algorithms. Algorithms are stateful objects that define the options for a given algorithm (e.g. the number of components used in PCA decomposition). The most direct usage of an algorithm is to use it to fit HCI data; that is, given a sample of pixels, apply the algorithm using the given options and return an object containing all the necessary information to reconstruct the input data.

See the extended help (??ADIAlgorithm) for interface details.

Extended help

Interface

To extend ADIAlgorithm you may implement the following

ADI.fit(::Alg, data::AbstractMatrix; kwargs...)

Fit the data (flattened into a matrix). To support RDI, ensure the ref keyword argument is usable (ref is also a flattened matrix). This is the only method you need to implement for a new ADIAlgorithm, along with a suitable ADIDesign.

ADI.jl automatically coaxes the cube input into a matrix for use with fit, appropriately handling the various geometries. When available, this input is a view, so if the algorithm requires dense arrays, make sure to call collect when appropriate. If a given algorithm doesn't support the default operations, all that needs to be done is override the default behavior (for an example, see the GreeDS implementation).

reconstruct(::Alg, cube; kwargs...)

Fit the data using the algorithm and return a cube with the estimate of the PSF. By default uses the reconstruction from the ADIDesign fit to the data.

subtract(::Alg, cube; kwargs...)

Fit the data using the algorithm and return a cube that has had the PSF estimate subtracted. By default, calls reconstruct and subtracts it from cube.

process(::ADIAlgorithm, cube; kwargs...)
(::ADIAlgorithm)(cube; kwargs...)

Fully process the data (estimate, subtract, collapse). By default, derotates and collapses output from subtract. You only need to define process, since the functor version is supplied automatically.

(::ADIAlgorithm)(cube, angles; [ref], kwargs...)

Fully process an ADI data cube using subtract and collapsing the residuals. This is a convenient alias for calling process Keyword arguments will be passed to ADI.fit.

ADI.ADIDesign

An abstract type used for holding the output of ADI.fit. The purpose of these types is to hold the minimum information required to reconstruct the PSF estimate (e.g. weights and principal components from PCA). ADI.jl includes two designs- ADI.ClassicDesign and ADI.LinearDesign.

Interface

When implementing a new algorithm, if your outputs do not fit into either of those designs, you will have to create your own ADIDesign with the following methods-

ADI.design(::Design)

accessor for the pertinent data to approximate the PSF

reconstruct(::Design)

return the approximate PSF estimate as a matrix

Classic(;method=median)
Classic(method)

Classic PSF subtraction using the median of entire data cube. If another statistic is desired, like the mean, it can be passed as an argument as long as it supports slicing multi-dimensional arrays.

References

1. Marois et al. 2006 Angular Differential Imaging: A Powerful High-Contrast Imaging Technique

ADI.ClassicDesign(n, frame)

Output for the Classic algorithm which contains the static frame unrolled into a vector (with size Npx). reconstruct will tile this vector in a non-allocating way n times to appear like a flattened cube. ADI.design will return the static frame .

DoubleSDI(alg)
DoubleSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses alg_spec to reduce each spectral cube slice in the SDI tensor. Then, the spectral residual frames will be reduced using alg_temp, which will include the derotation and final combination.

The difference between DoubleSDI and SliceSDI is that DoubleSDI does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. SliceSDI does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

Examples

julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = DoubleSDI(Classic(), GreeDS(15))(data, angles, scales)

Framewise(alg; limit=Inf, delta_rot=nothing)
Framewise(algs::AbstractVector; limit=Inf, delta_rot=nothing)

Wrap an algorithm such that the underlying data will be processed frame by frame. For each frame a reference library is created from the data. This reference can be filtered by rejecting frames which have not rotated a sufficient parallactic angle. delta_rot sets the required arc length for rotation in units of the FWHM. If delta_rot is nothing there will be now temporal filtering. The number of frames retained can be specified with limit, e.g. the 4 closest frames in time with the target frame.

The following keyword arguments must be provided to reconstruct or subtract

• angles - the measured parallactic angles for each frame

If delta_rot is provided, the following additional keyword arguments must be provided to reconstruct or subtract.

• fwhm - the FWHM of the instrument in pixels. Will be set to the width of a MultiAnnulusView
• r - The radius of the arc to calculate the parallactic angle threshold. Will be set automatically if using AnnulusView or MultiAnnulusView.

In addition, Framewise versions of algorithms do not implement ADI.fit and do not currently support RDI.

Annular reduction

In the case of reducing a MultiAnnulusView, a vector of algorithms can be used, each one corresponding to each annulus of the view. In this case, too, delta_rot can be given as a vector or as a tuple. If it is given as a tuple, delta_rot will increase linearly from the first value to the last value across each annulus.

Examples

julia> cube, angles = # load data

julia> alg = PCA(10) # the algorithm to use on each reference

julia> res = Framewise(alg)(cube, angles);

julia> mav = MultiAnnulusView(cube, 5; inner=5);

julia> res_ann = Framewise(alg, delta_rot=(0.1, 1))(mav, angles);

GreeDS(alg=PCA(); threshold=0)
GreeDS(ncomps; threshold=0, options...)

Performs the greedy disk subtraction (GreeDS) algorithm.

This method is an iterative approach to standard ADI reduction which seeks to minimize over-subtraction by constraining the low-rank matrix approximation from alg. By default, uses PCA. If ncomps or other PCA options are provided, they will be passed to the constructor.

Algorithms

The following algorithms work natively with GreeDS: PCA and NMF

References

1. Pairet et al. 2018 "Reference-less algorithm for circumstellar disks imaging"
2. Pairet et al. 2020 "MAYONNAISE: a morphological components analysis pipeline for circumstellar disks and exoplanets imaging in the near infrared"

LOCI(; dist_threshold=nothing, metric=Cityblock())

Local optimal combination of images (LOCI).

If provided, the frames used for the reference are filtered to only include the frames whose pairwise distances are within the dist_threshold quantile. The distances are measured using Distances.jl and the metric used for measuring the distacnes can be specified with metric (by default uses the Manhattan/cityblock metric).

# only use the most-similar 90th-percentile frames
alg = Framewise(LOCI(dist_threshold=0.90))

References

1. Lafreniere et al. (2007) "A New Algorithm for Point-Spread Function Subtraction in High-Contrast Imaging: A Demonstration with Angular Differential Imaging"

ADI.LinearDesign(coeffs, basis)

A "linear" design implies the use of some linear basis for reconstructing data along with a set of coefficients or weights. The reconstruction will be the matrix product of the basis and the weights, Z * w.

ADI.design will return (basis, ceoffs), and you can also extract them via iteration, like Z, w = design.

NMF(;ncomps=nothing)
NMF(ncomps)

Use non-negative matrix factorization (NMF, NNMF) to form a non-negative, low-rank, and orthonormal basis of the input. The implementation of the underlying NMF is provided by NMF.jl. The implementation uses a non-negative SVD for initialization and a coordinate-descent solver to fit.

If ncomps is nothing, it will be set to the number of frames in the reference cube when processed.

target = cube .- minimum(cube)
S = reconstruct(NMF(), target, angles)

References

1. Ren et al. 2018 Non-negative Matrix Factorization: Robust Extraction of Extended Structures

PCA(ncomps; options...)
PCA(;ncomps=nothing, options...)

Use principal components analysis (PCA) to form a low-rank orthonormal basis of the input. Uses deterministic singular-value decomposition (SVD) to decompose data.

If ncomps is nothing, the basis will not be truncated (i.e. ncomps is equal to the number of frames). ncomps can be set to :noise or :pratio to automatically choose the number of components using the residual frame noise or principal ratio, respectively. For more information, see the extended help.

References

1. Soummer, Pueyo, and Larkin (2012) "Detection and Characterization of Exoplanets and Disks Using Projections on Karhunen-Loève Eigenimages"

Extended help

Optimizing ncomps

There are a few ways to optimize ncomps using the input data. Additional options for the optimization are listed below

1. ncomps=:noise - residual noise optimization
2. ncomps=:pratio - principal ratio optimization

Residual noise optimization

This technique progressively increases ncomps at each step measuring the pixel-to-pixel noise (standard deviation) in the residual data. Iteration will stop when the noise is not improving beyond a threshold. This is suited for data with similar statistical characteristics, such as an annulus more so than a full-frame cube.

• collapse=false - if true, the temporal median of the residual data will be used for measuring the noise.
• noise_error=1e-3 - the threshold for the minimal noise improvement looking back 2 iterations

Principal ratio optimization

This technique chooses the number of components required to explain some ratio of the total variance in the data. This is known as the principal ratio or the explained variance ratio. The explained variance is measured by transforming the singular values of the SVD decomposition (Λ = @. S^2 / (n - 1)).

• pratio=0.9 - the target principal ratio (between 0 and 1)

ADI.SDIAlgorithm <: ADI.ADIAlgorithm

Spectral differential imaging (SDI) algorithms. These work on 4-D SDI tensors. To use these algorithms, simply treat them like functions (or call process)

(::SDIAlgorithm)(data::AbstractArray{T,4}, angles, scales; [ref] kwargs...)

The data is expected to be laid out in (nx, ny, nλ, nf) format, so you may need to permutedims before processing the data. The scales correspond to the relative wavelength scales for each spectrum, and can be retrieved with HCIToolbox.scale_list.

Algorithms

The current SDI implementations are

• SingleSDI
• DoubleSDI
• SliceSDI

SingleSDI(alg)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in a single pass. This means that each channel will be scaled and then concatenated together, so an SDI tensor (nx, ny, nλ, nf) becomes a stack (nx, ny, nλ * nf) which is processed by the underlying alg like ADI data.

SingleSDI(PCA(15))(data, angles, scales)

PCA(15)(data, angles, scales)

SliceSDI(alg)
SliceSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses alg_temp to reduce each temporal cube slice in the SDI tensor. These residuals will be rescaled and stacked into a new cube. Then, the temporal residual frames will be reduced using alg_spec, which will include the derotation and final combination.

The difference between SliceSDI and DoubleSDI is that DoubleSDI does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. SliceSDI does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

Examples

julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = SliceSDI(Classic(), GreeDS(15))(data, angles, scales)

ADI.design(::ADIDesign)

Return the pertinent data required to form the PSF approximation. For example, weights and components for PCA/NMF or the median frame for classic ADI.

expand_rotate(frame, angles, threshold; kwargs...)

Takes a frame, expands it into a cube, rotates it clockwise by angles, and min-clips at threshold. Keyword arguments will be passed to HCIToolbox.derotate.

ADI.fit(::ADIAlgorithm, cube; [ref], kwargs...)

Given the description of an algorithm and the appropriate options, take the pixels from cube and fit them, returning an (ADIDesign) containing the necessary information from the fit (e.g. the principal components from PCA decomposition).

If the algorithm supports reference differential imaging (RDI), the reference cube can be passed by the keyword argument ref.

process(alg, cube, angles; [ref], kwargs...)

Fully process an ADI data cube using subtract and collapsing the residuals. Keyword arguments will be passed to ADI.fit.

reconstruct(alg, cube; [ref], kwargs...)

Reconstruct the PSF approximation for the given algorithm, using ref as the reference cube if given and supported by the algorithm.

Examples

julia> cube, angles = # load data

julia> S = reconstruct(PCA(10), cube);

julia> size(S) == size(cube)
true

julia> flat_res = collapse(cube .- S, angles); # form resid, derotate, and combine

subtract(alg, cube; [ref], kwargs...)

Reconstruct the PSF approximation for the given algorithm and subtract it from cube, using ref as the reference cube if given and supported by the algorithm.

Examples

julia> cube, angles = # load data

julia> R = subtract(PCA(10), cube);

julia> size(R) == size(cube)
true

julia> flat_res = collapse(R, angles); # derotate, and combine