Statistical connection of binomial photon counting and photon averaging in high dynamic range beam-scanning microscopy.

Data from photomultiplier tubes are typically analyzed using either counting or averaging techniques, which are most accurate in the dim and bright signal limits, respectively. A statistical means of adjoining these two techniques is presented by recovering the Poisson parameter from averaged data and relating it to the statistics of binomial counting from Kissick et al. [Anal. Chem. 82, 10129 (2010)]. The point at which binomial photon counting and averaging have equal signal to noise ratios is derived. Adjoining these two techniques generates signal to noise ratios at 87% to approaching 100% of theoretical maximum across the full dynamic range of the photomultiplier tube used. The technique is demonstrated in a second harmonic generation microscope.

1. Introduction

Light intensity is often the core dependent variable in spectroscopic and microscopic measurements [1,2] and is governed by stochastic processes under low-light conditions. The number of photons generated in any one measurement is well-modeled by a Poisson distribution [3]. When one or more photons perturb a photonic detector, the corresponding voltage/current generated per photon is also governed by stochastic processes [4]. Furthermore, the photon arrival times will generally also be random unless using pulsed excitation. Precisely determining the mean of the underlying Poisson distribution can be challenging by nature of these collective random processes. This problem is notably apparent in quantitative beam-scanning microscopy, in which the mean of the Poisson distribution (λ) is desired for every pixel in the image with both high precision and routinely over the entire intrinsic dynamic range of the detector (e.g., photomultipler tube, or PMT).

Traditionally, there have been two major strategies in low-light microscopy for recovering the mean of the underlying Poisson distribution or a value proportional to it. The first common strategy is to average the detector voltage over many trials. While this does not directly generate λ, it does produce representative values for the light intensity that are proportional to λ. Of course, the precision of the measurements depends on both the signal and the noise. This strategy can routinely recover values for the signal to noise ratio (SNR) that approach the theoretical limit in the bright signal limit. However, the instrument dark noise, the Johnson noise (electronic white noise), and the intrinsic variability in the detector gain and in the corresponding voltage peak height distribution (PHD) can be significant or dominant in the low-light limit [5]. In cases such as these, photon counting strategies have the advantage of reducing contributions from Johnson noise and, more importantly, from the intrinsic variability arising from the detector PHD. Photons are counted by recording the number of times the PMT signal exceeds a minimum threshold value per trial, where the threshold is set to a value slightly larger than the Johnson noise floor [5]. Provided that there is no more than one photon arriving at the detector per measurement event, the number of counts divided by the number of trials directly gives λ. However, the two values depart when the probability of observing more than one photon per window becomes significant, in which case conventional photon counting introduces bias underestimating the value of λ.

Several attempts have been made to extend the usable range of photon counting out of the conventional Poisson counting limit. The most common approaches usually rely on diminishing the light that reaches the detector, either by attenuating the total signal or by dividing it across multiple detectors [6,7]. Unfortunately, all such methods suffer SNR losses in the high photon flux regime. The SNR in counting measurements scales with the square root of the number of counts, such that reducing the signal reduces the SNR accordingly.

Soukka et al. [8] took a different approach by using a high performance channel photomultiplier that had exceptionally low variance in the output voltage PHD, which allowed three counting thresholds to be carefully placed between each voltage distribution and extending the dynamic range of photon counting by ~three-fold. However, even with three thresholds, the increased variance in the PHDs with increasing number of photons n quickly resulted in overlap between the PHDs of multiple photon events through convolution, such that in practice only a relatively small number of thresholds can be used reliably with this approach and only for detectors with exceptionally low variance in the one-photon PHD.

