diff --git a/dev/algorithms/api/index.html b/dev/algorithms/api/index.html index 32ca4ba..9446057 100644 --- a/dev/algorithms/api/index.html +++ b/dev/algorithms/api/index.html @@ -1,17 +1,17 @@ API/Reference · ADI.jl
Edit on GitHub

API/Reference

    ADI.ADIAlgorithm — Type
    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.

    source
    ADI.reconstruct — Function
    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
+(::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.

    source
    ADI.reconstruct — Function
    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
    source
    ADI.subtract — Function
    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> flat_res = collapse(cube .- S, angles); # form resid, derotate, and combine
    source
    ADI.subtract — Function
    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
    source
    ADI.process — Function
    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.

    source
    ADI.fit — Function
    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.

    source
    ADI.design — Function
    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.

    source
    ADI.ADIDesign — Type
    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

    source
    ADI.LinearDesign — Type
    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.

    source
    +julia> flat_res = collapse(R, angles); # derotate, and combinesource
    ADI.process — Function
    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.

    source
    ADI.fit — Function
    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.

    source
    ADI.design — Function
    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.

    source
    ADI.ADIDesign — Type
    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

    source
    ADI.LinearDesign — Type
    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.

    source
    diff --git a/dev/algorithms/classic/index.html b/dev/algorithms/classic/index.html index c75ffef..2a25907 100644 --- a/dev/algorithms/classic/index.html +++ b/dev/algorithms/classic/index.html @@ -1,3 +1,3 @@ Classic · ADI.jl
    Edit on GitHub

    Classic

    API/Reference

    ADI.Classic — Type
    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
    source
    ADI.ClassicDesign — Type
    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 .

    source
    +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
    source
    ADI.ClassicDesign — Type
    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 .

    source
    diff --git a/dev/algorithms/greeds/index.html b/dev/algorithms/greeds/index.html index 4c8bdbb..f0d2b37 100644 --- a/dev/algorithms/greeds/index.html +++ b/dev/algorithms/greeds/index.html @@ -1,3 +1,3 @@ GreeDS · ADI.jl
    Edit on GitHub

    GreeDS

    API/Reference

    ADI.GreeDS — Type
    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.

    Note

    The GreeDS algorithm requires fully reconstructing a cube at each iteration, which requires knowing the geometry of the input (full-frame, annulus, etc.) and the corresponding parallactic angles. These angles must be passed as a keyword argument angles. In the case of reducing data, e.g. GreeDS()(cube, angles) the angles will be passed automatically. It is important to clarify, these angles should correspond to the reference data in the case of RDI, e.g. GreeDS()(cube, angles; ref=ref_cube, angles=ref_angles)

    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"
    source
    +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.

    Note

    The GreeDS algorithm requires fully reconstructing a cube at each iteration, which requires knowing the geometry of the input (full-frame, annulus, etc.) and the corresponding parallactic angles. These angles must be passed as a keyword argument angles. In the case of reducing data, e.g. GreeDS()(cube, angles) the angles will be passed automatically. It is important to clarify, these angles should correspond to the reference data in the case of RDI, e.g. GreeDS()(cube, angles; ref=ref_cube, angles=ref_angles)

    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"
    source diff --git a/dev/algorithms/loci/index.html b/dev/algorithms/loci/index.html index df16f88..e925011 100644 --- a/dev/algorithms/loci/index.html +++ b/dev/algorithms/loci/index.html @@ -1,3 +1,3 @@ LOCI · ADI.jl
    Edit on GitHub

    LOCI

    API/Reference

    ADI.LOCI — Type
    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).

    Traditional LOCI

    Unless LOCI is applied Framewise, no distance thresholding will occur. In order to recreate the traditional LOCI algorithm, consider constructing an algorithm like

    # 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"
    source
    +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"
    source diff --git a/dev/algorithms/nmf/index.html b/dev/algorithms/nmf/index.html index bf3d847..30b5662 100644 --- a/dev/algorithms/nmf/index.html +++ b/dev/algorithms/nmf/index.html @@ -1,4 +1,4 @@ NMF · ADI.jl
    Edit on GitHub

    NMF

    API/Reference

    ADI.NMF — Type
    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.

    Non-negativity constraint

    NMF is not designed to fit negative values. This algorithm will warn you (but will not error) if a target or reference cube contains negative values. The output may seem reasonable, but it is not well-defined with respect to the NMF algorithm. To overcome this, rescaling the data by its minimum before processing is necessary

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

    When doing full-frame reduction (e.g. NMF()(cube, angles)) this is handled automatically, so this constraint only applies to the lower-level API and methods which rely on those, such as GreeDS. In general, if you see warnings, heed them.

    References

    1. Ren et al. 2018 Non-negative Matrix Factorization: Robust Extraction of Extended Structures
    source
    +S = reconstruct(NMF(), target, angles)

    When doing full-frame reduction (e.g. NMF()(cube, angles)) this is handled automatically, so this constraint only applies to the lower-level API and methods which rely on those, such as GreeDS. In general, if you see warnings, heed them.

    References

    1. Ren et al. 2018 Non-negative Matrix Factorization: Robust Extraction of Extended Structures
    source diff --git a/dev/algorithms/pca/index.html b/dev/algorithms/pca/index.html index 686657e..d84100d 100644 --- a/dev/algorithms/pca/index.html +++ b/dev/algorithms/pca/index.html @@ -1,3 +1,3 @@ PCA · ADI.jl
    Edit on GitHub

    PCA

    API/Reference

    ADI.PCA — Type
    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)
    source
    +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.

    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)).

    source diff --git a/dev/api/index.html b/dev/api/index.html index 658212e..bc30d20 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,2 +1,2 @@ -Index · ADI.jl
    Edit on GitHub

    Index

    +Index · ADI.jl
    Edit on GitHub

    Index

    diff --git a/dev/benchmarks/index.html b/dev/benchmarks/index.html index 124c03b..1395955 100644 --- a/dev/benchmarks/index.html +++ b/dev/benchmarks/index.html @@ -22,225 +22,225 @@ ) - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - + + + - - - - - - - - - + + + - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - + + - - - - - - + + - - - - - - +

    Please note the log-scale for the left figure.

    Detection Maps

    This benchmark measures the duration to produce a signal-to-noise ratio (S/N) map. Rather than test exact cubes, these benchmarks test randomly generated frames of various sizes. The FWHM is fixed at 5.

    snrmap_data = CSV.File(benchdir("snrmap_benchmarks.csv")) |> DataFrame |> sort!
 snrmap_groups = groupby(snrmap_data, :framework)

    GroupedDataFrame with 2 groups based on key: framework

    First Group (5 rows): framework = InlineStrings.String7("ADI.jl")

    frameworkNtime
    String7Int64Float64
    1ADI.jl26010.0131519
    2ADI.jl102010.135207
    3ADI.jl404011.24764
    4ADI.jl906014.30251
    5ADI.jl16080110.2164

    ⋮

    Last Group (3 rows): framework = InlineStrings.String7("VIP")

    frameworkNtime
    String7Int64Float64
    1VIP26010.950529
    2VIP102016.9494
    3VIP4040147.5746
    @df snrmap_data scatter(
     :N,
@@ -252,103 +252,103 @@
 )
    - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + - - - - + + +

    Contrast Curves

    Finally, this benchmark measures the duration to generate a contrast curve for analyzing the algorithmic throughput of an ADI algorithm. For both benchmarks 3 azimuthal branches are used for throughput injections and a FWHM of 8. A Gaussian PSF function is evaluated in a (21, 21) grid for the injections. The data used are $\beta$ Pictoris and HR8799 from HCIDatasets.jl.

    contrast_data = CSV.File(benchdir("contrast_benchmarks.csv")) |> DataFrame |> sort!
 cube_labels = @. ifelse(contrast_data.N == 622261, "Beta Pictoris", "HR8799")
 insertcols!(contrast_data, 4, :cube => cube_labels)
