Bernoulli generalized likelihood ratio test for signal detection from photon counting images

. Because exoplanets are extremely dim, an electron multiplying charge-coupled device operating in photon counting (PC) mode is necessary to reduce the detector noise level and enable their detection. Typically, PC images are added together as a co-added image before processing. We present a signal detection and estimation technique that works directly with individual PC images. The method is based on the generalized likelihood ratio test (GLRT) and uses a Bernoulli distribution between PC images. The Bernoulli distribution is derived from a stochastic model for the detector, which accurately represents its noise characteristics. We show that our technique outperforms a previously used GLRT method that relies on co-added images under a Gaussian noise assumption and two detection algorithms based on signal-to-noise ratio. Furthermore, our method provides the maximum likelihood estimate of exoplanet intensity and background intensity while doing detection. It can be applied online, so it is possible to stop observations once a specified threshold is reached, providing confidence for the existence (or absence) of planets. As a result, the observation time is efficiently used. In addition to the observation time, the analysis of detection performance introduced in the paper also gives quantitative guidance on the choice of imaging parameters, such as the threshold. Lastly, though our work focuses on the example of detecting point source, the framework is widely applicable. © The


Introduction
Direct imaging of exoplanets is challenging; the flux ratio between an Earth-like exoplanet and its Sun-like host star is around 10 −10 in reflected light at visible wavelengths. 1 A starshade or internal coronagraph can suppress the starlight and leave only the planet's light to be detected; however, the planets are extremely faint and detecting them is still a challenge. An Earth-like planet ranges from 28th to 30th magnitude or fainter. As a result, the signal can be smaller than the camera read noise. An electron multiplying charge-coupled device (EMCCD) can alleviate this problem by amplifying the signal in an electron-multiplication (EM) register, thus reducing the effective readout noise to less than 1 electron. 2 Unfortunately, at the same time, a new noise is introduced-the multiplicative noise from the amplification process. This can be overcome by operating in photon counting (PC) mode. PC mode reports a value of 1 or 0 in each pixel for each integration time by thresholding the value at the final stage. The value reported in the pixel is 1 if the number of electrons in a pixel is bigger than a chosen threshold and 0 otherwise. The binary value only indicates the existence of photons in the pixels during the integration time but does *Address all correspondence to Mengya (Mia) Hu, mengyah@princeton.edu not reveal the exact number of those photons, so we need to choose the exposure time such that the expected photon count in any pixel is much less than 1. 3 Examples of simulated PC images are shown in Fig. 1 (all the simulations mentioned in this work use a 1-s integration time) and details on the generation of simulated images can be found in Ref. 4.
Operation of an EMCCD in PC mode is fairly new and thus techniques are still developing. Available literature on the design and characteristics of EMCCD detectors can be found in Refs. 2 and 5-11. Previous work on image processing focuses mainly on image stacking and Bayesian estimation methods, 12,13 which are applied to co-added PC images rather than designed for individual PC images. 14,15 As those methods generally require a high signal-to-noise ratio (SNR), large numbers of PC images are needed, which take a long observation time. In our previous work, we presented an alternative methodology for co-added PC images 16 which can efficiently detect even weak signals automatically. However, all these methods do not provide theoretical guidance on how to choose the integration times and the total number of PC images to combine into one co-added image.
In this work, we utilize a statistical model for the EMCCD to obtain a relationship between the detection probability and the photon rate. We use this distribution to formulate a Bernoulli distribution for the values of the same pixel on different PC images. Then, detection and estimation is performed. Consequently, detection and accurate intensity estimation can be achieved at low signal levels. We begin this paper with a description of the detector model, followed by the methodology of the Bernoulli generalized likelihood ratio test (BGLRT), and an application example. We finish with a comparison of the performance of our previous method, GLRT using co-added images, 16 and the new BGLRT method using PC images directly. This work builds on the previous work introduced in Ref. 17, where the method is called the sequential generalized likelihood ratio test (GLRT). We change the method name to Bernoulli generalized likelihood ratio test to better reflect that the improvement of performance is from the usage of an accurate model based on the Bernoulli distribution for the detector.

