6. PSF modeling and subtraction with reference star differential imaging
Author: Valentin Christiaens
Suitable for VIP v1.0.3 onwards
Last update: 2024/06/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, and leverage the RDI strategy;
how to use iterative (greedy) algorithms in VIP to correct for geometric biases.
Most PSF subtraction routines in VIP (and some other functions) have been implemented in VIP to allow for multiprocessing, in order to optimally harness the power of machines equipped with multiple CPUs. Any function where the nproc parameter is available in its call can be run in multi-processing, with the value of nproc setting the requested number of CPUs to use. Instead of an integer, one can set nproc=None to use half of all available CPUs. For optimal results in multiprocessing, set the following environment parameters in your terminal BEFORE launching your Jupyter notebook:
export MKL_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1
In the case of PCA, singular value decomposition can also be done on GPU by setting svd_mode to an appropriate value (see Sec. 3.5.1. and docstrings of vip_hci.psfsub.pca for details).
Let’s first import a couple of external packages needed in this tutorial:
%matplotlib inline
from hciplot import plot_frames, plot_cubes # plotting routines
from matplotlib import pyplot as plt
from multiprocessing import cpu_count
import numpy as np
from packaging import version
In the following box we check that your version of VIP is recent enough.
import vip_hci as vip
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) < version.parse("1.0.3"):
msg = "Please upgrade your version of VIP"
msg+= "It should be 1.0.3 or above to run this notebook."
raise ValueError(msg)
VIP version: 1.6.2
6.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:
from vip_hci.fits import open_fits
from astropy.utils.data import download_file
url_d = "https://github.com/vortex-exoplanet/VIP_extras/raw/master/datasets"
f1 = download_file("{}/naco_betapic_cube_cen.fits".format(url_d), cache=True)
f2 = download_file("{}/naco_betapic_psf.fits".format(url_d), cache=True)
f3 = download_file("{}/naco_betapic_derot_angles.fits".format(url_d), cache=True)
# alternatively, for local files simply provide their full or relative path. E.g.:
#f1 = '../datasets/naco_betapic_cube_cen.fits'
#f2 = '../datasets/naco_betapic_psf.fits'
#f3 = '../datasets/naco_betapic_derot_angles.fits'
cube = open_fits(f1)
psf = open_fits(f2)
angs = open_fits(f3)
FITS HDU-0 data successfully loaded. Data shape: (61, 101, 101)
FITS HDU-0 data successfully loaded. Data shape: (39, 39)
FITS HDU-0 data successfully loaded. Data shape: (61,)
For the purpose of illustrating RDI capabilities in VIP, let’s just define a fiducial reference cube as the science cube flipped along both the x and y axes. Normally you should load it the same way as your science cube.
cube_ref = np.flip(cube, axis=1) # just for the example
cube_ref = np.flip(cube_ref, axis=2)
Let’s fit the PSF with a 2D Gaussian to infer the FWHM, the flux in a 1-FWHM size aperture, and get a flux-normalized PSF:
%matplotlib inline
from vip_hci.fm import normalize_psf
psfn, flux, fwhm_naco = normalize_psf(psf, size=19, debug=True, full_output=True)
FWHM_y = 4.926059872957138
FWHM_x = 4.67577889500593
centroid y = 9.010992107833063
centroid x = 9.01917912265807
centroid y subim = 9.010992107833063
centroid x subim = 9.01917912265807
amplitude = 0.10032285220380605
theta = -38.446187060503924
Mean FWHM: 4.801
Flux in 1xFWHM aperture: 1.307
print(fwhm_naco)
4.800919383981534
The normalize_psf function performs internally a fine centering of the PSF. The input PSF should nevertheless already be centered within a few pixel accuracy for the fine centering to work.
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.
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:
from vip_hci.config import VLT_NACO
pxscale_naco = VLT_NACO['plsc']
print(pxscale_naco, "arcsec/px")
0.02719 arcsec/px
6.2. Classical RDI
6.2.1. median-RDI
The most straightforward way to leverage a stack of reference PSF images is to consider its median image as PSF model, and subtract it to each image of the science cube. This can be done with the median_sub function:
from vip_hci.psfsub import median_sub
med_rdi = median_sub(cube, angs, cube_ref=cube_ref)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:08:46
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time: 0:00:04.439478
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let’s visualize the results, after masking the central 13 px of the image.
from vip_hci.var import mask_circle
plot_frames(mask_circle(med_rdi, 13), grid=True)
As can be seen above, this approach is sufficient to detect exoplanet beta Pic b towards the bottom right of the image (after masking the central part of the image)! This is not too bad, considering we only considered as fake reference cube the science images flipped horizontally and vertically.
The default behaviour will subtract the median reference image of the reference stack however other options can be used. Although the median is more robust, the mean of the reference stack can be requested instead:
mean_rdi = median_sub(cube, angs, cube_ref=cube_ref, collapse_ref='mean')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:08:51
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time: 0:00:04.334821
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Alternatively, the mean/median reference image can be scaled to match in flux each science image before subtraction - this can be useful if the reference image is known to be brighter or fainter than the science image:
med_rdi_sc = median_sub(cube, angs, cube_ref=cube_ref, collapse_ref='sc_median')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:08:55
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time: 0:00:04.381679
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
By default when scaling is requested, the scaling factor is calculating from integrating the flux on the whole science and reference images. But an inner and outer radius can be provided, which can be relevant if a coronagraph is used or contamination is present at large radius, respectively. This is done by adding 2 numbers directly at the end of the collapse_ref string, separated by ‘-’. For example, if we want to consider an annulus going from 1 to 3 FWHM for this coronagraphic NaCo dataset (FWHM ~ 5pixels), we can set:
med_rdi_sc515 = median_sub(cube, angs, cube_ref=cube_ref, collapse_ref='sc_median5-15')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:09:00
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time: 0:00:04.414767
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Alternatively, median-RDI can be applied in concentric annuli with a scaling factor calculated per annulus. This is controlled with the mode argument (set to ‘fullfr’ by default). The width of the concentric annuli is set in pixels using argument asize, and the radius of the first annuli is set with radius_int. Let’s set both to the FWHM.
med_rdi_ann = median_sub(cube, angs, mode='annular', asize=fwhm_naco, radius_int=fwhm_naco,
cube_ref=cube_ref, collapse_ref='sc_median')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:09:04
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 9, FWHM = 4
Warning: No valid output stream.
Optimized median psf reference subtracted
Done derotating and combining
Running time: 0:00:04.485246
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let’s now visualize all the results obtained above, without any numerical mask:
plot_frames((mean_rdi, med_rdi_sc, med_rdi_sc515, med_rdi_ann), vmin=-200, vmax=200, grid=True)
As expected the results obtained with the different collapse_ref options are very similar – since here we used as fake reference cube the flipped science cube. We also notice strong residuals near the center of the image which shows the limitation of assuming a perfect centro-symmetry stellar halo. Nonetheless the signature of beta Pic b can be seen towards the bottom right of the image.
6.2.2. PCA-RDI
Using PCA leveraging the RDI strategy is as simple as calling the pca function in VIP with the cube_ref argument set to your cube of reference images. For example, to run PCA using 42 principal components:
from vip_hci.psfsub import pca
pca_rdi_fr = pca(cube, angs, ncomp=42, cube_ref=cube_ref)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:09:10
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 8.101 GB
System available memory = 6.418 GB
Done vectorizing the frames. Matrix shape: (61, 10201)
Done vectorizing the frames. Matrix shape: (61, 10201)
Done SVD/PCA with numpy SVD (LAPACK)
Running time: 0:00:00.059065
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time: 0:00:04.512246
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let’s visualize the results, after masking the central 13 px of the image.
from vip_hci.var import mask_circle
plot_frames(mask_circle(pca_rdi_fr, 13), grid=True)
We now see Beta Pic b a bit clearer than with the median-subtraction approach.
A lot of options are available with the PCA function, and are also relevant with a reference cube. We refer to Tutorial 3 (section 3.5), for all available options. Below we just illustrate the use of a convenient function providing the SNR of a candidate as a function of number of principal components.
Let’s first define the approximate coordinates of the companion based on the image above (remember to subtract 1 to the coordinates read in the image since Python uses a zero-based index system):
xy_b = (58, 35)
from vip_hci.psfsub import pca_grid
pca_rdi_opt = pca_grid(cube, angs, cube_ref=cube_ref, fwhm=fwhm_naco, range_pcs=(1,51,1),
source_xy=xy_b, mode='fullfr', full_output=True, plot=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:09:15
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time: 0:00:00.055566
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[19], line 2
1 from vip_hci.psfsub import pca_grid
----> 2 pca_rdi_opt = pca_grid(cube, angs, cube_ref=cube_ref, fwhm=fwhm_naco, range_pcs=(1,51,1),
3 source_xy=xy_b, mode='fullfr', full_output=True, plot=True)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/psfsub/utils_pca.py:349, in pca_grid(cube, angle_list, fwhm, range_pcs, source_xy, cube_ref, mode, annulus_width, svd_mode, scaling, mask_center_px, fmerit, collapse, ifs_collapse_range, verbose, full_output, debug, plot, save_plot, start_time, scale_list, initial_4dshape, weights, exclude_negative_lobes, **rot_options)
347 for pc in pclist:
348 if mode == 'fullfr':
--> 349 frame = truncate_svd_get_finframe(matrix, angle_list, pc, V)
350 elif mode == 'annular':
351 frame = truncate_svd_get_finframe_ann(matrix, annind,
352 angle_list, pc, V)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/psfsub/utils_pca.py:224, in pca_grid.<locals>.truncate_svd_get_finframe(matrix, angle_list, ncomp, V)
221 else:
222 residuals_reshaped = residuals_res
--> 224 residuals_res_der = cube_derotate(residuals_reshaped, angle_list,
225 **rot_options)
226 res_frame = cube_collapse(residuals_res_der, mode=collapse, w=weights)
227 return res_frame
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:369, in cube_derotate(array, angle_list, imlib, interpolation, cxy, nproc, border_mode, mask_val, edge_blend, interp_zeros, ker)
367 array_der = np.zeros_like(array)
368 for i in range(n_frames):
--> 369 array_der[i] = frame_rotate(array[i], -angle_list[i], imlib=imlib,
370 interpolation=interpolation, cxy=cxy,
371 border_mode=border_mode,
372 mask_val=mask_val,
373 edge_blend=edge_blend,
374 interp_zeros=interp_zeros, ker=ker)
375 elif nproc > 1:
377 res = pool_map(nproc, _frame_rotate_mp, iterable(array),
378 iterable(-angle_list), imlib, interpolation, cxy,
379 border_mode, mask_val, edge_blend, interp_zeros, ker)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:228, in frame_rotate(array, angle, imlib, interpolation, cxy, border_mode, mask_val, edge_blend, interp_zeros, ker)
225 raise ValueError(msg)
227 if imlib == 'vip-fft':
--> 228 array_out = rotate_fft(array_prep, angle)
230 elif imlib == 'skimage':
231 if interpolation == 'nearneig':
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:596, in rotate_fft(array, angle)
594 # TODO: make FFT padding work for other option than '0'.
595 s_x = _fft_shear(array_in, arr_x, a, ax=1, pad=0)
--> 596 s_xy = _fft_shear(s_x, arr_y, b, ax=0, pad=0)
597 s_xyx = _fft_shear(s_xy, arr_x, a, ax=1, pad=0)
599 if y_ori % 2 or x_ori % 2:
600 # set it back to original dimensions
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:621, in _fft_shear(arr, arr_ori, c, ax, pad, shift_ini)
619 s_x = np.exp(-2j*np.pi*c*arr_u*arr_ori)*s_x
620 s_x = fftshift(s_x)
--> 621 s_x = ifft(s_x, axis=ax)
622 s_x = fftshift(s_x)
624 return s_x
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/numpy/fft/_pocketfft.py:318, in ifft(a, n, axis, norm, out)
316 if n is None:
317 n = a.shape[axis]
--> 318 output = _raw_fft(a, n, axis, False, False, norm, out=out)
319 return output
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/numpy/fft/_pocketfft.py:94, in _raw_fft(a, n, axis, is_real, is_forward, norm, out)
90 elif ((shape := getattr(out, "shape", None)) is not None
91 and (len(shape) != a.ndim or shape[axis] != n_out)):
92 raise ValueError("output array has wrong shape.")
---> 94 return ufunc(a, fct, axes=[(axis,), (), (axis,)], out=out)
KeyboardInterrupt:
pca?
Signature: pca(*all_args: List, **all_kwargs: dict)
Docstring:
Full-frame PCA algorithm applied to PSF substraction.
The reference PSF and the quasi-static speckle pattern are modeled using
Principal Component Analysis. Depending on the input parameters this PCA
function can work in ADI, RDI or mSDI (IFS data) mode.
ADI: the target ``cube`` itself is used to learn the PCs and to obtain a
low-rank approximation model PSF (star + speckles). Both `cube_ref`` and
``scale_list`` must be None. The full-frame ADI-PCA implementation is based
on [AMA12]_ and [SOU12]_. If ``batch`` is provided then the cube is
processed with incremental PCA as described in [GOM17]_.
(ADI+)RDI: if a reference cube is provided (``cube_ref``), its PCs are used
to reconstruct the target frames to obtain the model PSF (star + speckles).
(ADI+)mSDI (IFS data): if a scaling vector is provided (``scale_list``) and
the cube is a 4d array [# channels, # adi-frames, Y, X], it's assumed it
contains several multi-spectral frames acquired in pupil-stabilized mode.
A single or two stages PCA can be performed, depending on ``adimsdi``, as
explained in [CHR19]_.
Parameters
----------
all_args: list, optional
Positionnal arguments for the PCA algorithm. Full list of parameters
below.
all_kwargs: dictionary, optional
Mix of keyword arguments that can initialize a PCA_Params and the
optional 'rot_options' dictionnary (with keyword values ``border_mode``,
``mask_val``, ``edge_blend``, ``interp_zeros``, ``ker``; see docstring
of ``vip_hci.preproc.frame_rotate``). Can also contain a PCA_Params
dictionary named `algo_params`.
PCA parameters
--------------
cube : str or numpy ndarray, 3d or 4d
Input cube (ADI or ADI+mSDI). If 4D, the first dimension should be
spectral. If a string is given, it must correspond to the path to the
fits file to be opened in memmap mode (incremental PCA-ADI of 3D cubes
only).
angle_list : numpy ndarray, 1d
Vector of derotation angles to align North up in your images.
cube_ref : 3d or 4d numpy ndarray, or list of 3D numpy ndarray, optional
Reference library cube for Reference Star Differential Imaging. Should
be 3D, except if input cube is 4D and no scale_list is provided,
reference cube can then either be 4D or a list of 3D cubes (i.e.
providing the reference cube for each individual spectral cube).
scale_list : numpy ndarray, 1d, optional
If provided, triggers mSDI reduction. These should be the scaling
factors used to re-scale the spectral channels and align the speckles
in case of IFS data (ADI+mSDI cube). Usually, the
scaling factors are the last channel wavelength divided by the
other wavelengths in the cube (more thorough approaches can be used
to get the scaling factors, e.g. with
``vip_hci.preproc.find_scal_vector``).
ncomp : int, float, tuple of int/None, or list, optional
How many PCs are used as a lower-dimensional subspace to project the
target frames.
* ADI (``cube`` is a 3d array): if an int is provided, ``ncomp`` is the
number of PCs extracted from ``cube`` itself. If ``ncomp`` is a float
in the interval [0, 1] then it corresponds to the desired cumulative
explained variance ratio (the corresponding number of components is
estimated). If ``ncomp`` is a tuple of two integers, then it
corresponds to an interval of PCs in which final residual frames are
computed (optionally, if a tuple of 3 integers is passed, the third
value is the step). If ``ncomp`` is a list of int, these will be used to
calculate residual frames. When ``ncomp`` is a tuple or list, and
``source_xy`` is not None, then the S/Ns (mean value in a 1xFWHM
circular aperture) of the given (X,Y) coordinates are computed.
* ADI+RDI (``cube`` and ``cube_ref`` are 3d arrays): ``ncomp`` is the
number of PCs obtained from ``cube_ref``. If ``ncomp`` is a tuple,
then it corresponds to an interval of PCs (obtained from ``cube_ref``)
in which final residual frames are computed. If ``ncomp`` is a list of
int, these will be used to calculate residual frames. When ``ncomp`` is
a tuple or list, and ``source_xy`` is not None, then the S/Ns (mean
value in a 1xFWHM circular aperture) of the given (X,Y) coordinates are
computed.
* ADI or ADI+RDI (``cube`` is a 4d array): same input format allowed as
above. If ``ncomp`` is a list with the same length as the number of
channels, each element of the list will be used as ``ncomp`` value
(be it int, float or tuple) for each spectral channel. If not a
list or a list with a different length as the number of spectral
channels, these will be tested for all spectral channels respectively.
* ADI+mSDI (``cube`` is a 4d array and ``adimsdi="single"``): ``ncomp``
is the number of PCs obtained from the whole set of frames
(n_channels * n_adiframes). If ``ncomp`` is a float in the interval
(0, 1] then it corresponds to the desired CEVR, and the corresponding
number of components will be estimated. If ``ncomp`` is a tuple, then
it corresponds to an interval of PCs in which final residual frames
are computed. If ``ncomp`` is a list of int, these will be used to
calculate residual frames. When ``ncomp`` is a tuple or list, and
``source_xy`` is not None, then the S/Ns (mean value in a 1xFWHM
circular aperture) of the given (X,Y) coordinates are computed.
* ADI+mSDI (``cube`` is a 4d array and ``adimsdi="double"``): ``ncomp``
must be a tuple, where the first value is the number of PCs obtained
from each multi-spectral frame (if None then this stage will be
skipped and the spectral channels will be combined without
subtraction); the second value sets the number of PCs used in the
second PCA stage, ADI-like using the residuals of the first stage (if
None then the second PCA stage is skipped and the residuals are
de-rotated and combined).
svd_mode : Enum, see `vip_hci.config.paramenum.SvdMode`
Switch for the SVD method/library to be used.
scaling : Enum, see `vip_hci.config.paramenum.Scaling`
Pixel-wise scaling mode using ``sklearn.preprocessing.scale``
function. If set to None, the input matrix is left untouched.
mask_center_px : None or int
If None, no masking is done. If an integer > 1 then this value is the
radius of the circular mask.
source_xy : tuple of int, optional
For ADI-PCA, this triggers a frame rejection in the PCA library, with
``source_xy`` as the coordinates X,Y of the center of the annulus where
the PA criterion is estimated. When ``ncomp`` is a tuple, a PCA grid is
computed and the S/Ns (mean value in a 1xFWHM circular aperture) of the
given (X,Y) coordinates are computed.
delta_rot : int, optional
Factor for tuning the parallactic angle threshold, expressed in FWHM.
Default is 1 (excludes 1xFWHM on each side of the considered frame).
fwhm : float, list or 1d numpy array, optional
Known size of the FWHM in pixels to be used. Default value is 4.
Can be a list or 1d numpy array for a 4d input cube with no scale_list.
adimsdi : Enum, see `vip_hci.config.paramenum.Adimsdi`
Changes the way the 4d cubes (ADI+mSDI) are processed. Basically it
determines whether a single or double pass PCA is going to be computed.
crop_ifs: bool, optional
[adimsdi='single'] If True cube is cropped at the moment of frame
rescaling in wavelength. This is recommended for large FOVs such as the
one of SPHERE, but can remove significant amount of information close to
the edge of small FOVs (e.g. SINFONI).
imlib : Enum, see `vip_hci.config.paramenum.Imlib`
See the documentation of ``vip_hci.preproc.frame_rotate``.
imlib2 : Enum, see `vip_hci.config.paramenum.Imlib`
See the documentation of ``vip_hci.preproc.cube_rescaling_wavelengths``.
interpolation : Enum, see `vip_hci.config.paramenum.Interpolation`
See the documentation of the ``vip_hci.preproc.frame_rotate`` function.
collapse : Enum, see `vip_hci.config.paramenum.Collapse`
Sets how temporal residual frames should be combined to produce an
ADI image.
collapse_ifs : Enum, see `vip_hci.config.paramenum.Collapse`
Sets how spectral residual frames should be combined to produce an
mSDI image.
ifs_collapse_range: str 'all' or tuple of 2 int
If a tuple, it should contain the first and last channels where the mSDI
residual channels will be collapsed (by default collapses all channels).
mask_rdi: tuple of two numpy array or one signle 2d numpy array, opt
If provided, binary mask(s) will be used either in RDI mode or in
ADI+mSDI (2 steps) mode. If two masks are provided, they will the anchor
and boat regions, respectively, following the denominations in [REN23]_.
If only one mask is provided, it will be used as the anchor, and the
boat images will not be masked (i.e., full frames used).
check_memory : bool, optional
If True, it checks that the input cube is smaller than the available
system memory.
batch : None, int or float, optional
When it is not None, it triggers the incremental PCA (for ADI and
ADI+mSDI cubes). If an int is given, it corresponds to the number of
frames in each sequential mini-batch. If a float (0, 1] is given, it
corresponds to the size of the batch is computed wrt the available
memory in the system.
nproc : None or int, optional
Number of processes for parallel computing. If None the number of
processes will be set to (cpu_count()/2). Defaults to ``nproc=1``.
full_output: bool, optional
Whether to return the final median combined image only or with other
intermediate arrays.
verbose : bool, optional
If True prints intermediate info and timing.
weights: 1d numpy array or list, optional
Weights to be applied for a weighted mean. Need to be provided if
collapse mode is 'wmean'.
left_eigv : bool, optional
Whether to use rather left or right singularvectors
This mode is not compatible with 'mask_rdi' and 'batch'
min_frames_pca : int, optional
Minimum number of frames required in the PCA library. An error is raised
if less than such number of frames can be found.
cube_sig: numpy ndarray, opt
Cube with estimate of significant authentic signals. If provided, this
will be subtracted before projecting considering the science cube as
reference cube.
Return
-------
final_residuals_cube : List of numpy ndarray
[(ncomp is tuple or list) & (source_xy=None or full_output=True)] List
of residual final PCA frames obtained for a grid of PC values.
frame : numpy ndarray
[ncomp is scalar or source_xy!=None] 2D array, median combination of the
de-rotated/re-scaled residuals cube.
pcs : numpy ndarray
[full_output=True, source_xy=None] Principal components. Valid for
ADI cubes 3D or 4D (i.e. ``scale_list=None``). This is also returned
when ``batch`` is not None (incremental PCA).
recon_cube, recon : numpy ndarray
[full_output=True] Reconstructed cube. Valid for ADI cubes 3D or 4D
(i.e. ``scale_list=None``)
residuals_cube : numpy ndarray
[full_output=True] Residuals cube. Valid for ADI cubes 3D or 4D
(i.e. ``scale_list=None``)
residuals_cube_ : numpy ndarray
[full_output=True] Derotated residuals cube. Valid for ADI cubes 3D or
4D (i.e. ``scale_list=None``)
residuals_cube_channels : numpy ndarray
[full_output=True, adimsdi='double'] Residuals for each multi-spectral
cube. Valid for ADI+mSDI (4D) cubes (when ``scale_list`` is provided)
residuals_cube_channels_ : numpy ndarray
[full_output=True, adimsdi='double'] Derotated final residuals. Valid
for ADI+mSDI (4D) cubes (when ``scale_list`` is provided)
cube_allfr_residuals : numpy ndarray
[full_output=True, adimsdi='single'] Residuals cube (of the big cube
with channels and time processed together). Valid for ADI+mSDI (4D)
cubes (when ``scale_list`` is provided)
cube_adi_residuals : numpy ndarray
[full_output=True, adimsdi='single'] Residuals cube (of the big cube
with channels and time processed together) after de-scaling the wls.
Valid for ADI+mSDI (4D) (when ``scale_list`` is provided).
ifs_adi_frames : numpy ndarray
[full_output=True, 4D input cube, ``scale_list=None``] This is the cube
of individual ADI reductions for each channel of the IFS cube.
medians : numpy ndarray
[full_output=True, source_xy=None, batch!=None] Median images of each
batch, in incremental PCA, for 3D input cubes only.
File: ~/GitHub/VIP/vip_hci/psfsub/pca_fullfr.py
Type: function
from vip_hci.config.paramenum import SvdMode
SvdMode?
Init signature:
SvdMode(
value,
names=None,
*,
module=None,
qualname=None,
type=None,
start=1,
)
Docstring:
Define the various modes to use with SVD in PCA as constant strings.
Modes
-----
* ``LAPACK``: uses the LAPACK linear algebra library through Numpy
and it is the most conventional way of computing the SVD
(deterministic result computed on CPU).
* ``ARPACK``: uses the ARPACK Fortran libraries accessible through
Scipy (computation on CPU).
* ``EIGEN``: computes the singular vectors through the
eigendecomposition of the covariance M.M' (computation on CPU).
* ``RANDSVD``: uses the randomized_svd algorithm implemented in
Sklearn (computation on CPU), proposed in [HAL09]_.
* ``CUPY``: uses the Cupy library for GPU computation of the SVD as in
the LAPACK version. `
* ``EIGENCUPY``: offers the same method as with the ``eigen`` option
but on GPU (through Cupy).
* ``RANDCUPY``: is an adaptation of the randomized_svd algorithm,
where all the computations are done on a GPU (through Cupy). `
* ``PYTORCH``: uses the Pytorch library for GPU computation of the SVD.
* ``EIGENPYTORCH``: offers the same method as with the ``eigen``
option but on GPU (through Pytorch).
* ``RANDPYTORCH``: is an adaptation of the randomized_svd algorithm,
where all the linear algebra computations are done on a GPU
(through Pytorch).
File: ~/GitHub/VIP/vip_hci/config/paramenum.py
Type: EnumMeta
Subclasses:
6.2.3. PCA-ARDI
One can also leverage both ADI and RDI in cases where there is enough field rotation in your science cube and an available cube of reference images - this is referred to as ARDI. PCA-ARDI is best done by calling the pca_annular function in VIP, with the cube_ref argument set to your cube of reference images – although a trick can be done for a faster usage in full-frame (see below).
Mind the different default behaviours of the pca and pca_annular functions when provided a cube_ref reference cube. pca will adopt by default a RDI-only strategy (i.e., not use the science images to calculate the principal components), while pca_annular will by default follow an ARDI stratey (i.e., use both the reference images and the science images that meet a certain rotation threshold to calculate the principal components associated to each image of the cube). The latter approach is thus slower.
First let’s consider two concentric annuli with respective width set to one quarter of the image size, using 13 principal components.
from vip_hci.psfsub import pca_annular
asize = cube.shape[-1]//4 - 1
pca_ardi_fr = pca_annular(cube, angs, ncomp=13, asize=asize, fwhm=fwhm_naco, cube_ref=cube_ref)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 12:45:20
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 2, FWHM = 4.801
PCA per annulus (or annular sectors):
Ann 1 PA thresh: 2.29 Ann center: 12 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:07:52.442157
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2 PA thresh: 7.85 Ann center: 35 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:15:56.485911
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time: 0:15:59.441470
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let’s visualize the results, after masking the central 1 FWHM in the image.
plot_frames(mask_circle(pca_ardi_fr, fwhm_naco), grid=True)
We now see Beta Pic b even clearer - with a result more similar to PCA-ADI (tutorial 3 (section 3.5)).
By default, pca_annular uses both science and reference images for the PCA library, but the PCA library is recalculated for each image of the cube and for each concentric annulus (which is why it is much slower than full-frame pca). Only science images meeting a certain rotation threshold with respect to the considered image are considered for the associated PCA library. By default, the threshold is a minimum rotation of 0.1xFWHM for the innermost annulus, linearly increasing up to 1xFWHM for the outermost annulus. This behaviour is controlled with the delta_rot argument (set by default to (0.1, 1)).
The delta_rot parameter is expressed in FWHM units, and corresponds to the amount of required linear motion (in pixels) at the middle of each concentric annulus. A single value, a tuple of 2 elements, or a list of n_annuli values, can be provided for that argument corresponding to: a single threshold in rotation to be used for all annuli, the lower and upper bounds of a linearly increasing rotation threshold, and the exact rotation threshold to be used for each of all concentric annuli, respectively.
Let’s consider a few example cases - the most relevant choice for you will depend on your data:
ARDI using all science and reference images in the PCA library, in 2 concentric annuli:
pca_ardi_all_fr = pca_annular(cube, angs, ncomp=13, asize=asize, cube_ref=cube_ref, delta_rot=0)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 14:07:13
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 2, FWHM = 4.000
PCA per annulus (or annular sectors):
Ann 1 Ann center: 12 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:08:07.345447
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2 Ann center: 35 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:16:46.945749
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time: 0:16:50.021504
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Pure RDI in 2 concentric annuli. This is done by setting
delta_rotto an arbitrarily high value for no science image to meet the criterion for inclusion in the PCA library.
pca_rdi_ann_fr = pca_annular(cube, angs, ncomp=13, asize=asize, cube_ref=cube_ref, delta_rot=100)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 13:26:21
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 2, FWHM = 4.000
PCA per annulus (or annular sectors):
Ann 1 PA thresh: 173.13 Ann center: 12 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:02:34.106366
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2 PA thresh: 160.15 Ann center: 35 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:05:16.320061
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time: 0:05:19.488118
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
ARDI in a single annulus, with a rotation threshold of 1 FWHM for inclusion of the science images to the PCA library:
asize1 = cube.shape[-1]//2 - 1
pca_ardi_thr1 = pca_annular(cube, angs, ncomp=13, asize=asize1, cube_ref=cube_ref, delta_rot=1)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 13:31:40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 1, FWHM = 4.000
PCA per annulus (or annular sectors):
Ann 1 PA thresh: 9.73 Ann center: 24 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:06:54.914208
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time: 0:06:57.950685
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
ARDI in 4 concentric annuli, with a rotation threshold linearly increasing from 0.2 FWHM for the innermost annulus up to 1 FWHM for the outermost annulus:
asize4 = cube.shape[-1]//8 - 1
pca_ardi_ann4 = pca_annular(cube, angs, ncomp=13, asize=asize4, cube_ref=cube_ref, delta_rot=[0.1, 0.25, 0.5, 1])
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 13:38:38
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 4, FWHM = 4.000
PCA per annulus (or annular sectors):
Ann 1 PA thresh: 4.17 Ann center: 6 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:06:43.615960
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2 PA thresh: 3.47 Ann center: 16 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:13:38.191523
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3 PA thresh: 4.17 Ann center: 28 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:21:10.653297
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4 PA thresh: 6.11 Ann center: 38 N segments: 1
Done PCA with lapack for current annulus
Running time: 0:28:18.134590
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time: 0:28:21.212663
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
ARDI in full-frame considering all science and reference images (the fastest of all available options for ARDI). Here the trick is to redefine the reference cube to be a concatenation of the original science and reference cubes before passing it to the
pcafunction:
cube_ref_tmp = np.concatenate((cube, cube_ref))
pca_ardi_full = pca(cube, angs, ncomp=13, cube_ref=cube_ref_tmp)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-09-11 14:06:59
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 3.515 GB
Done vectorizing the frames. Matrix shape: (61, 10201)
Done vectorizing the frames. Matrix shape: (122, 10201)
Done SVD/PCA with numpy SVD (LAPACK)
Running time: 0:00:09.125319
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time: 0:00:12.217991
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
If you run ARDI without any threshold in rotation on the science images, you may get very similar results to either ADI or RDI depending on the relative number of science and reference images in each cube, as the larger the stack the more it will dominate the principal components.
Let’s visualize the results. For this, let’s first place all the processed images in a numpy array, and mask the inner 10px of each image:
pca_ardi_res = np.array((pca_ardi_fr, pca_ardi_all_fr, pca_rdi_ann_fr,
pca_ardi_thr1, pca_ardi_ann4, pca_ardi_full))
pca_ardi_res_mask = mask_circle(pca_ardi_res, 12)
plot_frames(pca_ardi_res_mask, rows=2, grid=True)
Remember that the optimal strategy and relative weighting of science vs reference images will depend on the data at hand… Is a lot of field rotation available? Are the reference images well correlated to the science images? Do the radial profile of the reference and science PSF match?
6.2.4. Contrast curves
Let’s now calculate and compare contrast curves obtained with PCA-ADI, PCA-RDI and PCA-ARDI for this dataset. For simplicity, we’ll just consider a single set-up for each case in terms of number of principal components and rotation threshold.
Note that the correct way to do this is to first estimate the parameters of any companion (see Tutorial 5A) and remove the companion from the data. We quickly do this with the cube_planet_free function in the next cell, considering the parameters of beta Pic b reported in Absil et al. (2013).
from vip_hci.fm import cube_planet_free
r_b = 0.452/pxscale_naco # Absil et al. (2013)
theta_b = 211.2+90 # Absil et al. (2013)
f_b = 648.2
cube_emp = cube_planet_free([(r_b, theta_b, f_b)], cube, angs, psfn=psfn)
Among the parameters of the contrast_curve function, one needs to provide starphot which corresponds to the (expected) flux of the star in the coronagraphic images. The latter can be obtained from the non-coronagraphic PSF (e.g. through the normalize_psf function), and after rescaling this flux to the integration time ratio used for the coronagraphic vs non-coronagraphic observations. Let’s assume here that this task was done separately:
starphot = flux # 764939.6
Another relevant parameter of the contrast curve function is algo, which takes any function in VIP for model PSF subtraction, or a user-defined imported function. Optional parameters of the algo can also be passed as regular parameters of contrast_curve.
Let’s now calculate the contrast curves:
from vip_hci.metrics import contrast_curve
cc_adi = contrast_curve(cube_emp, angs, psfn, fwhm=fwhm_naco, pxscale=pxscale_naco, starphot=starphot,
sigma=5, nbranch=3, algo=pca, ncomp=25)
cc_rdi = contrast_curve(cube_emp, angs, psfn, fwhm=fwhm_naco, pxscale=pxscale_naco, starphot=starphot,
sigma=5, nbranch=3, algo=pca, ncomp=10, cube_ref=cube_ref)
cc_ardi = contrast_curve(cube_emp, angs, psfn, fwhm=fwhm_naco, pxscale=pxscale_naco, starphot=starphot,
sigma=5, nbranch=3, algo=pca_annular, ncomp=20, cube_ref=cube_ref, asize=asize, delta_rot=0,
verbose=False)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:23:35
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
ALGO : pca, FWHM = 4.800919383981533, # BRANCHES = 3, SIGMA = 5, STARPHOT = 764939.6
Finished the throughput calculation
Running time: 0:00:19.781555
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:23:55
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
ALGO : pca, FWHM = 4.800919383981533, # BRANCHES = 3, SIGMA = 5, STARPHOT = 764939.6
Finished the throughput calculation
Running time: 0:00:20.264567
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let’s now compare the 3 contrast curves in a single plot, and also place Beta Pic b in the plot:
plt.figure(figsize=(8,5))
plt.plot(cc_adi['distance']*pxscale_naco,
-2.5*np.log10(cc_adi['sensitivity_student']),
'b-', label='5-sigma contrast (PCA-ADI)', alpha=0.5)
plt.plot(cc_rdi['distance']*pxscale_naco,
-2.5*np.log10(cc_rdi['sensitivity_student']),
'r-', label='5-sigma contrast (PCA-RDI)', alpha=0.5)
plt.plot(cc_ardi['distance']*pxscale_naco,
-2.5*np.log10(cc_ardi['sensitivity_student']),
'm-', label='5-sigma contrast (PCA-ARDI)', alpha=0.5)
plt.plot(r_b*pxscale_naco,
-2.5*np.log10(f_b/starphot), 'gs', label='Beta Pic b')
plt.gca().invert_yaxis()
plt.ylabel('Contrast (mag)')
plt.xlabel('Separation (arcsec)')
_ = plt.legend(loc='best')
plt.show()
Note that for well-correlated reference images (not the case here with our fake reference cube), the contrast achieved with PCA-RDI can become deeper than PCA-ADI at short separation, where little linear motion for field rotation is available.
Check out tutorial 4 for more examples of contrast curve calculations.
6.2.5. NEGFC
Tutorial 5A showcases in details how to use the NEGFC technique (combined with the ADI strategy) to estimate the position and flux of a directly imaged companion in your data. The typical workflow involves obtaining a first guess with a simplex algorithm, followed by using a MCMC algorithm for better estimates on the uncertainties on the parameters of the companion.
Below we just briefly illustrate how to use the NEGFC technique with the simplex (firstguess) algorithm for either RDI or ARDI strategies, but the same applies to the NEGFC algorithm.
First, let’s inject a fake companion in the data:
from vip_hci.fm import cube_inject_companions
rad_fc = 30.5
theta_fc = 240
flux_fc = 400.
gt = [rad_fc, theta_fc, flux_fc]
cubefc = cube_inject_companions(cube, psf_template=psfn, angle_list=angs, flevel=flux_fc, plsc=pxscale_naco,
rad_dists=[rad_fc], theta=theta_fc, n_branches=1)
from vip_hci.var import frame_center
cy, cx = frame_center(cube[0])
xy_test = cx+rad_fc*np.cos(np.deg2rad(theta_fc)), cy+rad_fc*np.sin(np.deg2rad(theta_fc))
xy_test
(34.749999999999986, 23.58622518457463)
6.2.5.1. PCA-RDI
By default, the algorithm used by NEGFC is PCA in a single annulus (pca_annulus) encompassing the location of the companion. This function behaves in the same way as full-frame PCA when provided with a reference cube (i.e., it triggers pure RDI). The reference cube is provided through the algo_options argument, which allows to provide additional arguments associated to the algorithm:
from vip_hci.fm import firstguess
r_rdi, theta_rdi, f_rdi = firstguess(cubefc, angs, psfn, ncomp=13, planets_xy_coord=[xy_test],
fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
aperture_radius=2, simplex=True, plot=True,
algo_options={'cube_ref':cube_ref}, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:25:11
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux | chi2r
1/30 0.100 1.025
2/30 0.149 1.025
3/30 0.221 1.025
4/30 0.329 1.024
5/30 0.489 1.024
6/30 0.728 1.023
7/30 1.083 1.022
8/30 1.610 1.021
9/30 2.395 1.018
10/30 3.562 1.014
11/30 5.298 1.007
12/30 7.880 0.999
13/30 11.721 0.984
14/30 17.433 0.962
15/30 25.929 0.932
16/30 38.566 0.892
17/30 57.362 0.833
18/30 85.317 0.752
19/30 126.896 0.631
20/30 188.739 0.426
21/30 280.722 0.200
22/30 417.532 0.053
23/30 621.017 0.332
24/30 923.671 1.098
25/30 1373.824 2.331
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 105, nfev: 229, chi2r: 0.0510203564208559
message: Optimization terminated successfully.
Planet 0 simplex result: (r, theta, f)=(30.558, 239.944, 402.470) at
(X,Y)=(34.70, 23.55)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time: 0:00:49.165751
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
print("Ground truth injected parameters:", gt)
print("Estimated parameters by NEGFC (simplex):", [r_rdi[0], theta_rdi[0], f_rdi[0]])
Ground truth injected parameters: [30.5, 240, 400.0]
Estimated parameters by NEGFC (simplex): [30.558242070331474, 239.943979957012, 402.4698017374459]
6.2.5.2. PCA-ARDI
To activate the ARDI strategy, one needs to use the pca_annular function (note the minor but important difference with pca_annulus). In this case, we also request a rotation threshold (delta_rot) to be 1 in the algorithm options dictionary.
r_ardi, theta_ardi, f_ardi = firstguess(cubefc, angs, psfn, ncomp=13, planets_xy_coord=[xy_test],
fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
aperture_radius=2, simplex=True, plot=True, algo=pca_annular,
algo_options={'cube_ref':cube_ref, 'verbose':False,
'delta_rot':1},
verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:26:00
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux | chi2r
1/30 0.100 1.890
2/30 0.149 1.890
3/30 0.221 1.889
4/30 0.329 1.888
5/30 0.489 1.887
6/30 0.728 1.885
7/30 1.083 1.882
8/30 1.610 1.877
9/30 2.395 1.871
10/30 3.562 1.861
11/30 5.298 1.846
12/30 7.880 1.824
13/30 11.721 1.791
14/30 17.433 1.744
15/30 25.929 1.674
16/30 38.566 1.570
17/30 57.362 1.426
18/30 85.317 1.212
19/30 126.896 0.928
20/30 188.739 0.592
21/30 280.722 0.233
22/30 417.532 0.035
23/30 621.017 0.474
24/30 923.671 2.723
25/30 1373.824 9.889
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 85, nfev: 195, chi2r: 0.033955191901336564
message: Optimization terminated successfully.
Planet 0 simplex result: (r, theta, f)=(30.534, 240.161, 417.392) at
(X,Y)=(34.81, 23.51)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time: 0:01:39.080224
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
print("Ground truth injected parameters:", gt)
print("Estimated parameters by NEGFC (simplex):", [r_ardi[0], theta_ardi[0], f_ardi[0]])
Ground truth injected parameters: [30.5, 240, 400.0]
Estimated parameters by NEGFC (simplex): [30.534119831048265, 240.16124692799315, 417.39168439970615]
6.3. RDI for disk imaging
6.3.1. PCA with data imputation
Now let’s try the Principal Component Analysis (PCA)-based algorithms (Amara & Quanz 2012, Soummer et al. 2012) in vip.pca.
We will more specifically use PCA-RDI with data imputation, as presented in Ren (2023). We refer to that paper for the maths and nomenclature used in this section, including the concepts of anchor and boat masks.
Let’s first define a binary mask, with zeros covering the area of the image which could host circumstellar signals. Let’s assume it covers the inner 20px radius:
from vip_hci.var import mask_circle
anchor_mask = np.ones_like(cube[0])
anchor_mask = mask_circle(anchor_mask, 20)
plot_frames(anchor_mask)
Let’s set to 30 the number of principal components ncomp considered for model creation, and let’s run the PCA algorithm (feel free to test other values):
from vip_hci.psfsub import pca
ncomp = 30
mask_rdi1 = anchor_mask # can be a 2D bool array
fr_pca_di1 = pca(cube, angs, ncomp=ncomp, cube_ref=cube_ref, mask_rdi=mask_rdi1, nproc=5)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:27:39
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 3.793 GB
Done de-rotating and combining
Running time: 0:00:00.530680
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_pca_di1, grid=True)
There are strong residuals because the reference cube is not well correlated in this example. A better approach would be to define 2 masks, anchor and boat masks, as in Ren et al. (2023), where the boat mask only covers the relevant part of the field of view you are interested in (e.g. mask the IWA of the coronagraph or the areas too far out in the field).
In the case of a single mask as above, it is assumed to correspond to the anchor mask, and that the boat mask simply covers the full image.
Let’s consider more sophisticated annular masks than above, and leverage the get_annulus_segments utility function in VIP:
from vip_hci.var import get_annulus_segments
get_annulus_segments?
Signature:
get_annulus_segments(
data,
inner_radius,
width,
nsegm=1,
theta_init=0,
optim_scale_fact=1,
mode='ind',
out=False,
)
Docstring:
Return indices or values in segments of a centered annulus.
The annulus is defined by ``inner_radius <= annulus < inner_radius+width``.
Parameters
----------
data : 2d numpy ndarray or tuple
Input 2d array (image) ot tuple with its shape.
inner_radius : float
The inner radius of the donut region.
width : float
The size of the annulus.
nsegm : int
Number of segments of annulus to be extracted.
theta_init : int
Initial azimuth [degrees] of the first segment, counting from the
positive x-axis counterclockwise.
optim_scale_fact : float
To enlarge the width of the segments, which can then be used as
optimization segments (e.g. in LOCI).
mode : {'ind', 'val', 'mask'}, optional
Controls what is returned: indices of selected pixels, values of
selected pixels, or a boolean mask.
out : bool; optional
Return all indices or values outside the centered annulus.
Returns
-------
indices : list of ndarrays
[mode='ind'] Coordinates of pixels for each annulus segment.
values : list of ndarrays
[mode='val'] Pixel values.
masked : list of ndarrays
[mode='mask'] Copy of ``data`` with masked out regions.
Note
----
Moving from ``get_annulus`` to ``get_annulus_segments``:
.. code::python
# get_annulus handles one single segment only, so note the ``[0]`` after
the call to get_annulus_segments if you want to work with one single
segment only.
get_annulus(arr, 2, 3, output_indices=True)
# is the same as
get_annulus_segments(arr, 2, 3)[0]
get_annulus(arr, inr, w, output_values=True)
# is the same as
get_annulus_segments(arr, inr, w, mode="val")[0]
get_annulus(arr, inr, w)
# is the same as
get_annulus_segments(arr, inr, w, mode="mask")[0]
# the only difference is the handling of the border values:
# get_annulus_segments is `in <= ann < out`, while get_annulus is
# `in <= ann <= out`. But that should make no difference in practice.
File: ~/GitHub/VIP/vip_hci/var/shapes.py
Type: function
coro = 12 # discard signal within ~IWA of coronagraphic mask
ring_in = 24
ring_out = 40
ones = np.ones_like(cube[0])
anchor_mask = get_annulus_segments(ones, ring_in, ring_out-ring_in, mode="mask")[0]
boat_mask = get_annulus_segments(ones, coro, ring_out-coro, mode="mask")[0]
plot_frames((anchor_mask, boat_mask))
from vip_hci.psfsub import pca
ncomp = 30
mask_rdi2 = (anchor_mask, boat_mask) # can be a tuple of 2 2D bool array
fr_pca_di2 = pca(cube, angs, ncomp=ncomp, cube_ref=cube_ref, mask_rdi=mask_rdi2)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 18:27:40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 3.964 GB
Done de-rotating and combining
Running time: 0:00:02.099145
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_pca_di2, grid=True)
6.3.2. Iterative PCA
IPCA-ADI
from vip_hci.greedy import ipca
fr_ipca_ADI = ipca(cube=cube, angle_list=angs, ncomp=15, nit=20)
Iterating...
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:02:28
%matplotlib inline
plot_frames(fr_ipca_ADI, grid=True)
To check more outputs, including the cube of images obtained at each iteration, it’s a matter of setting full_output=True:
from vip_hci.greedy import ipca
res_ipca_ADI = ipca(cube=cube, angle_list=angs, ncomp=15, nit=20, full_output=True)
Iterating...
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:03:01
fr_ipca, it_cube, sig_images, residuals_cube, residuals_cube_, stim_cube, it_cube_nd = res_ipca_ADI
plot_cubes(it_cube, grid=True)
:Dataset [x,y,time] (flux)
:Cube_shape [101, 101, 20]
IPCA-ARDI (recommended)
For all cases above, and this one below, the number of PCs is constant, and the number of iterations set by nit.
from vip_hci.greedy import ipca
res_ipca_ARDI = ipca(cube=cube, angle_list=angs, ncomp=15, nit=20, cube_ref=cube_ref, strategy='ARDI', full_output=True)
Iterating...
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:05:00
fr_ipca_ARDI, it_cube_ARDI, sig_images, residuals_cube, residuals_cube_, stim_cube, it_cube_nd = res_ipca_ARDI
%matplotlib inline
plot_frames(fr_ipca_ARDI)
plot_cubes(it_cube_ARDI, grid=True)
:Dataset [x,y,time] (flux)
:Cube_shape [101, 101, 20]
It is possible to increase the number of PCs with the iterations by setting mode='Pairet21', in this case ncomp is the final number of PCs at the last iteration, and nit is the number of iterations with each tested number of PCs:
from vip_hci.greedy import ipca
res_ipca_ARDI_incr = ipca(cube=cube, angle_list=angs, ncomp=30, nit=5, cube_ref=cube_ref, strategy='ARDI',
mode='Pairet21', full_output=True)
Iterating...
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:31:21
fr_ipca_ARDI_incr, it_cube_ARDI_incr, sig_images, residuals_cube, residuals_cube_, stim_cube, it_cube_nd = res_ipca_ARDI_incr
%matplotlib inline
plot_frames(fr_ipca_ARDI_incr, grid=True)
plot_cubes(it_cube_ARDI_incr, grid=True, vmin=-20, vmax=50)
:Dataset [x,y,time] (flux)
:Cube_shape [101, 101, 150]
IPCA-ARDI with data imputation at the first iteration (not yet much tested, but possibly even better than regular IPCA-ARDI)
from vip_hci.greedy import ipca
res_ipca_ARDI_di = ipca(cube=cube, angle_list=angs, ncomp=30, nit=42, cube_ref=cube_ref,
strategy='ARDI', mask_rdi=mask_rdi2, full_output=True)
Iterating...
0% [########################## ] 100% | ETA: 00:00:48Final Convergence criterion met after 37 iterations
An internal criterion of convergence exists, and if met the algorithm can stop earlier than the requested number of iterations. The convergence criterion can be adapted with arguments rtol and atol, the relative and absolute tolerance, respectively, which concerns the difference between the images of 2 consecutive iterations.
fr_ipca_ARDI_di, it_cube_ARDI_di, sig_images, residuals_cube, residuals_cube_, stim_cube, it_cube_nd = res_ipca_ARDI_di
%matplotlib inline
plot_frames(fr_ipca_ARDI_di, grid=True)
plot_cubes(it_cube_ARDI_di, grid=True)
:Dataset [x,y,time] (flux)
:Cube_shape [101, 101, 38]