@@ -362,119 +362,119 @@ 
    
 )
    - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + + - - - - - - + + - - - - - - -

    Please note the log-scale.

     + +

    Please note the log-scale.

    diff --git a/dev/examples/betapictoris/index.html b/dev/examples/betapictoris/index.html index 67f4b01..58deae5 100644 --- a/dev/examples/betapictoris/index.html +++ b/dev/examples/betapictoris/index.html @@ -16,93 +16,93 @@ imshow(reduced) - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - @@ -4192,93 +4197,93 @@

    - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - + - - + - - - + - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - +

    The contrast uses a robust estimator for the noise, which means the bright companion doesn't overly bias the contrast measurement. Nonetheless, it is good form to remove the companion signal in a maximum likelihood framework. For convenience here, let's use the :cube_empty entry for BetaPictoris, which already has the companion removed.

    cube_empty = BetaPictoris[:cube_empty]
 reduced_empty = alg(cube_empty, angles)
 imshow(reduced_empty)
    - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - -

    This page was generated using Literate.jl.

    + +

    This page was generated using Literate.jl.

    diff --git a/dev/examples/geometries/index.html b/dev/examples/geometries/index.html index bcbd551..2685f6d 100644 --- a/dev/examples/geometries/index.html +++ b/dev/examples/geometries/index.html @@ -19,93 +19,93 @@ imshow(res, title=name) - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - - - - - - - - - - - - - - - - - + - + - + - - + - - + - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + - - - + - - - + - + - + -S/N maps - + - - + - - + - - - + S/N maps - - + - - - + S/N maps - - + - - - + S/N maps - - + - - - + S/N maps - - + - - - + S/N maps - - + - - - + S/N maps -

    This page was generated using Literate.jl.

    +

    This page was generated using Literate.jl.

    diff --git a/dev/framewise/index.html b/dev/framewise/index.html index 85344f9..b63565d 100644 --- a/dev/framewise/index.html +++ b/dev/framewise/index.html @@ -8,7 +8,7 @@ julia> mav = MultiAnnulusView(cube, 5; inner=5); -julia> res_ann = Framewise(alg, delta_rot=(0.1, 1))(mav, angles);source

    API/Reference

    HCIToolbox.AnnulusView — Type
    AnnulusView(cube::AbstractArray{T,3};
+julia> res_ann = Framewise(alg, delta_rot=(0.1, 1))(mav, angles);
    source

    API/Reference

    HCIToolbox.AnnulusView — Type
    AnnulusView(cube::AbstractArray{T,3};
             inner=0, outer=first(size(parent))/2 + 0.5,
             fill=0)

    Cut out an annulus with inner radius inner and outer radius outer. Values that fall outside of this region will be replaced with fill. This does not copy any data, it is merely a view into the data.

    (::AnnulusView)(asview=false)

    Return the pixels that fall within the annulus as a matrix. This matrix is equivalent to unrolling each frame and then spatially filtering the pixels outside the annulus. If asview is true, the returned values will be a view of the parent array instead of a copy.

    Examples

    julia> ann = AnnulusView(ones(101, 101, 10); inner=5, outer=20);
 
@@ -56,4 +56,4 @@
 
 julia> out ≈ -ann
 true
    HCIToolbox.inverse! — Function
    inverse!(::AnnulusView, out, mat)

    In-place version of inverse that fills out in-place.

    inverse!(::MultiAnnulusView, out, idx, mat)
-inverse!(::MultiAnnulusView, out, mats...)

    In-place version of inverse that fills out with annuli defined by the geometry of the view.

    +inverse!(::MultiAnnulusView, out, mats...)

    In-place version of inverse that fills out with annuli defined by the geometry of the view.

    diff --git a/dev/gettingstarted/index.html b/dev/gettingstarted/index.html index 319e736..d6a2950 100644 --- a/dev/gettingstarted/index.html +++ b/dev/gettingstarted/index.html @@ -24,4 +24,4 @@ # use different algorithms for each annulus N_ann = length(anns.indices) algs = [PCA(10), PCA(9), PCA(8), ...] -res = process(algs, anns, angles)

    Comparison to VIP

    ADI.jl took a lot of ideas from VIP and expanded them using the power of Julia. To begin with, Julia typically has smaller and more self-contained packages, so most of the basic image-processing that is used here is actually written in the HCIToolbox.jl package. In the future, I have plans to incorporate forward-modeling distributions in Firefly.jl, which currently is an active research project.

    Some technical distinctions to VIP

    • Julia is 1-indexed. This means all the positions for apertures, bounds, images, etc. start at 1. This is distinct from 0-based indexing in python, but is equivalent to the indexing in DS9 and IRAF.
    • Julia is column-major. This means all data, when loaded from disk, will appear "transposed" from what is loaded by astropy.io.fits.
    • Julia's std uses the sample statistic (n-1 degrees of freedom) while numpy's std uses the population statistic (n degrees of freedom). This may cause very slight differences in measurements that rely on this.
    • Aperture mapping - many of the Metrics are derived by measuring statistics in an annulus of apertures. In VIP, this ring is not equally distributed- the angle between apertures is based on the exact number of apertures rather than the integral number of apertures that are actually measured. In ADI.jl the angle between apertures is evenly distributed. The same number of pixels are discarded in both packages, but in VIP they all end up in the same region of the image (see this figure).
    • Collapsing - by default VIP collapses a cube by derotating it then finding the median frame. In ADI.jl, the default collapse method is a weighted sum using the inverse of the temporal variance for weighting. This is documented in HCIToolbox.collapse and can be overridden by passing the keyword argument method=median or whichever statistical function you want to use.
    • Image interpolation - by default VIP uses a lanczos4 interpolator from opencv, by default ADI.jl uses a bilinear b-spline interpolator through Interpolations.jl
    • Annular and framewise processing - some of the VIP algorithms allow you to go annulus-by-annulus and optionally filter the frames using parallactic angle thresholds. ADI.jl does not bake these options in using keyword arguments; instead, the geometric filtering is achieved through AnnulusView and MultiAnnulusView. Parallactic angle thresholds are implemented in the Framewise algorithm wrapper. I've separated these techniques because they are fundamentally independent and because it greatly increases the composability of the algorithms.

    The biggest difference, though, is Julia's multiple-dispatch system and how that allows ADI.jl to do more with less code. For example, the GreeDS algorithm was designed explicitly for PCA, but the formalism of it is more generic than that. Rather than hard-coding in PCA, the GreeDS algorithm was written generically, and Julia's multiple-dispatch allows the use of, say, NMF instead of PCA. By making the code generic and modular, ADI.jl enables rapid experimentation with different post-processing algorithms and techniques as well as minimizing the code required to implement a new algorithm and be able to fully use the ADI.jl API.

    Feature comparison

    Some notable libraries for HCI tasks include VIP, pyKLIP, and PynPoint. A table of the feature sets of these packages alongside ADI.jl is presented below. In general VIP offers the most diversity in algorithms and their applications, but not all algorithms are as feature-complete as the PCA implementation. VIP also contains many useful utilities for pre-processing and a pipeline framework. pyKLIP primarily uses the PCA (KLIP) algorithm, but offers many forward modeling implementations. PynPoint has a highly modular pre-processing module that is focused on pipelines.

    -Pre.Algs.Techs.D.M.F.M.
    ADI.jl✗median, LOCI, PCA, NMF, fixed-point GreeDSFull-frame ADI/RDI, SDI (experimental), annular ADI*detection maps, STIM, SLIMask, contrast curve✗
    VIP✓median, LOCI, PCA, NMF, LLSG, ANDROMEDA, pairwise frame differencingFull-frame ADI/RDI, SDI, annular ADI/RDI*detection maps, blob detection, STIM, ROC, contrast curveNegFC
    pyKLIP✗PCA, NMF, weighted PCAFull-frame ADI/RDI, SDI, annular ADI/RDIdetection maps, blob detection, contrast curve, cross-correlationKLIP-FM, Planet Evidence, matched filter (FMMF), spectrum fitting, DiskFM
    PynPoint✓median, PCAFull-frame ADI/RDI, SDIdetection maps, contrast curve✗

    Column labels: Pre-processing, Algorithms, Techniques, Detection Metrics, Forward Modeling.

    Techniques marked with * indicate partial support, meaning that not all algorithms are supported.

    +res = process(algs, anns, angles)

    Comparison to VIP

    ADI.jl took a lot of ideas from VIP and expanded them using the power of Julia. To begin with, Julia typically has smaller and more self-contained packages, so most of the basic image-processing that is used here is actually written in the HCIToolbox.jl package. In the future, I have plans to incorporate forward-modeling distributions in Firefly.jl, which currently is an active research project.

    Some technical distinctions to VIP

    • Julia is 1-indexed. This means all the positions for apertures, bounds, images, etc. start at 1. This is distinct from 0-based indexing in python, but is equivalent to the indexing in DS9 and IRAF.
    • Julia is column-major. This means all data, when loaded from disk, will appear "transposed" from what is loaded by astropy.io.fits.
    • Julia's std uses the sample statistic (n-1 degrees of freedom) while numpy's std uses the population statistic (n degrees of freedom). This may cause very slight differences in measurements that rely on this.
    • Aperture mapping - many of the Metrics are derived by measuring statistics in an annulus of apertures. In VIP, this ring is not equally distributed- the angle between apertures is based on the exact number of apertures rather than the integral number of apertures that are actually measured. In ADI.jl the angle between apertures is evenly distributed. The same number of pixels are discarded in both packages, but in VIP they all end up in the same region of the image (see this figure).
    • Collapsing - by default VIP collapses a cube by derotating it then finding the median frame. In ADI.jl, the default collapse method is a weighted sum using the inverse of the temporal variance for weighting. This is documented in HCIToolbox.collapse and can be overridden by passing the keyword argument method=median or whichever statistical function you want to use.
    • Image interpolation - by default VIP uses a lanczos4 interpolator from opencv, by default ADI.jl uses a bilinear b-spline interpolator through Interpolations.jl
    • Annular and framewise processing - some of the VIP algorithms allow you to go annulus-by-annulus and optionally filter the frames using parallactic angle thresholds. ADI.jl does not bake these options in using keyword arguments; instead, the geometric filtering is achieved through AnnulusView and MultiAnnulusView. Parallactic angle thresholds are implemented in the Framewise algorithm wrapper. I've separated these techniques because they are fundamentally independent and because it greatly increases the composability of the algorithms.

    The biggest difference, though, is Julia's multiple-dispatch system and how that allows ADI.jl to do more with less code. For example, the GreeDS algorithm was designed explicitly for PCA, but the formalism of it is more generic than that. Rather than hard-coding in PCA, the GreeDS algorithm was written generically, and Julia's multiple-dispatch allows the use of, say, NMF instead of PCA. By making the code generic and modular, ADI.jl enables rapid experimentation with different post-processing algorithms and techniques as well as minimizing the code required to implement a new algorithm and be able to fully use the ADI.jl API.

    Feature comparison

    Some notable libraries for HCI tasks include VIP, pyKLIP, and PynPoint. A table of the feature sets of these packages alongside ADI.jl is presented below. In general VIP offers the most diversity in algorithms and their applications, but not all algorithms are as feature-complete as the PCA implementation. VIP also contains many useful utilities for pre-processing and a pipeline framework. pyKLIP primarily uses the PCA (KLIP) algorithm, but offers many forward modeling implementations. PynPoint has a highly modular pre-processing module that is focused on pipelines.

    -Pre.Algs.Techs.D.M.F.M.
    ADI.jl✗median, LOCI, PCA, NMF, fixed-point GreeDSFull-frame ADI/RDI, SDI (experimental), annular ADI*detection maps, STIM, SLIMask, contrast curve✗
    VIP✓median, LOCI, PCA, NMF, LLSG, ANDROMEDA, pairwise frame differencingFull-frame ADI/RDI, SDI, annular ADI/RDI*detection maps, blob detection, STIM, ROC, contrast curveNegFC
    pyKLIP✗PCA, NMF, weighted PCAFull-frame ADI/RDI, SDI, annular ADI/RDIdetection maps, blob detection, contrast curve, cross-correlationKLIP-FM, Planet Evidence, matched filter (FMMF), spectrum fitting, DiskFM
    PynPoint✓median, PCAFull-frame ADI/RDI, SDIdetection maps, contrast curve✗

    Column labels: Pre-processing, Algorithms, Techniques, Detection Metrics, Forward Modeling.

    Techniques marked with * indicate partial support, meaning that not all algorithms are supported.

    diff --git a/dev/index.html b/dev/index.html index 73f1780..1f9c4cb 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,4 +1,4 @@ Home · ADI.jl
    Edit on GitHub

    ADI.jl

    GitHub Build Status Coverage License

    JOSS DOI

    A package for angular differential imaging (ADI) along with its variants, such as reference differential imaging (RDI) and spectral differential imaging (SDI).

    Installation and Setup

    ADI.jl is a registered package and can be installed using the Julia package manager. From the Julia REPL, enter Pkg mode (by pressing ])

    julia>]
 
