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
import vip_hci
print(vip_hci.__file__)

/Users/nasedkin/python-packages/VIP/vip_hci/__init__.py

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.invprob.paco import FastPACO
    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 coronagraphic (VORTEX AGPM) dataset of beta Pictoris published in Absil et al. (2013). 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: throughout this notebook, questions will be raised in orange. Corresponding answers will be provided in green 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 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):