Skip to content

Commit

Permalink
Merge pull request #395 from timothygebhard/refactor__false_alarm
Browse files Browse the repository at this point in the history 
Slightly refactor false_alarm() function
  • Loading branch information
Tomas Stolker committed Nov 5, 2019
2 parents da966c7 + 1f628bd commit b11a354
Showing 1 changed file with 82 additions and 34 deletions.
116 changes: 82 additions & 34 deletions pynpoint/util/analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,72 +25,120 @@ def false_alarm(image: np.ndarray,
size: float,
ignore: bool) -> Tuple[float, float, float, float]:
"""
Function for the formal t-test for high-contrast imaging at small working angles and the
related false positive fraction (Mawet et al. 2014).
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
Input image (2D).
The input image as a 2D numpy array. For example, this could be a residual frame returned by
a :class:`.PcaPsfSubtractionModule`.
x_pos : float
Position (pix) along the horizontal axis. The pixel coordinates of the bottom-left
corner of the image are (-0.5, -0.5).
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
Position (pix) along the vertical axis. The pixel coordinates of the bottom-left corner
of the image are (-0.5, -0.5).
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
Aperture radius (pix).
The radius of the references apertures (in pixels). Usually, this values should be chosen
close to lambda over D, that is, the typical FWHM of the PSF.
ignore : bool
Ignore the neighboring apertures for the noise estimate.
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
-------
float
Signal.
float
Noise level.
float
Signal-to-noise ratio.
float
False positive fraction (FPF).
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).
"""

# 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.)
radius = math.sqrt((center[0] - y_pos)**2 + (center[1] - x_pos)**2)

num_ap = int(math.pi*radius/size)
ap_theta = np.linspace(0, 2.*math.pi, num_ap, endpoint=False)
# 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])
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):
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)
# 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 this 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']

# Note: ddof=1 is a necessary argument in order to compute the *unbiased* estimate of the
# standard deviation, as suggested by eq. 8 of Mawet et al. (2014).
noise = np.std(ap_phot[1:], ddof=1) * math.sqrt(1.+1./float(num_ap-1))
t_test = (ap_phot[0] - np.mean(ap_phot[1:])) / noise

# 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.
return ap_phot[0], noise, t_test, 1.-t.cdf(t_test, num_ap-2)
# 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(ap_phot[1:], ddof=1) * math.sqrt(1 + 1 / (num_ap - 1))

# 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=(num_ap - 2))

return signal_aperture, noise, snr, fpf


@typechecked
Expand Down

0 comments on commit b11a354

Please sign in to comment.