-(@v1.5) pkg> add ADI

    To exit Pkg mode, just backspace. Once the package is installed it can be imported with

    julia> using ADI

    For more information, see the Pkg documentation.

    Citations

    If you use ADI.jl or derivatives in your work, please consider citing both the JOSS paper and the code record. The JOSS paper citation can be found in CITATION.bib. The code will have a unique reference for each released version, so visit the Zenodo record to grab the BibTeX for whichever version you used.

    Contributing and Support

    ColPrac: Contributor's Guide on Collaborative Practices for Community Packages

    In general contributions should follow ColPrac. If you are interested in extending/improving ADI.jl, head to the discussions to reach out. For support with using ADI.jl, please open an issue describing the problem and steps to reproduce it.

    License

    This work is distributed under the MIT "expat" license. See LICENSE for more information.

    +(@v1.5) pkg> add ADI

    To exit Pkg mode, just backspace. Once the package is installed it can be imported with

    julia> using ADI

    For more information, see the Pkg documentation.

    Citations

    If you use ADI.jl or derivatives in your work, please consider citing both the JOSS paper and the code record. The JOSS paper citation can be found in CITATION.bib. The code will have a unique reference for each released version, so visit the Zenodo record to grab the BibTeX for whichever version you used.

    Contributing and Support

    ColPrac: Contributor's Guide on Collaborative Practices for Community Packages

    In general contributions should follow ColPrac. If you are interested in extending/improving ADI.jl, head to the discussions to reach out. For support with using ADI.jl, please open an issue describing the problem and steps to reproduce it.

    License

    This work is distributed under the MIT "expat" license. See LICENSE for more information.

    diff --git a/dev/introduction/index.html b/dev/introduction/index.html index bf3460a..8ffbdf5 100644 --- a/dev/introduction/index.html +++ b/dev/introduction/index.html @@ -1,2 +1,2 @@ -Introduction to HCI · ADI.jl
    Edit on GitHub

    Introduction to High-Contrast Imaging

    This will serve as a brief primer to high-contrast imaging (HCI) to illustrate the concepts and themes prevalent within this package. This is not meant to be an exhaustive lecture note on the topic, but more of a gentle introduction.

    If you are comfortable with HCI topics, consider skipping to the Getting Started section to get introduced to ADI.jl, browse the Examples to see sample workflows, or browse through the API to start exploring the capabilities of ADI.jl.

    What is HCI

    HCI is an advanced imaging technique comprising modern instrumentation, clever observational techniques, and post-processing algorithms. The goal of HCI is to probe the circumstellar regions of a star in the search for companions, debris disks, and more. The use of large aperture telescopes, adaptive optics (AO), and coronagraphs are the basis of HCI. An example of an image taken with these instruments is shown below.

    What is notable in this image is that there is a lot of structured noise in the image that is overwhelming any potential companion signal. The center of the image is particularly noisy, which is precisely where we are most interested in searching for exoplanets. This noise is the effect of quasi-static speckles in the focal-plane. These speckles occur from non-common path aberrations in the AO system and are a fundamental part of the data received by these instruments. Improving the quality of instrumentation is an active topic of HCI research, but it is beyond the scope of this introduction.

    Angular Differential Imaging (ADI)

    Because this noise is fundamental to the data, post-processing must be performed in order to see any circumstellar signal. This post-processing requires us to fit the PSF of the speckles and then remove it. An example of the above frame with the speckles removed is shown below.

    Unfortunately, there is no companion evident; but the speckles have been removed, so what is left? The exoplanet is still sitting below the statistical noise in this frame, but the noise can be averaged out by combining many frames together. Since we are concerned with subtracting the speckles, we need to be careful and consider how do we fit and subtract the speckles without removing potential companion signal?

    This is where angular differential imaging (ADI) comes in. ADI is an observational technique pioneered in the early 2000s as an extension of roll deconvolution for ground-based telescopes. The core of this process is that the quasi-static speckles are a function of the optical system, not the astrophysical source. Throughout a night of observing we can leverage the rotation of the Earth to make the field-of-view (FOV) appear to rotate (on an Alt-Az mounted telescope with the field rotator disabled). Even though the sky appears to rotate, because the speckles are due to the telescope optics they will not appear to rotate. The animation below shows a cube of data with a bright fake companion that illustrates the sky rotation typical of ADI.

    By taking this sequence of images (commonly referred to as a cube) we can more easily model and fit the speckle signal separate from any companions. If you median combine the cube as-is, the non-stationary companion signal will attenuate leaving just the speckles. If we derotate the sequence according to the parallactic angles for each frame we align the sky to a common heading. Now we can collapse the derotated sequence and the planet will constructively interfere while the now-rotating speckles will attenuate. The figure below shows these two competing reductions.

    Post-Processing Algorithms

    Using data cubes (as described in the ADI section), we are tasked with fitting the speckles without capturing the rotating companion signal. Quite a few algorithms have been proposed and a thorough discussion of them is beyond the scope of this introduction. For now, let's assume the algorithms are a black-box that produce speckle approximation cubes.

    If we have this cube, all we need to post-process the data is

    1. Retrieve a speckle estimate cube
    2. Subtract the speckle estimate from the target cube and form a residual cube
    3. Derotate the residual cube according to the parallactic angles of the target
    4. Collapse the derotated residual cube

    Steps 2-4 are shown in the following figure

    After all this processing, finally the substellar companion HR8799e is evident. Hopefully this shows the difficulty of HCI and builds up part of the process that occurs outside of the reduction you'll be doing with ADI.jl.

    References

    Here is a selection of further reading for information about high-contrast imaging, ADI, and similar techniques

    +Introduction to HCI · ADI.jl
    Edit on GitHub

    Introduction to High-Contrast Imaging

    This will serve as a brief primer to high-contrast imaging (HCI) to illustrate the concepts and themes prevalent within this package. This is not meant to be an exhaustive lecture note on the topic, but more of a gentle introduction.

    If you are comfortable with HCI topics, consider skipping to the Getting Started section to get introduced to ADI.jl, browse the Examples to see sample workflows, or browse through the API to start exploring the capabilities of ADI.jl.

    What is HCI

    HCI is an advanced imaging technique comprising modern instrumentation, clever observational techniques, and post-processing algorithms. The goal of HCI is to probe the circumstellar regions of a star in the search for companions, debris disks, and more. The use of large aperture telescopes, adaptive optics (AO), and coronagraphs are the basis of HCI. An example of an image taken with these instruments is shown below.

    What is notable in this image is that there is a lot of structured noise in the image that is overwhelming any potential companion signal. The center of the image is particularly noisy, which is precisely where we are most interested in searching for exoplanets. This noise is the effect of quasi-static speckles in the focal-plane. These speckles occur from non-common path aberrations in the AO system and are a fundamental part of the data received by these instruments. Improving the quality of instrumentation is an active topic of HCI research, but it is beyond the scope of this introduction.

    Angular Differential Imaging (ADI)

    Because this noise is fundamental to the data, post-processing must be performed in order to see any circumstellar signal. This post-processing requires us to fit the PSF of the speckles and then remove it. An example of the above frame with the speckles removed is shown below.

    Unfortunately, there is no companion evident; but the speckles have been removed, so what is left? The exoplanet is still sitting below the statistical noise in this frame, but the noise can be averaged out by combining many frames together. Since we are concerned with subtracting the speckles, we need to be careful and consider how do we fit and subtract the speckles without removing potential companion signal?

    This is where angular differential imaging (ADI) comes in. ADI is an observational technique pioneered in the early 2000s as an extension of roll deconvolution for ground-based telescopes. The core of this process is that the quasi-static speckles are a function of the optical system, not the astrophysical source. Throughout a night of observing we can leverage the rotation of the Earth to make the field-of-view (FOV) appear to rotate (on an Alt-Az mounted telescope with the field rotator disabled). Even though the sky appears to rotate, because the speckles are due to the telescope optics they will not appear to rotate. The animation below shows a cube of data with a bright fake companion that illustrates the sky rotation typical of ADI.

    By taking this sequence of images (commonly referred to as a cube) we can more easily model and fit the speckle signal separate from any companions. If you median combine the cube as-is, the non-stationary companion signal will attenuate leaving just the speckles. If we derotate the sequence according to the parallactic angles for each frame we align the sky to a common heading. Now we can collapse the derotated sequence and the planet will constructively interfere while the now-rotating speckles will attenuate. The figure below shows these two competing reductions.

    Post-Processing Algorithms

    Using data cubes (as described in the ADI section), we are tasked with fitting the speckles without capturing the rotating companion signal. Quite a few algorithms have been proposed and a thorough discussion of them is beyond the scope of this introduction. For now, let's assume the algorithms are a black-box that produce speckle approximation cubes.

    If we have this cube, all we need to post-process the data is

    1. Retrieve a speckle estimate cube
    2. Subtract the speckle estimate from the target cube and form a residual cube
    3. Derotate the residual cube according to the parallactic angles of the target
    4. Collapse the derotated residual cube

    Steps 2-4 are shown in the following figure

    After all this processing, finally the substellar companion HR8799e is evident. Hopefully this shows the difficulty of HCI and builds up part of the process that occurs outside of the reduction you'll be doing with ADI.jl.

    References

    Here is a selection of further reading for information about high-contrast imaging, ADI, and similar techniques

    diff --git a/dev/metrics/index.html b/dev/metrics/index.html index 1eec444..4ee8c64 100644 --- a/dev/metrics/index.html +++ b/dev/metrics/index.html @@ -1,15 +1,15 @@ -Metrics · ADI.jl
    Edit on GitHub

    Metrics

    ADI.Metrics — Module
    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.

    source

    Detection maps

    ADI.Metrics.detectionmap — Function
    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.
    Tip

    This code is automatically multi-threaded, so be sure to set JULIA_NUM_THREADS before loading your runtime to take advantage of it!

    source
    ADI.Metrics.snr — Function
    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.

    Note

    SNR is not equivalent to significance, use significance instead

    source
    ADI.Metrics.significance — Function
    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

    source
    ADI.Metrics.noise — Function
    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.

    source
    ADI.Metrics.stimmap — Function
    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