Kissick et al. [9] showed that photon counting can be extended beyond the low-signal limit for time-coincident photon detection by embracing and accounting for the multiple convolutions present in the net PHD as the photon flux increased. This feat was accomplished by recognizing that voltage discrimination conforms to a binomial distribution at both high and low photon fluxes. The counts observed for an arbitrary threshold setting were shown to be related to the cumulative density function (CDF) of the total net Poisson-weighted PHD, which in turn can be analytically generated from the one-photon PHD and λ. Using binomial statistics, the probability of success p or of failure q describes the likelihood of a voltage transient exceeding the discriminator threshold or not, respectively, for each pulse of the laser. Therefore, provided the one-photon PHD is known, the measured counts for successes (or failures) from any threshold value can in principle be related back to λ in both the low-light and high-light regimes. The SNR in binomial statistics is generally optimized when half of the events (e.g., laser pulses) produce a positive outcome (i.e., a voltage transient exceeding a threshold setting). Importantly, this approach removes the measurement bias introduced with conventional photon counting performed outside the Poisson limit. With this method, the full ~8 decade intrinsic dynamic range of the PMT could be accessed with SNR approaching the theoretical limit by using multiple threshold values. Practical implementation of this multiple-threshold approach in microscopy measurements is complicated in part by the need to split the amplified electrical signal into multiple closely calibrated discriminators. Furthermore, the number of calculations required to iteratively perform the multiple numerical integrations over the Poisson-weighted PHD in order to extract λ is arguably prohibitively time-consuming for microscopy applications, where such calculations would be required for each pixel. Therefore, the binomial counting technique is most conveniently utilized with a single threshold or a small number of thresholds [10,11] while measuring dim to moderate signals. However, bright signals are challenging to reliably quantify in a single-threshold approach as the probability of a laser pulse producing zero counts approaches zero. In this limit, the uncertainty in λ approaches infinity, and the SNR approaches zero.

One alternative way of spanning the entire dynamic range of the detector is to simultaneously perform conventional counting and averaging of the signal, and then selectively preserve the answer that has the highest SNR. Nau and Nieman [12] have created an electronic device capable of such a task in hardware. The device has a manually adjustable output gain for normalization of the methods as well as a manually adjustable crossover point. A hardware switch toggles between saving either conventional photon counting data or signal averaging data, depending on the intensity of the detected signal. Staton et al. [13] patented a similar device. The drawbacks of these systems are that the calibration can be defeated by fluctuation in the gain or DC offset of the PMT [14]. Further, averaging by RC integration requires time on the order of microseconds to discharge to baseline after every integration even for modern integrating circuits, which is too slow to be compatible with modern ~100 MHz ultrafast pulsed lasers [15]. Finally and most importantly, the implemented approach can only be successful if the linear dynamic ranges of both conventional photon counting and signal averaging significantly overlap. In practice, this criterion can be challenging to meet due to the bias introduced in conventional photon counting at high photon flux.

In the present study, the statistical analysis of binomial photon counting developed by Kissick et al. [9] was combined with signal averaging in second harmonic generation (SHG) microscopy through entirely digital means post-acquisition. The peak of every PMT signal event was directly flash analog-to-digital-converted (ADC) synchronously with an ultrafast pulsed laser. SHG is particularly suited to characterize this technique since all detected signal photons are synchronous with the laser pulse, which filters out many of the temporally random dark counts of the detector and largely removes the need to account for convolution with a temporal response function associated with random photon arrival times. Further, the high peak powers in ultrafast pulsed lasers generate prominent and discrete SHG signal events at repetition rates approaching and exceeding 100 MHz, allowing testing of the data throughput and processing requirements for this technique. The captured data were analyzed with a personal computer. Binomial counting and signal averaging were joined by mathematically extracting the mean and experimental SNR of the underlying Poisson distribution from the averaged signal, effectively allowing precise extrapolation of photon counting through the entire dynamic range of the detector. A distinction is made between ‘signal averaging’ and ‘photon averaging,’ since the most likely value of λ is recovered from this averaging technique, instead of a value proportional to λ but with an unknown proportionality constant as in conventional averaging of the detector signal. The preferential crossover point of the data analysis between photon counting and photon averaging is derived analytically based upon SNR maximization. Preliminary efforts for this technique were centered in SHG microscopy measurements, where each pixel of the image was analyzed to optimize the SNR over ranges spanning photon averaging and photon counting. SHG is particularly suited to flash ADC of the signal transients since all detected signal photons are synchronous with the laser pulse.

2. Adjoining binomial photon counting and photon averaging: theoretical foundation

