High Flux Passive Imaging with Single-Photon Sensors

TL;DR

This paper introduces passive free-running SPAD imaging, enabling high dynamic range 2D intensity images under ambient light by estimating scene brightness from photon timing, surpassing conventional sensors by over 100 times.

Contribution

The paper presents a novel PF-SPAD imaging method that uses photon timing to estimate scene brightness, achieving unprecedented dynamic range without active illumination.

Findings

Achieves over 2 orders of magnitude improvement in dynamic range.

Can image scenes with a 1,000,000:1 brightness ratio.

Provides a theoretical model and brightness estimator for PF-SPADs.

Abstract

Single-photon avalanche diodes (SPADs) are an emerging technology with a unique capability of capturing individual photons with high timing precision. SPADs are being used in several active imaging systems (e.g., fluorescence lifetime microscopy and LiDAR), albeit mostly limited to low photon flux settings. We propose passive free-running SPAD (PF-SPAD) imaging, an imaging modality that uses SPADs for capturing 2D intensity images with unprecedented dynamic range under ambient lighting, without any active light source. Our key observation is that the precise inter-photon timing measured by a SPAD can be used for estimating scene brightness under ambient lighting conditions, even for very bright scenes. We develop a theoretical model for PF-SPAD imaging, and derive a scene brightness estimator based on the average time of darkness between successive photons detected by a PF-SPAD pixel. Our key insight is that due to the stochastic nature of photon arrivals, this estimator does not suffer from a hard saturation limit. Coupled with high sensitivity at low flux, this enables a PF-SPAD pixel to measure a wide range of scene brightness, from very low to very high, thereby achieving extreme dynamic range. We demonstrate an improvement of over 2 orders of magnitude over conventional sensors by imaging scenes spanning a dynamic range of 1,000,000:1.

Equations

$$\hat{\Phi }=\frac{1}{q\left(\bar{X}-{\tau }_d\right)},$$

$$\underset{\text{PF-SPAD Flux Estimator}}{\underbrace{\boxed{\hat{\Phi }=\frac{N_T}{q\left(T-N_T{\tau }_d\right)}}}}.$$

$$\mathbf{E}[N_T]=\frac{q\Phi T}{1+q\Phi {\tau }_d}.$$

$$\mathrm{Var}[N_T]=\frac{q\Phi T}{(1+q\Phi {\tau }_d{)}^3}.$$

$$\mathsf{RMSE}(\hat{\Phi })=\sqrt{\mathbf{E}[(\hat{\Phi }-\Phi {)}^2]},$$

$$\mathsf{RMSE}(\hat{\Phi })=\sqrt{({\Phi }_{\mathrm{dark}}+B_{\mathrm{ap}}{)}^2+V_{\text{shot}}+V_{\mathrm{quantization}}}.$$

$$V_{\mathrm{shot}}=\frac{\Phi (1+q\Phi {\tau }_d)}{qT}.$$

$$V_{\mathrm{quantization}}=\frac{(1+q\Phi {\tau }_d{)}^4}{12q^2{T}^2}.$$

$$B_{\mathrm{ap}}=p_{\mathrm{ap}}q\Phi (1+\Phi {\tau }_d)e^{-q\Phi {\tau }_d}.$$

$$\mathsf{SNR}(\Phi )=20{\log }_{10}\left(\frac{\Phi }{\mathsf{RMSE}(\hat{\Phi })}\right).$$

$$\boxed{\mathsf{SNR}(\Phi )=-10{\log }_{10}\left[{\left(\frac{{\Phi }_{\mathrm{dark}}}{\Phi }+q(1+\Phi {\tau }_d)p_{\mathrm{ap}}{e}^{-q\Phi {\tau }_d}\right)}^2+\frac{(1+q\Phi {\tau }_d)}{q\Phi T}+\frac{(1+q\Phi {\tau }_d{)}^4}{12q^2{\Phi }^2{T}^2}\right]}.$$