Stochastic Model for EMCCDs in Photon Counting Mode
Hirsch et al. 2 present a detailed imaging process for an EMCCD and build a statistical model for each stage of the imaging process. Their statistical model is built as follows. 2 First, incident photons follow a Poisson process. Second, the EM register is represented by a gamma distribution. Third, readout noise is added using a Gaussian distribution. Fourth, a threshold is applied to give a binary result.
Hirsch et al. use some approximations to simplify and arrive at the equation for the probability distribution for the number of electrons in a pixel, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 1 1 6 where λ is the expected number of input electrons for the EM register, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 2 3 where s is the photon rate; q is the quantum efficiency; t is the integration time; d is the dark current; CIC is the clock-induced charge; n ic is the number of counts after the EM register, including readout noise; pðn ic Þ is the probability of counts n ic ; σ is the standard deviation of the readout noise; g is the EM gain; B PC is the bias in PC mode; F χ ð2λ; 4; 2n ic g Þ is the non-central χ 2 distribution for 2λ with 4 degrees of freedom and the non-centrality parameter 2n ic g . Units for these parameters are given in Table 1.
An EMCCD decreases the read-out noise to sub-electron level by amplifying the signal before readout. However, the amplification process introduces an excess noise factor (ENF) that reaches 2 1 2 at high gains. 18 The effect on the SNR is the same as if the quantum efficiency of the EMCCD was halved. 19 This can be overcome completely by operating in PC mode. PC mode applies a single threshold at the final stage and reports a binary result. The pixel registers a detected photon if its count level is higher than the threshold, and reports no photon otherwise. This solution alleviates the ENF without making any assumption on the signal's stability across multiple images. 19 Combining models for all the stages described above produces the final model, which calculates the probability of getting a value of 1 in a detector pixel given a certain incident intensity. Thus, with a threshold of B PC þ T × σ, set by the threshold parameter T, the probability of detecting a value of 1 is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 4 6 6 fðsÞ ¼ where fðsÞ is shortened for fðs; q; t; CIC; d; g; B PC ; σ; TÞ (the detector parameters are not shown explicitly). The value in each pixel in each image only has two outcomes: either 1 or 0. The probability of getting a value of 1, fðsÞ, is decided by the ground truth of the flux, s. As the flux in each pixel is fixed for stable observations, the probability is fixed. That is to say, the measurement in each pixel in each image satisfies a Bernoulli distribution [the Bernoulli distribution is simply the probabilities of "heads" (p) and "tails" (1 − p) in a weighted coin flip]. In this paper, we assume that all the detector parameters, such as integration time, are fixed and the values are given in Table 1; so for each pixel, the imaging result follows a Bernoulli distribution whose probability of value 1 is only related to the flux value in the pixel. Figure 2 shows results for the detection probability given different flux levels and the calculated derivative of this probability. We use the detector parameters values given in Table 1, which are similar to the parameter values for the Nancy Grace Roman Space Telescope 20 and assume a starshade for starlight Threshold parameter T 5.5 -suppression. The detection probability increases as the flux increases, where a very small flux or a very large flux tend to give deterministic 0 or 1 measurements. The choice of imaging parameters (integration time, PC threshold, etc.) helps efficiently detect signals. The binary value only indicates the existence of photons in a pixel during the integration time but does not reveal the exact number of the photons. To utilize the photons efficiently, the exposure time is chosen so that the expected photon count in any pixel is less than 1. To decide a threshold for PC, we need to also pay attention to signal intensity and integration time, because they have similar influence on the result. Suppose that we have decided on the integration time and the expected range of signal intensity. To determine the imaging parameters, we can directly utilize the area under the receiver operating characteristic (ROC) curve (AUC) (details of ROC and AUC are in Sec. 3.4), which is a number reflecting the overall detection performance. For one threshold, we can calculate the average of the AUCs for different signal intensities in the expected intensity range. The threshold with the largest average AUC should be chosen.