The underlying Poisson distribution in any pixel can be described entirely by its average value, λ, a good estimate of which is the value that is sought to describe the spectrochemical measurement. The probability of n photons being generated is The mean and variance of the Poisson distribution are both equal to λ.

When these photons reach a PMT detector and discharge electrons, the number of electrons that are ejected off of each dynode in the chain per detected photoelectron is again probabilistic. Courtesy of the central limit theorem, products of many random events converge to lognormal probability density functions [4], which describe the distribution of voltage expected per detected photoelectron. For n photoelectrons generated by the detector, the n-photon probability density function (PDF) is described as the single photon lognormal PDF convolved with itself n-1 times. Although no closed form solution exists for the convolution of two lognormal PDFs, the resulting distribution is often well-approximated by another lognormal PDF [16]. The overall PDF of detected voltage for any point in the sample is a combination of these Poisson and lognormal processes, and can be intuited as the linear combination of each n-photon lognormal PDF, weighted by the Poisson probability of event n. In the binomial counting technique, the probability of observing a count, p, is the probability of a signal event exceeding a user-defined threshold: In Eqs. (3) and (4), μ and σn correspond to the mean and standard deviation, respectively, of the n-photon voltage lognormal PDF describing the peak height distribution, while M and S are the standard lognormal parameters, corresponding to the mean and standard deviation, respectively, of ln(V).

In the case where the detection threshold is set low enough to observe essentially every photoelectron event, the normalized lognormal distribution integrates to one. In this limit, Kissick et al. [9] derived the binomial counting equation [8,9,17] by integrating the Poisson distribution over all values greater than zero, and solving for λ: where the variance in the recovery of λ was found to be [9] In Eq. (6), N is the number of events capable of producing photons, in this case equal to the number of laser pulses.

The SNR of this technique can be found by dividing the mean of the signal by the standard deviation in the recovery of λ: In the limit of small p, , , , and , consistent with conventional photon counting in the Poisson limit.

Similarly, the parameter λ can also be extracted from signal averaged analyses. The mean voltage of the pooled samples μ is found in signal averaging, which analytically corresponds to the expectation value of the underlying distribution. From the standard definition, the expectation value is the integral of a random variable multiplied by its probability density function (PDF). For time-coincident photon detection, the peak PMT voltage, V, is the measured random variable, and its expectation value µ is given in Eq. (8): Note that within this equation, the integral of V multiplied by the lognormal distribution of V will lead to the expectation value of that lognormal distribution for all n.

It was briefly described in Eqs. (3) and (4) that the mean and variance scale linearly with the number of convolutions, n. Using Eq. (3) as the expectation value for the lognormal distribution, Eq. (8) reduces to A further simplification can be applied to Eq. (9) by recognizing that multiplying the Poisson distribution by n and summing over all values results in yet another expectation value. The mean of the Poisson distribution is λ, so Eq. (9) reduces to the result we could have expected from the outset: Thus, dividing the measured mean voltage by the mean of the one photon lognormal PDF will recover the mean value of the underlying Poisson distribution, effectively providing a “calibration”to relate the averaged signal directly back to the average number of photons. This result is consistent with the expectation value for a sum of a random amount of random numbers [18], where a Poisson random number of lognormal random variables are summed.

Although this explicit analysis recovers the fairly trivial result in Eq. (10) for the mean signal, its utility is clearer in comparisons of the noise. The variance for a Poisson weighting of lognormal PDFs is easily derived through the standard definition of variance for any PDF, : This equation can be simplified by substituting the definition of the second moment of the n-photon lognormal distribution, , where is defined in Eq. (4). Equation (11) reduces to The n in the first term will combine with the infinite sum and reduce to λ. The n2 in the second term will combine with the infinite sum and reduce to . After substitution and combining terms, Eq. (12) reduces to This result is also consistent with the variance for a sum of a random amount of random numbers [18], where a Poisson random number of lognormal random variables are summed. The variance of the Johnson noise mixed in the signal will add linearly to this variance. In many signal averaging techniques, the signal is integrated over a number of trials, and the relative variance of an integrated signal increases linearly with the number of samples taken. Therefore, the expression for the sample variance is The SNR for signal averaging is obtained by dividing the sample mean from Eq. (10) for N laser pulses by the standard deviation from Eq. (14): The maximum SNR theoretically achievable is defined by SNR of the underlying Poisson distribution, . In the high signal limit or when is negligible, the ratio SNR goes to Therefore, the performance of signal averaging in comparison to the theoretical maximum is defined entirely by the performance of the detector. This result is in agreement with simulation, as shown in Fig. 1 .

