For VIP v0.9.11
By: Carlos Alberto Gomez Gonzalez
Date: 2019/07/26
Let's assume VIP
and HCIplot
are properly installed. We import the packages along with other useful libraries:
%matplotlib inline
from matplotlib.pyplot import *
import numpy as np
import vip_hci as vip
from hciplot import plot_frames, plot_cubes
vip.__version__
'0.9.10'
The HCIplot
package aims to be the "Swiss army" solution for plotting and visualizing multi-dimensional high-contrast imaging datacubes on JupyterLab. It replaces the vip_hci.var.pp_subplots
function present in older version of VIP
. VIP
also allows working with DS9
provided that it's installed on your system along with pyds9
(again, see VIP
's documentation). A DS9
session (window) can be started with:
>>> ds9 = vip.Ds9Window()```
Then the ds9 methods will allow interact with the DS9
window. For example, for sending an in-memory array to the DS9
window you can use:
>>> ds9.display(2d_or_3d_array)```
In the same folder where this notebook is located you can find a toy ADI (Angular Differential Imaging) cube and a NACO PSF to demostrate the capabilities of VIP
. This is a L-band NACO coronagraphic (VORTEX AGPM) dataset of beta Pictoris. The sequence is heavily sub-sampled temporarily to make it smaller. Also the frames were 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 inspect the fits file:
vip.fits.info_fits('../datasets/naco_betapic_cube.fits')
Filename: ../datasets/naco_betapic_cube.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 7 (101, 101, 61) float32
It contains a cube and its correspoding parallactic angles. An ADI sequence is acquired in pupil tracking mode, which means that the instrument field derotator remains off, thereby keeping the instrument and telescope optics aligned while the image rotates with time.
Now we load the data "in" memory with the functions below:
psfnaco = '../datasets/naco_betapic_psf.fits'
cube = '../datasets/naco_betapic_cube.fits'
angs = '../datasets/naco_betapic_pa.fits'
cube_orig = vip.fits.open_fits(cube)
psf = vip.fits.open_fits(psfnaco)
angs = vip.fits.open_fits(angs)
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, 1)
The file is opened using astropy
functionlity. The fits file is now a ndarray or numpy
array on memory. If you open the file (or send it to the DS9
window) and adjust the cuts you will see the beta Pic b planet movin in a circular trajectory.
Let's measure the FWHM by fitting a 2d Gaussian to the core of the PSF:
_ = vip.var.fit_2dgaussian(psf, crop=True, cropsize=9, debug=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([4.73, 4.47])
print(fwhm_naco)
4.6
Let's normalize the flux to one in a 1xFWHM aperture and crop the PSF array:
psfn = vip.metrics.normalize_psf(psf, fwhm_naco, size=19)
Flux in 1xFWHM aperture: 1.227
plot_frames(psfn, grid=True, size_factor=4)
Let's finally define the pixel-scale for NACO (L band), which we get from a dictionaries stored in the conf
subpackage:
pxscale_naco = vip.conf.VLT_NACO['plsc']
print(pxscale_naco)
0.02719
The subpackage preproc
contains all sort of low-level image transformation routines: translation, rotation, cropping, pixel resamplling. It also contains functions dedicated for processing 3D or 4D arrays. In this tutorial we'll focuse on ADI sequences, which come in the form of 3D arrays.
For ADI sequences, it's critical to have the star aligned at the center of the frames, a task that can be particularly difficult for coronagraphic observations. Image registration is an important step for any ADI post-processing algorithm, maximizing the S/N of any potential companion. preproc
includes several re-centering functions.
If we take a look at one single image we will see the dark hole clearly off-center (here we naively assume it corresponds to the location of the star):
plot_frames(vip.preproc.frame_crop(cube_orig[0], 15), grid=True, size_factor=4)
New shape: (15, 15)
For this ADI dataset we could use either a negative 2d Gaussian fit (assuming the dark hole corresponds with the position of the star) or a Fourier cross-correlation based method (using more information including the bright donut and the speckles).
from vip_hci.preproc import cube_recenter_2dfit, cube_recenter_dft_upsampling
For the 2d Gaussian fit it's important to run the algorithm once with the debug
parameter set to True in order to check the results of the 2d Gaussian fit. You may need to adjust the position of the center of the structure to be fitted and the subi_size
(in pixels). When using the debug mode set nproc
to 1 or the plots won't show up.
cube1, shy1, shx1 = cube_recenter_2dfit(cube_orig, xy=(51, 51), fwhm=fwhm_naco, nproc=1, subi_size=5,
model='gauss', negative=True, full_output=True, debug=False)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 14:59:18 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― 2d Gauss-fitting frames 0% [##############################] 100% | ETA: 00:00:00 Total time elapsed: 00:00:00 Shifting 0% [##############################] 100% | ETA: 00:00:00 Total time elapsed: 00:00:00 Running time: 0:00:00.670850 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
cube2, shy2, shx2 = cube_recenter_dft_upsampling(cube_orig, center_fr1=(51, 51), fwhm=fwhm_naco, subi_size=5, full_output=True,
debug=False)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:00:07 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Shift for first frame X,Y=(1.188, 12.770) The rest of the frames will be shifted by cross-correlation wrt the 1st frames 0% [##############################] 100% | ETA: 00:00:00 Total time elapsed: 00:00:00 Running time: 0:00:00.761144 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Frames re-centering is a case dependent procedure. For the purpose of this tutorial let's proceed with the results of the negative 2d Gaussian fit. Let's check the first frame again after the shift was applied:
cube = cube1
im1 = vip.preproc.frame_crop(cube_orig[0], 15, verbose=False)
im2 = vip.preproc.frame_crop(cube[0], 15, verbose=False)
plot_frames((im1, im2), grid=True, label=('1st frame', '1st fr. after re-centering'))
For working with VIP
the cubes must have been calibrated previously (dark/sky subtracted and flat-fielded).
For saving any results (1, 2 or 3d arrays) in FITS format use:
>>> vip.fits.write_fits?```
preproc
also includes a few bad-pixel correction routines, several bad-frames detection functions as well as pixel upscaling routines. In vip.var
are contained high and low pass filtering functions. We won't do any further pre-processing to our toy cube.
Let's take a small detour for defining the S/N in high-contrast imaging (in the case of point-like sources). By default we adopt the definition of S/N given in Mawet el al. 2014, where it's proposed to use a two samples t-test for the problem of planet detection in the case of small sample statistics, which arise in high contrast imaging at small angles.
The main idea is to test a given speckle against the background resolution elements (at the same angular separation or radial distance from the center):
$$S/N≡ \frac{\overline{x}_1 - \overline{x}_2}{s_2\sqrt{1+\frac{1}{n_2}}},$$where $\overline{x}_1$ is the flux of the tested resolution element (blue dot in the figure below), $\overline{x}_2$ and $s_2$ are the mean and empirical standard deviation of the fluxes of the background resolution elements (red dots in the figure below) and $n_2$ the number of such background resolution elements.
vip.metrics.snr(vip.pca.pca(cube, angs, ncomp=20, verbose=False),
source_xy=(62,62), fwhm=fwhm_naco, plot=True)
11.05809676192747
VIP
also contains a function for estimating the S/N in a more classical way (without taking into account resolution elements but considering pixels as independent). In this case the S/N is taken as the ratio of pixel value of the planet (test speckle) and the noise computed as the standard deviation of the pixels in an annulus at the same radial distance from the center of the frame. The diameter of the signal aperture and the annulus width is in both cases 1 FWHM ~ 1 lambda/D.
In the metrics
subpackage we have also implemented a function for S/N map generation, by computing the S/N for each pixel of a 2D array. It has a parameter nproc
for exploiting multi-core systems which by default is None, meaning that it will get the number of physical cores available in the system.
In case we really need the S/N map of a big frame (which can take quite some time to compute depending on your hardware), a good option is to use the function snrmap_fast
. It uses an approximated S/N definition that yields close results to the one mentioned earlier.
In the framework of the differential imaging techniques, we ultimately rely on modelling and subtracting the post-coronagraphic star PSF and associated speckle noise pattern. Algorithms of different complexities and performances have been proposed in the last ~13 years. Several of those algorithms are implemented in VIP
.
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):
fr_adi = vip.medsub.median_sub(cube, angs, mode='fullfr')
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:02:03 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Median psf reference subtracted Done derotating and combining Running time: 0:00:00.074116 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
For the annular mode, using a parallactic angle threshold:
fr_adi_an = vip.medsub.median_sub(cube, angs, fwhm_naco, asize=1, mode='annular', delta_rot=0.5, nframes=4)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:02:05 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― N annuli = 50, FWHM = 4 Processing annuli: 0% [##############################] 100% | ETA: 00:00:00 Total time elapsed: 00:00:01 Optimized median psf reference subtracted Done derotating and combining Running time: 0:00:01.402177 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
b is clearly visible in the final frames:
plot_frames((fr_adi, fr_adi_an), vmin=-50, vmax=None)
Another simple approach is to perform pairwise frame subtraction, taking into account a rotation threshold. This is similar to the algorithm proposed by Hegelber et al 2015. This implementation can also impose a distance threshold (similarity)
fr_fdiff = vip.frdiff.frame_diff(cube, angs, fwhm_naco, metric='l1', dist_threshold=90,
delta_rot=0.2, radius_int=4, asize=4, nproc=1, verbose=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:02:24 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― 11 annuli. Performing pair-wise subtraction: Ann 1 PA thresh: 8.77 Ann center: 6 N segments: 1 Ann 2 PA thresh: 5.27 Ann center: 10 N segments: 1 Ann 3 PA thresh: 3.76 Ann center: 14 N segments: 1 Ann 4 PA thresh: 2.93 Ann center: 18 N segments: 1 Ann 5 PA thresh: 2.40 Ann center: 22 N segments: 1 Ann 6 PA thresh: 2.03 Ann center: 26 N segments: 1 Ann 7 PA thresh: 1.76 Ann center: 30 N segments: 1 Ann 8 PA thresh: 1.55 Ann center: 34 N segments: 1 Ann 9 PA thresh: 1.39 Ann center: 38 N segments: 1 Ann 10 PA thresh: 1.26 Ann center: 42 N segments: 1 Ann 11 PA thresh: 1.17 Ann center: 45 N segments: 1 Done processing annuli Running time: 0:00:02.102389 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_fdiff, vmin=-5)
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 linear combination of the reamining frames (after discarding the close ones).
fr_lstsq = vip.leastsq.xloci(cube, angs, fwhm=fwhm_naco, asize=4, n_segments='auto', nproc=1,
metric='correlation', dist_threshold=90, delta_rot=0.2, optim_scale_fact=2,
solver='lstsq', tol=0.01)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:07:16 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Building 12 annuli: Ann 1 PA thresh: 25.91 Ann center: 2 N segments: 2 Ann 2 PA thresh: 8.77 Ann center: 6 N segments: 3 Ann 3 PA thresh: 5.27 Ann center: 10 N segments: 5 Ann 4 PA thresh: 3.76 Ann center: 14 N segments: 7 Ann 5 PA thresh: 2.93 Ann center: 18 N segments: 9 Ann 6 PA thresh: 2.40 Ann center: 22 N segments: 11 Ann 7 PA thresh: 2.03 Ann center: 26 N segments: 12 Ann 8 PA thresh: 1.76 Ann center: 30 N segments: 14 Ann 9 PA thresh: 1.55 Ann center: 34 N segments: 16 Ann 10 PA thresh: 1.39 Ann center: 38 N segments: 18 Ann 11 PA thresh: 1.26 Ann center: 42 N segments: 20 Ann 12 PA thresh: 1.17 Ann center: 45 N segments: 22 Patch-wise least-square combination and subtraction: 0% [##############################] 100% | ETA: 00:00:00 Total time elapsed: 00:00:15 Done processing annuli Running time: 0:00:15.740182 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_lstsq, vmin=-5)
Now let's try the PCA based algorithms (Amara & Quanz 2012, Soummer et al. 2012) in vip.pca
:
fr_pca1 = vip.pca.pca(cube, angs, ncomp=4, mask_center_px=None)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:07:41 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― System total memory = 17.180 GB System available memory = 5.663 GB Done vectorizing the frames. Matrix shape: (61, 10201) Done SVD/PCA with numpy SVD (LAPACK) Running time: 0:00:00.062889 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Done de-rotating and combining Running time: 0:00:00.109022 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Here we compare the median subtraction with the PCA result (on the right). The improvement is clear:
plot_frames(fr_pca1, vmin=-5)
Clearly modeling the background (leaked starlight + static and quasi-static speckles) with PCA is better than the subtraction of the median frame. We chose 4 principal componenets (PCs) and we already improved the residual noise.
With the function below we can find the optimal number of PCs for this blob at this specific location on the image:
fr_pca2 = vip.pca.pca(cube, angs, fwhm=fwhm_naco, source_xy=(62,62),
mask_center_px=None, ncomp=(1, 61, 2))
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:12:13 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― System total memory = 17.180 GB System available memory = 7.124 GB Done SVD/PCA with numpy SVD (LAPACK) Running time: 0:00:00.061726 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Number of steps 31 Optimal number of PCs = 19, for S/N=9.289 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Coords of chosen px (X,Y) = 62.0, 62.0 Flux in a centered 1xFWHM circular aperture = 57.090 Central pixel S/N = 11.809 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Inside a centered 1xFWHM circular aperture: Mean S/N (shifting the aperture center) = 9.289 Max S/N (shifting the aperture center) = 11.809 stddev S/N (shifting the aperture center) = 2.286
Now let's look at the image computed using 19 PCs and compute the S/N maps:
vip.metrics.snr(fr_pca2, source_xy=(62, 62), fwhm=fwhm_naco)
11.809138382671737
snrmap1 = vip.metrics.snrmap(fr_pca1, fwhm_naco, plot=False)
snrmap2 = vip.metrics.snrmap(fr_pca2, fwhm_naco, plot=False)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:12:56 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― S/N map created using 2 processes Running time: 0:00:16.472500 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:13:13 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― S/N map created using 2 processes Running time: 0:00:15.091563 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames((fr_pca1, fr_pca2, snrmap1, snrmap2), rows=2, dpi=100, colorbar=True, axis=False, versp=0.01, horsp=0.2,
label=('4 PCs', '19 PCs', 'S/N map with 4 PCs', 'S/N map with 19 PCs'), label_pad=8)
With 19 PCs we've increased the S/N (although the recovered flux is smaller).
This dataset has a good range of rotation:
figure(figsize=(9,3))
plot(angs, 'o-', alpha=0.5)
ylabel('Parallactic angles')
Text(0, 0.5, 'Parallactic angles')
We can partially avoid the companion self-subtraction with full-frame PCA by applying a PA threshold for a given distance from the center (exploiting the ~80 degrees of rotation). The PCA function in VIP
accepts a parameter source_xy
for defining this location and the corresponding PA threshold. The parameter delta_rot
sets the amount of rotation in terms of the FWHM for rejecting frames from the library.
fr_pca_optlib = vip.pca.pca(cube, angs, ncomp=30, source_xy=(62,62), delta_rot=1, fwhm=fwhm_naco)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:20:27 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― System total memory = 17.180 GB System available memory = 6.608 GB Size LIB: min=30.0 / 1st QU=30.0 / med=30.0 / 3rd QU=30.0 / max=30.0 Done de-rotating and combining Running time: 0:00:00.755030 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_pca_optlib)
We can have a quick report on the frames by choosing a position with source_xy
:
vip.metrics.snr(fr_pca2, (62,62), fwhm_naco, verbose=True)
S/N for the given pixel = 11.809 Integrated flux in FWHM test aperture = 57.090 Mean of background apertures integrated fluxes = -3.348 Std-dev of background apertures integrated fluxes = 5.005
11.809138382671737
vip.metrics.snr(fr_pca_optlib, (62,62), fwhm_naco, verbose=True)
S/N for the given pixel = 6.257 Integrated flux in FWHM test aperture = 494.039 Mean of background apertures integrated fluxes = -37.807 Std-dev of background apertures integrated fluxes = 83.135
6.256747464375844
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). Although, the S/N has decreased due to the less smooth background at the same separation. The S/N is a very tricky metric!
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. Let's compute the incremental PCA on our toy cube (even if it's size does not require it):
fr_increm = vip.pca.pca(cube, angs, ncomp=19, fwhm=fwhm_naco, batch=31)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:31:08 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― System total memory = 17.180 GB System available memory = 5.354 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.109331 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Reconstructing and obtaining residuals Running time: 0:00:00.169448 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_increm)
We can exploit the range of rotation by using annular PCA approximations and applying a PA threshold for different annuli. This is the same idea used before but the PA threshold will be adjusted depending on the distance from the star. PCA can be computed in the whole annulus or in separate segments (n_segments
). The computational cost increases accordingly.
The function vip.pca.pca_annular
processes the cube in annular fashion in a serial or parallel way. The improvement in speed is noticeable in multi-core machines.
from vip_hci.pca import pca_annular
fr_pca_an15 = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp=15, asize=4,
delta_rot=0.5, nproc=None, svd_mode='lapack')
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:36:48 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― N annuli = 12, FWHM = 4.600 PCA per annulus (or annular sectors): PA threshold 59.80 is too big, will be set to 36.61 Ann 1 PA thresh: 36.61 Ann center: 2 N segments: 1 Ann 2 PA thresh: 21.70 Ann center: 6 N segments: 1 Ann 3 PA thresh: 13.12 Ann center: 10 N segments: 1 Ann 4 PA thresh: 9.39 Ann center: 14 N segments: 1 Ann 5 PA thresh: 7.31 Ann center: 18 N segments: 1 Ann 6 PA thresh: 5.98 Ann center: 22 N segments: 1 Ann 7 PA thresh: 5.07 Ann center: 26 N segments: 1 Ann 8 PA thresh: 4.39 Ann center: 30 N segments: 1 Ann 9 PA thresh: 3.87 Ann center: 34 N segments: 1 Ann 10 PA thresh: 3.47 Ann center: 38 N segments: 1 Ann 11 PA thresh: 3.14 Ann center: 42 N segments: 1 Ann 12 PA thresh: 2.93 Ann center: 45 N segments: 1 Done derotating and combining. Running time: 0:00:02.724960 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The parameter ncomp
can be set to auto for letting the algorithm define automatically the number of PCs for each patch (instead of fixing it for all of them) by minimizing the change in the standard deviation of the residuals after the subtraction of the PCA approximation, when increasing the PCs. This procedure is similar to the Least-squares/LOCI. Let's try this with 3 segments per annulus:
fr_pca_auto = pca_annular(cube, angs, fwhm=fwhm_naco, ncomp='auto', svd_mode='lapack',
asize=8, delta_rot=0.5, n_segments=3, tol=0.1, nproc=1, verbose=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:36:58 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― N annuli = 6, FWHM = 4.600 PCA per annulus (or annular sectors): Ann 1 PA thresh: 32.08 Ann center: 4 N segments: 3 Ann 2 PA thresh: 10.95 Ann center: 12 N segments: 3 Ann 3 PA thresh: 6.58 Ann center: 20 N segments: 3 Ann 4 PA thresh: 4.70 Ann center: 28 N segments: 3 Ann 5 PA thresh: 3.66 Ann center: 36 N segments: 3 Ann 6 PA thresh: 3.06 Ann center: 43 N segments: 3 Done derotating and combining. Running time: 0:00:08.516464 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames((fr_pca_an15, fr_pca_auto), dpi=100, vmin=-5, colorbar=True)
_ = vip.metrics.snr(fr_pca_an15, fwhm=fwhm_naco, source_xy=(62,62), verbose=True)
S/N for the given pixel = 8.817 Integrated flux in FWHM test aperture = 291.805 Mean of background apertures integrated fluxes = -12.014 Std-dev of background apertures integrated fluxes = 33.702
_ = vip.metrics.snr(fr_pca_auto, fwhm=fwhm_naco, source_xy=(62,62), verbose=True)
S/N for the given pixel = 12.195 Integrated flux in FWHM test aperture = 154.842 Mean of background apertures integrated fluxes = -5.335 Std-dev of background apertures integrated fluxes = 12.846
A PSF reference can be modelled using a low-rank approximation different than the PCA one. Non-negative matrix factorization serves for such purpose. In VIP's
paper there will be more details about it.
fr_nmf = vip.nmf.nmf(cube, angs, ncomp=22, max_iter=100, random_state=0, mask_center_px=None)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:40:49 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Done vectorizing the frames. Matrix shape: (61, 10201) Done NMF with sklearn.NMF. Running time: 0:00:01.262731 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Done derotating and combining. Running time: 0:00:01.452533 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
plot_frames(fr_nmf, dpi=100, colorbar=True)
_ = vip.metrics.snr(fr_nmf, fwhm=fwhm_naco, source_xy=(62,62), verbose=True)
S/N for the given pixel = 8.165 Integrated flux in FWHM test aperture = 63.638 Mean of background apertures integrated fluxes = -4.745 Std-dev of background apertures integrated fluxes = 8.191
Local Low-rank plus Sparse plus Gaussian-noise decomposition (Gomez Gonzalez et al. 2016) proposes a three term decomposition for improving the detectability of point-like sources in ADI data. This first version of our algorithm 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) allow us to increase the S/N of potential planets.
Let's try it out:
fr_llsg = vip.llsg.llsg(cube, angs, fwhm_naco, rank=6, thresh=1, max_iter=20, random_seed=10)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:40:56 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― 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:02.016635 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let's plot it along the full-frame ADI PCA result, after the optimization of the number of PCs:
plot_frames(fr_llsg)
_ = vip.metrics.snr(fr_llsg, fwhm=fwhm_naco, source_xy=(62,62), verbose=True)
S/N for the given pixel = 25.193 Integrated flux in FWHM test aperture = 177.370 Mean of background apertures integrated fluxes = -1.983 Std-dev of background apertures integrated fluxes = 6.963
The compution of the S/N map when the noise has been almost totally supressed becomes problematic (both from a theoretical and computational points of view), but nevertheless let's compute it:
snrmap_llsg = vip.metrics.snrmap(fr_llsg, fwhm_naco, plot=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:42:33 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
S/N map created using 2 processes Running time: 0:00:14.640871 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Let's take a look at the final residual frames produced by different algorithms:
plot_frames((fr_adi, fr_adi_an, fr_fdiff, fr_lstsq, fr_pca2, fr_pca_optlib, fr_pca_auto, fr_nmf, fr_llsg),
rows=3, colorbar=True, label_size=12, label_pad=8, axis=False, versp=0.05, horsp=0.15,
label=('median subtraction', 'annular median subtraction', 'pairwise frdiff', 'least-squares approx',
'ffPCA 17 PCs','ffPCA 17 PCs + PA thresh', 'annPCA auto PCs + PA thresh', 'ffNMF 20 comps', 'LLSG'))
Let's try the detection module, providing the processed frame, the FWHM in pixels and a PSF (what the planet should look like). Take a look at the help/docstring for an detailed explanation of the function. Depending on the mode
the results can be different. As a final step a S/N minimum criterion can be checked, with snr_thresh
.
We'll try a computer vision blob-detection method, the Laplacian of Gaussian filter (http://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_log) on the median subtracted frame and the annular PCA.
vip.metrics.detection(fr_adi, fwhm=fwhm_naco, psf=psfn, debug=False, mode='log', snr_thresh=5,
bkg_sigma=5, matched_filter=False)
Blobs found: 4 ycen xcen ------ ------ 71.943 71.869 63.370 51.213 54.679 60.082 55.026 52.975 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (61.9,61.9) S/N constraint NOT fulfilled (S/N = 2.266) ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (41.2,53.4) S/N constraint NOT fulfilled (S/N = 1.886) ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (50.1,44.7) S/N constraint NOT fulfilled (S/N = 3.121) ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (43.0,45.0) S/N constraint NOT fulfilled (S/N = 0.983) ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
y | x |
---|
Planet b is highlighted but with rather small S/N (~2). We can try the frame obtained with the annular PCA:
vip.metrics.detection(fr_pca_auto, fwhm=fwhm_naco, psf=psfn, bkg_sigma=5, debug=False, mode='log',
snr_thresh=5, plot=True, verbose=True)
Blobs found: 2 ycen xcen ------ ------ 71.397 71.443 55.973 66.035 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (61.4,61.4) ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Coords of chosen px (X,Y) = 61.4, 61.4 Flux in a centered 1xFWHM circular aperture = 171.785 Central pixel S/N = 13.797 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Inside a centered 1xFWHM circular aperture: Mean S/N (shifting the aperture center) = 8.734 Max S/N (shifting the aperture center) = 12.195 stddev S/N (shifting the aperture center) = 2.604 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― X,Y = (56.0,46.0) S/N constraint NOT fulfilled (S/N = 1.119) ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
y | x | |
---|---|---|
0 | 61.397496 | 61.44332 |
With VIP
we can we can inject fake companions with an arbitrary flux using the normalized PSF. Other procedures, e.g. the negative fake companion technique and the contrast curve generation, rely on the injection of fake companions.
from vip_hci.metrics import cube_inject_companions
cubefc = cube_inject_companions(cube, psf_template=psfn, angle_list=angs, flevel=400, plsc=pxscale_naco,
rad_dists=[30.5], theta=240, n_branches=1)
Branch 1: (X,Y)=(34.75, 23.59) at 0.83 arcsec (30.50 pxs from center)
Let's optimize the S/N with full-frame ADI PCA at the position of the injected companion:
fr_pca3 = vip.pca.pca(cubefc, angs, fwhm=fwhm_naco, source_xy=(35,23),
mask_center_px=None, ncomp=(1, 61, 2))
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 15:46:44 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― System total memory = 17.180 GB System available memory = 5.268 GB Done SVD/PCA with numpy SVD (LAPACK) Running time: 0:00:00.062613 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Number of steps 31 Optimal number of PCs = 15, for S/N=12.413 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Coords of chosen px (X,Y) = 35.0, 23.0 Flux in a centered 1xFWHM circular aperture = 89.794 Central pixel S/N = 14.439 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Inside a centered 1xFWHM circular aperture: Mean S/N (shifting the aperture center) = 12.413 Max S/N (shifting the aperture center) = 18.244 stddev S/N (shifting the aperture center) = 4.443
plot_frames(fr_pca3, vmin=-2, colorbar=True, circle=(35, 23))
_ = vip.metrics.snr(fr_pca3, fwhm=fwhm_naco, source_xy=(35,23), verbose=True)
S/N for the given pixel = 14.439 Integrated flux in FWHM test aperture = 89.794 Mean of background apertures integrated fluxes = -2.638 Std-dev of background apertures integrated fluxes = 6.325
VIP
implements the Negative fake companion technique for point-like sources position and flux estimation. It follows this steps:
$\chi^2 = sum(|I_j|),$ where $j \in {1,...,N}$ and $N$ the total number of pixels contained in a circular aperture. We will take 3xFWHM for the aperture radius to make sure it fits inside the image.
With the funciton vip.negfc.firstguess
we can obtain a first estimation of the flux and position by running a grid naive minimization and a Nelder-mead based minimization. Initially a naive minimization of the $\chi^2$ by trying a grid of values for the flux (parameter f_range
). If the parameter simplex
is set to True then a Nelder-Mead (simplex) minimization is done based on the preliminary guess of the naive minimization. The maximum number of iterations and error can be set with the parameter simplex_options
as a dicitionary (see scipy.minimize
function for the Nelder-Mead options).
from vip_hci.negfc import firstguess
Fisrt we define the position of the sources by examining a flux frame or S/N map. source_xy
takes a list or array of X,Y pairs like ((x1,y1),(x2,y2)...(x_n,y_n)). Let's take the (approximate) coordinates of the previously injected companion. Also, let's set the # of PCs to 20, which is higher than the optimal value (for full-frame PCA) and leads to a more aggresive reduction:
r_0, theta_0, f_0 = firstguess(cubefc, angs, psfn, ncomp=30, plsc=pxscale_naco,
planets_xy_coord=[(35, 23)], fwhm=fwhm_naco,
f_range=None, annulus_width=3, aperture_radius=3,
simplex=True, plot=True, verbose=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 16:09:05 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Planet 0 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Planet 0: flux estimation at the position [35,23], running ... Step | flux | chi2r 1/100 0.000 0.127 2/100 50.505 0.127 3/100 101.010 0.127 4/100 151.515 0.122 5/100 202.020 0.120 6/100 252.525 0.119 7/100 303.030 0.112 8/100 353.535 0.110 9/100 404.040 0.112 10/100 454.545 0.116
Planet 0: preliminary guess: (r, theta, f)=(30.9, 240.9, 353.5) Planet 0: Simplex Nelder-Mead minimization, running ... Planet 0: Success: True, nit: 77, nfev: 191, chi2r: 0.0923788054915499 message: Optimization terminated successfully. Planet 0: simplex result: (r, theta, f)=(30.622, 240.133, 371.483) at (X,Y)=(34.75, 23.45) ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― DONE ! ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Running time: 0:00:22.650508 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The Nelder-Mead optimization has found a (close-to-optimal) minimum. The estimated parameters are unfortunately not equal to the original values (radius=30.5, theta=240, flux=400), which reflects the limitations of this heuristic minization procedure (depending on the initial guess the minimization can get trapped in a different local minimum). Moreover, these are point estimates without uncertainties (error bars). We will come back to this questions later on.
Let's use the values obtained with the simplex optimization at try to subtract the planet with the function vip.negfc.cube_planet_free
.
from vip_hci.negfc import cube_planet_free
First we define a list with the parameters (r, theta, flux) is each companion that we obtained via the NegFC, in this case one:
plpar_bpicb = [(r_0, theta_0, f_0)]
cube_emp = cube_planet_free(plpar_bpicb, cubefc, angs, psfn, pxscale_naco)
fr_pca_emp = vip.pca.pca(cube_emp, angs, ncomp=19, verbose=False)
Let's take a look at the PSF of the planet in the full-frame PCA final image and the same PSF in the frame resulting of processing the planet-subtracted cube:
cropped_frame1 = vip.preproc.frame_crop(fr_pca3, cenxy=(35, 23), size=15)
New shape: (15, 15)
cropped_frame2 = vip.preproc.frame_crop(fr_pca_emp, cenxy=(35, 23), size=15)
New shape: (15, 15)
plot_frames((cropped_frame1, cropped_frame2), mode='surface', vmax=8)
plot_frames((fr_pca3, fr_pca_emp), axis=False)
Inspite of the understimated flux, the subtraction does not leave a noticeable artifact/defect.
MCMC is a more robust way of obtaining the flux and position. It samples the posterior distributions of the parameters and from them we can generate proper error bars. With parameter initial_state
we set the initial first estimation of the planets parameters (r, theta, flux). We can grabe the result of the simplex optimization.
from vip_hci.negfc import mcmc_negfc_sampling
init = np.array([30.619, 240.056, 358.064])
MCMC is a very CPU intensive procedure. We use the affine invariante sampler from emcee
which can be run in parallel. At least 100 walkers are needed for convergence on our MCMC chain, but it's a case dependent parameter.
nwalkers, itermin, itermax = (100, 200, 500)
chain = mcmc_negfc_sampling(cubefc, angs, psfn, ncomp=ncomp_fc, plsc=pxscale_naco,
fwhm=fwhm_naco, svd_mode='lapack', annulus_width=3,
aperture_radius=3, initial_state=init, nwalkers=nwalkers,
bounds=None, niteration_min=itermin, rhat_count_threshold=1,
niteration_limit=itermax, check_maxgap=50, nproc=20,
display=True, verbosity=1, save=False)
---------------------------------------------------------------------- Starting time: 2018-03-14 23:15:54 ---------------------------------------------------------------------- MCMC sampler for the NEGFC technique ----------------------------------------------------------------------
r_hat = [1.00415119 1.00374326 1.00461943] r_hat <= threshold = [ True True True] ... ==> convergence reached We break the loop because we have reached convergence Running time: 0:21:37.007073 ----------------------------------------------------------------------
Saved results, can be loaded from disk like this:
Due to the small number of walkers and iterations convergence was not achieved.
from vip_hci.negfc import show_corner_plot, show_walk_plot
show_corner_plot(chain)
show_walk_plot(chain)
For determining the most highly probable value for each model parameter, as well as the 1-sigma confidence interval:
from vip_hci.negfc import confidence
burnin = 0.3
isamples_flat = chain[:, int(chain.shape[1]//(1/burnin)):, :].reshape((-1,3))
val_max, conf = confidence(isamples_flat, cfd=68, gaussian_fit=True,
verbose=True, save=False, title='fake planet')
percentage for r: 69.48226950354612% percentage for theta: 69.70921985815603% percentage for flux: 69.68794326241135% Confidence intervals: r: 30.67055933292905 [-0.2821393746500398,0.1307475150817261] theta: 239.99866533229192 [-0.1601776469700269,0.1880346290517707] flux: 386.10415866046395 [-30.302214874376773,27.667239667909257] Gaussian fit results: r: 30.594695369471484 +-0.18326400602191925 theta: 240.0397534824626 +-0.1612041360058565 f: 382.9752887205708 +-30.15118978071996
The mean and sigma of the posterior distributions of the estimated parameters give us the parameter and its error bar. We can see that, with the MCMC approach, the mean flux value is closer to the true value. Also the true value lies within the 1-sigma error bars. Perhaps with a longer chain (for example setting rhat_count_threshold
to 3) we could obtain even better results.
VIP
allows to measure the throughput of its algorithms by injectin fake companions. That will give us an idea of how much they subtract or bias the signal from companions as a function of the distance from the center. We'll use vip.metrics.throughput
. This function requires an empty cube, for avoiding messing up with the throughput curve (due to the presence of a companion).
First, let's create a cube without the real companion (beta Pic b):
cube_emp = cube_planet_free([(16.583, 44.721, 700)], cube, angs, psfn=psfn, plsc=pxscale_naco)
plot_frames(vip.pca.pca(cube_emp, angs, ncomp=20, verbose=False), axis=False)
res_thr = vip.metrics.throughput(cube_emp, angs, psfn, fwhm_naco, pxscale_naco,
ncomp=15, algo=vip.pca.pca, nbranch=1, full_output=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 16:10:37 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Cube without fake companions processed with pca Running time: 0:00:00.100743 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Measured annulus-wise noise in resulting frame Running time: 0:00:00.134676 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Flux in 1xFWHM aperture: 1.000 Fake companions injected in branch 1 (pattern 1/3) Running time: 0:00:00.325750 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Cube with fake companions processed with pca Measuring its annulus-wise throughput Running time: 0:00:00.409095 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Fake companions injected in branch 1 (pattern 2/3) Running time: 0:00:00.613222 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Cube with fake companions processed with pca Measuring its annulus-wise throughput Running time: 0:00:00.701397 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Fake companions injected in branch 1 (pattern 3/3) Running time: 0:00:00.982051 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Cube with fake companions processed with pca Measuring its annulus-wise throughput Running time: 0:00:01.214620 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Finished measuring the throughput in 1 branches Running time: 0:00:01.231270 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
figure(figsize=(6,3))
plot(res_thr[2], res_thr[0][0,:], 'o-', alpha=0.5)
ylabel('throughput')
xlabel('separation in pixels')
Text(0.5, 0, 'separation in pixels')
Let's compare this with the annular PCA result:
res_thr2 = vip.metrics.throughput(cube_emp, angs, psfn, fwhm_naco, pxscale_naco,
algo=vip.pca.pca_annular, nbranch=1, verbose=False,
full_output=True, ncomp=15, delta_rot=0.5,
svd_mode='lapack', asize=2)
figure(figsize=(6,3))
plot(res_thr[2], res_thr[0][0,:], 'o-', label='PCA', alpha=0.5)
plot(res_thr2[2], res_thr2[0][0,:], 'o-', label='annular PCA', alpha=0.5)
ylabel('throughput')
xlabel('separation in pixels')
_ = legend(loc='best')
We clearly see the gain in throughput by applying a parallactic angle rejection in our annular PCA processing. Here the delta_rot
or amount of minimum rotation was decreased due the small number of frames, but increasing it to 1 or 1.5 (with a non temporally sub-sampled cube) will lead to better throughput.
Now let's see how to generate 5-sigma contrast curves for ADI data using vip.metrics.contrast_curve
. algo
takes any function in VIP
for model PSF subtraction. In the example below, we'll use full-frame PCA. Parameter starphot
sets the flux of the star: aperture photometry of the off-axis or non-coronagraphic star PSF (we measured it on the unnormalized PSF image).
starphot = 764939.6
_ = vip.metrics.contrast_curve(cube_emp, angs, psfn, fwhm=fwhm_naco, pxscale=pxscale_naco, starphot=starphot,
sigma=5, nbranch=1, algo=vip.pca.pca, ncomp=19, debug=True)
―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Starting time: 2019-07-26 16:11:57 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― ALGO : pca, FWHM = 4.6, # BRANCHES = 1, SIGMA = 5, STARPHOT = 764939.6 ―――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――― Finished the throughput calculation Running time: 0:00:00.982519 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――