Bernoulli Generalized Likelihood Ratio Test
The GLRT is a powerful tool for detection. The detection problem can be seen as deciding between the null hypothesis, H 0 , that there is no signal, and the alternative hypothesis, H 1 , that there is a signal. GLRT uses the ratio of the likelihood of observing the data under H 1 to the likelihood of observing the data under H 0 as the test statistic. If the ratio is large enough, we will decide H 1 . Compared to other likelihood ratio tests, the GLRT can include hypothesis tests which have unknown parameters in the two hypothesis. GLRT use maximum likelihood estimation (MLE) to estimate the unknowns and then calculates the maximum likelihood ratio. In other words, the MLEs of the unknown parameters such as signal intensity and background intensity, under both positive and negative hypotheses, are first required to be calculated. Then, a likelihood ratio test using the MLEs can be calculated and applied to decide whether a signal has been detected. Our previous work approximated the noise in co-added images by a Gaussian distribution. 16 In this work, we directly use the detector's probability function derived in the previous section, which does not rely on co-added images with a large number of PC images. As a result, measurements in each pixel in each image follow a Bernoulli distribution.

Signal Estimation
In this section and Sec. 3.2, we examine an image area that has the size of PSF's central core (the main part of a signal). We explain how to apply the method on the whole image, which is normally larger than the PSF core, in Sec. 3.3. We begin with the simple model of the signal received by the detector,  Table 1, which are similar to the parameter values for the Roman Telescope 20 and assume a starshade for starlight suppression.
where x is a reference point spread function (PSF) of the telescope; α is the intensity of the planet relative to the source intensity of the reference PSF; β is the background light (such as dust and light leakage from starshade defects); and 1 is a column vector where each element equals one. We stack pixel values in the target area into a column vector for easier mathematical manipulation. That is to say, s and x are column vectors.
In real life, due to the randomness of photons and the detector noise, we would not be able to observe s directly, but rather collect the noisy image output from the detector. Assume that N PC images have been collected, which is denoted as y ¼ fy 1 ; · · · ; y N g. Given the PSF shape, x ¼ ðx 1 ; · · · ; x K Þ T , where the subscripts of x represent the pixel indices and K is the total number of pixels in a PSF, the conditional probability of an image y n ¼ ðy n;1 ; : : : ; y n;K Þ T is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 9 2 pðy n jα; βÞ ¼ where the function fð·Þ is the probability function of detecting a one value in Eq. (3). The measurement in each pixel satisfies a Bernoulli distribution, and ½1 − fðαx k þ βÞ 1−y n;k fðαx k þ βÞ y n;k is the probability for the measurement in the k'th pixel (for PC images, y n;k is either 1 or 0). The probability for the whole image is simply the product of the probability of every measurement. The interesting part is that the probabilities in different pixels are correlated because of the underlying signal model in Eq. (4). Figure 3 illustrates the meaning of the indices in the equation. For this work, as a demonstration of the detection method, we use a PSF template as the underlying signal model. However, different signal models can be used and the later discussion of the detection method is still valid; the only difference would be detecting signal with different templates. This provides the Bernoulli GLRT with a broad applicability to various observation scenarios.
Furthermore, the conditional probability of N images, y ¼ fy 1 ; · · · ; y N g, can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 9 7 where θ is the simpler vector representation for the two unknown parameters α β ; N 0;k ¼ N − P N n¼1 y n;k and N 1;k ¼ P N n¼1 y n;k represent the frequency of 0 and 1 occurring in N measurements of the k'th pixel, respectively. The equality N ¼ N 0;k þ N 1;k holds for every pixel. As the data y is known and the parameters θ are unknown, this probability function is a likelihood function for the unknown parameters [so it is denoted as LðθÞ above].
Taking the natural logarithm of both sides of Eq. (6), the log-likelihood of the series of measurements is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 3 8 To estimate α and β, we can conduct an MLE based on Eq. (7) by taking the derivatives with respect to α and β and setting them equal to 0, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 5 6 4 ∂lðα; βÞ where f 0 ðsÞ ¼ dfðsÞ ds . These equations can be solved for estimates of α and β using a gradient descent method. We denote the solution, i.e., the MLE, with a hat, such asα andβ. As can be seen, the derivatives are weighted summations of the difference between the measurements and the distribution means, where Nfð·Þ and Nfð·Þ½1 − fð·Þ are, respectively, the mean and variance (each measurement satisfies a Bernoulli distribution). According to Fig. 2(b), the function f 0 ðsÞ is zero when s → þ∞. That means that the method neglects the difference between different intensities that are too large in the estimation since the PC detector always gives 1 in these cases.
In the problem of parameter estimation, we obtain information about the unknown parameters from the observed data of the random variables from the probability distribution governed by the parameters. The Fisher information matrix is a way to quantify the amount of information that the observable random variables carry about the unknown parameters. In our case, the Fisher information matrix is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 3 3 0 where the notation Var means the variance, E means expectation, and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 2 6 8 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 2 1 0 (11) For simpler notation, here, f is fðαx k þ βÞ; f 0 is f 0 ðαx k þ βÞ; f 00 is f 00 ðαx k þ βÞ; and f 00 ðsÞ ¼ d 2 fðsÞ d 2 s . With the Fisher information matrix, we can also derive the confidence intervals for the MLE, 21 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 9 7α and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 7 2 3β where z is the appropriate critical value (for example, 1.96 for 95% confidence), and the notation ðIðθÞ −1 Þ ii means that we invert the Fisher information matrix first and then take the ii'th component of the inverted matrix.

