3. PSF modeling and subtraction¶
Authors: Carlos Alberto Gomez Gonzalez and Valentin ChristiaensSuitable for VIP v1.0.0 onwardsLast update: 2022/05/04
Table of contents
- 3.1. Loading ADI data
- 3.2. median-ADI
- 3.3. Pairwise frame difference
- 3.4. Least-squares approximation - LOCI
- 3.5. PCA
- 3.6. NMF
- 3.7. LLSG
- 3.8. ANDROMEDA
- 3.9. FMMF
- 3.10. PACO
- 3.11. Summary mosaic
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:
%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__)
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
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:
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.
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:
%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
fwhm_naco = np.mean([DF_fit['fwhm_x'],DF_fit['fwhm_y']]) print(fwhm_naco)
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:
psfn = normalize_psf(psf, fwhm_naco, size=19, imlib='ndimage-fourier')
Flux in 1xFWHM aperture: 1.228
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
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
pxscale_naco = VLT_NACO['plsc'] print(pxscale_naco, "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.
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
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).
imlib_rot = 'vip-fft' interpolation=None