1B. Quick-start with object-oriented support

Authors: Thomas Bédrine and Valentin Christiaens
Suitable for VIP v1.4.2 and Python v3.10 onwards
Last update: 2023/07/07

Table of contents

This quick-start tutorial shows:

  • how to load Angular Differential Imaging (ADI) datacubes;
  • how to use some of the algorithms implemented in VIP to model and subtract the stellar halo in order to produce final post-processed images.

Let’s assume VIP and hciplot are properly installed. We import the packages along with other useful libraries:

[1]:
%matplotlib inline
%load_ext autoreload
from matplotlib.pyplot import *
from matplotlib import pyplot as plt
import numpy as np

import vip_hci as vip
from hciplot import plot_frames, plot_cubes
from packaging import version
from dataclass_builder import update

In the following box we import all the routines that will be used in this tutorial.

Parts of this notebook will NOT run for VIP versions below 1.4.2, also make sure the Python version used is 3.10 or higher.

[2]:
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) < version.parse("1.4.2"):
    msg = "Please upgrade your version of VIP"
    msg+= "It should be 1.4.2 or above to run this notebook."
    raise ValueError(msg)

from platform import python_version
vpy = python_version()
print("Python version: ", vpy)
if version.parse(vpy) < version.parse("3.10.0"):
    msg = "To use this notebook, you'll need to upgrade your Python version to 3.10 or above."
    raise ValueError(msg)

from vip_hci.config import VLT_NACO

# common to all versions:
from vip_hci.fits import info_fits
from vip_hci.config import VLT_NACO
from vip_hci.var import frame_center
from vip_hci.objects import (MedianBuilder, PCABuilder, Dataset)
VIP version:  1.5.0
Python version:  3.11.4

The hciplot package aims to be the “Swiss army” solution for plotting and visualizing multi-dimensional high-contrast imaging datacubes. It replaces the vip_hci.var.pp_subplots function present in older version of VIP. There are some object-oriented options to display your results, but hciplot allows lots of customization and interaction with the outputs. Interactive windows can also be generated with DS9 - more information on that in 1A - Quick start.

1.1. Loading ADI data

In order to demonstrate the capabilities of VIP, you can find in the ‘dataset’ folder of the VIP_extras repository: a toy Angular Differential Imaging (ADI) datacube of images and a point spread function (PSF) associated to the NACO instrument of the Very Large Telescope (VLT). 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.

Let’s first inspect the fits file:

[3]:
info_fits('../datasets/naco_betapic_cube_cen.fits')
Filename: ../datasets/naco_betapic_cube_cen.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       7   (101, 101, 61)   float32

The fits file contains an ADI datacube of 61 images 101x101 px across.

The parallactic angle vector for this observation is saved in the naco_betapic_pa.fits file.

Now we load all the data in memory with the Dataset object:

[4]:
psfname = '../datasets/naco_betapic_psf.fits'
cubename = '../datasets/naco_betapic_cube_cen.fits'
angname = '../datasets/naco_betapic_pa.fits'

betapic = Dataset(cube=cubename, psf=psfname, angles=angname)
Cube array shape: (61, 101, 101)
Angles array shape: (61,)
PSF array shape: (39, 39)

Dataset uses the open_fits function, which is a wrapper of the astropy fits functionality. The fits file is now a ndarray or numpy array in memory that are stored as attributes of betapic.

If you open the datacube with ds9 (or send it to the DS9 window) and adjust the cuts you will see the beta Pic b planet moving on a circular trajectory, among similarly bright quasi-static speckle. Alternatively, you can use hciplot.plot_cubes:

[5]:
plot_cubes(betapic.cube, vmin=0, vmax=400)
:Dataset   [x,y,time]   (flux)
:Cube_shape     [101, 101, 61]
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/panel/pane/plot.py:347: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "orientation" which is no longer supported as of 3.3 and will become an error two minor releases later
  obj.canvas.print_figure(
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/panel/pane/plot.py:347: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "facecolor" which is no longer supported as of 3.3 and will become an error two minor releases later
  obj.canvas.print_figure(
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/panel/pane/plot.py:347: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "edgecolor" which is no longer supported as of 3.3 and will become an error two minor releases later
  obj.canvas.print_figure(
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/panel/pane/plot.py:347: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "bbox_inches_restore" which is no longer supported as of 3.3 and will become an error two minor releases later
  obj.canvas.print_figure(
[5]:

Note that the parallactic angles are not the final derotation angles to be applied to the images to match north up and east left. The following offsets, taken from Absil et al. (2013), are also added:

[6]:
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)
betapic.angles -= derot_off+TN

The Full-Width Half-Maximum (FWHM) will automatically be computed and stored as the attribute fwhm while we do another operation : normalizing the PSF. This can be done with the function normalize_psf of the Dataset object :

[7]:
betapic.normalize_psf(size=19, imlib='ndimage-fourier')
betapic.psf = betapic.psfn

Mean FWHM: 4.801
Flux in 1xFWHM aperture: 1.307
Normalized PSF array shape: (19, 19)
The attribute `psfn` contains the normalized PSF
`fwhm` attribute set to
4.801

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' allows paper-quality figures with annotations which can be saved (default), while 'bokeh' enables interactive visualization.

[8]:
%matplotlib inline
plot_frames(betapic.psf, grid=True, size_factor=4)
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "orientation" which is no longer supported as of 3.3 and will become an error two minor releases later
  fig.canvas.print_figure(bytes_io, **kw)
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "facecolor" which is no longer supported as of 3.3 and will become an error two minor releases later
  fig.canvas.print_figure(bytes_io, **kw)
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "edgecolor" which is no longer supported as of 3.3 and will become an error two minor releases later
  fig.canvas.print_figure(bytes_io, **kw)
/Users/valentin/opt/miniconda3/envs/testenv311/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "bbox_inches_restore" which is no longer supported as of 3.3 and will become an error two minor releases later
  fig.canvas.print_figure(bytes_io, **kw)
../_images/tutorials_01B_quickstart_with_objects_23_1.png

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

[9]:
betapic.px_scale = VLT_NACO["plsc"]
print(betapic.px_scale)
0.02719

1.2. Model PSF subtraction for ADI

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 usage of VIP pre-processing 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 illustrated in Tutorial 7.

[10]:
imlib = 'vip-fft' #'opencv'
interpolation=None # 'lanczos4'

1.2.1. median-ADI

The most simple approach is to model the PSF with the median of the ADI sequence, subtract it, then derotate and median-combine all residual images (Marois et al. 2006). We can do this through the PPMedianSub object, obtained with a MedianBuilder:

[11]:
medsub_obj = MedianBuilder(dataset=betapic, imlib=imlib, interpolation=interpolation).build()
medsub_obj.run()
ff_med = medsub_obj.frame_final
No changes were made to the dataset.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 00:53:15
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Median psf reference subtracted
Done derotating and combining
Running time:  0:00:01.346295
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

This can also be done in concentric annuli of annular width asize, considering a parallactic angle threshold delta_rot and the nframes closest images complying with that criterion as model for subtraction to each image of the ADI cube:

[12]:
update(medsub_obj, MedianBuilder(asize=betapic.fwhm, mode='annular', delta_rot=0.5, radius_int=4, nframes=4))
medsub_obj.run()
ann_med = medsub_obj.frame_final
No changes were made to the dataset.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 00:53:16
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 9, FWHM = 4.800919383981534
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.480336
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

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 visualize the final images:

[13]:
plot_frames((ff_med, ann_med))
../_images/tutorials_01B_quickstart_with_objects_38_0.png

Strong residuals can be seen near the center of the image, with a faint potential point source in the south-west of the final image.

Go to the top

1.2.2. Full-frame PCA

Now let’s try the Principal Component Analysis (PCA) algorithm (Amara & Quanz 2012, Soummer et al. 2012), and set to 5 the number of principal components ncomp considered for model creation.

[14]:
pca_obj = PCABuilder(dataset=betapic, ncomp=5, mask_center_px=None, imlib=imlib, interpolation=interpolation, svd_mode='lapack').build()
pca_obj.run()
No changes were made to the dataset.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 00:53:18
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 0.786 GB
Done vectorizing the frames. Matrix shape: (61, 10201)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:03.141101
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:04.407638
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
These attributes were just calculated:
        frame_final
        cube_reconstructed
        cube_residuals
        cube_residuals_der
        pcs
        cube_residuals_per_channel
        cube_residuals_per_channel_der
        cube_residuals_resc
        final_residuals_cube
        medians
        dataframe
        opt_number_pc
The following attributes can be calculated now:
        detection_map   with .make_snrmap()

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

[15]:
plot_frames(pca_obj.frame_final, grid=True)
../_images/tutorials_01B_quickstart_with_objects_45_0.png

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.

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:

[16]:
xy_b = (58.5,35.5)

Now with this point and the previously obtained frame, we are able to compute the significance of that candidate, by using the built-in function compute_significance. This function needs at least the coordinates of the candidate and a corrected frame. If no detection map was generated through make_snrmap, it will be created inside compute_significance.

[17]:
pca_obj.compute_significance(source_xy=xy_b)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 00:53:23
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
S/N map created using 5 processes
Running time:  0:00:02.738463
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
5.2 sigma detection

Congratulations! The detection is significant! You obtained a conspicuous direct image of an exoplanet! Although from this image alone (at a single wavelength) 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).

You’re now ready for more advanced tutorials. If you wish to know more about:

  • pre-processing functionalities available in VIP (star centering, bad frame trimming, bad pixel correction…) –> go to Tutorial 2;
  • more advanced post-processing algorithms –> go to Tutorial 3A or 3B (the latter for OO example usage);
  • metrics to evaluate the significance of a detection or the achieved contrast in your image –> go to Tutorial 4;
  • how to characterize a directly imaged companion –> go to Tutorial 5;
  • how to characterize a circumstellar disk –> go to Tutorial 6;
  • the difference between FFT- and interpolation-based image operations (rotation, shifts, scaling) –> go to Tutorial 7;
  • the object-oriented implementation of VIP –> go to Tutorial 8.