Signal Detection
The detection of an exoplanet signal can be modeled as a composite hypotheses testing problem. The null and alternative hypotheses are as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 5 8 5 Here, since we have no prior information about the exoplanet intensity α and the background light β, the posterior probability of the observation under the two hypotheses cannot be calculated to decide which hypothesis is more likely. One solution for this challenge is to use the maximum likelihood. We can conduct a GLRT (or also called maximum-likelihood test), i.e., we determine whether there is an exoplanet signal based on the ratio, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 4 9 4 R ¼ max α;β p 1 ðyjα; βÞ max α;β p 0 ðyjα; βÞ where p 0 ð·Þ and p 1 ð·Þ are the probability under hypotheses H 0 and H 1 , respectively, and the two maximal probabilities and the corresponding parameters (α and β) are calculated using the estimation algorithm in Sec. 3.1. R is the generalized likelihood ratio. For easier calculation, we usually take the log of R, that is, the log-likelihood ratio. In a real space mission, we can conduct detection while sequentially collecting images. After receiving every new image, we update the test ratio in Eq. (15); when R ≥ π U or R ≤ π L , we conclude H 1 or H 0 is accepted and stop taking new images, otherwise we take a new image and repeat the detection procedure, where π U and π L are thresholds chosen beforehand. The thresholds can be determined according to users' desired true positive and false alarm rates. According to Wilks' theorem, the probability distribution of 2R under H 0 , i.e., twice the ratio, is approximately a chi-squared distribution with 1 degree of freedom. 22 There is no simple closed-form theoretical solution for the true positive rate and false alarm rate given a threshold for this model, so we numerically calculate them via Monte Carlo simulation. Given the planet intensity that users want to detect and the number of PC images that they plan to take, they can numerically compute the true positive and false alarm rate for each threshold using this detection algorithm. They can choose a desirable true positive and false alarm rate value pair and thus its corresponding threshold.
We simulate a signal (signal only) with intensity 1e-8 Jy (for reference, Venus is 2.99e-8 Jy and Earth is 4.85e-9 Jy if viewed from 10 pc at 0.552 μm) and run multiple trials to get the statistics of the method's performance, shown in Fig. 4. The Bernoulli GLRT method correctly estimates the intensity of the signal quickly. The value of the log-likelihood ratio is high and increases quickly with the increasing number of observations. This means that it is possible to confidently detect the existence of a signal with just a few images and the more observations the larger the gain in confidence for the detection.
The following is a note on using Bernoulli (rather than binomial) in the method name. Each count in each pixel follows a Bernoulli distribution, while the sum of photon counts in each pixel follows a binomial distribution. Equation (6) describes the probability for a specific count sequence rather than their sum. Furthermore, the binomial distribution and Bernoulli distribution only differ by a binomial coefficient, which will not change the MLE theoretically. In our likelihood ratio in Eq. (15), the coefficient will also be cancelled out from the numerator and denominator, so Bernoulli and binomial will have the same detection result theoretically. In our algorithm, we include the PC images sequentially rather than work on a co-added image. (The MLE calculated using the first k PC images would be the initial value for the first k þ 1 PC images and the process repeats until all the N images are included.) Thus, "Bernoulli distribution" is the more accurate name for our algorithm.
Theoretically, a cold start on a co-added image using BGLRT (i.e., we directly calculate the MLE for the co-added image from all the N PC images with a chosen set of initial values) should have the same result as calculating sequentially. However, calculating sequentially converges to better results in practice.