Fig. 1

(Left) SNR of photon averaging and binary counting. (Right) SNR of photon averaging and binary counting as a ratio of the theoretical maximum SNR (the SNR of the underlying Poisson distribution). SNR at the crossover point is at ~87% of the theoretical limit. μ1 = 7.2 mV, σ1 = 4 mV, σJ = 0.3 mV. Data were simulated per sample by summing a Poisson distributed random number of lognormal random numbers, with an additional normally distributed random number also added to represent the Johnson noise (between 2 × 106 values at low λ to 5000 at high λ). Binary counting offers higher SNR at low λ, and photon counting offers higher SNR at high λ.

Using Eqs. (5) and (10), λ can be directly extracted from counting data as well as from photon averaged data. Photon averaging will provide a higher SNR in the bright signal regime, where the binomial counting algorithm in Eq. (5) poorly resolves λ as p approaches one. For the dim signal regime (λ< 0.5), the binomial counting algorithm offers a considerable SNR advantage. For such a preferential analysis, there is some intermediate signal value where the SNR for both binomial counting and signal averaging are equal, and that value of λ will be the preferential crossover point for the analysis. This crossover point is pictorially described in Fig. 1. It can be derived by comparing the SNR for photon averaging to the SNR for binomial photon counting.

Setting the SNR for photon averaging in Eq. (15) equal to the SNR for photon counting in Eq. (7) will provide the preferential analysis crossover point: For , can be approximated as being negligible, and the preferential crossover point for the PMT one photon parameters can be described ratiometrically as a function of λ: The crossover value of λ can be solved iteratively, or the expression can be approximated with an exponential function, and the crossover point can be easily calculated over a wide range of lognormal distribution parameters as described in Eq. (19) and shown in Fig. 2 : where 0.75 < 𝜇1/𝜎1 < 10. Therefore, if either counting or averaging yields a value of λ that is less than that obtained through Eq. (19), binomial photon counting will produce a larger SNR. Conversely, for signals greater than this threshold value, photon averaging offers a SNR advantage. In the limit of negligible Johnson noise, the preferential crossover point is defined entirely by the mean and variance per photon in the response of the PMT.

Fig. 2

Approximation of crossover point from a power function, where 0.75 < 𝜇1/𝜎1 < 10

3. Experimental

The ~800 nm output of an 80 MHz Ti:Sapphire laser (Spectra-Physics, Mountain View, CA) was directed into an in-house built microscope which scanned the beam across a crystalline urea sample surface by a 4 kHz resonant mirror/galvanometer mirror pair. The SHG was collected in transmission (Hamamatsu H10721-01), with band-pass filters (HQ 400/20m-2p; Chroma Technology, Bellows Falls, VT) to reject the incident beam. Urea samples were prepared by solvent evaporation of concentrated urea (Sigma-Aldrich U5128-100G). A flash digitizer (AlazarTech ATS-9462) was externally clocked by the laser photodiode to flash ADC the peak of every SHG signal event from the PMTs. A pulse generator module (Quantum Composer 9530) was phase locked to the laser and used to clock the resonant mirror, as well as provide a trigger for the digitizer card. The galvanometer mirror was controlled with a National Instruments USB DAC (NI 9263). All computation and computer interfacing was accomplished using in-house developed software in Python (www.python.org). PMT calibration, λ extraction from the raw data, and λ image compilation was performed in MATLAB. Minor post-processing image cropping and contrast adjustments were accomplished with ImageJ [19].

4. Adjoining binomial photon counting and photon averaging: application in SHG microscopy