$$X_n\stackrel{iid}{\sim }{f}_{X_n}(t)=\{\begin{array}{ll}q\Phi e^{-q\Phi (t-{\tau }_d)}& \text{mbox}fort\ge {\tau }_d\\ 0& \text{mbox}otherwise.\end{array}$$

$$\log l(q\Phi ;X_1,\dots ,X_{N_T})$$

$$=-q\Phi \left(\sum _{n=1}^{N_T}{X}_n-{\tau }_d{N}_T\right)+N_T\log q\Phi$$

$$=-q\Phi N_T\left(\bar{X}-{\tau }_d\right)+N_T\log q\Phi$$

$$\frac{N_T}{q\hat{\Phi }}-N_T(\bar{X}-{\tau }_d)=0$$

$$\hat{\Phi }=\frac{1}{q}\frac{1}{\bar{X}-{\tau }_d}.$$

$$\mathbf{E}[N_T]=\frac{q\Phi T}{1+q\Phi {\tau }_d}$$

$$\mathrm{Var}[N_T]=\frac{q\Phi T}{(1+q\Phi {\tau }_d{)}^3}.$$

$$F_{\hat{\Phi }}(x)$$

$$=\mathsf{Pr}\left(\frac{1}{q}\frac{N_T}{T-{\tau }_d{N}_T}\le x\right)$$

$$\approx \mathsf{Pr}\left(N_T\le \frac{qxT}{1+qx{\tau }_d}\right)$$

$$=\frac{1}{2}\left(1+\mathrm{erf}\left(\frac{\frac{qxT}{1+qx{\tau }_d}-\frac{q\Phi T}{1+q\Phi {\tau }_d}}{\sqrt{2}\sqrt{\frac{q\Phi T}{(1+q\Phi {\tau }_d{)}^3}+\frac{1}{12}}}\right)\right)$$

$$\approx \frac{1}{2}\left(1+\mathrm{erf}\left(\frac{x-\Phi }{\sqrt{2}\sqrt{\frac{\Phi (1+q\Phi {\tau }_d)}{qT}+\frac{(1+q\Phi {\tau }_d{)}^4}{12q^2{T}^2}}}\right)\right)$$

$$\Delta \hat{\Phi }=\hat{\Phi }\frac{T}{T-N_T{\tau }_d}\frac{\Delta N_T}{N_T}=\frac{T}{T-N_T{\tau }_d}{p}_{\mathrm{ap}}{e}^{-q\Phi {\tau }_d}=q\Phi (1+\Phi {\tau }_d)p_{\mathrm{ap}}{e}^{-q\Phi {\tau }_d}.$$

$$\mathsf{RMSE}(\hat{\Phi })=\sqrt{({\Phi }_{\mathrm{dark}}+q\Phi (1+\Phi {\tau }_d)p_{\mathrm{ap}}{e}^{-q\Phi {\tau }_d}{)}^2+\frac{\Phi (1+q\Phi {\tau }_d)}{qT}+\frac{(1+q\Phi {\tau }_d{)}^4}{12q^2{T}^2}}$$

$$p_n:=\mathsf{Pr}(N_T=n)=F_{S_n}(T)-F_{S_{n+1}}(T)$$

$$F_{S_n}(T)=1-\sum _{k=0}^{n-1}\frac{(T-n{\tau }_d{)}^k(q\Phi {)}^k}{k!}{e}^{-(T-n{\tau }_d)q\Phi }=1-Q(n-1,(T-n{\tau }_d)q\Phi )$$

$$Q_{q,\Phi ,T,{\tau }_d}(k):=Q(k,(T-k{\tau }_d)q\Phi ).$$

$$p_n=\{\begin{array}{ll}{Q}_{q,\Phi ,T,{\tau }_d}(n)-Q_{q,\Phi ,T,{\tau }_d}(n-1)& \text{mbox}for1\le n\le \lfloor \frac{T}{{\tau }_d}\rfloor \\ 0& \text{mbox}otherwise.\end{array}$$

$$\mathsf{RMSE}(\hat{\Phi })=\sqrt{({\Phi }_{\mathrm{dark}}+q\Phi (1+\Phi {\tau }_d)p_{\mathrm{ap}}{e}^{-q\Phi {\tau }_d}{)}^2+\sum _{n=1}^{\lfloor \frac{T}{{\tau }_d}\rfloor }{p}_n{\left(\frac{1}{q}\frac{n}{T-n{\tau }_d}-\Phi \right)}^2}$$

$${\hat{\Phi }}_{\text{CCD}}=\{\begin{array}{ll}\frac{N_T}{qT},& N_T<N_{\mathrm{FWC}}\\ \infty ,& N_T=N_{\mathrm{FWC}}.\end{array}$$

$$\mathsf{RMSE}({\hat{\Phi }}_{\text{CCD}})=\sqrt{\mathbf{E}[({\hat{\Phi }}_{\text{CCD}}-\Phi {)}^2]}=\{\begin{array}{ll}\frac{\sqrt{q\Phi T+{\sigma }_r^2}}{qT},& \Phi <\frac{N_{\mathrm{FWC}}}{qT}\\ \infty ,& \Phi \ge \frac{N_{\mathrm{FWC}}}{qT}.\end{array}$$

Full text

\dagger$$\daggerfootnotetext: Equal contribution.

This research was supported in part by ONR grants N00014-15-1-2652 and N00014-16-1-2995 and DARPA grant HR0011-16-C-0025.

1 Introduction

Single-photon avalanche diodes (SPADs) can count individual photons and capture their temporal arrival statistics with very high precision [1]. Due to this capability, SPADs are widely used in low light scenarios [2, 3, 4], LiDAR [5, 6] and non-line of sight imaging [7, 8, 9]. In these applications, SPADs are used in synchronization with an active light source (e.g., a pulsed laser). In this paper, we propose passive free-running SPAD (PF-SPAD) imaging, where SPADs are used in a free-running mode, with the goal of capturing 2D intensity images of scenes under passive lighting, without an actively controlled light source. Although SPADs have so far been limited to low flux settings, using the timing statistics of photon arrivals, PF-SPAD imaging can successfully capture much higher flux levels than previously thought possible.

We build a detailed theoretical model and derive a scene brightness estimator for PF-SPAD imaging that, unlike a conventional sensor pixel, does not suffer from full well capacity limits [10] and can measure high incident flux. Therefore, a PF-SPAD remains sensitive to incident light throughout the exposure time, even under very strong incident flux. This enables imaging scenes with large brightness variations, from extreme dark to very bright. Imagine an autonomous car driving out of a dark tunnel on a bright sunny day, or a robot inspecting critical machine parts made of metal with strong specular reflections. These scenarios require handling large illumination changes, that are often beyond the capabilities of conventional sensors.

Intriguing Characteristics of PF-SPAD Imaging: Unlike conventional sensor pixels that have a linear input-output response (except past saturation), a PF-SPAD pixel has a non-linear response curve with an asymptotic saturation limit as illustrated in Figure 1. After each photon detection event, the SPAD enters a fixed dead time interval where it cannot detect additional photons. The non-linear response is a consequence of the PF-SPAD adaptively missing a fraction of the incident photons as the incident flux increases (see Figure 1 top-right). Theoretically, a PF-SPAD sensor does not saturate even at extremely high brightness values. Instead, it reaches a soft saturation limit beyond which it still stays sensitive, albeit with a lower signal-to-noise ratio (SNR). This soft saturation point is reached considerably past the saturation limits of conventional sensors, thus, enabling PF-SPADs to reliably measure high flux values.

Various noise sources in PF-SPAD imaging also exhibit counter-intuitive behavior. For example, while in conventional imaging, photon noise increases monotonically (as square-root) with the incident flux, in PF-SPAD imaging, the photon noise first increases with incident flux, and then decreases after reaching a maximum value, until eventually, it becomes even lower than the quantization noise. Quantization noise dominates at very high flux levels. In contrast, for conventional sensors, quantization noise affects SNR only at very low flux; and when operating in realistic flux levels, photon noise dominates other sources of noise.

Extreme Dynamic Range Imaging with PF-SPADs: Due to their ability to measure high flux levels, combined with single-photon sensitivity, PF-SPADs can simultaneously capture a large range of brightness values in a single exposure, making them well suited as high dynamic range (HDR) imaging sensors. We provide theoretical justification for the HDR capability of PF-SPAD imaging by modeling its photon detection statistics. We build a hardware prototype and demonstrate single-exposure imaging of scenes with an extreme dynamic range of $10^{6}\!\!:\!\!1$, over 2 orders of magnitude higher than conventional sensors. We envision that the proposed approach and analysis will expand the applicability of SPADs as general-purpose, all-lighting-condition, passive imaging sensors, not limited to specialized applications involving low flux conditions or active illumination, and play a key role in applications that witness extreme variations in flux levels, including astronomy, microscopy, photography, and computer vision systems.

Scope and Limitations: The goal of this paper is to present the concept of adaptive temporal binning for passive flux sensing and related theoretical analysis using a single-pixel PF-SPAD implementation. Current SPAD technology is still in a nascent stage, not mature enough to replace conventional CCD and CMOS image sensors. Megapixel PF-SPAD arrays have not been realized yet. Various technical design challenges that must be resolved to enable high resolution PF-SPAD arrays are beyond the scope of this paper.

2 Related Work

HDR Imaging using Conventional Sensors: The key idea behind HDR imaging with digital CMOS or CCD sensors is similar to combination printing [11] — capture more light from darker parts of the scene to mitigate sensor noise and less light from brighter parts of the scene to avoid saturation. A widely used computational method called exposure bracketing [12, 13] captures multiple images of the scene using different exposure times and blends the pixel values to generate an HDR image.Exposure bracketing algorithms can be adapted to the PF-SPAD image formation model to further increase their dynamic range.

Hardware Modifications to Conventional Sensors: Spatially varying exposure technique modulates the amount of light reaching the sensor pixels using fixed [14] or adaptive [15] light absorbing neutral density filters. Another method [16] involves the use of beam-splitters to relay the scene onto multiple imaging sensors with different exposure settings. In contrast, our method can provide improved dynamic range without having to trade off spatial resolution.

Sensors with Non-Linear Response: Logarithmic image sensors [17] use additional hardware in each pixel that applies logarithmic non-linearity to obtain dynamic range compression. Quanta image sensors (QIS) obtain logarithmic dynamic range compression by exploiting fine-grained (sub-diffraction-limit) spatial statistics, through spatial oversampling [18, 19, 20]. We take a different approach of treating a SPAD as an adaptive temporal binary sensor which subdivides the total exposure time into random non-equispaced time bins at least as long as the dead time of the SPAD. Experimental results in recent work [21] have shown the potential of this method for improved dynamic range over the QIS approach. Here we provide a comprehensive theoretical justification by deriving the SNR from first principles and also show simulated and experimental imaging results demonstrating dynamic range improvements of over two orders of magnitude.

3 Passive Imaging with a Free-Running SPAD

In this section we present an image formation model for a PF-SPAD and derive a photon flux estimator that relies on inter-photon detection times and photon counts. This provides formal justification for the notion of adaptive photon rejection and the asymptotic response curve of a PF-SPAD.

Each PF-SPAD pixel passively measures the photon flux from a scene point by detecting incident photons over a fixed exposure time. The time intervals between consecutive incident photons vary randomly according to a Poisson process [22]. If the difference in the arrival times of two consecutive photons is less than the SPAD dead time, the later photon is not detected. The free-running operating mode means that the PF-SPAD pixel is ready to capture the next available photon as soon as the dead time interval from the previous photon detection event elapses111In contrast, conventionally, SPADs are triggered at fixed intervals, for example, synchronized with a laser pulse in a LiDAR application, and the SPAD detects at most one photon for each laser pulse.. In this free-running, passive-capture mode the PF-SPAD pixel acts as a temporal binary sensor that divides the total exposure time into random, non-uniformly spaced time intervals, each at least as long as the dead time. As shown in Figure 1, the PF-SPAD pixel detects at most one photon within each interval; additional incident photons during the dead time interval are not detected. The same figure also shows that as the average number of photons incident on a SPAD increases, the fraction of the number of detected photons decreases.

PF-SPAD Image Formation Model: Suppose the PF-SPAD pixel is exposed to a constant photon flux of $\Phi$ photons per unit time over a fixed exposure time $T$. Let $N_{T}$ denote the total number of photons detected in time $T$, and $\{X_{1},X_{2},\ldots,X_{N_{T}-1}\}$ denote the inter-detection time intervals. We define the average time of darkness as $\bar{X}=\frac{1}{N_{T}-1}\sum_{i=1}^{N_{T}-1}{X_{i}}$. Intuitively, a larger incident flux should correspond to a lower average time of darkness, and vice versa. Based on this intuition, we derive the following estimator of the incident flux as a function of $\bar{X}$ (see \NoHyperSupplementary Note 1\endNoHyper for derivation):

[TABLE]

where $\hat{\Phi}$ denotes the estimated photon flux, $0<q<1$ is the photon detection probability of the SPAD pixel, and $\tau_{d}$ is the dead time. Note that since $X_{i}\geq\tau_{d}\,\,\forall\,i$, the estimator in Equation (1) is positive and finite. In a practical implementation, it is often more efficient to use fast counting circuits that only provide a count of the total number of SPAD detection events in the exposure time interval, instead of storing timestamps for individual detection events. In this case, the average time of darkness can be approximated as $\bar{X}\approx T/N_{T}$. The flux estimator that uses only photon counts is given by:

[TABLE]

Interpreting the PF-SPAD Flux Estimator: The photon flux estimator in Equation (2) is a function of the number of photons detected by a dead time-limited SPAD pixel and is valid at all incident flux levels. The image formation procedure applies this inverse non-linear mapping to the photon counts from each PF-SPAD pixel to recover flux values, even for bright parts of the scene. The relationship between the estimated flux, $\hat{\Phi}$, and the number of photons detected, $N_{T}$, is non-linear, and is similar to the well-known non-paralyzable detector model used to describe certain radioactive particle detectors [23, 24].

To obtain further insight into the non-linear behaviour of a SPAD pixel in the free-running mode, it is instructive to analyze the average number of detected photons as a function of $\Phi$ for a fixed $T$. Using the theory of renewal processes [24] we can show that:

[TABLE]

This non-linear SPAD response curve is shown in Figure 1. The non-linear behavior is a consequence of the ability of a SPAD to perform adaptive photon rejection during the exposure time. The shape of the response curve is similar to a gamma-correction or tone-mapping curve used for displaying an HDR image. As a result, the SPAD response curve provides dynamic range compression, gratis, with no additional hardware modifications. The key observation about Equation (3) is that it has an asymptotic saturation limit given by $\lim_{\Phi\rightarrow\infty}\mathbf{E}[N_{T}]=T/\tau_{d}$. Therefore, in theory, the photon counts never saturate because this asymptotic limit can only be achieved with an infinitely bright light source. In practice, as we discuss in the following sections, due to the inherent quantized nature of photon counts, the estimator in Equation (2) suffers from a soft saturation phenomenon at high flux levels and limits the SNR.

4 Peculiar Noise Characteristics of PF-SPADs

In this section, we list various noise sources that affect a PF-SPAD pixel, derive mathematical expressions for the bias and variance they introduce in the total photon counts, and provide intuition on their surprising, counter-intuitive characteristics as compared to a conventional pixel. Ultimately, the flux estimation performance limits will be determined by the cumulative effect of these sources of noise as a function of the incident photon flux.

Shot Noise: For a conventional image sensor, due to Poisson distribution of photon arrivals, the variance of shot noise is proportional to the incident photon flux [22], as shown in Figure 2. A PF-SPAD, however, adaptively rejects a fraction of the incident photons during the dead time. Therefore, although the incident photons follow Poisson statistics, the photon counts (number of detected photons) do not. We define shot noise for PF-SPADs as the variance in the detected number of photon counts. This is approximately given by (see \NoHyperSupplementary Note 2\endNoHyper):

[TABLE]

As shown in Figure 2, the variance first increases as a function of incident flux, reaches a maximum and then decreases at very high flux levels. This peculiar behavior can be understood intuitively from the PF-SPAD photon detection timelines in Figure 1 and observing how the dead time intervals are spread within the exposure time. At low flux, when $\Phi\ll 1/\tau_{d}$, the dead time windows, on average, have large intervening time gaps. So the detected photon count statistics behave approximately like a conventional image sensor with Poisson statistics: $\mathrm{Var}[N_{T}]\approx q\Phi T$. This explains the monotonically increasing trend in variance at low flux. However, for large incident flux $\Phi\gg 1/\tau_{d}$ the time of darkness between consecutive dead time windows becomes sufficiently small that the PF-SPAD detects a photon soon after the preceding dead time interval expires. This causes a decrease in randomness which manifests as a monotonically decreasing photon count variance. In theory, as $\Phi\rightarrow\infty$ the process becomes deterministic with zero variance: the PF-SPAD detects exactly one photon per dead time window.

Quantization Noise and Saturation: For a PF-SPAD, since the photon counts are always integer valued, the source of quantization noise is inherent in the measurement process. As a first order approximation, this can be modeled as being uniformly distributed in the interval $[0,1]$ which has a variance of $\nicefrac{{1}}{{12}}$ for all incident flux levels.222For exact theoretical analysis refer to \NoHyperSupplementary Note 3\endNoHyper. A surprising consequence of the monotonically decreasing behavior of PF-SPAD shot noise is that at sufficiently high photon flux, quantization noise exceeds shot noise and becomes the dominant source of noise. This shown in Figure 2 (zoomed inset). We refer to this phenomenon as soft saturation, and discuss this in more detail in the next section.

In contrast, for a conventional imaging sensor, quantization noise is often ignored at high incident flux levels because state of the art CMOS and CCD sensors have analog-to-digital conversion (ADC) with sufficient bit depths. However, these sensors suffer from full well capacity limits beyond which they can no longer detect incident photons. As shown in Figure 2, we incorporate this hard saturation limit into quantization noise by allowing the quantization variance to jump to infinity when the full well capacity is reached.

Dark Count and Afterpulsing Noise: Dark counts are spurious counts caused by thermally generated electrons and can be modeled as a Poisson process with rate $\Phi_{\mathrm{dark}}$, independent of the true photon arrivals. Afterpulsing noise refers to spurious counts caused due to charged carriers that remain trapped in the SPAD from preceding photon detections. In most modern SPAD detectors dark counts and afterpulsing effects are usually negligible and can be ignored.

Effect of Noise on Scene Brightness Estimation: Since the output of a conventional sensor pixel is linear in the incident brightness, the variance in estimated brightness is simply equal (up to a constant scaling factor) to the noise variance. This is not the case for a PF-SPAD pixel due to its non-linear response curve — the variance in photon counts due to different sources of noise must be converted to a variance in brightness estimates, by accounting for the non-linear dependence of $\hat{\Phi}$ on $N_{T}$ in Equation (2). This raises a natural question: Given the various noise sources that affect the photon counts obtained from a PF-SPAD pixel, how reliable is the estimated scene brightness?

5 Extreme Dynamic Range of PF-SPADs

The various sources of noise in a PF-SPAD pixel described in the previous section cause the estimated photon flux $\hat{\Phi}$ to deviate from the true value $\Phi$. In this section we derive mathematical expressions for the bias and variance introduced by these different sources of noise in the PF-SPAD flux estimate. The cumulative effect of these errors is captured in the root-mean-squared error (RMSE) metric:

[TABLE]

where the expectation operation averages over all the sources of noise in the SPAD pixel. Using the bias-variance decomposition, the RMSE of the PF-SPAD flux estimator can be decomposed as a sum of flux estimation errors from the different sources of noise:

[TABLE]

The variance in the estimated flux due to shot noise (Equation (4)) is given by:

[TABLE]

The variance in estimated flux due to quantization is:

[TABLE]

The dark count bias $\Phi_{\mathrm{dark}}$ depends on the operating temperature. Finally, the afterpulsing bias $B_{\mathrm{ap}}$ can be expressed in terms of the afterpulsing probability $p_{\mathrm{ap}}$:

[TABLE]

See \NoHyperSupplementary Note 2\endNoHyper and \NoHyperSupplementary Note 3\endNoHyper for detailed derivations of Equations (6–8).

Figure 3(a) shows the flux estimation errors introduced by the various noise sources as a function of the incident flux levels for a conventional and a PF-SPAD pixel.333The effects of dark counts and afterpulsing noise are usually negligible and are discussed in \NoHyperSupplementary Note 4\endNoHyper and shown in Supplementary Figure 1. The performance of the PF-SPAD flux estimator can be expressed in terms of its SNR, formally defined as the ratio of the true photon flux to the RMSE of the estimated flux [18]:

[TABLE]

By substituting the expressions for various noise sources from Equations (5-7) into Equation (9), we get an expression for the SNR of the SPAD-based flux estimator shown in Equation (10). Figure 3(b) shows the theoretical SNR as a function of incident flux for the PF-SPAD flux estimator, and a conventional sensor. A conventional sensor suffers from an abrupt drop in SNR due to hard saturation (see \NoHyperSupplementary Note 5\endNoHyper). In contrast, the SNR achieved by a SPAD sensor degrades gracefully, even beyond the soft saturation point.

The Soft Saturation Phenomenon: It is particularly instructive to observe the behavior of quantization noise for the SPAD pixel. Although the quantization noise in the detected photon counts remains small and constant at all flux levels, the variance in the estimated flux due to quantization increases monotonically with incident flux. This is due to the non-linear nature of the estimator in Equation (2). At high incident flux levels, a single additional detected photon maps to a large range of estimated flux values, resulting in large errors in estimated flux. We call this phenomenon soft saturation. Beyond the soft saturation flux level, quantization dominates all other noise sources, including shot noise. The soft saturation limit, however, is reached at considerably higher flux levels as compared to the hard saturation limit of conventional sensors, thus, enabling PF-SPADs to reliably estimate very high flux levels.

Effect of Varying Exposure Time: For conventional imaging sensors, increasing the exposure time causes the sensor pixel to saturate at a lower value of the incident flux level. This is equivalent to a horizontal translation of the conventional sensor’s SNR curve in Figure 3(b). This does not affect its dynamic range. However, for a PF-SPAD pixel, the asymptotic saturation limit increases linearly with the exposure time, hence increasing the SNR at all flux levels. This leads to a remarkable behavior of increasing the dynamic range of a PF-SPAD pixel with increasing exposure time. See \NoHyperSupplementary Note 6\endNoHyper and Supplementary Figure \NoHyper2\endNoHyper.

Simulated Megapixel PF-SPAD Imaging System: Figure 1 (bottom row) shows simulated images for a conventional megapixel image sensor array and a hypothetical megapixel PF-SPAD array. The ground truth photon flux image was obtained from an exposure bracketed HDR image captured using a Canon EOS Rebel T5 DSLR camera with 10 stops rescaled to cover a dynamic range of $10^{6}:1.$ An exposure time of $T=$5\text{,}\mathrm{ms} was used to simulate both images. For fair comparison, the SPAD dead time was set to $149.7\text{\,}\mathrm{ns}$, which corresponds to an asymptotic saturation limit of $\nicefrac{{T}}{{\tau_{d}}}=$34\,000, equal to the conventional sensor full well capacity. The quantum efficiencies of the conventional sensor and PF-SPAD were set to 90% and 40%. Observe that the PF-SPAD can simultaneously capture details in the dark regions of the scene (e.g. the text in the shadow) and bright regions in the sun-lit sky. The conventional sensor array exhibits saturation artifacts in the bright regions of the scene. (See \NoHyperSupplementary Note 7\endNoHyper).

The human eye has a unique ability to adapt to a wide range of brightness levels ranging from a bright sunny day down to single photon levels [25, 26]. Conventional sensors cannot simultaneously reliably capture very dark and very bright regions in many natural scenes. In contrast, a PF-SPAD can simultaneously image dark and bright regions of the scene in a single exposure. Additional simulation results are shown in Supplementary Figures \NoHyper7–9\endNoHyper.

6 Experimental Results

SNR and Dynamic Range of a Single-Pixel PF-SPAD:

Figure 3(c) shows experimental SNR measurements using our prototype single-pixel SPAD sensor together with the SNR predicted by our theoretical model. Our hardware prototype has an additional $6\text{\,}\mathrm{ns}$ jitter introduced by the digital electronics that control the dead time window duration. This is not included in the SNR curve of Figure 3(b) but is accounted for in the theoretical SNR curve shown in Figure 3(c). See \NoHyperSupplementary Note 8\endNoHyper for details. We define dynamic range as the ratio of largest to smallest photon flux values that can be measured above a specified minimum SNR. Assuming a minimum acceptable SNR of 30 dB, the SPAD pixel achieves a dynamic range improvement of over 2 orders of magnitude compared to a conventional sensor.

Point-Scanning Setup: The imaging setup shown in Figure 4 consists of a SPAD module mounted on a pair of micro-translation stages (VT-21L Micronix USA) to raster-scan the image plane of a variable focal length lens (Fujifilm DV3.4x3.8SA-1). Photon counts were recorded using a single-photon counting module (PicoQuant HydraHarp 400), with the SPAD in the free-running mode. A monochrome machine vision camera (FLIR GS3-U3-23S6M-C) was used for qualitative comparisons with the images acquired using the SPAD setup. The machine vision camera uses the same variable focal length lens with identical field of view as the scene imaged by the SPAD point-scanning setup. This ensures a comparable effective incident flux on a per-pixel basis for both the SPAD and the machine vision camera. The sensor pixel parameters are identical to those used in simulations. Images captured with the machine vision camera were downsampled to match the resolution of the raster-scanned PF-SPAD images.

Extreme HDR: Results of single-shot HDR images from our raster-scanning PF-SPAD prototype are shown in Figure 5 and Supplementary Figure 10, for different scenes spanning a wide dynamic range ($\geq 10^{6}:1$) of brightness values. To reliably visualize the wide range of brightnesses in these scenes, three different tone-mapping algorithms were used to tone-map the main figures, the dark zoomed insets and the bright zoomed insets, respectively. The machine vision camera fails to capture bright text outside the tunnel (Fig. 5(a)) and dark text in the tunnel (Fig. 5(b)) in a single exposure interval. The PF-SPAD successfully captures the entire dynamic range (Fig. 5(c)). In Fig. 5(f), the PF-SPAD even captures the bright filament of an incandescent bulb simultaneously with dark text in the shadow. The halo artifacts in Figure 5(d–f) are due to a local adaptation-based non-invertible tone-map that was used to simultaneously visualize the bright filament and the dark text. This ability of the PF-SPAD flux estimator to capture a wide range of flux from very low to high in a single capture can have implications in many applications [27, 28, 29, 30]. that require extreme dynamic range.

7 Discussion

Quanta Image Sensor: An alternative realization [31] of a SPAD-based imaging sensor divides the total exposure time $T$ into uniformly spaced intervals of duration $\tau_{b}\geq\tau_{d}$. This “uniform-binning” method leads to a different image formation model which is known in literature as the oversampled binary image sensor [18] or quanta image sensor (QIS) [19, 20]. In \NoHyperSupplementary Note 9\endNoHyper, we show that in theory, this uniform-binning implementation has a smaller dynamic range as compared to a PF-SPAD that allows the dead time windows to shift adaptively [21]. Note, however, that state of the art QIS technology provides much higher resolution and fill factor with high quantum efficiencies, and lower read noise than current SPAD arrays.

Limitations and Future Outlook: Our proof-of-concept imaging system uses a SPAD that is not optimized for operating in the free-running mode. The duration of the dead time window, which is a crucial parameter in our flux estimator, is not stable in current SPAD implementations (such as silicon photo-multipliers) as it is not crucial for active time-of-flight applications. Various research and engineering challenges must be met to realize a high resolution SPAD-based passive image sensor. State of the art SPAD pixel arrays that are commercially available today consist of thousands of pixels with row or column multiplexed readout capabilities and do not support fully parallel readout. Current SPAD arrays also have very low fill factors due to the need of integrating counting and storage electronics within each pixel [32, 33]. Our method and results make a case for developing high resolution fabrication and 3D stacking techniques that will enable high fill-factor SPAD arrays, which can be used as general purpose, passive sensors for applications requiring extreme dynamic range imaging.

Supplementary Note 1 Image Formation Model and Flux Estimator for a PF-SPAD

Pixel

A PF-SPAD sensor pixel and a time-correlated photon counting module are used to obtain total photon counts over a fixed exposure time together with picosecond resolution measurements of the time elapsed between successive photon detection events. We will assume that the PF-SPAD pixel is exposed to a true photon flux of $\Phi$ photons/second for an exposure time of $T$ seconds and it records $N_{T}$ photons in exposure time interval $(0,T]$. For mathematical convenience, we assume that the exposure interval starts with a photon detection event at time $t=0$.

Photons arrive at the SPAD according to a Poisson process. Accounting for an imperfect photon detection efficiency of $0<q<1$, the time between consecutive incident photons follows an exponential distribution with rate $q\Phi.$ After each detection event, the SPAD enters a dead time window of duration $\tau_{d}$. Due to the memoryless property of Poisson processes [6], the time interval between the end of a dead time window and the next photon arrival is also exponentially distributed and has the same rate $q\Phi$ as the incident Poisson process. Let $X_{1}$ be the time of the first photon detection after $t=0$ and $X_{n}$ be the time between the $n-1^{st}$ and $n^{th}$ detection event for $n\geq 2$. Then the inter-detection time duration $X_{n}$ follows a shifted exponential distribution given by:

[TABLE]

This provides a probabilistic model of the photon inter-detection times. We now derive a flux estimator from a sequence of observed inter-detection times captured by a PF-SPAD pixel.

Estimating Flux from Inter-Detection Time

Intervals

The log-likelihood function for the observed inter-detection times is given by

[TABLE]

where $\bar{X}:=\frac{1}{N_{T}}\sum_{n=1}^{N_{T}}X_{n}$ is the mean time between photon detection events. The maximum likelihood estimate $\hat{\Phi}$ of the true photon flux is computed by setting the derivative of Equation (S2) to zero:

[TABLE]

which implies

[TABLE]

Supplementary Note 2 Approximate Closed Form Formula for SNR of a PF-SPAD

pixel

We first derive an approximate formula for the SNR of a PF-SPAD pixel using a continuous Gaussian distribution approximation for the number of counts $N_{T}$. The effective incident photon flux for a quantum efficiency of $0<q<1$ is equal to $q\Phi$ photons/second.

The random process describing the detections of this PF-SPAD pixel is not a Poisson process, but can be modeled as a renewal process [6] with a shifted exponential inter-arrival distribution which has a mean $\tau_{d}+\frac{1}{q\Phi}$ and variance $\frac{1}{q^{2}\Phi^{2}}$. Using the central limit theorem for renewal processes, $N_{T}$ is approximately Gaussian distributed with mean:

[TABLE]

and variance:

[TABLE]

Quantization Noise: An additional source of variance arises due to quantization noise which we can treat as uniformly distributed between [math] and $1$ with variance $1/12$. The c.d.f. of the estimated photon flux $\hat{\Phi}$ can be computed using the delta method [3]:

[TABLE]

where Equation (S6) follows from the fact that in practice the denominator is always non-negative since $N_{T}\leq\lfloor T/\tau_{d}\rfloor$, Equation (S7) follows from the formula for the Gaussian c.d.f. of $N_{T}$ with $\mathrm{erf}$ denoting the error function [1], and Equation (S8) is follows from a first order Taylor series approximation. This result shows that $\hat{\Phi}$ is approximately normally distributed with mean equal to the true photon flux and variance given by the denominator in Equation (S8).

Dark Count and Afterpulsing Bias: In addition to quantization and shot noise that introduce variance in the estimated photon flux, PF-SPADs also suffer from dark counts and afterpulsing noise that introduce a bias in the estimated flux. The dark count rate $\Phi_{\mathrm{dark}}$ is often given in published datasheets and can be used as the bias term. Afterpulsing noise is quoted in datasheets as afterpulsing probability which denotes the probability of observing a spurious afterpulse after the dead time $\tau_{d}$ has elapsed. Due to an exponentially distributed waiting time, the probability of observing a gap between true photon-induced avalanches is equal to $e^{-q\Phi\tau_{d}}$. A fraction $p_{\mathrm{ap}}e^{-q\Phi\tau_{d}}$ of these gaps will contain afterpluses, on average. The bias $\Delta\hat{\Phi}$ in the estimated flux is given by:

[TABLE]

Using the bias-variance decomposition of mean-squared error, we have

[TABLE]

and the approximate closed from SNR is obtained by plugging Equation (S10) into Equation (\NoHyper9\endNoHyper) in the main text.

Supplementary Note 3 Exact Formula for Numerical Computation of SNR of a SPAD

Pixel

It is possible to model the exact discrete distribution of the number of counts $N_{T}$ for a PF-SPAD pixel using non-asymptotic renewal theory. The times between consecutive counts for a PF-SPAD pixel can be modeled as a shifted exponential distribution as before. Let $X_{n}$ be the time between when the SPAD detects the $(n-1)^{st}$ and $n^{th}$ photons ($n\geq 1$). For mathematical convenience, we assume $X_{0}=0.$ Let $F_{S_{n}}$ be the c.d.f. of the sum $S_{n}:=\sum_{i=1}^{n}X_{n}.$ Then, by definition, $F_{S_{n}}(T)=\mathsf{Pr}(N_{T}\geq n).$ Therefore we can write

[TABLE]

where

[TABLE]

and $Q(\cdot,\mu)$ is the c.d.f. of a Poisson random variable with rate $\mu,$ also known as the regularized gamma function [1]. For convenience, define:

[TABLE]

The following formula can now be used to numerically compute the probability mass function of $N_{T}$:

[TABLE]

Using the bias-variance decomposition, the RMSE can be written as:

[TABLE]

and SNR can be computed by plugging Equation (S11) and Equation (S12) in Equation (\NoHyper9\endNoHyper) in the main text. Note that although this formula is exact, it does not lend itself to an intuitive interpretation as the approximate formula in Equation (S10) which decomposes the sources of variance into shot noise and quantization noise.

Supplementary Note 4 Various Sources of Noise Affecting the PF-SPAD Flux Estimate

Various unique properties of the shot noise and quantization noise were discussed in the main text. Another surprising result is that the effect of afterpulsing bias first increases and then decreases with incident photon flux. Recall that afterpulses are correlated with past avalanche events. At very low incident photon flux there are very few photon-induced avalanches which implies that there are even fewer afterpulsing avalanches. At very high photon flux values, the afterpulsing noise is overwhelmed by the large number of true photon-induced avalanches that leave negligible temporal gaps between consecutive dead time windows. However, for most modern SPAD pixels, afterpulsing noise is so small that it can often be ignored. The plot in Supplementary Figure 1 shows an afterpulsing error curve using an unrealistically high afterpulsing rate to accentuate the trend as a function of incident flux.

Supplementary Note 5 SNR of a Conventional Sensor Pixel

A conventional CCD or CMOS pixel suffers from a hard saturation limit due to its full well capacity, $N_{\mathrm{FWC}}$. Assuming a quantum efficiency of $0<q<1$, an incident photon flux of $\Phi$ photons/second and an exposure time $T$ seconds, the photon counts $N_{T}$ follow a Gaussian distribution with mean $q\Phi T$ and variance $q\Phi T+\sigma_{r}^{2}$ where $\sigma_{r}$ is the read noise of the pixel. The estimated flux is given by [7]:

[TABLE]

The RMSE of the estimated flux is given by:

[TABLE]

which leads to the following formula for SNR of conventional pixel:

[TABLE]

This formula does not account for dark current noise because it is only relevent at extremely low incident photon flux values with very long exposure times of many minutes or longer.

Supplementary Note 6 Effect of Varying Exposure Time

The notions of quantum efficiency and exposure time are interchangeable in case of conventional image sensors; Equation (S13) remains unchanged if the symbols $q$ and $T$ were to be swapped. This is not true for a PF-SPAD sensor where changing $q$ and changing $T$ has different effects on the overall SNR. This is because the SPAD pixel has an asymptotic saturation limit of $T/\tau_{d}$ counts which is a function of exposure time, unlike a conventional sensor whose full well capacity is a fixed constant independent of exposure time. As shown in Supplementary Figure 2, decreasing exposure time decreases the maximum achievable SNR value of a SPAD sensor. Experimental results were obtained from our hardware prototype using a dead time of $300\text{\,}\mathrm{ns}$ and capturing photon counts with two different exposure times of $0.5\text{\,}\mathrm{ms}$ and $$5\text{,}\mathrm{ms}$.$

Supplementary Note 7 Details of SPAD Simulation Model and Experimental Setup

We implemented a time-domain simulation model for a PF-SPAD pixel to validate our theoretical formulas for the PF-SPAD response curve and SNR. Photons impinge the simulated PF-SPAD pixel according to a Poisson process; a fraction of these photons are missed due to limited quantum efficiency. The PF-SPAD pixel counts an incident photon when it arrives outside a dead time window. The simulation model also accounts for spurious detection events to dark counts and afterpulsing. The pseudo-code is shown in Supplementary Figure 3.

PF-SPAD and Conventional Sensor Specifications

Each pixel in the simulated PF-SPAD array mimics the specifications of the single-pixel hardware prototype. Each pixel in our simulated conventional sensor array uses slightly higher specifications than the one we used in our experiments. It has a full well capacity of 33,400 electrons, quantum efficiency of 90% and read noise of 5 electrons.

The single-pixel SPAD simulator was used for generating synthetic color images from a hypothetical megapixel SPAD array camera. The ground truth photon flux values were obtained from an exposure bracketed HDR image that covered over 10 orders of magnitude in dynamic range. Unlike regular digital images that use 8 bit integers for each pixel value, an HDR image is represented using floating point values that represent the true scene radiance at each pixel. These floating point values were appropriately scaled and used as the ground-truth photon flux to generate a sequence of photon arrival times following Poisson process statistics. Red, green and blue color channels were simulated independently. Results of simulated HDR images are shown in Supplementary Figures 7, 8 and 9.

Experimental Setup Details

The single-pixel SPAD from our hardware prototype has a pitch of $25\text{\,}\mathrm{\SIUnitSymbolMicro m}$, quantum efficiency of 40%, dark count rate of 100 photons/second and 1% afterpulsing probability. The dead time is programmable and was set to $149.7\text{\,}\mathrm{ns}$ and exposure time to $5\text{\,}\mathrm{ms}$. This corresponds to an asymptotic saturation limit of 33,400 photons.

Each pixel in our machine vision camera (Point Grey GS3-U3-23S6M-C) has a full well capacity of 32,513 electrons, a peak quantum efficiency of 80% and a Gaussian-distributed read noise with a standard deviation of 6.83 electrons. Note that the asymptotic saturation limit of the PF-SPAD pixel is similar to the full well capacity of this machine vision camera to enable fair comparison.

Supplementary Note 8 Effect of Dead Time Jitter

In practice the dead time window is controlled using digital timer circuits that have a limited precision dictated by the clock speed. The hardware used in our experiments has a clock speed of 167 MHz which introduces a variance of $6\text{\,}\mathrm{ns}$ in the duration of the dead time window. As a result the dead time $\tau_{d}$ can no longer be treated as a constant but must be treated as a random variable $T_{d}$ with mean $\mu_{d}$ and variance $\sigma_{d}.$ The inter-arrival distribution in Equation (S1) must be understood as a conditional distribution, conditioned on $T_{d}=\tau_{d}$. The mean and variance of the time between photon detections can be computed using the law of iterated expectation [6]:

[TABLE]

and

[TABLE]

Using similar computations as those leading to Equation (S8), we can derive a modified shot noise variance term equal to $\frac{\Phi(1+q^{2}\Phi^{2}\sigma_{d}^{2})(1+q\Phi\mu_{d})}{qT}$ that must be used in Equation (S10) to account for dead time variance. All instances of $\tau_{d}$ in Equation (S10) must be replaced by its mean value $\mu_{d}$. Supplementary Figure 5 shows theoretical SNR curves for a PF-SPAD pixel with a nominal dead time duration of $149.7\text{\,}\mathrm{ns}$. Observe that the 30 dB dynamic range degrades by almost 3 orders of magnitude when the dead time jitter increases from $0.01\text{\,}\mathrm{ns}$ to $50\text{\,}\mathrm{ns}$. For reference, our hardware prototype has a dead time jitter of $6\text{\,}\mathrm{ns}$ RMS.

Supplementary Note 9 Comparison with Quanta Image Sensors

A quanta image sensor (QIS) [4,5,8] improves dynamic range by spatially oversampling the 2D scene intensities using sub-diffraction limit sized pixels called jots. Each jot has a limited full well capacity, usually just one photo-electron. The PF-SPAD imaging modality presented in this paper is different from these methods. Instead of using a SPAD as a binary pixel [4] and relying on spatial oversampling, the PF-SPAD achieves dynamic range compression by allowing the dead time windows to shift randomly based on the most recent photon detection time and performing adaptive photon rejection. This is equivalent to the “event-driven recharge” method described in [2].

We now derive the image formation model and an expression for SNR for a QIS and other related methods that use equi-spaced time bins [8], and show that their dynamic range is lower than what can be achieved using a PF-SPAD.

QIS Image Formation and Flux Estimator

The output response of a QIS is logarithmic and mimics silver halide photographic film [5]. Each jot has a binary output, and the final image is formed by spatio-temporally combining groups of jots called a “jot-cube”. Let $\tau_{b}$ be the temporal bin width and for mathematical convenience, assume that the exposure time $T$ is an integer multiple of $\tau_{b}$, so that there are $N=T/\tau_{b}$ uniformly spaced time bins that split the total exposure duration. Suppose the jot-cube is exposed to a constant photon flux of $\Phi$ photons/second, and each jot has a quantum efficiency $0<q<1$. The number of photons received by each jot in a time interval $\tau_{b}$ follows a Poisson distribution with mean $q\Phi\tau_{b}$. Therefore the probability that the binary output of a jot is [math] is given by:

[TABLE]

and the probability that the binary output of a jot is $1$ is given by the probability of observing $1$ or more photons:

[TABLE]

Let $N_{T}$ denote the total photon counts output from a jot-cube with $N$ jots. Then $N_{T}$ follows a binomial distribution given by:

[TABLE]

for $0\leq k\leq N$. The maximum-likelihood estimate of the photon flux is given by:

[TABLE]

Our PF-SPAD flux extimator has a higher dynamic range than this uniform binning method. This can be intuitively understood by noting that in the limiting case of $\tau_{b}=\tau_{d}$ both schemes have an upper limit on photon counts given by $N_{T}\leq T/\tau_{d},$ but the QIS estimator in Equation (S14) saturates and flattens out more rapidly than the PF-SPAD estimator in Equation (\NoHyper2\endNoHyper) from the main text:

[TABLE]

SNR of a QIS

The variance of the QIS flux estimator can be computed numerically using the binomial probability mass function of $N_{T}$. For convenience, a closed form expression can be derived using a Gaussian approximation, similar to the approximation techniques used for deriving the SNR of a PF-SPAD pixel in Equation (S8). The Gaussian approximation to a binomial distribution suggests $N_{T}$ has a normal distribution with mean $N\,(1-e^{-q\Phi\tau_{b}})$ and variance $N\,e^{-q\Phi\tau_{b}}(1-e^{-q\Phi\tau_{b}})$. Next, the c.d.f. of the estimated flux is given by:

[TABLE]

where Equation (S15) follows from the formula for the c.d.f. of a Gaussian distribution and Equation (S16) is obtained after making a Taylor series approximation. The final form of Equation (S17) suggests that the estimated photon flux follows a normal distribution with variance $\frac{(1-e^{-\Phi q\tau_{b}})}{q^{2}T\tau_{b}e^{-\Phi q\tau_{b}}}.$

The read noise of each jot affects the RMSE of the QIS sensor at low incident flux. We note that at low flux values there are, on average, $q\Phi T$ bins already filled by true photon counts leaving $N-q\Phi T$ bins empty. Read noise will cause some of these empty bins to contain false positives and introduce additional noisy counts equal to $\frac{1}{2}(N-q\Phi T)\left(1+\mathrm{erf}\left(\frac{1}{2\sqrt{2}\sigma_{r}}\right)\right)$, where $\sigma_{r}$ is the read noise standard deviation. This corresponds to a bias of $\frac{1}{2}(\frac{1}{q\tau_{b}}-\Phi)\left(1-\mathrm{erf}\left(\frac{1}{2\sqrt{2}\sigma_{r}}\right)\right)$ in the estimated photon flux. Incorporating this bias term together with the variance associated with the Gaussian distribution of the estimated photon flux, the RMSE is given by

[TABLE]

A single jot only generates a binary output and must be combined into a jot-cube to generate the final image. One way to obtain a fair comparison between a PF-SPAD pixel and a jot-cube is computing the SNR for a fixed image pixel size and fixed exposure time. We use a square grid of jots that spatially occupy the same area as our single PF-SPAD pixel that has a pitch of $25\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Supplementary Figure 6 shows the SNR curves obtained using our theoretical derivations for a QIS jot-cube and a single PF-SPAD pixel. State of the art jot arrays are limited to a pixel size of around $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and frame rates of a few kHz. These SNR curves show that we will require large spatio-temporal oversampling factors and extremely small jots to obtain similar dynamic range as a single PF-SPAD pixel. For example, a $150\text{\,}\mathrm{nm}$ jot size can accommodate almost 30,000 jots in a 25\times$$25\text{\,}\mathrm{\SIUnitSymbolMicro m}2 area occupied by the PF-SPAD pixel and can provide similar dynamic range and higher SNR than our PF-SPAD pixel when operated at a frame rate of 1 kHz. QIS technology will require an order of magnitude increase in frame readout rate or an order of magnitude reduction in jot size to bring it closer to the dynamic range achievable with our PF-SPAD prototype.

Supplementary Note 10 Additional Simulated and Experimental Results

Supplementary References

1. [1]

Abramowitz, M. & Stegun, I. A. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, vol. 55 (Dover American Nurses Association Publications, 1964), 9 edn. 2. [2]

Antolovic, I. M., Bruschini, C. & Charbon, E. Dynamic range extension for photon counting arrays. Optics Express 26, 22234–22248 (2018). 3. [3]

Casella, G. & Berger, R. L. Statistical Inference (Pacific Grove, CA: Duxbury/Thomson Learning, 2002), 2nd edn. Sec. 5.5.4. 4. [4]

Dutton, N. A. W. et al. A SPAD-based QVGA image sensor for single-photon counting and quanta imaging. IEEE Transactions on Electron Devices 63, 189–196 (2016). 5. [5]

Fossum, E., Ma, J., Masoodian, S., Anzagira, L. & Zizza, R. The quanta image sensor: Every photon counts. Sensors 16, 1260 (2016). 6. [6]

Grimmett, G. R. & Stirzaker, D. R. Probability and Random Processes (Oxford University Press, 2001), 3rd edn. 7. [7]

Hasinoff, S. W. et al. Noise-Optimal Capture for High Dynamic Range Photography. Proc. 23rd IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 553–560 (2010). 8. [8]

Itzler, M. A. Apparatus comprising a high dynamic range single-photon passive 2D imager and methods therefor (2017).