3. PSF modeling and subtraction¶

Authors: Carlos Alberto Gomez Gonzalez and Valentin Christiaens

Suitable for VIP v1.0.0 onwards

Last update: 2022/05/04

Table of contents

This tutorial shows:

how to load ADI-ready datacubes;
how to use the stellar PSF subtraction algorithms implemented in VIP to produce final post-processed images (more details and higher completeness than the quick-start tutorial).

Let’s first import a couple of external packages needed in this tutorial:

[1]:

%matplotlib inline
from hciplot import plot_frames, plot_cubes
from matplotlib.pyplot import *
from matplotlib import pyplot as plt
from multiprocessing import cpu_count
import numpy as np
from packaging import version

In the following box we import all the VIP routines that will be used in this tutorial. The path to some routines has changed between versions 1.0.3 and 1.1.0, which saw a major revamp of the modular architecture, hence the if statements.

[2]:

import vip_hci as vip
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) < version.parse("1.0.0"):
    msg = "Please upgrade your version of VIP"
    msg+= "It should be 1.0.0 or above to run this notebook."
    raise ValueError(msg)
elif version.parse(vvip) <= version.parse("1.0.3"):
    from vip_hci.andromeda import andromeda
    from vip_hci.conf import VLT_NACO
    from vip_hci.frdiff import frame_diff
    from vip_hci.leastsq import xloci
    from vip_hci.llsg import llsg
    from vip_hci.medsub import median_sub
    from vip_hci.metrics import normalize_psf
    from vip_hci.metrics import compute_stim_map as stim_map
    from vip_hci.metrics import compute_inverse_stim_map as inverse_stim_map
    from vip_hci.nmf import nmf, nmf_annular
    from vip_hci.pca import pca, pca_annular, pca_annulus, pca_grid
else:
    if version.parse(vvip) >= version.parse("1.2.4"):
        from vip_hci.invprob import fmmf, FastPACO
    from vip_hci.config import VLT_NACO
    from vip_hci.fm import normalize_psf
    from vip_hci.invprob import andromeda
    from vip_hci.metrics import inverse_stim_map, stim_map
    from vip_hci.psfsub import (frame_diff, llsg, median_sub, nmf, nmf_annular,
                                pca, pca_annular, pca_annulus, pca_grid, xloci)

# common to all versions:
from vip_hci.fits import open_fits, write_fits, info_fits
from vip_hci.metrics import significance, snr, snrmap
from vip_hci.var import fit_2dgaussian, frame_center

VIP version:  1.2.4

3.1. Loading ADI data¶

In the ‘dataset’ folder of the VIP_extras repository you can find a toy ADI (Angular Differential Imaging) cube and a NACO point spread function (PSF) to demonstrate the capabilities of VIP. This is an L’-band VLT/NACO dataset of beta Pictoris published in Absil et al. (2013) obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells.

More info on this dataset, and on opening and visualizing fits files with VIP in general, is available in Tutorial 1. Quick start.

Let’s load the datacube, associated parallactic angles and non-coronagraphic PSF:

[3]:

cubename = '../datasets/naco_betapic_cube_cen.fits'
angname = '../datasets/naco_betapic_pa.fits'
psfnaco = '../datasets/naco_betapic_psf.fits'

cube = open_fits(cubename)
pa = open_fits(angname)
psf = open_fits(psfnaco)

Fits HDU-0 data successfully loaded. Data shape: (61, 101, 101)
Fits HDU-0 data successfully loaded. Data shape: (61,)
Fits HDU-0 data successfully loaded. Data shape: (39, 39)

Question 1.1: When observing a celestial object with a given telecope, what 3 parameters do the values of parallactic angle depend on?

Note: different questions will be raised throughout this notebook. Corresponding answers will be provided at the end of the respective (sub)sections.

[4]:

derot_off = 104.84 # NACO derotator offset for this observation (Absil et al. 2013)
TN = -0.45         # Position angle of true north for NACO at the epoch of observation (Absil et al. 2013)

angs = pa-derot_off-TN

Let’s measure the Full-Width at Half Maximum (FWHM) by fitting a 2D Gaussian to the core of the unsaturated non-coronagraphic PSF:

[5]:

%matplotlib inline
DF_fit = fit_2dgaussian(psf, crop=True, cropsize=9, debug=True, full_output=True)

FWHM_y = 4.733218722257407
FWHM_x = 4.473682405059958

centroid y = 19.006680059041216
centroid x = 18.999424475165455
centroid y subim = 4.006680059041214
centroid x subim = 3.9994244751654535

amplitude = 0.10413004853269707
theta = -34.08563676836685

[6]:

fwhm_naco = np.mean([DF_fit['fwhm_x'],DF_fit['fwhm_y']])
print(fwhm_naco)

4.603450563658683

Question 3.2: The FWHM measured above can slightly vary depending on observing conditions and quality of the adaptive optics (AO) correction. In excellent conditions, the PSF gets narrower and can approach the theoretical diffraction limit. What is that angular resolution (in arcsec and in pixels) at L’ band, for a 8.1m diameter telescope such as the VLT?

Let’s normalize the flux to one in a 1xFWHM aperture and crop the PSF array:

[7]:

psfn = normalize_psf(psf, fwhm_naco, size=19, imlib='ndimage-fourier')

Flux in 1xFWHM aperture: 1.228

The normalize_psf function performs internally a fine centering of the PSF (note the input PSF should already be centered within a few pixel accuracy for this fine centering to work). The imlib argument sets the library to be used for sub-px shifts (more details in Tutorial 7).

Let’s visualize the normalized PSF with hciplot.plot_frames. Feel free to adapt the backend argument throughout the notebook: 'matplotlib' (default) allows paper-quality figures with annotations which can be saved (default), while 'bokeh' enables interactive visualization.

[8]:

plot_frames(psfn, grid=True, size_factor=4)

Let’s finally define the pixel scale for NACO (L’ band), which we get from a dictionary stored in the conf subpackage:

[9]:

pxscale_naco = VLT_NACO['plsc']
print(pxscale_naco, "arcsec/px")

0.02719 arcsec/px

Question 3.3: When building an instrument, one wants to make sure that even in the best observing conditions the PSF is always sampled at least at the Nyquist frequency - what would then be the coarsest pixel scale that would still allow a Nyquist sampling of the PSF at L’ at the VLT?

Answer 3.1: The parallactic angle depends on the celestial coordinates of the target, the latitude of the observatory and the observation time (commonly expressed in local hour angle for the target).

Answer 3.2: The theoretical diffraction limit is $\sim \lambda/D$, or $\sim$ 0.10 arcsec after conversion from radian to arcsec. This is $:nbsphinx-math:`sim`$3.9 pixels at the NACO pixel scale.

Answer 3.3: A maximum pixel scale of 0.05 arcsec/px is required to ensure Nyquist sampling at L’ band even in the best observing conditions.

3.2. median-ADI¶