+Metrics · ADI.jl
    Edit on GitHub
    Metrics
    ADI.Metrics — Module
    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.
    source
    Detection maps
    ADI.Metrics.detectionmap — Function
    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.
    Tip
    This code is automatically multi-threaded, so be sure to set JULIA_NUM_THREADS before loading your runtime to take advantage of it!
    source
    ADI.Metrics.snr — Function
    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.
    Note
    SNR is not equivalent to significance, use significance instead
    source
    ADI.Metrics.significance — Function
    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
    source
    ADI.Metrics.noise — Function
    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.
    source
    ADI.Metrics.stimmap — Function
    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
    source
    ADI.Metrics.stim_threshold — Function
    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> 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
    source
    ADI.Metrics.stim_threshold — Function
    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
    source
    ADI.Metrics.stim — Function
    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.
    source
    Ensemble methods
    The following methods utilize the results of multiple ADI reductions, in some form.
    ADI.Metrics.slimmap — Function
    slimmap(residuals, angles; N)
+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
    source
    ADI.Metrics.stim — Function
    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.
    source
    Ensemble methods
    The following methods utilize the results of multiple ADI reductions, in some form.
    ADI.Metrics.slimmap — Function
    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);
@@ -18,19 +18,19 @@
 
 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
    source
    Throughput
    Throughput
    ADI.Metrics.throughput — Function
    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.
    source
    throughput(alg, cube, angles, psf, position;
+           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.
    source
    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
    source
    throughput(alg, sdi_cube, angles, scales, psf, position;
+           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
    source
    throughput(alg, sdi_cube, angles, scales, psf, position;
            fwhm, snr=100, flux_scale=1,
            reduced_empty=nothing,
-           verbose=true, kwargs...)
    Calculate throughput for an SDI cube with radially scaling list scales and flux scaling list flux_scale.
    source
    Contrast curve
    ADI.Metrics.contrast_curve — Function
    contrast_curve(alg, cube, angles, psf;
+           verbose=true, kwargs...)
    Calculate throughput for an SDI cube with radially scaling list scales and flux scaling list flux_scale.
    source
    Contrast curve
    ADI.Metrics.contrast_curve — Function
    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.
    Tip
    If you prefer a tabular format, simply pipe the output of this function into any type supporting the Tables.jl interface, e.g.
    contrast_curve(alg, cube, angles, psf; fwhm=fwhm) |> DataFrame
    source
    ADI.Metrics.subsample_contrast — Function
    Metrics.subsample_contrast(empty_frame, distance, throughput;
+               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.
    Tip
    If you prefer a tabular format, simply pipe the output of this function into any type supporting the Tables.jl interface, e.g.
    contrast_curve(alg, cube, angles, psf; fwhm=fwhm) |> DataFrame
    source
    ADI.Metrics.subsample_contrast — Function
    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
 
@@ -39,5 +39,5 @@
 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)
    source
    ADI.Metrics.estimate_starphot — Function
    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.
    source
    
+                                    fwhm=8.4, starphot=sp)
    source
    ADI.Metrics.estimate_starphot — Function
    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.

    source
    diff --git a/dev/sdi/index.html b/dev/sdi/index.html index e6fd7f5..8650d20 100644 --- a/dev/sdi/index.html +++ b/dev/sdi/index.html @@ -1,16 +1,16 @@ -SDI · ADI.jl
    Edit on GitHub

    SDI

    Warning

    SDI should be considered very experimental.

    ADI.SDIAlgorithm — Type
    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

    source

    API/Reference

    ADI.SingleSDI — Type
    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.

    Tip

    SingleSDI is the default SDI mode. This means instead of writing

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

    you can just write

    PCA(15)(data, angles, scales)
    source
    ADI.DoubleSDI — Type
    DoubleSDI(alg)
+SDI · ADI.jl
    Edit on GitHub
    SDI
    Warning
    SDI should be considered very experimental.
    ADI.SDIAlgorithm — Type
    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
    source
    API/Reference
    ADI.SingleSDI — Type
    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.
    Tip
    SingleSDI is the default SDI mode. This means instead of writing
    SingleSDI(PCA(15))(data, angles, scales)
    you can just write
    PCA(15)(data, angles, scales)
    source
    ADI.DoubleSDI — Type
    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)
    source
    ADI.SliceSDI — Type
    SliceSDI(alg)
+julia> res = DoubleSDI(Classic(), GreeDS(15))(data, angles, scales)
    source
    ADI.SliceSDI — Type
    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)
    source
    HCIToolbox.scale_list — Function
    scale_list(wavelengths)
    Returns a list of scaling factors for aligning SDI tensors from a list of wavelengths.
    Examples
    julia> scale_list([0.5, 2, 4])
+julia> res = SliceSDI(Classic(), GreeDS(15))(data, angles, scales)
    source
    HCIToolbox.scale_list — Function
    scale_list(wavelengths)
    Returns a list of scaling factors for aligning SDI tensors from a list of wavelengths.
    Examples
    julia> scale_list([0.5, 2, 4])
 3-element Vector{Float64}:
  8.0
  2.0
- 1.0
    
+ 1.0
    diff --git a/dev/search/index.html b/dev/search/index.html index 52c8ceb..8b6cd00 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -Search · ADI.jl

    Loading search...

      +Search · ADI.jl

      Loading search...