Binomial counting and photon averaging techniques were mathematically adjoined to maximize SNR in every pixel of SHG microscopy images. A 512 × 512 pixel image of the urea sample was captured with 500 samples per pixel, amounting to 250MB of raw data acquired and analyzed in 6.4 seconds. Binomial counting and averaging analyses were used to extract the λ values from every pixel in the image in realtime, and a final image was then compiled through the preferential analysis. Calibration of the one-photon PMT lognormal PHD was performed within the image data by recovering λ with binomial counting for pixels where SNR is maximized (p = 0.8), and then using Eq. (10) to determine the mean of the one photon distribution.

Images of λ values recovered through the photon averaging analysis are shown in Fig. 3A and 3B. Photon averaging accesses the upper end of the dynamic range up to the point of PMT saturation, particularly when the background Johnson noise is small. By only performing the digitization at the peak of the voltage distribution following the firing of the laser, the background noise remains least susceptible to noise from integration over a drifting DC background. Even parts of the sample that had only one or two detected photons out of 500 trials are evident under these optimized photon averaging conditions. Nevertheless, instrument noise (primarily from comparatively slow DC drift) is clearly observable in the image as horizontal streaks over what should ideally be a completely black background.

Fig. 3

SHG images of crystalline urea. (Left Column) Full contrast image, (Right Column) Contrast adjusted to λmax = 0.02. (A,B) Analysis with photon averaging only. The majority of the image is silhouetted in comparison to the brightest pixels in (A), up to λ = 74. The dimmest pixels are evident in (B), but are as prominent as the horizontal streaks and noise in the image. (C,D) Analysis with binomial counting only. The largest recoverable value was λ = 6.23; pixels brighter than this were clipped to this value in (C). The dimmest pixels are easily identifiable in (D), and the instrument noise is not evident. (E,F) Preferential crossover analysis incorporating photon averaging and binomial counting. Pixels brighter than λ = 0.48 were analyzed by photon averaging, so the full upper range of detection is preserved in (E). Pixels dimmer than λ = 0.48 were analyzed by binomial counting, so the lower range of detection is preserved in (F). Selection of a crossover point defined by Eq. (18) maximizes SNR across the entire range of detection.

The same set of raw data was used to generate an image of λ values using binomial photon counting with a threshold set to a value slightly larger than the noise floor, shown in Fig. 3C and 3D. The largest recoverable value of λ was 6.23, after which the binomial counting algorithm returned divergent answers (i.e., every laser pulse produced at least one photon in that pixel). However, counting removed the background noise in the weak signal limit.

A final image was compiled to maximize SNR on a pixel by pixel basis by preferentially using binomial counting for dim pixels and photon averaging for bright pixels shown in Fig. 3E and 3F. The PMTs used in this experiment were found to have a single-photon mean response of 7.2 mV with a variance of 16 mV, and a Johnson noise floor standard deviation of 0.4 mV. From Eq. (19), the preferential crossover point where SNR is maximized across the entire dynamic range was calculated to be λ = 0.48 photons per laser pulse for these PMT characteristics.

5. Summary

We have developed a method for mathematically bridging binomial counting and photon averaging by analytically recovering the mean and standard deviation of the underlying Poisson distribution for signal averaging measurements, making the results of the two methods directly comparable. The Poisson parameter can be recovered simply by dividing the signal average by the mean of the per-photon voltage output of the detector. From statistical characterization of the theoretical SNR for each approach, an analytical expression and mathematical approximation were derived to determine the crossover point for which each approach offers a SNR advantage in terms of the measured characteristics of the detector response. A flash digitizer was used to capture the peak height of every SHG signal event from a PMT at 80 MHz, which filters out many of the temporally random dark counts. The simplicity of the calculations enables rapid analysis and calibration of the large image data set. From the histograms of voltages stored for every laser pulse within each pixel, photon counting and signal averaging could be performed post-acquisition, with selection of the approach yielding the highest SNR assessed on a per-pixel basis and preserving optimized quantitation across the full intrinsic dynamic range of the detector, in this case spanning nearly 5 decades in intensity.