In the framework of differential imaging techniques, we ultimately rely on modelling and subtracting the stellar PSF and associated speckle noise pattern. Algorithms of different complexities and performances have been proposed since 2006. Several of those algorithms are implemented in VIP and showcased in the next subsections.

When loading the ADI cube above, we assumed all calibration and preprocessing steps were already performed. In particular, the star is assumed to be already centered in the images. If your data require calibration/preprocessing, feel free to check tutorial 2. Pre-processing for example uses of VIP preprocessing routines on NACO and SPHERE/IFS data.

IMPORTANT - VIP’s convention regarding centering is:

for odd number of pixels along the x and y directions: the star is placed on the central pixel;
for even number of pixels: the star is placed on coordinates (dim/2,dim/2) in 0-based indexing.

All ADI-based algorithms involve image rotation. We therefore first set the preferred method to be used for image rotation. This is controlled with the imlib and interpolation optional arguments (see tutorial 7. FFT- vs. interpolation-based image operations for more details). Feel free to adapt the following box to test other methods (keeping in mind that FFT-based rotation is only available in versions >= 1.0.0).

[10]:

imlib_rot = 'vip-fft'
interpolation=None

3.2.1. Full-frame median-ADI¶

The most simple approach is to model the PSF with the median of the ADI sequence. The median subtraction algorithm is the original post-processing approach proposed for ADI data (see Marois et al. 2006):

Figure 3.1 - Median-ADI concept. A sequence of images $A_i$ is obtained in pupil-tracking mode. Any putative planet would rotate throughout the images, along with the field of view, while the stellar halo and speckle pattern remains quasi-static. The median $B_i$ of the sequence does not compare and is subtracted to $A_i$ to create a residual cube $C_i$. The residual cube is derotated ($D_i$) to align the field-of-view. The derotated residual cube is then stacked (typically with a median) to create a finale image $E$.

Question 3.4: In practice, what are the main limitations of median-ADI?

[11]:

fr_adi = median_sub(cube, angs, mode='fullfr', imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:28:50
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time:  0:00:01.102045
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

3.2.2. Smart median-ADI¶

A “smart median-ADI” can be performed in concentric annuli and using a parallactic angle threshold. For each image of the input cube the subtracted model consists of the median of the nframes closest images in time that have rotated by at least delta_rotx FWHM (linear motion) - i.e. it is a different model PSF for each image. This model is constructed in concentric annuli of width asizex FWHM, starting at radius_int pixels radius (i.e. the inner radius_int pixels are masked).

Question 3.5: What are the pros and cons of using a parallactic angle threshold?

[12]:

fr_adi_an = median_sub(cube, angs, fwhm_naco, asize=fwhm_naco, mode='annular', delta_rot=0.5,
                       radius_int=4, nframes=4, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:28:51
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 10, FWHM = 4
Processing annuli:
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:00:00
Optimized median psf reference subtracted
Done derotating and combining
Running time:  0:00:01.274225
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Note that the parallactic angle threshold corresponding to 0.5 FWHM linear motion could not be met for the innermost annulus, and was therefore automatically relaxed to find nframes eligible frames to be used for the model.

Let’s vizualize the final images:

[13]:

plot_frames((fr_adi, fr_adi_an))

Strong residuals can be seen near the center of the image, with a faint potential point source in the north-east of the final image:

Answer 3.4 Quasi-static speckles (i.e. not exactly static) are one of the main limitations. These stem from slowly varying aberrations in the instrument and optics while tracking the object during the observation (due e.g. to different mechanical flexures, imperfections in the mirrors, non-common path aberrations in the instrument). If observing conditions are not stable (e.g. sudden bursts in seeing), the adaptive optics performance may also fluctuate during the observing sequence, which

can also lead to a decorrelation of the speckle pattern. The latter should be mitigated with bad frame removal (done in preprocessing).

Answer 3.5 A non-zero parallactic angle threshold helps minimizing self-subtraction of putative planets, by only bulding a model for subtraction using images where a putative planet has sufficiently rotated. On the other hand, a large parallactic angle threshold can force the use of images taken far away in time from the one considered for modeling. The decorrelation of quasi-static speckle can then lead to a poor PSF model for subraction. While a value of delta_rot=1 will lead to low

self-subtraction and is a good default value for observations obtained in good and stable conditions with sufficient field rotation, keep in mind that the optimal value of delta_rot for a specific dataset and given planet radial separation may be different than 1.

3.3. Pairwise frame difference¶

Another simple approach is to perform pairwise frame subtraction, taking into account both a rotation and distance (similarity) threshold. For this example, we select the L1-distance metric, with a distance threshold set to the 90th percentile of all distances, a 0.5xFWHM linear rotation threshold, an innermost radius starting at 2 px, and an annular size set to fwhm_naco pixels.

[14]:

fr_fdiff = frame_diff(cube, angs, fwhm_naco, metric='l1', dist_threshold=90, delta_rot=0.5,
                      radius_int=4, asize=fwhm_naco, nproc=None, imlib=imlib_rot,
                      interpolation=interpolation, verbose=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:28:52
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
10 annuli. Performing pair-wise subtraction:
Ann 2    PA thresh: 12.05    Ann center:  11    N segments: 1 Ann 3    PA thresh:  8.49    Ann center:  16    N segments: 1 Ann 4    PA thresh:  6.55    Ann center:  20    N segments: 1 Ann 1    PA thresh: 20.70    Ann center:   6    N segments: 1 Ann 5    PA thresh:  5.33    Ann center:  25    N segments: 1




Ann 6    PA thresh:  4.50    Ann center:  29    N segments: 1
Ann 7    PA thresh:  3.89    Ann center:  34    N segments: 1
Ann 8    PA thresh:  3.42    Ann center:  39    N segments: 1
Ann 9    PA thresh:  3.06    Ann center:  43    N segments: 1
Ann 10    PA thresh:  2.82    Ann center:  47    N segments: 1
Done processing annuli
Running time:  0:00:05.354684
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s vizualize the final image:

[15]:

plot_frames(fr_fdiff, vmin=-5)

As for the median-ADI results, there seems to be a faint point-like source towards the north-east, although strong residuals dominate the center of the image.

3.4. Least-squares approximation (LOCI)¶

A continuation of the idea of imposing a rotation/PA threshold was proposed by Lafreniere et al. (2007), but this time working in annular patches and using a more sophisticated model PSF. In this case the model PSF is built as a Locally Optimized linear Combination of Images (LOCI; after discarding the ones that have not rotated enough). The coefficients of the linear combination are found to minimize local residuals in each annular subsection (i.e. locally optimized).

[16]:

fr_lstsq = xloci(cube, angs, fwhm=fwhm_naco, asize=fwhm_naco, n_segments='auto', nproc=None,
                 metric='correlation', dist_threshold=90, delta_rot=0.1, optim_scale_fact=2,
                 solver='lstsq', tol=0.01, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:28:58
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Building 10 annuli:
Ann 1    PA thresh: 11.42    Ann center:   2    N segments: 2
Ann 2    PA thresh:  3.82    Ann center:   7    N segments: 3
Ann 3    PA thresh:  2.29    Ann center:  12    N segments: 5
Ann 4    PA thresh:  1.64    Ann center:  16    N segments: 7
Ann 5    PA thresh:  1.27    Ann center:  21    N segments: 9
Ann 6    PA thresh:  1.04    Ann center:  25    N segments: 11
Ann 7    PA thresh:  0.88    Ann center:  30    N segments: 12
Ann 8    PA thresh:  0.76    Ann center:  35    N segments: 14
Ann 9    PA thresh:  0.67    Ann center:  39    N segments: 16
Ann 10    PA thresh:  0.62    Ann center:  43    N segments: 18
Patch-wise least-square combination and subtraction: with 5 processes
Done processing annuli
Running time:  0:00:12.250467
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s vizualize the final image:

[17]:

plot_frames(fr_lstsq)

This time the residuals are better subtracted. But the companion candidate also suffers significant self-subtraction. Note as well the comparatively longer processing time.

3.5. Principal Component Analysis (PCA)¶

3.5.1. Full-frame PCA¶

Now let’s try the Principal Component Analysis (PCA)-based algorithms (Amara & Quanz 2012, Soummer et al. 2012) in vip.pca. As a reminder, below is a schematics summarizing how full-frame PCA-ADI works:

Figure 3.2 - PCA-ADI concept. Compared to median-ADI (Fig. 3.1), only the model used for subtraction is refined. It is constructed from the projection of each image on the first $n_{\rm pc}$ principal components. The principal components are typically found by singular value decomposition, after converting the observed sequence of images $A_i$ into a 2D matrix (rows: time dimension, columns: linearized version of the images).

Let’s set to 5 the number of principal components ncomp considered for model creation, and let’s run the PCA algorithm (feel free to test other values):

[18]:

ncomp = 5
fr_pca1 = pca(cube, angs, ncomp=ncomp, mask_center_px=None, imlib=imlib_rot, interpolation=interpolation,
              svd_mode='lapack')

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:10
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 1.139 GB
Done vectorizing the frames. Matrix shape: (61, 10201)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.039818
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:01.136650
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

The svd_mode argument sets how the Singular Value Decomposition should be calculated. ‘lapack’ is the default mode (exact SVD calculation by CPU through numpy), although other options may get faster results depending on input matrix and your machine (see description of svd_mode argument).

Let’s visualize the final image. This time let’s set the grid option on, to better read coordinates from the image:

[19]:

plot_frames(fr_pca1, grid=True)

The improvement is clear compared to previous reductions. Very low residuals are seen near the star, and the companion candidate appears as a bright point source surrounded by negative side-lobes. The procedure is also 15-20x faster than LOCI.

Question 3.6: What is the origin of the negative side lobes?

Question 3.7: Are they always expected to be symmetric?

Now let’s see whether the candidate is significant. For that let’s first set the x,y coordinates of the test point source based on the above image:

[20]:

xy_b = (58.5, 35.5)

Now let’s compute the signal-to-noise ratio at that location, following the definition in Mawet et al. (2014).

[21]:

snr1 = snr(fr_pca1, source_xy=xy_b, fwhm=fwhm_naco)
print(r"S/N = {:.1f}".format(snr1))

S/N = 8.1

One can also calculate a S/N ratio map over the whole image (this may take a couple of seconds depending on the size of your image):

[22]:

snrmap1 = snrmap(fr_pca1, fwhm_naco, plot=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:11
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

S/N map created using 5 processes
Running time:  0:00:02.721794
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Remember that S/N ratio is NOT the same as significance Mawet et al. 2014. Let’s convert the measured S/N (Students statistics) into a Gaussian “sigma” with equivalent false alarm probability. This first involves calculating the radial separation of the candidate:

[23]:

cy, cx = frame_center(snrmap1)
rad = np.sqrt((cy-xy_b[1])**2+(cx-xy_b[0])**2)

Now let’s use the significance routine in VIP to operate the conversion:

[24]:

sig1 = significance(snr1, rad, fwhm_naco, student_to_gauss=True)
print(r"{:.1f} sigma detection".format(sig1))

5.4 sigma detection

Congratulations! The detection is significant! You obtained a conspicuous direct image of an exoplanet! Although from this image (at a single wavelength) alone one cannot disentangle a true physically bound companion from a background star, this source has now been extensively studied: it corresponds to the giant planet Beta Pic b (e.g. Lagrange et al. 2009, Absil et al. 2013).

3.5.2. Optimizing the number of PCs for full-frame PCA-ADI¶

Clearly, modeling the background (leaked starlight + static and quasi-static speckles) with PCA leads to better images than subtracting a median frame. We chose 5 principal components (PCs), which greatly reduced the residual noise. A different number of PCs may lead to an even better detection.

For a given point source, the pca function can search for the optimal number of principal components which maximizes the S/N ratio of a given point source. When calling pca, this is done by setting ncomp to a tuple of 3 integers (initial value, last value, step) setting the range of explored npc values, and providing the source_xy coordinates of the point source of interest:

[25]:

fr_pca2 = pca(cube, angs, fwhm=fwhm_naco, source_xy=xy_b, mask_center_px=None,
              ncomp=(1, 61, 2), imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:14
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 1.112 GB
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026285
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 31
Optimal number of PCs = 13, for S/N=9.028
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 58.5, 35.5
Flux in a centered 1xFWHM circular aperture = 114.501
Central pixel S/N = 11.321
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 9.028
Max S/N (shifting the aperture center) = 11.526
stddev S/N (shifting the aperture center) = 2.167

The optimal number of principal components is found to be 13 (note this value may change slightly if using a different imlib/interpolation, but should be in that ballpark).

[26]:

ncomp_opt=13

Question 3.8: Why is the shape of the S/N vs Principal components curve as seen above (i.e. increasing and then decreasing again)?

Now let’s look at the images computed using 5 PCs and 13 PCs (optimal npc found above), and compute the S/N maps:

[27]:

snrmap1 = snrmap(fr_pca1, fwhm_naco, plot=False)
snrmap2 = snrmap(fr_pca2, fwhm_naco, plot=False)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:50
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
S/N map created using 5 processes
Running time:  0:00:02.830185
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:53
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
S/N map created using 5 processes
Running time:  0:00:02.761397
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[28]:

plot_frames((fr_pca1, fr_pca2, snrmap1, snrmap2), rows=2, dpi=100, colorbar=True, axis=False,
            versp=0.01, horsp=0.2,
            label=('{} PCs'.format(ncomp), '{} PCs'.format(ncomp_opt),
                   'S/N map with {} PCs'.format(ncomp), 'S/N map with {} PCs'.format(ncomp_opt)),
            label_pad=8)

With 13 PCs we’ve increased the S/N of the point source (although the recovered flux is smaller). Let’s compare the difference it makes on the significance:

[29]:

snr2 = snr(fr_pca2, source_xy=xy_b, fwhm=fwhm_naco)
sig2 = significance(snr2, rad, fwhm_naco, student_to_gauss=True)

print(r"S/N = {:.1f} ({} PCs) VS. S/N = {:.1f} ({} PCs)".format(snr1, ncomp, snr2, ncomp_opt))
print(r"{:.1f} sigma detection ({} PCs) VS. {:.1f} sigma detection ({} PCs)".format(sig1, ncomp, sig2, ncomp_opt))

S/N = 8.1 (5 PCs) VS. S/N = 11.3 (13 PCs)
5.4 sigma detection (5 PCs) VS. 6.3 sigma detection (13 PCs)

Question 3.9: The result above was obtained for the full-frame PCA algorithm. For a given point source (i.e. at a given radius), is the optimal npc expected to be the same for PCA applied in a single annulus (encompassing the point source) and PCA on full frame? If not, would you intuitively expect it to be lower or higher for PCA on a single annulus?

The answer to Q5.6 will come in Section 3.5.6.

Answer 3.6: The negative side lobes stem from the presence of the point source in the constructed model used for subtraction. Since the field of view is rotating, the point source will span a range of position angles, and the stellar halo model built for subtraction for each frame of the cube (step B) will typically include some signal from the rotating point-source.

Answer 3.7: The negative side lobes are not necessarily symmetric, as the field rotation rate (hence the rate at which the point source rotates) varies during the observation. The parallactic angle variation is the fastest when the target crosses the meridian. That being said, well-designed ADI observations typically aim for maximal field rotation, which is the case for an observation that is centered around meridian crossing, hence typically resulting in relatively symmetric negative

signatures.

Answer 3.8: The S/N first increases with $n_{\rm pc}$ due to a better modeling of the speckle pattern by including more PCs. After reaching the optimal $n_{\rm pc}$, the planet self-subtraction starts to dominate (i.e. a large $n_{\rm pc}$ will systematically lead to overfitting of the stellar PSF).

3.5.3. Full-frame PCA-ADI with a parallactic angle threshold¶

We can partially avoid the companion self-subtraction with full-frame PCA by applying a PA threshold for a given distance from the center. This is more efficient if a lot of rotation is present in the data.

This dataset happens to have a good range of rotation (about ~80 deg), and is therefore suitable for that:

[30]:

figure(figsize=(9,3))
plot(angs, 'o-', alpha=0.5)
ylabel('Derotation angles')

[30]:

Text(0, 0.5, 'Derotation angles')

../_images/tutorials_03_psfsub_102_1.png

Question 3.10 For a given observatory, where are the most ideal targets (in terms of celestial coordinates) to apply the ADI strategy (i.e. the ones leading to an optimal parallactic angle variation)?

The pca function in VIP accepts a parameter source_xy for defining a location of interest for which the rotation threshold will be applied. The parameter delta_rot sets the minimum amount of linear motion in terms of the FWHM to include other frames when building the PCA library for each image of the sequence. In the following example we set the threshold to 1 FWHM, which will significantly reduce companion self-subtraction. We also increase the npc to 30.

[31]:

fr_pca_optlib = pca(cube, angs, ncomp=30, source_xy=xy_b, delta_rot=1, fwhm=fwhm_naco,
                    imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:29:56
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 1.075 GB
Size LIB: min=30.0 / 1st QU=33.0 / med=39.0 / 3rd QU=44.0 / max=52.0
Done de-rotating and combining
Running time:  0:00:04.147441
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[32]:

plot_frames(fr_pca_optlib)#, backend='bokeh')

../_images/tutorials_03_psfsub_106_0.png

Question 3.11: As can be seen above, a PA threshold reduces the self-subtraction (notice the maximum in the scale), but also leads to a less efficient modeling and subtraction of speckle in the innermost part of the image. Why is that so?

We can have a quick report on the last two frames, obtained without and with parallactic angle threshold, and in particular inspect the effect on the flux and S/N ratio of the companion:

[33]:

snr(fr_pca2, xy_b, fwhm_naco, verbose=True)

S/N for the given pixel = 11.321
Integrated flux in FWHM test aperture = 114.501
Mean of background apertures integrated fluxes = -4.382
Std-dev of background apertures integrated fluxes = 10.260

[33]:

11.320734961541241

[34]:

snr(fr_pca_optlib, xy_b, fwhm_naco, verbose=True)

S/N for the given pixel = 7.454
Integrated flux in FWHM test aperture = 538.744
Mean of background apertures integrated fluxes = -26.475
Std-dev of background apertures integrated fluxes = 74.088

[34]:

7.453634382886601

The companion flux on fr_pca_optlib is much larger than on fr_pca2 (best full-frame PCA frame without the application of a rotation threshold) because we’ve reduced self-subtraction (companion over subtraction). On the other hand, the S/N has decreased due to the larger residual speckle noise.

Answer 3.10 Targets with a declination $\delta$ very different to the latitude $L$ of the observatory will have a very slow parallactic angle variation rate, and are hence not ideal for ADI. On the other extreme, an object whose $\delta$ matches almost exactly $L$ (within a few degrees) will have all the field rotation happening in a matter of minutes, around the time when the object crosses the local meridian (north-south imaginary line) - called transit. This is

not ideal either. The best targets for ADI are thus the ones where $|\delta-L|$ is neither too small ($>5$deg) nor too large ($\lesssim 35$deg), as is the case for beta Pic as observed from the Very Large Telescope ($|\delta-L| \approx 27$deg)

Answer 3.11 Using a PA threshold implies that for a given image to be modeled, the PCA library will be built using exclusively images taken sufficiently far away in time (hence with less correlated speckle). This in turn can lead to a less efficient speckle subtraction. This effect can be particularly stronger in the innermost parts of the image, where a given threshold in terms of linear motion implies a larger delay than for the outer parts of the image.

3.5.4. PCA for big datacubes¶

Computing the PCA model can be a CPU and memory intensive procedure, when working with thousands of large frames, let’s say cubes with shapes [10000, 600, 600]. Such cube won’t fit on the memory of a normal laptop or desktop computer. With VIP we just need to provide a positive integer value to the batch parameter of the pca function to alleviate this issue. Let’s compute the incremental PCA on our toy cube (even if its size does not require it):

[35]:

fr_increm = pca(cube, angs, ncomp=19, fwhm=fwhm_naco, batch=31, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:00
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 1.032 GB
Cube size = 0.002 GB (61 frames)
Batch size = 31 frames (0.001 GB)

Batch 1/2       shape: (31, 101, 101)   size: 1.3 MB
Batch 2/2       shape: (30, 101, 101)   size: 1.2 MB
Running time:  0:00:00.201156
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Reconstructing and obtaining residuals
Running time:  0:00:01.368411
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[36]:

plot_frames(fr_increm)#, backend='bokeh')

../_images/tutorials_03_psfsub_117_0.png

3.5.5. Annular PCA¶

PCA can also be performed on concentric annuli, and including a parallactic angle threshold when building the PCA library associated to each image. This is the same idea used before in full-frame when the position of a source of interest was defined except that the PA threshold here will be adjusted depending on the radial distance of each annulus from the star. PCA can be computed in full annuli or in separate annular segments (n_segments). The computational cost increases accordingly.

The function pca_annular processes the cube in annular fashion in a serial or parallel way. Setting nproc to an integer will use that number of CPUs in multiprocessing, while setting it to None will automatically set it to half the number of CPUs available on the machine. The improvement in speed is noticeable in multi-core machines.

Let’s try pca_annular with 3 segments per annulus, 6 principal components, a 0.2 FWHM rotation threshold, and in multi-processing mode:

[37]:

fr_pca_an6 = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp=6, asize=fwhm_naco, delta_rot=0.2, nproc=None,
                          n_segments=3, svd_mode='lapack', imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:01
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 10, FWHM = 4.603
PCA per annulus (or annular sectors):
Ann 1    PA thresh: 22.62    Ann center:   2    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:00.154152
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.63    Ann center:   7    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:00.316621
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  4.58    Ann center:  12    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:00.494537
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh:  3.27    Ann center:  16    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:00.679275
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh:  2.55    Ann center:  21    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:00.861368
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  2.08    Ann center:  25    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:01.052701
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 7    PA thresh:  1.76    Ann center:  30    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:01.249771
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 8    PA thresh:  1.53    Ann center:  35    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:01.448816
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 9    PA thresh:  1.35    Ann center:  39    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:01.648350
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 10    PA thresh:  1.23    Ann center:  43    N segments: 3
Done PCA with lapack for current annulus
Running time:  0:00:01.856600
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:02.253299
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

The parameter ncomp can be set to auto for letting the algorithm define automatically the number of PCs for each annular patch (instead of fixing it for all of them). The optimal value is found when the standard deviation of the residuals after the subtraction of the PCA approximation drops below a given (absolute) tolerance tol while progressively increasing the number of principal components. Note that this is different to finding the optimal $n_{\rm pc}$ that maximizes the S/N ratio of a given companion companion (as the former does not require any companion).

[38]:

fr_pca_auto = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp='auto', svd_mode='lapack', asize=fwhm_naco,
                          delta_rot=0.2, tol=0.1, nproc=None, imlib=imlib_rot, interpolation=interpolation,
                          verbose=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:04
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 10, FWHM = 4.603
PCA per annulus (or annular sectors):
Ann 1    PA thresh: 22.62    Ann center:   2    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.075258
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.63    Ann center:   7    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.173753
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  4.58    Ann center:  12    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.289932
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh:  3.27    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.400212
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh:  2.55    Ann center:  21    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.505738
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  2.08    Ann center:  25    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.627375
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 7    PA thresh:  1.76    Ann center:  30    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.764332
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 8    PA thresh:  1.53    Ann center:  35    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.903481
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 9    PA thresh:  1.35    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:01.043284
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 10    PA thresh:  1.23    Ann center:  43    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:01.192695
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:01.615928
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Now let’s compare both images:

[39]:

plot_frames((fr_pca_an6, fr_pca_auto), dpi=100, vmin=-5, colorbar=True)

../_images/tutorials_03_psfsub_124_0.png

And let’s compare the S/N ratio measured in each case:

[40]:

snr1 = snr(fr_pca_an6, fwhm=fwhm_naco, source_xy=xy_b, verbose=True)

S/N for the given pixel = 9.401
Integrated flux in FWHM test aperture = 195.983
Mean of background apertures integrated fluxes = -5.573
Std-dev of background apertures integrated fluxes = 20.946

[41]:

snr2 = snr(fr_pca_auto, fwhm=fwhm_naco, source_xy=xy_b, verbose=True)

S/N for the given pixel = 10.533
Integrated flux in FWHM test aperture = 84.835
Mean of background apertures integrated fluxes = -3.679
Std-dev of background apertures integrated fluxes = 8.211

Despite the relatively high S/N ratio, the planet appears significantly self-subtracted. This is because delta_rot was set to a low value. Setting it to a value larger than 0.3 leads to a bug (feel free to test it). This is because for a given image to be modelled, not enough images are found in the rest of the datacube complying with the PA threshold condition for the innermost annuli. To solve this issue, one can increase the radius radius_int of the innermost annulus where annular PCA starts.

[42]:

fr_pca_ann_mask = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp='auto', asize=fwhm_naco, radius_int=5,
                              delta_rot=1, nproc=None, svd_mode='lapack', imlib=imlib_rot,
                              interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:06
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 9, FWHM = 4.603
PCA per annulus (or annular sectors):
Ann 1    PA thresh: 34.99    Ann center:   7    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.072635
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh: 21.88    Ann center:  12    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.165516
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh: 15.87    Ann center:  17    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.269982
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh: 12.44    Ann center:  21    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.372420
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh: 10.23    Ann center:  26    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.479060
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  8.68    Ann center:  30    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.602953
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 7    PA thresh:  7.54    Ann center:  35    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.728285
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 8    PA thresh:  6.67    Ann center:  40    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.856712
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 9    PA thresh:  6.11    Ann center:  43    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.990339
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:01.381445
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

When using PCA with a numerical mask covering the inner part of the image, it is recommended to explicitly set the optional argument mask_val=0 (the value used for missing values not calculated by PCA). This will tell the image (de)rotation algorithm which pixels are part of a mask. These values will then be interpolated before rotation (using a Gaussian kernel size set by ker), and masked again after rotation. This will prevent a sharp transition between the mask and the signals near the mask. This sharp transition in the inner part of the image would otherwise lead to either Gibbs-related artefacts (when using the FFT-based rotation method), or ringing artefacts (when using an interpolation-based method). See the difference below. Similar tests are also shown in Tutorial 7. FFT- vs interpolation-based image operations.

[43]:

fr_pca_ann_mask0 = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp='auto', asize=fwhm_naco, radius_int=5,
                              delta_rot=1, nproc=None, svd_mode='lapack', imlib=imlib_rot,
                              interpolation=interpolation, mask_val=0, ker=1)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:07
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 9, FWHM = 4.603
PCA per annulus (or annular sectors):
Ann 1    PA thresh: 34.99    Ann center:   7    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.082763
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh: 21.88    Ann center:  12    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.188862
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh: 15.87    Ann center:  17    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.297967
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh: 12.44    Ann center:  21    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.399461
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh: 10.23    Ann center:  26    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.506975
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  8.68    Ann center:  30    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.630319
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 7    PA thresh:  7.54    Ann center:  35    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.757159
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 8    PA thresh:  6.67    Ann center:  40    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.885795
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 9    PA thresh:  6.11    Ann center:  43    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:01.024896
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:01.396278
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[44]:

plot_frames((fr_pca_ann_mask, fr_pca_ann_mask0), vmin=-15, vmax=15)

../_images/tutorials_03_psfsub_132_0.png

For more flexibility, delta_rot can also be set to a tuple of 2 values instead of a scalar. If a tuple, it corresponds to the PA thresholds for the innermost and outermost annuli in the image (linearly interpolated in between).

[45]:

fr_pca_ddrot = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp='auto', asize=fwhm_naco, delta_rot=(0.2,1), nproc=None,
                          svd_mode='lapack', imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:09
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 10, FWHM = 4.603
PCA per annulus (or annular sectors):
Ann 1    PA thresh: 22.62    Ann center:   2    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.069465
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh: 11.00    Ann center:   7    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.162790
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.64    Ann center:  12    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.275536
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh:  7.63    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.382676
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh:  7.06    Ann center:  21    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.486519
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  6.71    Ann center:  25    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.594576
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 7    PA thresh:  6.46    Ann center:  30    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.717737
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 8    PA thresh:  6.28    Ann center:  35    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.841542
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 9    PA thresh:  6.14    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.969251
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 10    PA thresh:  6.17    Ann center:  43    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:01.101693
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:01.483061
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s compare with the image obtained with a single PA threshold at all radii.

[46]:

plot_frames((fr_pca_an6, fr_pca_ddrot), dpi=100, vmin=-5, colorbar=True)

../_images/tutorials_03_psfsub_136_0.png

The residual speckle noise is now lower in the innermost part of the image.

3.5.6. PCA in a single annulus¶

Once a companion candidate has been identified, the fastest way to recover (and hence characterize) it is to apply PCA on a single annulus encompassing the companion candidate. This is used in particular for the characterization of companion candidates with the negative fake companion technique (see Tutorial 5. Planet forward modeling). The relevant function is called pca_annulus, which requires the radius r_guess of the annulus to be provided, and its annulus_width.

[47]:

ncomp_test = 5
pca_ann_test = pca_annulus(cube, angs, ncomp=ncomp_test, annulus_width=3*fwhm_naco, r_guess=rad,
                           imlib=imlib_rot, interpolation=interpolation)

[48]:

plot_frames(pca_ann_test, label='PCA annulus (npc={:.0f})'.format(ncomp_test),
            dpi=100, vmin=-5, colorbar=True)

../_images/tutorials_03_psfsub_141_0.png

Let’s now use pca_grid to search for the optimal $n_{\rm pc}$ that maximizes the S/N ratio of a given point source. pca_grid encompasses both the pca and pca_annulus function, depending on whether mode='fullfr' or mode='annular', respectively. This will allow us to answer Q3.9 we posed earlier: Is the optimal number of PCs the same for full-frame PCA and PCA in an annulus?

[49]:

# full frame
res_ff_opt = pca_grid(cube, angs, fwhm=fwhm_naco, range_pcs=(1,31,1), source_xy=xy_b, mode='fullfr',
                      imlib=imlib_rot, interpolation=interpolation, full_output=True, plot=True)

# single annulus with no PA threshold
res_ann_opt = pca_grid(cube, angs, fwhm=fwhm_naco, range_pcs=(1,31,1), source_xy=xy_b, mode='annular',
                       annulus_width=3*fwhm_naco, imlib=imlib_rot, interpolation=interpolation,
                       full_output=True, plot=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.114386
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 31
Optimal number of PCs = 13, for S/N=9.028
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 58.5, 35.5
Flux in a centered 1xFWHM circular aperture = 114.501
Central pixel S/N = 11.321
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 9.028
Max S/N (shifting the aperture center) = 11.526
stddev S/N (shifting the aperture center) = 2.167

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:30:47
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.013453
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 31
Optimal number of PCs = 7, for S/N=10.032
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 58.5, 35.5
Flux in a centered 1xFWHM circular aperture = 199.955
Central pixel S/N = 11.180
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.032
Max S/N (shifting the aperture center) = 14.606
stddev S/N (shifting the aperture center) = 2.763

../_images/tutorials_03_psfsub_143_1.png

../_images/tutorials_03_psfsub_143_2.png

For full-frame PCA the optimal npc is found to be 13 (as already seen in Sec. 3.5.2).

For PCA on a single annulus, the optimal npc is found to be 7, hence lower.

Let’s plot the final images below. Note that since we ran pca_grid with the full_output=True option, both final images and optimal number of pcs are returned (among other outputs).

[50]:

_, final_ff_opt, _, opt_npc_ff = res_ff_opt
_, final_ann_opt, _, opt_npc_ann = res_ann_opt

[51]:

plot_frames((final_ff_opt, final_ann_opt), label=('PCA-opt full-frame (npc={:.0f})'.format(opt_npc_ff),
                                                  'PCA-opt annulus (npc={:.0f})'.format(opt_npc_ann)),
            dpi=100, vmin=-5, colorbar=True)

../_images/tutorials_03_psfsub_146_0.png

Answer 5.6: The optimal $n_{\rm pc}$ is in general lower for PCA on a single annulus than on full frame. Intuitively, one could expect this result due to the smaller number of pixels in the annular case compared to the full frame case. With less pixels (i.e. a smaller matrix for singular value decomposition), a lower $n_{\rm pc}$ is sufficient to achieve similar (or better) modeling of speckle. Conversely a lower $n_{\rm pc}$ would be enough to cause overfitting.

3.6. Non-negative Matrix Factorization (NMF)¶

A PSF reference can be modelled using a low-rank approximation different than PCA (which relies on singular value decomposition). Non-negative matrix factorization can be used instead (Gomez Gonzalez et al. 2017; Ren et al. 2018). Instead of $M = U\Sigma V$ (SVD), non-negative matrix factorization can be written as $M = \mathcal{W} \mathcal{H}$, where both matrices are non-negative, and the columns of $\mathcal{H}$ are the non-zero components onto which the images are projected (note that the residuals after model subtraction can be negative though).

The two matrices are typically found iteratively. Currently, the default method in VIP is based on the scikit-learn implementation of NMF, which uses the principal components obtained by SVD as first guess of the iterative process. Therefore, the results can look fairly similar to PCA, in particular if the PCs are mostly non-negative.

3.6.1. Full-frame NMF¶

Full-frame NMF is simply done by calling the nmf function. Input arguments are mostly similar to pca, apart from init_svd which sets the method used to estimate the initial SVD from which the solution is iteratively found.

[52]:

fr_nmf = nmf(cube, angs, ncomp=14, max_iter=10000, init_svd='nndsvdar', mask_center_px=None, imlib=imlib_rot,
             interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:22
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:05.078957
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done NMF with sklearn.NMF.
Running time:  0:00:05.081559
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:06.183213
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[53]:

plot_frames(fr_nmf, dpi=100, colorbar=True)#, backend='bokeh')

../_images/tutorials_03_psfsub_153_0.png

[54]:

snr_nmf = snr(fr_nmf, fwhm=fwhm_naco, source_xy=xy_b, verbose=True)

S/N for the given pixel = 6.454
Integrated flux in FWHM test aperture = 120.744
Mean of background apertures integrated fluxes = -2.919
Std-dev of background apertures integrated fluxes = 18.722

3.6.2. Annular NMF¶

The same as above can be performed in concentric annuli (new since v1.0.0):

[55]:

fr_nmf_ann = nmf_annular(cube, angs, ncomp=9, max_iter=10000, init_svd='nndsvdar', radius_int=0, nproc=None,
                         fwhm=fwhm_naco, asize=fwhm_naco, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:28
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 10, FWHM = 4.603
NMF per annulus (or annular sectors):
Ann 1    PA thresh: 11.42    Ann center:   2    N segments: 1
Ann 2    PA thresh:  7.63    Ann center:   7    N segments: 1
Ann 3    PA thresh:  6.87    Ann center:  12    N segments: 1
Ann 4    PA thresh:  6.54    Ann center:  16    N segments: 1
Ann 5    PA thresh:  6.36    Ann center:  21    N segments: 1
Ann 6    PA thresh:  6.24    Ann center:  25    N segments: 1
Ann 7    PA thresh:  6.16    Ann center:  30    N segments: 1
Ann 8    PA thresh:  6.11    Ann center:  35    N segments: 1
Ann 9    PA thresh:  6.06    Ann center:  39    N segments: 1
Ann 10    PA thresh:  6.17    Ann center:  43    N segments: 1
Done derotating and combining.
Running time:  0:00:15.829312
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[56]:

plot_frames(fr_nmf_ann, dpi=100, colorbar=True)

../_images/tutorials_03_psfsub_158_0.png

[57]:

snr_nmf_ann = snr(fr_nmf_ann, fwhm=fwhm_naco, source_xy=xy_b, verbose=True)

S/N for the given pixel = 8.379
Integrated flux in FWHM test aperture = 363.671
Mean of background apertures integrated fluxes = -13.543
Std-dev of background apertures integrated fluxes = 43.986

[58]:

snrmap_nmf_ann = snrmap(fr_nmf_ann, fwhm_naco, plot=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:44
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

../_images/tutorials_03_psfsub_160_1.png

S/N map created using 5 processes
Running time:  0:00:02.770179
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

3.7. Local Low-rank plus Sparse plus Gaussian-noise decomposition (LLSG)¶

Local Low-rank plus Sparse plus Gaussian-noise decomposition (LLSG; Gomez Gonzalez et al. 2016) proposes a three terms decomposition to improve the detectability of point-like sources in ADI data. It aims at decomposing ADI cubes into L+S+G (low-rank, sparse and Gaussian noise) terms. Separating the noise from the S component (where the moving planet should stay) allows us to increase the S/N of potential planets.

Let’s try it out:

[59]:

fr_llsg = llsg(cube, angs, fwhm_naco, rank=5, thresh=1, max_iter=20, random_seed=10,
               imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:47
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Annuli = 5
Processing annulus:
1 : in_rad=0, n_segm=4
2 : in_rad=10, n_segm=4
3 : in_rad=20, n_segm=4
4 : in_rad=30, n_segm=4
5 : in_rad=40, n_segm=4

Running time:  0:00:03.388970
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s plot it along the full-frame PCA-ADI result, after the optimization of the number of PCs:

[60]:

plot_frames((final_ff_opt, fr_llsg))#, backend='bokeh')

../_images/tutorials_03_psfsub_165_0.png

The computation of S/N ratio and S/N maps when the noise has been almost totally supressed becomes problematic (both from a theoretical and computational points of view). Let’s nevertheless compute them for reference:

[61]:

snr_llsg = snr(fr_llsg, fwhm=fwhm_naco, source_xy=xy_b, verbose=True)

S/N for the given pixel = 28.185
Integrated flux in FWHM test aperture = 225.660
Mean of background apertures integrated fluxes = -1.612
Std-dev of background apertures integrated fluxes = 7.878

[62]:

snrmap_llsg = snrmap(fr_llsg, fwhm_naco, plot=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:50
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

../_images/tutorials_03_psfsub_168_1.png

S/N map created using 5 processes
Running time:  0:00:02.766376
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

3.8. ANgular Differential OptiMal Exoplanet Detection Algorithm (ANDROMEDA)¶

Another approach to find point sources in ADI sequences consists in i) subtracting pairs of images separated by a known amount of rotation (typically 0.5 FWHM to produce a specific positive-negative signature); ii) find this residual signature using a maximum match-filter; iii) evaluate the likelihood that a point source is present at each pixel and the corresponding photometry. This is what ANDROMEDA does (Mugnier et al. (2009), Cantalloube et al. (2015)). The principle is summarized in this figure:

ANDROMEDA follows an inverse-problem approach since it involves a model for the expected positive-negative signature of a point source, and finding the optimal coefficients which minimize residuals between model and observed residual images.

The andromeda function in VIP requires the calculation of the oversampling factor (i.e. ratio between Nyquist sampling and actual pixel sampling):

[63]:

lbda = VLT_NACO['lambdal']
diam = VLT_NACO['diam']
resel = (lbda/diam)*206265 # lambda/D in arcsec
nyquist_samp = resel/2.
oversamp_fac = nyquist_samp/pxscale_naco
oversamp_fac

[63]:

1.7577458534791308

Let’s now try it:

[64]:

res = andromeda(cube, oversamp_fac, angs, psf, filtering_fraction=.25,
              min_sep=.5, annuli_width=1., roa=2, opt_method='lsq',
              nsmooth_snr=18, iwa=2, owa=None, precision=50, fast=False,
              homogeneous_variance=True, ditimg=1.0, ditpsf=None, tnd=1.0,
              total=False, multiply_gamma=True, nproc=1, verbose=False)

[65]:

flux, snr, snr_norm, stdcontrast, stdcontrast_norm, likelihood, ext_radius = res

One can then plot the SNR and flux of a putative companion at each pixel:

[66]:

plot_frames((snr_norm, flux))#, backend='bokeh')

../_images/tutorials_03_psfsub_178_0.png

3.9. Forward-Model Matched Filter (FMMF)¶

Let’s now test the forward-model matched filter (FMMF) with both LOCI and KLIP (Pueyo 2016; Ruffio et al. 2017; Dahlqvist et al. 2021):

[67]:

if version.parse(vvip) >= version.parse("1.2.4"):
    loci_fmmf_f, loci_fmmf_snr = fmmf(cube, angs, psf, fwhm_naco, model='LOCI', var='FR', nproc=None,
                                      min_r = int(2*fwhm_naco), max_r = int(6*fwhm_naco),
                                      param={'ncomp': 10, 'tolerance': 0.005, 'delta_rot': 0.5},
                                      crop=5, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:31:58
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Radial distance 9 done!
Radial distance 11 done!
Radial distance 10 done!
Radial distance 12 done!
Radial distance 13 done!
Radial distance 14 done!
Radial distance 15 done!
Radial distance 16 done!
Radial distance 17 done!
Radial distance 18 done!
Radial distance 19 done!
Radial distance 20 done!
Radial distance 21 done!
Radial distance 22 done!
Radial distance 23 done!
Radial distance 24 done!
Radial distance 25 done!
Radial distance 26 done!
Running time:  0:09:43.078826
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[68]:

if version.parse(vvip) >= version.parse("1.2.4"):
    klip_fmmf_f, klip_fmmf_snr = fmmf(cube, angs, psf, fwhm_naco, model='KLIP', var='FR', nproc=None,
                                      min_r = int(2*fwhm_naco), max_r = int(6*fwhm_naco),
                                      param={'ncomp': 10, 'tolerance': 0.005, 'delta_rot': 0.5},
                                      crop=5, imlib=imlib_rot, interpolation=interpolation)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-10 17:41:42
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Radial distance 9 done!
Radial distance 11 done!
Radial distance 10 done!
Radial distance 12 done!
Radial distance 13 done!
Radial distance 14 done!
Radial distance 15 done!
Radial distance 16 done!
Radial distance 17 done!
Radial distance 18 done!
Radial distance 19 done!
Radial distance 20 done!
Radial distance 21 done!
Radial distance 22 done!
Radial distance 23 done!
Radial distance 24 done!
Radial distance 25 done!
Radial distance 26 done!
Running time:  0:10:14.201550
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[69]:

if version.parse(vvip) >= version.parse("1.2.4"):
    plot_frames((loci_fmmf_f, loci_fmmf_snr, klip_fmmf_f, klip_fmmf_snr), rows=2,
                 label=('LOCI-FMMF (flux)', 'LOCI-FMMF (SNR)','KLIP-FMMF (flux)','KLIP-FMMF (SNR)'))#, backend='bokeh')

../_images/tutorials_03_psfsub_183_0.png

3.10. Patch Covariances (PACO)¶

The PAtch COvariances (PACO) algorithm was introduced by Flasseur et al. (2009). It provides a method for extracting the faint signal of the exoplanet from the background noise by building up an empirical model of the background noise. The motion of the planet through the ADI frame is used to identify statistically significant point sources within the field of view. Finally, the flux of the planet can be estimated using an iterative fitting procedure.

Two versions of the PACO algorithm are implemented. While the FullPACO method is more statistically robust, FastPACO produces adequate results in a much shorter time, and has been parallelised using the multiprocessing library.

[70]:

if version.parse(vvip) >= version.parse("1.2.4"):
    paco = FastPACO(cube=cube, angles=angs, psf = psf, dit_psf = 1.0, dit_science = 1.0, nd_transmission = 1.0,
                    fwhm=pxscale_naco*fwhm_naco,#PACO takes the fwhm in arcsec
                    pixscale = pxscale_naco, rescaling_factor = 1.0, verbose = True)
    paco_snr, paco_flux = paco.run(cpu = cpu_count()//2)

----------------------
Summary of PACO setup:

Image Cube shape = (61, 101, 101)
PIXSCALE = 0.02719
PSF |  Area  |  Rad   |  Width |
    |   49   |  04    |  039   |
Patch width: 11
----------------------

Precomputing Statistics using 5 Processes...
FWHM: 4.876
Flux in 1xFWHM aperture: 1.338
Running Fast PACO...
Done

[71]:

if version.parse(vvip) >= version.parse("1.2.4"):
    plot_frames((paco_snr, paco_flux), label=("PACO SNR", "PACO flux"))

../_images/tutorials_03_psfsub_188_0.png

We can now estimate the flux from the snr map.

[72]:

if version.parse(vvip) >= version.parse("1.2.4"):
    sources = paco.subpixel_threshold_detect(paco_snr, threshold = 5.0) #Threshold in sigma
    posns = [(sources[entry]['x'], sources[entry]['y']) for entry in sources]
    initial_guesses = [paco_flux[int(s[0]),int(s[1])] for s in posns]
    flux_estimates = paco.flux_estimate(posns,eps = 0.1, initial_est = initial_guesses)

Blobs found: 1
 ycen   xcen
------ ------
35.809   59.507

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
X,Y = (59.5,35.8)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 59.5, 35.8
Flux in a centered 1xFWHM circular aperture = 77.149
Central pixel S/N = 8.785
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 6.204
Max S/N (shifting the aperture center) = 8.809
stddev S/N (shifting the aperture center) = 1.639

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Computing unbiased flux estimate...
Initial guesses:
Positions:  [(59.50743857609105, 35.80875467973141)]
Contrasts:  [1.7155835246740885]
FWHM: 4.876
Flux in 1xFWHM aperture: 1.338
Precomputing Statistics...
114.35862566624024
Contrast estimate: [109.97867044]
Contrast estimate: [107.82367406]
Extracted contrasts
-------------------
x: 59.50743857609105, y: 35.80875467973141, contrast: [107.82367406]±[4.11762806]

3.11. Summary mosaic¶

Let’s take a look at the final residual frames produced by the different algorithms:

[73]:

if version.parse(vvip) <= version.parse("1.2.3"):
    plot_frames((fr_adi, fr_adi_an, fr_fdiff, fr_lstsq, fr_pca1, fr_pca_an6,
                 fr_pca_auto, final_ff_opt, final_ann_opt, fr_nmf, fr_nmf_ann,
                 fr_llsg, snr_norm),
                rows=3, colorbar=True, label_size=12, label_pad=8, axis=False, versp=0.05, horsp=0.15,
                label=('median ADI', 'ann. median ADI', 'pairwise frame diff.', 'LOCI',
                       'full-frame PCA (5 PCs)', 'ann. PCA (3 seg., 6 PCs)',
                       'ann. PCA auto PCs (noise tol.)', 'optimal full-frame PCA (max SNR)',
                       'optimal ann. PCA (max SNR)', 'full-frame NMF (13 comps)', 'ann. NMF (9 comps)',
                       'LLSG', 'ANDROMEDA (SNR)'))
else:
    plot_frames((fr_adi, fr_adi_an, fr_fdiff, fr_lstsq, fr_pca1, fr_pca_an6,
                 fr_pca_auto, final_ff_opt, final_ann_opt, fr_nmf, fr_nmf_ann,
                 fr_llsg, snr_norm, loci_fmmf_snr, klip_fmmf_snr, paco_snr),
                rows=4, colorbar=True, label_size=12, label_pad=8, axis=False, versp=0.05, horsp=0.15,
                label=('median ADI', 'ann. median ADI', 'pairwise frame diff.', 'LOCI',
                       'full-frame PCA (5 PCs)', 'ann. PCA (3 seg., 6 PCs)',
                       'ann. PCA auto PCs (noise tol.)', 'optimal full-frame PCA (max SNR)',
                       'optimal ann. PCA (max SNR)', 'full-frame NMF (13 comps)', 'ann. NMF (9 comps)',
                       'LLSG', 'ANDROMEDA (SNR)', 'LOCI-FMMF (SNR)', 'KLIP-FMMF (SNR)', 'PACO (SNR)'))

../_images/tutorials_03_psfsub_193_0.png

Check out tutorial 4. Metrics for more details on evaluating i) the significance of a potential point source, ii) the achieved contrast in your image, in the form of contrast curves.