"""
Functions for point source analysis.
"""
import math
from typing import Optional, Tuple
import numpy as np
from typeguard import typechecked
from scipy.stats import t
from scipy.ndimage import gaussian_filter
from skimage.feature import hessian_matrix
from photutils.aperture import aperture_photometry, CircularAperture
from pynpoint.util.image import shift_image, center_subpixel, pixel_distance, select_annulus, \
cartesian_to_polar, create_mask
from pynpoint.util.psf import pca_psf_subtraction
from pynpoint.util.residuals import combine_residuals
[docs]def compute_aperture_flux_elements(image: np.ndarray,
x_pos: float,
y_pos: float,
size: float,
ignore: bool):
"""
Computes the average fluxes inside apertures with the same separation from the center.
This function can be used to to estimate the residual flux of a planet at position
(x_pos, y_pos) and the respective noise elements with same separation (see function false_alarm)
It can also be used to compute the noise apertures is if no planet is present
(needed for contrast curves).
Parameters
----------
image : numpy.ndarray
The input image as a 2D numpy array. For example, this could be a residual frame returned by
a :class:`.PcaPsfSubtractionModule`.
x_pos : float
The planet position (in pixels) along the horizontal axis. The pixel coordinates of the
bottom-left corner of the image are (-0.5, -0.5). If no planet is present x_pos and y_pos
determine the separation from the center.
y_pos : float
The planet position (pix) along the vertical axis. The pixel coordinates of the bottom-left
corner of the image are (-0.5, -0.5). If no planet is present x_pos and y_pos
determine the separation from the center.
size : float
The radius of the reference apertures (in pixels). Usually, this value is chosen close to
one half of the typical FWHM of the PSF (0.514 lambda over D for a perfect Airy pattern; in
practice, however, the FWHM is often larger than this).
ignore : bool
Whether or not to ignore the immediate neighboring apertures for the noise estimate. This is
desirable in case there are "self-subtraction wings" left and right of the planet which
would bias the estimation of the noise level at the separation of the planet if not ignored.
Returns
-------
ap_phot :
A list of aperture photometry values. If a planet was present ap_phot[0] contains the flux
of the planet and ap_phot[1:] contains the noise. If not planet was present ap_phot[...]
gives the aperture photometry of the noise elements.
"""
# Compute the center of the current frame (with subpixel precision) and use it to compute the
# radius of the given position in polar coordinates (with the origin at the center of the frame)
center = center_subpixel(image)
radius = math.sqrt((center[0] - y_pos)**2 + (center[1] - x_pos)**2)
# Compute the number of apertures which we can place at the separation of the given position
num_ap = int(math.pi * radius / size)
# Compute the angles at which to place the reference apertures
ap_theta = np.linspace(0, 2 * math.pi, num_ap, endpoint=False)
# If ignore is True, delete the apertures immediately right and left of the aperture placed on
# the planet signal. These apertures often contain "self-subtraction wings", which means they
# cannot be considered to originate from the same distribution. In accordance with section 3.2
# of Mawet et al. (2014), such apertures are ignored to prevent bias.
if ignore:
num_ap -= 2
ap_theta = np.delete(ap_theta, [1, np.size(ap_theta) - 1])
# If the number of apertures is 2 or less, we cannot compute the false positive fraction
if num_ap < 3:
raise ValueError(
f'Number of apertures (num_ap={num_ap}) is too small to calculate the '
'false positive fraction.')
# Initialize a numpy array in which we will store the integrated flux of all reference apertures
ap_phot = np.zeros(num_ap)
# Loop over all reference apertures and measure the integrated flux
for i, theta in enumerate(ap_theta):
# Compute the position of the current aperture in polar coordinates and convert to Cartesian
x_tmp = center[1] + (x_pos - center[1]) * math.cos(theta) - \
(y_pos - center[0]) * math.sin(theta)
y_tmp = center[0] + (x_pos - center[1]) * math.sin(theta) + \
(y_pos - center[0]) * math.cos(theta)
# Place a circular aperture at a position and sum up the flux inside the aperture
aperture = CircularAperture((x_tmp, y_tmp), size)
phot_table = aperture_photometry(image, aperture, method='exact')
ap_phot[i] = phot_table['aperture_sum']
return ap_phot
[docs]@typechecked
def false_alarm(image: np.ndarray,
x_pos: float,
y_pos: float,
size: float,
ignore: bool) -> Tuple[float, float, float, float]:
"""
Compute the signal-to-noise ratio (SNR), which is formally defined as the test statistic of a
two-sample t-test, and related quantities (such as the FPF) at a given position in an image.
For more detailed information about the definition of the signal-to-noise ratio and the
motivation behind it, please see the following paper:
Mawet, D. et al. (2014): "Fundamental limitations of high contrast imaging set by small
sample statistics". *The Astrophysical Journal*, 792(2), 97.
DOI: `10.1088/0004-637X/792/2/97 <https://dx.doi.org/10.1088/0004-637X/792/2/97>`_.
Parameters
----------
image : numpy.ndarray
The input image as a 2D numpy array. For example, this could be a residual frame returned by
a :class:`.PcaPsfSubtractionModule`.
x_pos : float
The planet position (in pixels) along the horizontal axis. The pixel coordinates of the
bottom-left corner of the image are (-0.5, -0.5).
y_pos : float
The planet position (pix) along the vertical axis. The pixel coordinates of the bottom-left
corner of the image are (-0.5, -0.5).
size : float
The radius of the reference apertures (in pixels). Usually, this value is chosen close to
one half of the typical FWHM of the PSF (0.514 lambda over D for a perfect Airy pattern; in
practice, however, the FWHM is often larger than this).
ignore : bool
Whether or not to ignore the immediate neighboring apertures for the noise estimate. This is
desirable in case there are "self-subtraction wings" left and right of the planet which
would bias the estimation of the noise level at the separation of the planet if not ignored.
Returns
-------
signal_sum :
The integrated (summed up) flux inside the signal aperture.
Please note that this is **not** identical to the numerator of the fraction defining the SNR
(which is given by the `signal_sum` minus the mean of the noise apertures).
noise :
The denominator of the SNR, i.e., the standard deviation of the integrated flux of the noise
apertures, times a correction factor that accounts for small sample statistics.
snr :
The signal-to-noise ratio (SNR) as defined by Mawet et al. (2014) in eq. (8).
fpf :
The false positive fraction (FPF) as defined by Mawet et al. (2014) in eq. (10).
"""
ap_phot = compute_aperture_flux_elements(image=image,
x_pos=x_pos,
y_pos=y_pos,
size=size,
ignore=ignore)
# Define shortcuts to the signal and the noise aperture sums
signal_aperture = ap_phot[0]
noise_apertures = ap_phot[1:]
# Compute the "signal", that is, the numerator of the signal-to-noise ratio: According to
# eq. (8) in Mawet et al. (2014), this is given by the difference between the integrated flux
# in the signal aperture and the mean of the integrated flux in the noise apertures
signal = signal_aperture - np.mean(noise_apertures)
# Compute the "noise", that is, the denominator of the signal-to-noise-ratio: According to
# eq. (8) in Mawet et al. (2014), this is given by the standard deviation of the integrated flux
# in the noise apertures times a correction factor to account for the small sample statistics.
# NOTE: `ddof=1` is a necessary argument for np.std() in order to compute the *unbiased*
# estimate (i.e., including Bessel's corrections) of the standard deviation.
noise = np.std(noise_apertures, ddof=1) *\
math.sqrt(1 + 1 / (noise_apertures.shape[0]))
# Compute the signal-to-noise ratio by dividing the "signal" through the "noise"
snr = signal / noise
# Compute the false positive fraction (FPF). According to eq. (10) in Mawet et al. (2014), the
# FPF is given by 1 - F_nu(SNR), where F_nu is the cumulative distribution function (CDF) of a
# t-distribution with `nu = n-1` degrees of freedom (see Section 3 of Mawet et al. (2014) for
# more details on the Student's t distribution).
# For numerical reasons, we use the survival function (SF), which is defined precisely as 1-CDF,
# but may give more accurate results according to the scipy documentation.
fpf = t.sf(snr, df=(noise_apertures.shape[0] - 1))
return signal_aperture, noise, snr, fpf
[docs]@typechecked
def student_t(t_input: Tuple[str, float],
radius: float,
size: float,
ignore: bool) -> float:
"""
Function to calculate the false positive fraction for a given sigma level (Mawet et al. 2014).
Parameters
----------
t_input : tuple(str, float)
Tuple with the input type ('sigma' or 'fpf') and the input value.
radius : float
Aperture radius (in pixels).
size : float
Separation of the aperture center from the center of the frame (in pixels).
ignore : bool
Whether or not to ignore the immediate neighboring apertures of the point source to exclude
potential self-subtraction lobes.
Returns
-------
float
False positive fraction (FPF).
"""
num_ap = int(math.pi * radius / size)
if ignore:
num_ap -= 2
# Note that the number of degrees of freedom is given by nu = n-1 with n the number of samples.
# The number of samples is equal to the number of apertures minus 1 (i.e. the planet aperture).
# See Section 3 of Mawet et al. (2014) for more details on the Student's t distribution.
if t_input[0] == 'sigma':
t_result = t.sf(t_input[1], num_ap-2, loc=0., scale=1.)
elif t_input[0] == 'fpf':
t_result = t.ppf(1. - t_input[1], num_ap-2, loc=0., scale=1.)
else:
raise ValueError('First element of t_input needs to be "sigma" or "fpf"!')
return t_result
[docs]@typechecked
def fake_planet(images: np.ndarray,
psf: np.ndarray,
parang: np.ndarray,
position: Tuple[float, float],
magnitude: float,
psf_scaling: float,
interpolation: str = 'spline') -> np.ndarray:
"""
Function to inject artificial planets in a dataset.
Parameters
----------
images : numpy.ndarray
Input images (3D).
psf : numpy.ndarray
PSF template (3D).
parang : numpy.ndarray
Parallactic angles (deg).
position : tuple(float, float)
Separation (pix) and position angle (deg) measured in counterclockwise with respect to the
upward direction.
magnitude : float
Magnitude difference used to scale input PSF.
psf_scaling : float
Extra factor used to scale input PSF.
interpolation : str
Interpolation type ('spline', 'bilinear', or 'fft').
Returns
-------
numpy.ndarray
Images with artificial planet injected.
"""
sep = position[0]
ang = np.radians(position[1] + 90. - parang)
flux_ratio = 10. ** (-magnitude / 2.5)
psf = psf*psf_scaling*flux_ratio
x_shift = sep*np.cos(ang)
y_shift = sep*np.sin(ang)
im_shift = np.zeros(images.shape)
for i in range(images.shape[0]):
if psf.shape[0] == 1:
im_shift[i, ] = shift_image(psf[0, ],
(float(y_shift[i]), float(x_shift[i])),
interpolation,
mode='reflect')
else:
im_shift[i, ] = shift_image(psf[i, ],
(float(y_shift[i]), float(x_shift[i])),
interpolation,
mode='reflect')
return images + im_shift
[docs]@typechecked
def merit_function(residuals: np.ndarray,
merit: str,
aperture: Tuple[int, int, float],
sigma: float,
var_noise: Optional[float]) -> float:
"""
Function to calculate the figure of merit at a given position in the image residuals.
Parameters
----------
residuals : numpy.ndarray
Residuals of the PSF subtraction (2D).
merit : str
Figure of merit for the chi-square function ('hessian', 'poisson', or 'gaussian').
aperture : tuple(int, int, float)
Position (y, x) of the aperture center (pix) and aperture radius (pix).
sigma : float
Standard deviation (pix) of the Gaussian kernel which is used to smooth the residuals
before the chi-square is calculated.
var_noise : float, None
Variance of the noise which is required when `merit` is set to 'gaussian' or 'hessian'.
Returns
-------
float
Chi-square value.
"""
rr_grid, _, _ = pixel_distance(residuals.shape, position=(aperture[0], aperture[1]))
indices = np.where(rr_grid <= aperture[2])
if merit == 'hessian':
hessian_rr, hessian_rc, hessian_cc = hessian_matrix(image=residuals,
sigma=sigma,
mode='constant',
cval=0.,
order='rc',
use_gaussian_derivatives=False)
hes_det = (hessian_rr*hessian_cc) - (hessian_rc*hessian_rc)
chi_square = np.sum(hes_det[indices]**2)/var_noise
elif merit == 'poisson':
if sigma > 0.:
residuals = gaussian_filter(input=residuals, sigma=sigma)
chi_square = np.sum(np.abs(residuals[indices]))
elif merit == 'gaussian':
chi_square = np.sum(residuals[indices]**2)/var_noise
else:
raise ValueError('Figure of merit not recognized. Please use \'hessian\', \'poisson\' '
'or \'gaussian\'. Previous use of \'sum\' should now be set as '
'\'poisson\'.')
return chi_square
[docs]@typechecked
def pixel_variance(var_type: str,
images: np.ndarray,
parang: np.ndarray,
cent_size: Optional[float],
edge_size: Optional[float],
pca_number: int,
residuals: str,
aperture: Tuple[int, int, float],
sigma: float) -> float:
"""
Function to calculate the variance of the noise. After the PSF subtraction, images are rotated
in opposite direction of the regular derotation, therefore dispersing any companion or disk
signal. The noise is measured within an annulus.
Parameters
----------
var_type : str
Variance type ('gaussian' or 'hessian').
images : numpy.ndarray
Input images (3D).
parang : numpy.ndarray
Parallactic angles.
cent_size : float, None
Radius of the central mask (pix). No mask is used when set to None.
edge_size : float, None
Outer radius (pix) beyond which pixels are masked. No outer mask is used when set to
None.
pca_number : int
Number of principal components (PCs) used for the PSF subtraction.
residuals : str
Method for combining the residuals ('mean', 'median', 'weighted', or 'clipped').
aperture : tuple(int, int, float)
Aperture position (y, x) and radius (pix).
sigma : float, None
Standard deviation (pix) of the Gaussian kernel which is used to smooth the images.
Returns
-------
float
Variance of the pixel values. Either the variance of the pixel values ('gaussian') or
the variance of the determinant of the Hessian ('hessian').
"""
mask = create_mask(images.shape[-2:], (cent_size, edge_size))
_, im_res_derot = pca_psf_subtraction(images*mask, parang, pca_number)
res_noise = combine_residuals(residuals, im_res_derot)
sep_ang = cartesian_to_polar(center_subpixel(res_noise), aperture[0], aperture[1])
if var_type == 'gaussian':
selected = select_annulus(res_noise[0, ], sep_ang[0]-aperture[2], sep_ang[0]+aperture[2])
elif var_type == 'hessian':
hessian_rr, hessian_rc, hessian_cc = hessian_matrix(image=res_noise[0, ],
sigma=sigma,
mode='constant',
cval=0.,
order='rc',
use_gaussian_derivatives=False)
hes_det = (hessian_rr*hessian_cc) - (hessian_rc*hessian_rc)
selected = select_annulus(hes_det, sep_ang[0]-aperture[2], sep_ang[0]+aperture[2])
return float(np.var(selected))