Multi-Signal Detection in an Image
The hypothesis testing process described in the previous section assumes the planet PSF is in the center of the set of K test pixels (the size of PSF's central core), which will be called an image window from here on. In reality, an image much bigger than the size of a PSF core will be examined for an unknown number of planets in an unspecified location. Thus, the detection procedure is repeated for each set of K pixels across the whole image. This generates a loglikelihood ratio map for the image.
In this section, we demonstrate the performance of the Bernoulli GLRT detection method on simulated starshade images. The imaging system in all simulations in this paper is the Nancy Grace Roman Space Telescope with a starshade. The pixel size of the detector is 0.021 arc sec × 0.021 arc sec. The optical model to calculate the starshade diffraction uses Fresnel diffraction theory 23 to propagate light past the starshade. A detailed description of the simulation process can be found in Ref. 4.
We apply the Bernoulli GLRT method to simulated starshade images to detect possible exoplanet signals. As shown in Fig. 5(a), we observe our solar system from 10 pc away with a Starshade-Roman Telescope system, where two observable planets are in view, Venus and Earth. The Sun is present but is not recognizable because the starshade successfully suppressed its light. For this demonstration, we assume no local or extrasolar zodiacal dust. Here, we use a fixed integration time (1 s) for each single exposure and wavelength at λ ¼ 0.552 μm with a 0.128-μm bandwidth. More details about the simulation are available in Ref. 4. After including the PC mode detector described in Sec. 2, PC images are produced, an example of which is given in Fig. 5(b). As the image size is much bigger than the size of a PSF core, the detection procedure is repeated for each image window with the size of a PSF core across the whole image. Examples of an image window are also shown in Fig. 5(b). The white box is of the size of PSF. It is the area used to calculate the likelihood ratio and then decide whether there is a planet at the position of the white asterisk, i.e., pixel (7,12). The magenta box is for the magenta asterisk, i.e., pixel (8,14). Repeating the process, we will get the likelihood ratio at each pixel and thus generates a log-likelihood ratio map for the whole image. Though the method is applied to each PC image separately rather than on a co-added image, there is insufficient space to show all PC images one-by-one here. Thus, to give a sense of what the data look like, the co-added images are shown in Figs. 6(a)-6(c). Their corresponding loglikelihood ratio maps given by Bernoulli GLRT are shown in Figs. 6(d)-6(f). Venus is located at pixel (7,12) and Earth at (14,9). However, as we see in Figs. 6(d)-6(f), not only the pixel (7,12), and the pixel (14,9), where the planets are located, have high log-likelihood ratio values, but the pixels around them also have high log-likelihood ratio values compared to the background. This Examples of an image window are also shown. The white box is of the size of PSF core. It is the area used to detect whether there is a planet at the position of the white asterisk, i.e., pixel (7,12). The magenta box is for the magenta asterisk, i.e., pixel (8,14). makes sense as an image window centered on these pixels close to the signal position includes part of the signal. For example, for pixel (8,14), the magenta asterisk, its image window, the magenta box in Fig. 5(b), covers most part of Venus. Our methods only consider two cases. The first is that a signal is at the center of the window. The second one is that no signal exists in the window (only constant background). Therefore, the bright partial signal will have a high log-likelihood ratio. However, the image window only includes part of the signal, so the light distribution does not fit the PSF template perfectly. That is to say, it does not follow hypothesis 1, so the log-likelihood ratio is lower than that in the real signal pixel. To decide the exact position of the signal, after thresholding the log-likelihood map, we will just choose the center of the detected area's circumscribed circle as the position of the planet. An example is shown in Fig. 7. This planet position estimation should be more robust, especially when only a partial signal is in the image. Our method detects the existence of the signals quickly and separates them from the background successfully. The more observations, the larger the gain in confidence for the detection.
The changes of log-likelihood ratio at Venus, Earth, and a background pixel at (5, 5) after each observation are shown in Fig. 8. As the figures demonstrate, our method can detect the existence of the signals quickly and separate them from the background successfully. The corresponding false alarm rate for the log-likelihood value of Venus and Earth for the three cases are in Table 2. Details on finding the false alarm rate are given in Sec. 3.4. The false alarm rate is how likely a background pixel will have a log-likelihood equal or bigger than the chosen threshold value and thus get wrongly considered as a detection. The smaller the false alarm rate is, the more confidence we have for the detection. As shown by the example in this section, the more observations, the more confidence we gain for the detection.  When two signals overlap, the light distribution in the search area centered on the signals will be distorted. This violates the hypothesis that the image area contains a PSF-shaped signal and constant background. Therefore, the estimation and detection performance may be influenced. If there is not too much overlap, the intensity estimate for both signals will be higher than the true values and the position estimate will be biased toward each other. If the signals are close enough, the algorithm may consider them as one signal. Future work can look into solving this problem by expanding the image area and modeling overlapping signals. Moreover, astronomers can take another set of images after a while, because the signals' relative positions are likely to change and will be separate. In this work, we will not further consider the problem.
In the case of non-uniform background, the algorithm's performance depends on how well we know the non-uniform background. The current algorithm finds the signal looking like the PSF template against the uniform background in an area of the size as PSF core. We traverse the whole image by checking the PSF-core sized area one by one. Thus, the assumption of the algorithm is that the background is locally constant but can be non-uniform in the whole image. If the non-uniform background changes slowly spatially, a constant background may still be a good approximation and the detection will not be influenced much. If we have zero knowledge about the distribution of the bright background close to the planets, the results are affected by the distorted signal. However, if we have some knowledge about the non-uniform background, we can accommodate it in the model. For example, it is a reasonable approximation to assume the face-on exozodiacal light being axisymmetric. We have developed an iterative GLRT to do detection and estimation, 24 which is essentially an EM algorithm. We iteratively estimate either the planets' signal or the exozodiacal dust first, and then use the estimation as a known factor to estimate the other until both estimates converge. 24

Comparison with Other Methods
In Refs. 16 and 24, we introduced a GLRT method that was applied to co-added images generated from the combination of many frames. That method assumed the values in each pixel followed a Gaussian distribution, so the summation of frames was necessary in order to build enough signal in each pixel for the Gaussian assumption to be valid. Alternatively, the Bernoulli distribution in the new BGLRT method works for single PC images, which makes it easier to use online as it can process even a few PC images.
The underlying relationship between the probability in the Bernoulli distribution and flux intensity can be derived theoretically based on the detector model or directly measured in experiment. No approximation is used and thus the method is accurate and efficient at extracting information. In Figs. 9(a)-9(c), we show the T map using the Gaussian GLRT for the co-added images of Fig. 6. 16 The T value is a proxy for the log-likelihood ratio that is easier to compute. The T value in a pixel is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 1 6 3 where L G ðx k Þ is the log-likehood ratio of the co-added image x k under a Gaussian assumption.
The T value consists of basic operations (matrix multiplication, summation, and division) of the image and parameters, so it is easier to calculate than L G . 16,24 The T value is used to compared with a threshold to get a detection result in the Gaussian GLRT. Figure 9 shows that the T map is noisier than the log-likelihood map from the Bernoulli GLRT. The corresponding false alarm rates for the T value of Venus and Earth for the three cases are given in Table 2. We generally have higher confidence using Bernoulli GLRT. The above comparison indicates that the Bernoulli model for PC images contributes to the performance improvement. To further analyze the contribution from GLRT, we also compared the Bernoulli GLRT with the performance of detection based on the SNR from our Bernoulli model. With the MLE of signal intensity and estimated standard deviation of the MLE derived from the Fisher information matrix in Sec. 3.1, we define their ratio as the SNR from our model. We will refer to this SNR definition as Bernoulli SNR (BSNR). For each pixel, we utilize the image window centered at it and calculate the estimated signal intensity and the standard deviation as in Sec. 3.1. Then, we calculate the ratio. After repeating the process for all the pixels in the image, we obtain an SNR map. The SNR maps for the three example co-added images in Fig. 6 are shown in Fig. 10 to help visually compare the performance with the BGLRT method.
We also compared the Bernoulli GLRT with the performance of the detection method based on the SNR map implemented in pyKLIP. 25 For each pixel, the algorithm masks its surrounding pixels within the signal area in question (we chose the size of signal area as the size of PSF core) and then computes the standard deviation using the rest of pixels in concentric annuli. The ratio of the pixel value and this standard deviation is the SNR value at this pixel. The width of the annuli used is the diameter of the PSF core in this paper. Repeating the process for all the pixels in the image generates an SNR map. The SNR maps for the three example co-added images in Fig. 6 are shown in Fig. 11 to help visually compare the performance with the GLRT methods.
To further demonstrate the different properties of the methods, we compare their ROC curves for the detection of Venus and Earth. It is difficult to analytically derive a closed-form relationship between the false alarm rate and true positive rate and the threshold chosen in the new Bernoulli GLRT, unlike our previous GLRT, 16 which assumes Gaussian noise. Thus, we use 50 PC images 5 1 0 1 5 (a)  Fig. 9 Results of a Gaussian-based GLRT for images in Fig. 6. (a) T value of each pixel for the coadded image with 50 sequential PC images for Fig. 6(a). (b) T value of each pixel for the co-added image with 200 sequential PC images for Fig. 6(b). (c) T value of each pixel for the co-added image with 700 sequential PC images for Fig. 6(c).  Fig. 10 The SNR map based on the Bernoulli model for images in Fig. 6. (a) The SNR map of each pixel for the co-added image with 50 sequential PC images for Fig. 6(a). (b) The SNR map of each pixel for the co-added image with 200 sequential PC images for Fig. 6(b). (c) The SNR map of each pixel for the co-added image with 700 sequential PC images for Fig. 6(c).

PC images
Monte Carlo simulations, using multiple thresholds for each simulation, to calculate the ROC curves, shown in Fig. 12. The simulation is the same as described in Sec. 3.3, which is the starshade images with the Sun, Venus, and Earth. We compare the performance of the methods using different numbers of total images, which is denoted as n pc . We apply the Bernoulli GLRT on the set of images and obtain the log-likelihood ratio map; we also apply the Gaussian GLRT to get the false alarm rate map. We also apply BSNR and SNR from pyKLIP on the images. We apply a set of different thresholds with the resulting detection or missed detection of Earth, Venus, and a background pixel. We run a large number of trials, which is denoted as n trials , and record the ratio of detection of Earth and Venus as the true positive rate for Earth and Venus, and record the ratio of detection of the background pixel as the false positive rate. For each ROC curve, we also calculate the confidence interval shown as the shaded areas in Fig. 12. The confidence interval is standard Monte Carlo statistics, which is detailed in the Appendix. As can be seen in Fig. 12, the confidence interval is tight. For Fig. 12(a), we ran 5000 trials with 50 PC images generated for Venus and Earth in each trial. The Bernoulli GLRT and BSNR were applied to the 50 PC images sequentially; the Gaussian GLRT and SNR from pyKLIP were applied to the co-added image of the 50 PC images. Results for different decision thresholds were recorded and combined into the ROC curve. For Fig. 12(b), we ran 5000 trials and 200 PC images for each trial. For Fig. 12(c),  Fig. 11 The SNR map from the pyKLIP package 25 for images in Fig. 6. (a) SNR of each pixel for the co-added image with 50 sequential PC images as shown in Fig. 6(a). (b) SNR of each pixel for the co-added image with 200 sequential PC images as shown in Fig. 6(b). (c) SNR of each pixel for the co-added image with 700 sequential PC images as shown in Fig. 6(c).

Fig. 12
ROC curves with confidence interval for Earth and Venus detection with the four different methods. "BGLRT" means the result used the Bernoulli GLRT described in this paper, which processes each PC image sequentially using the GLRT based on a Bernoulli distribution. "BSNR" means the detection method based on our BSNR, which utilizes the final MLEs from our Bernoulli model and estimated standard deviation from Fisher information matrix. "GGLRT" in the legends means the result used the previous GLRT, 16 which processes co-added images and assumes Gaussian noise. "SNR" means the detection method is based on SNR implemented in pyKLIP, which is applied to co-added images. The shaded region behind each ROC curve is its 95% confidence interval. (a) ROC curves using 50 PC images calculated from 5000 trials. (b) ROC curves using 200 PC images calculated from 5000 trials. (c) ROC curves using 700 PC images calculated from 2000 trials.
we ran 2000 trials and 700 PC images for each trial. The performance for Venus is better than Earth when we use the same method and the same number of PC images. One reason is that Venus is brighter than Earth. Moreover, Venus is further away from the star at the center, so the influence from the residual starlight is smaller compared to that for Earth. The performance for both Venus and Earth is better with more images. It is easy to understand that the method performs better with more information, i.e., more images. Generally speaking, the Bernoulli GLRT outperforms the other methods. For easier comparison, we also list the AUC for all the curves in Table 3. AUC is an aggregate measure of performance across all possible thresholds. It can be interpreted as the probability that the model ranks a random positive example higher than a random negative example. AUC is 1 if the model's decisions are all correct and 0 if all wrong. More comparisons of GLRT with Gaussian assumption and the SNR method can be found in our other work. 24 Overall, Bernoulli GLRT outperforms GLRT with Gaussian assumption and the SNR method.

Early Stopping for Observation When No Planet Exists
As shown in the previous examples, the Bernoulli GLRT successfully detects a planet even for short integration times. Another important problem is to know when to stop when no detection is made so time wasted on a system with no planets is minimized. As no closed-form analytical solution for Bernoulli GLRT can be derived among false alarm rate and intensity and number of PC images, which will define the total time needed, it is impossible to solve for the time needed to reach a specific false alarm rate given the planet intensity via inverting a function. Instead, we will rely on numerical calculation via Monte Carlo simulation.
First, three parameters are specified: the minimum planet intensity to be detected, the maximum false alarm rate that can be accepted, and the minimum true positive rate that is acceptable. Then, for a different number of images, the ROC is calculated via Monte Carlo simulation. Finally, the minimum number of PC images that can reach the requirements are chosen. For example, if we assume the dimmest planet has the same intensity as Earth, the maximum acceptable false alarm rate is 0.16 and the minimum acceptable true positive rate is 0.85. The acceptable false alarm rate and true positive rate pairs are in the shaded green region in Fig. 13.