Open Access
Add to BibSonomy
Abstract
Single-photon lidar (SPL) is a promising technology for depth measurement at long range or from weak reflectors because of the sensitivity to extremely low light levels. However, constraints on the timing resolution of existing arrays of single-photon avalanche diode (SPAD) detectors limit the precision of resulting depth estimates. In this work, we describe an implementation of subtractively-dithered SPL that can recover high-resolution depth estimates despite the coarse resolution of the detector. Subtractively-dithered measurement is achieved by adding programmable delays into the photon timing circuitry that introduce relative time shifts between the illumination and detection that are shorter than the time bin duration. Careful modeling of the temporal instrument response function leads to an estimator that outperforms the sample mean and results in depth estimates with up to 13 times lower root mean-squared error than if dither were not used. The simple implementation and estimation suggest that globally dithered SPAD arrays could be used for high spatial- and temporal-resolution depth sensing.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Fast and precise three-dimensional imaging is a topic of great interest to countless industries, including autonomous navigation, motion capture for video games, terrestrial mapping, and environmental monitoring, among others. Single-photon lidar (SPL) has become one of the top depth measurement modalities due to the extraordinary timing precision that can be achieved via time-correlated single photon counting (TCSPC), which time stamps individual photon detection times with picosecond resolution. Taking advantage of the sensitivity of SPL to extremely weak incident light, accurate depth imaging has been demonstrated for extremely small numbers of photon detections [1–5], over long ranges [6–8], through partially-transmissive surfaces [9,10], through fog [11,12], or even underwater [13]. Unfortunately, most of these approaches use single-pixel detectors such as single-photon avalanche diodes (SPADs) and thus rely on raster scanning of the laser illumination over the scene to provide spatial resolution. The serial process of sequentially scanning each point in the scene is limited by the mechanical scanning speed of the galvo mirrors used. Increasing scanning speed, especially by non-mechanical means, is an area of ongoing industrial development [14].
An alternative to raster scanning is the use of a detector array, which parallelizes the acquisition process. Instead of illuminating one point at a time, the use of a linear or 2D array allows for the illumination of a corresponding line or region of a scene, with spatial resolution then derived from the array of time-resolved pixels [15–19]. Combined with fast parallel processing, detector arrays enable real-time depth imaging [20], which is crucial for applications such as navigation in autonomous vehicles [21]. However, the combination of long measurement range and high-speed clocks found in single-pixel systems is difficult to replicate for detector arrays, which have competing demands that include the photon-sensitive area, timing resolution, power consumption, and memory [22]. Although some time-to-digital converter (TDC) designs enable minimum bin durations in tens of picoseconds [22–27], fixed TDC bit depths require coarsening the timing resolution to cover the full range for long-distance measurements, which can lead to the least significant bit (LSB) of the TDC having a duration of 100s of picoseconds or longer [15,28–30]. Crucially, laser illumination with pulse widths on the order of 10s of picoseconds is easily achievable, so the resolution of photon time of flight measurements is limited by the temporal quantization of the detector, not by the pulse width. Repeated illumination of a stationary scene patch—as is typical in SPL—may then cause photons to all be detected in the same time bin, so the uncertainty in a surface depth does not decrease, even as the number of photon detections increases.
Rather than modifying the array architecture to have higher bit-depth TDCs with shorter LSBs, which requires increased power consumption for faster clocks, several approaches aim to recover finer time resolution by combining adjustments to the illumination with additional signal processing. In particular, dither is a longstanding approach to reducing the effect of coarse quantization in analog-to-digital conversion by introducing a small amount of noise before quantization [31–37]. For SPL, Shin et al. [16] proposed increasing the laser pulse duration so that photons were more likely to fall in multiple time bins, a form of self-dithering that could allow depth resolution beyond the limits of the time bins via spatio-temporal averaging. In [38,39], we introduced subtractive dither in the time domain to resolve depth information smaller than the least significant bit of a TCSPC system. Subtractive dither means known noise values are added to the time-of-flight signal and then subtracted off after quantization, resulting in improved statistical properties over nonsubtractive dither. The subtractive dithering was implemented with a sequence of sub-bin delays in the TCSPC synchronization during acquisition. The generalized Gaussian (GG) distribution was proposed to model the error resulting from an approximately Gaussian instrument response function (IRF) plus uniform noise resulting from the dithering. In [40], we demonstrated that computationally simple order statistics-based estimators could outperform conventional averaging for general dithered data modeled as GG when the ratio of Gaussian pulse width to quantization bin size is small. Subtractive dither was also shown to always be at least as effective in reducing measurement error per detected photon as increasing the pulse duration. However, the GG-based estimator performed worse than the sample mean for dithered lidar in [39] because of the naïve Gaussian pulse shape assumption.
In this work, we introduce a complete subtractively-dithered lidar system and a new processing approach that can achieve high-resolution depth estimates despite coarse timing resolution. Key to the estimation algorithm’s success is an instrument response function model that better matches the physical behavior of SPAD detectors. Specifically, we use an exponentially-modified Gaussian (EMG) model that accounts for both an approximately Gaussian laser pulse shape and the exponential temporal decay common to SPADs. We modify the GG estimator to incorporate the exponential component of the pulse shape. Finally, we describe in detail a dithered lidar system implementation emulated with a single-pixel detector and present results using our new estimator for several experimental datasets. The estimator based on the EMG model outperforms the sample mean—which would be optimal under a naïve Gaussian temporal distribution assumption—and resulting depth estimates have up to 13 times lower root mean-squared error (RMSE) than if dither had not been used.
2. Related work
Similar ideas using shifted acquisition windows to acquire range information have been used in the past. In [41], Busck and Heiselberg implemented a “laser-gated viewing” (LGV) system based on early concepts from Gillespie [42], using a high-speed charge-coupled device (CCD) array with no inherent time resolution to form depth images. In their work, the scene was repeatedly illuminated by a pulsed laser, and a gate exposed the camera to light for a short duration on the order of a few hundred picoseconds. The light was integrated on the CCD over many illuminations, and once the result was stored, a short delay was added to the gate activation time, and the acquisition process was repeated with the new gate delay. By stepping through a sequence of delays, slices of the transient light response were collected from the scene, and a depth estimate was created by weighting the time of each slice by the grayscale intensity value of the image acquired for that delay. An LGV system was used by Christnacher et al. [43] to improve vision through scattering media such as smoke or fog. The same principle was recently used with a quanta image sensor (QIS) [44], a large-factor SPAD array [19], and with quantum parametric mode sorting (QPMS) [45] to also take advantage of single-photon sensitivity. Crucially, the LGV approach uses detectors with no inherent time resolution. Since light is only detected within the short gate time, all of the light arriving outside this window is discarded, which is an inefficient use of the reflected light and makes the process slow to acquire a full sequence of time slices. TCSPC is a more time- and photon-efficient approach to acquisition, as each incident photon (not arriving during a dead time) can be registered and resolved to a time bin.
Chen et al. recently proposed a method of dithered SPL using two optical paths with different lengths, such that the longer path length effectively delays laser transmission by a half-bin duration [46]. In simulations, they showed that a simple averaging of signals from different bin shifts reduced the effect of quantization on the depth estimation error. While we share the same basic idea of combining measurements from different sub-bin delays, our electronic delay implementation is far more flexible in terms of the bin sizes and numbers of delays that can be achieved in practice.
After the introduction of our work in [39], Raghuram et al. [47] used a similar electronic delay implementation and a sampling theory-based processing approach to avoid the limits of coarse temporal quantization. Our method assumes a single depth per pixel and negligible ambient light, which allows for an estimator tailored to the system IRF and the temporal quantization, enabling accurate estimates from far fewer photons than required for the deconvolution and denoising approach in [47]. In strong ambient light or when multiple depths per pixel are present, our current estimator would likely become less reliable, although robust variations on our approach could still be more photon-efficient under those conditions.
3. Photon detection model
3.1 Quantization and subtractive dither
We begin by describing the photon detection model introduced in [39] and expanded in [40]. In a standard SPL system, photon arrival times $X_i$, $i=1,\,2,\,\ldots$, are time stamped by a TDC producing measurements
where the mid-tread quantizer function with quantization step (time bin) size $\Delta$ is We assume all detected photons are due to illumination back-reflected from a single depth and ambient light is negligible, so if the temporal distribution of photon arrivals times has variance $\sigma _X \ll \Delta$, then it is likely for all $X_i$s after repeated illumination to be detected in the same bin $U$, which results in a depth estimate of the true depth limited by the quantization bin size $\Delta$. Furthermore, the estimation error is signal-dependent, depending on the true depth relative to the center of the bin.Instead, suppose the photon arrival times $X$ are measured by a quantizer that is subtractively dithered. If the dither signal $D\sim \mathcal {U}(-\Delta /2,\Delta /2)$ is uniformly distributed, then the dithered measurements
are equal in distribution to where is i.i.d. and independent of $X$ [48]. Crucially, this means that adding noise to the measurement (in the form of dither) means the error is no longer signal-dependent and has a known distribution, so we can reduce estimation error by collecting more samples and using an appropriate estimator.The distribution of the total noise can be determined as follows. Let $X$ be a random variable with probability density function (PDF) $f_X(x)$ and cumulative distribution function (CDF) $F_X(x)$. Let $W$ be a uniform random variable as defined in Eq. (5), with PDF
3.2 Inaccurate Gaussian modeling
We assume dark counts are negligible and that optical filtering, perhaps in conjunction with algorithmic methods of rejecting background detections (e.g., [4]), is sufficient to remove all detections due to ambient light, so that the photon detections are only due to the back-reflected illumination signal. If we assume the IRF yields a Gaussian distribution of arrival times $X = \mu _X+Z$, where $Z\sim \mathcal {N}(0,\sigma _Z^{2})$, then the measurements $Y_i = \mu _X+Z_i+W_i$ are the desired time-of-flight signal $\mu _X$ plus a combination of Gaussian noise from the IRF and uniform noise from the dither. We showed in [40] that such measurements are well-approximated by $Y_i = \mu _X+V_i$, where $V$ has a generalized Gaussian (GG) distribution $GG(\mu _V,\sigma _V,p_V)$, where $\mu _V = \mu _X+\mu _W$ and $\sigma _V^{2} = \sigma _Z^{2}+ \sigma _W^{2}$. The GG shape parameter $p_V$ is approximated by matching the kurtosis of the true noise distribution to the generalized Gaussian kurtosis.
Extensive simulation in [40] showed that this GG approximation is well suited to the combination of Gaussian and uniform noise, and estimators based on the order statistics of GG random variables required fewer samples than those based on Gaussian noise (the sample mean) or uniform noise (the sample midrange) alone to reach a specified mean-squared error (MSE). Both in principle and in practice [49], weighted combinations of photon detection order statistics result in estimators that outperform simple averaging, so long as the order statistics weights are correctly chosen based on the underlying detection time distribution. However, the GG-based estimators produced worse depth estimates than simply applying the sample mean when applied to actual SPL data [39]. The reason for this under-performance is that the Gaussian IRF model is not accurate for SPAD detections. In a SPAD, carriers generated as intended in the junction depletion layer are accelerated by the bias voltage, and fluctuations in the buildup of the resulting avalanche do cause a Gaussian response [50]. However, if a photon reaches the neutral region of the semiconductor before generating a carrier, diffusion of the carrier to the depletion layer causes an exponential delay in the avalanche timing [50], [51, Ch. 6]. Thus, while other single-photon detectors such as superconducting nanowire single-photon detectors (SNSPDs) tend to have instrument response functions (IRFs) well described by a Gaussian distribution [52], many SPADs additionally have a significant exponential decay [53].
3.3 Proposed exponentially-modified Gaussian response function
Since the experimental results suggest the GG-based estimator is not robust to this model mismatch, we here propose an updated pulse model that accounts for the exponential decay. The sum of Gaussian and exponential random variables, suggested by Hernandez et al. as a good approximation to the SPAD detection time distribution [53], has been long been used to describe the asymmetric timing responses in chromatography [54], for which the distribution became known as the exponentially-modified Gaussian (EMG) [55]. The same distribution was also found to be a good fit for reaction times in psychological studies, for which it was called the ex-Gaussian distribution [56], and for the IRF of a space-based lidar system [57].
Let $H$ be a random variable representing the exponential component of the SPAD response, with PDF
where $\eta _H = 1/ \tau _H$ with expected value $E[H] = \tau _H$ and variance $\sigma _H^{2} = \tau _H^{2}$. Then the PDF of the total SPAD response time $R=Z+H$ is an exponentially-modified Gaussian, with PDF given by
Fig. 1. The IRF of a lidar system is experimentally acquired by illuminating a single planar target, forming a histogram of back-reflected photon detection times and censoring background detections (see Sec. 5 for implementation details). A Gaussian fit, determined by matching the mean and variance of the empirical distribution, is a poor approximation to the true pulse shape, whereas the exponentially-modified Gaussian (EMG) distribution captures the exponential tail behavior much more accurately.
4. Dithered depth estimation
4.1 Generalized Gaussian order statistics-based estimators
One major benefit of approximating the sum of Gaussian and uniform noise with a GG random variable in [39,40] was that GG shape parameter $p_V$ could be used to parameterize simple, computationally-efficient estimators based on linear combinations of the measurement order statistics $Y_{(1)}, Y_{(2)}, \dots , Y_{(K)}$, which simply sorts $K$ samples $Y_{(1)} \leq Y_{(2)} \leq \dots \leq Y_{(K)}$ from least to greatest. The estimators weight the various order statistics based on the particular distribution as
where the $a_i$s are the order statistics coefficients.In [61], Beaulieu and Guo proposed selecting the coefficients in a data-dependent manner as
4.2 Estimator modification for the EMG
One of the main challenges in transferring the GG-based approach used in previous work [40] to experimental single-photon lidar data is that the Gaussian model of the pulse shape does not capture the asymmetry of the actual SPAD response. Symmetry has often been used to greatly simplify the form of estimators, as the mean, median, midrange, and other symmetric linear combinations of order statistics are all co-located [63]. However, it is still valuable to make estimates of location when that is not the case, although the parameter being estimated must more clearly be identified and the estimation done more carefully [64].
One approach to getting an unbiased estimator that is sensitive to the distribution is to use the best linear unbiased estimator (BLUE) of the order statistics as suggested by [65]. Unfortunately, this requires the ability to compute the expected values and covariance matrix of order statistics, which changes for every change in the number of samples $K$ and can be complicated depending on the parent distribution. A simplification suggested by [66] is to assume the covariance matrix is the identity. The location estimator then reweights the coefficients of the order statistics so that the estimate is unbiased. However, this still requires computation of the expected value of each order statistic, which for the dithered EMG requires repeated numerical integration using the substitution of Eq. (12) into Eq. (8), which is so unwieldy that we have omitted writing out the full expression. While it would be possible to precompute and tabulate reweighted coefficients for every reasonable number of detections and for fixed system parameters, we find such computation prohibitively expensive for a negligible change in the estimator, so we continue without exactly-unbiased estimators here.
Rather than choosing order statistics and their coefficients to exactly match the dithered detection distribution, we choose to ignore the effect of asymmetry and only modify the GG-based approximation to account for the longer tails incurred from the exponential component. First, we note that the dithered EMG model changes the expected value of the measurements to $E[Y] = \mu _X+ \tau _H$ and the variance to $\textrm {Var}[Y] = \sigma _Z^{2}+ \tau _H^{2}+\Delta ^{2}/12$. Furthermore, the kurtosis of the true distribution changes:
Fig. 2. A histogram of simulated dithered EMG data compared to the true distribution and generalized Gaussian approximation. The distribution parameters were chosen to resemble experimental values, with $\sigma _Z = 40$ ps, $\tau _H = 200$ ps, and $\Delta = 2000$ ps, which were then standardized for the overall distribution to have zero mean and unit variance. While the GG does not match the tail behavior exactly on either side, the value of $\widehat {p}$ determined through kurtosis matching does get the width and overall shape generally correct.
5. Implementation
The process of TCSPC can typically be described like a stopwatch. At the illumination of a laser pulse, the laser driver sends a control signal to the TDC to start timing to keep track of the pulse’s time of flight. If a photon is detected (or a dark count occurs) the avalanche of the SPAD sends a signal to stop the TDC. The time of flight is then recorded as the time difference between the stop and start signals at the resolution of a time bin (the least-significant bit of the TDC). While it is usually desired to have the timing start simultaneously with the illumination, controlling the two systems separately adds a degree of freedom to the measurement. In particular, we implemented dither by inserting relative delays of variable duration between the true illumination time and the start of timing. We note that the underlying theory for subtractive dither requires continuously-distributed delay values. Since the values must also be known (or measured separately) in order to perform the subtraction, it is far easier in practice to implement delays with a discrete set of values. For a large-enough set of values, the effect of dither does approach the theoretical performance of continuous systems [67].
Because we do not have access to a SPAD camera, we use the system shown in Fig. 3 to emulate array imaging by raster scanning a laser spot over the scene and detecting with a single-pixel SPAD. The laser driver (PicoQuant PDL 828 “Sepia II”) is connected to a 640-nm diode laser (PicoQuant LDH P-C-640B), which has a full-width at half-maximum (FWHM) pulse duration of less than 90 ns according to the manufacturer, and is repeatedly illuminated at 10 MHz. The laser power was kept low ($\approx$ 0.25 mW) to maintain an approximately Gaussian profile and avoid pulse shape distortions. The SPAD (Micro Photon Devices PDM-series) is connected to the TCSPC module (PicoQuant HydraHarp 400). A narrowband bandpass filter (BPF) centered at the laser wavelength is placed in front of the SPAD to reduce the amount of ambient light incident on the sensor. An electronic signal synchronizing the timing of the illumination and detection is transmitted via a cable connecting the sync output from the laser driver to the sync input of the HydraHarp. A digital acquisition (DAQ) board (a LabJack U3 with an LJTick-DAC card) controls the 2-axis galvo (ThorLabs GVS012) and connects to a Marker input of the HydraHarp to simultaneously indicate when the galvo position changes. The system IRF is shown in Fig. 1.
Fig. 3. A photograph of the bistatic, subtractively-dithered single-photon lidar implementation. The only difference from the conventional setup is the external timing control of a digital delay generator, which inserts 10-picosecond-resolution delays between the laser and SPAD triggers to implement dither.
To incorporate dither into the experiment, instead of using the internal clock of the laser driver, we externally control the experiment via a digital delay generator (DDG, Stanford Research Systems DG645). As shown in Fig. 4, we connect one channel of the DDG to the external trigger input of the laser driver, and another DDG channel connects to the sync input of the HydraHarp. A third channel could be used to control the SPAD gate and further reduce ambient light detection if some depth information is known a priori. Normal lidar operation without dither can proceed by setting no delay between the channels. However, by changing the relative delay between the DDG channels, the perceived time of flight of the illumination can be changed in a predictable and repeatable manner. The DG645 is remotely controlled via a GPIB connection. In order to keep track of the dither step for each detection, an independent DAQ channel is connected to a second Marker input of the HydraHarp.
Fig. 4. A diagram showing the conventional SPL setup modified for a subtractively-dithered implementation. Varying delays $d_i$ are inserted by the digital delay generator between the laser and TCSPC triggers to implement dither. The amount of delay is controlled by the computer, which allows the delay to be subtracted from the photon detection times after acquisition.
The experiment is controlled by a MATLAB script, with dithered acquisition proceeding as follows. A grid of laser points is chosen to cover the desired portion of the scene. A set of $n_\textrm {dv}$ dither delay values is chosen to evenly subdivide $[0,\Delta )$. For each point in the grid, MATLAB sets the galvo position and resets the channel delay to zero, recording each change with a marker at the HydraHarp. The laser is pulsed, a MATLAB command increases the delay in channel AB and sends a corresponding marker to the HydraHarp. This process is repeated until all $n_\textrm {dv}$ dither delays have been stepped through in sequence. Then the galvo position is moved, the delay is reset, and the process starts again until all grid points have been addressed. When acquisition is complete, the detection data is sorted into a 3D cell array, with the laser grid position as the first two dimensions and the delay step as the third dimension. Each cell in the array contains a random number of photon detection times.
Since the system IRF is relatively long, we demonstrate the effectiveness of dither by downsampling photon detection time stamps to emulate a larger quantization bin size of $\Delta =$ 2048 ps. The dither is implemented with discrete delay steps $d_i\in \{0,\dots ,204\}\times$10 ps to evenly subdivide the quantization bins. The delays are stepped through in sequence at each pixel, and a Poisson-distributed random number of photons is detected for each delay.
We note that although the proof of principle experiments are performed here with a single-pixel detector system, incorporating subtractive dither into an array-based system is possible. Calibration could be used to negate the effects of systematic delays between TDCs. Similarly, variations in individual pixel IRFs could be calibrated to extract the necessary parameters of the EMG model, leading to slightly different estimators being applied to the detection times at each pixel. The number of dither levels used would need to be considered along with the desired number of photon detections, the IRF, and tradeoffs in other aspects of the hardware implementation.
6. Results
Dithered lidar scans were performed for the test scenes in Fig. 5. The 3D-printed depth resolution chart has boxes of various heights (labelled in Fig. 5(a) in millimeters) to help quantify the precision achieved when dither is used. The “egghead” trophy scene shows greater contrast in scene reflectivity and a wider range of scene depths. The 3D-printed mask highlights how the various methods estimate smoothly varying surfaces and fine features. The results of dithered depth estimation for an 80$\times$80-pixel scan of the depth resolution chart, an 80$\times$80-pixel scan of the egghead trophy, and a 45$\times$45-pixel scan of the mask are shown in Fig. 6.
Fig. 5. Photographs showing the 3D-printed depth resolution chart, “egghead” trophy, and 3D-printed mask used for lidar experiments. Note for the depth chart that the heights in millimeters are displayed over each block but are not physically present in the scene.
Fig. 6. Point cloud estimates for the (first row) Depth Chart, (second row) Egghead, and (third row) Mask datasets. Additional quantitative results comparing each algorithm to the baseline can be found in Table 1. The use of subtractive dither yields a large decrease in the RMSE, and the trimmed-mean estimator $\widehat {\mu }_\alpha$ based on the EMG IRF modeling produces a small additional improvement.
Table 1. RMSE performance of the depth estimates in Fig. 6 compared to several other methods. The generalized Gaussian-based methods are distinguished by whether the kurtosis matching is performed using a Gaussian (G) or exponentially-modified Gaussian (EMG) pulse shape approximation. The trimmed mean using the EMG-based kurtosis matching, denoted $\widehat {\mu }_\alpha$ (EMG), often has the best performance and consistently has lower RMSE than the sample mean.
Baseline depth estimates in Fig. 6(a) are formed by averaging thousands of photon detections per pixel measured with 4-ps resolution. Datasets with coarsely-quantized detection times were acquired for 6.15 ms per pixel for each of the datasets, and the mean number of photons per pixel (summed over all delay steps in the case of dither) in each acquisition is shown in Fig. 7. The depth estimates $\widehat {\mu }_\textrm {Q}$, formed from the pixelwise mean of quantized detections, are strongly biased by the quantization, as shown in Fig. 7. In particular, the depth chart is estimated to be a planar surface with no fine structures, the egghead trophy is flattened while its backplane is shifted to the wrong depth, and the features of the mask are lost.
Fig. 7. Detected photons per pixel (ppp) for the point clouds in Fig. 6. The Depth Chart is fairly uniform, whereas the Egghead and Mask reveal differences in reflectivity and appear to show shadows due to the bistatic lidar setup at short range.
Dithered data was acquired for 30 $\mu$s at each of the 205 dither delay steps of 10-ps each for the same total acquisition time per pixel and mean photon count as the coarse measurements without dither. We then compared the pixelwise sample mean of detection times $\widehat {\mu }_\textrm {mean}$ (Fig. 6(c)) to the estimators $\widehat {\mu }_\textrm {BG}$ and $\widehat {\mu }_\alpha$ using both the Gaussian and EMG IRF models. The sample mean result $\widehat {\mu }_\textrm {mean}$ in Fig. 6(c) shows a dramatic improvement over the result without dither, recovering many of the fine features that had been lost with the coarse quantization. Finally, Fig. 6(d) shows the results of the trimmed-mean estimator $\widehat {\mu }_\alpha$ defined in Eqs. (13) and (15), with $p_V$ determined by kurtosis matching with the dithered EMG distribution as in Eq. (17). We show only the results for $\widehat {\mu }_\alpha$ (EMG) since all methods applied to dithered data are visually similar, but $\widehat {\mu }_\alpha$ consistently improves upon $\widehat {\mu }_\textrm {mean}$. Judging from the depth chart, small features on the order of 4 mm are now visible, despite the 2048-ps bin size corresponding to approximately 30-cm resolution.
Quantitative comparisons to the baseline estimates in Table 1 show the $\widehat {\mu }_\alpha$ estimator using the EMG distribution always outperforming the sample mean applied to the dithered data and producing the best results more consistently than $\widehat {\mu }_\textrm {BG}$. The differences in estimator performance is likely because the coefficients for $\widehat {\mu }_\textrm {BG}$ are data-dependent and thus sensitive to “outliers” from the long tail of the exponential distribution. For instance, while $\widehat {\mu }_\textrm {BG}$ shows the best performance for the depth chart, it actually performs worse than the sample mean for the mask. Quantitative improvement factors in RMSE over the estimates from coarsely-quantized data range from 6.5$\times$ to 13$\times$. Thus, we have demonstrated both that subtractive dither can improve depth estimation for SPL with coarse quantization and that modeling of the resulting detection time distribution leads to estimators that yield additional performance improvements.
7. Conclusion
In this work, we presented how single-photon lidar systems can be used to make precise depth estimates despite coarse hardware timing resolution. The hardware limitations were surpassed by introducing subtractive dither, i.e., inserting short delays between the illumination and timing circuitry such that the time stamps of photons returning from the same depth could vary over multiple time bins. The addition of the dither signal was paired with depth estimators that took not only the dither signal but also the instrument response function into account. Experimental acquisitions with a laboratory lidar system confirmed that tailoring the estimator to the system response indeed outperformed more naïve existing estimators. The use of subtractive dither resulted in depth error improvement factors of up to 13.
One potential application of this work is in new SPAD array designs. Rather than provide each pixel in an array with a high-resolution TDC, dither could be implemented globally for the array, thus enabling TCSPC in megapixel SPAD arrays. Photons at each pixel would be time stamped at lower resolution, and a common delay would be added to all pixels that would vary for each illumination period. Keeping track of the delay for each period would then allow depths to be reconstructed at finer resolution than could otherwise be achieved from the hardware. Future work should address both the hardware challenges of implementing subtractive dither in a detector array as well as robust algorithmic approaches to dealing with a broader array of experimental conditions, including multiple depths per pixel or strong ambient light, which would introduce a mismatch with the EMG distribution assumed by the estimator. Additional study could also explore the use of spatial regularization to further improve the photon efficiency of dithered depth imaging.
Funding
National Science Foundation (1422034, 1815896, 1955219); Defense Advanced Research Projects Agency (HR0011-16-C-0030); Charles Stark Draper Laboratory (Draper Fellowship).
Acknowledgments
The authors would like to thank Dr. J. Hollman of Draper for assistance with the optical setup and Dr. D. Shin for helpful discussions on the role of nonsubtractive dither in improving estimation accuracy.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. A. Kirmani, D. Venkatraman, D. Shin, A. Colaço, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, “First-photon imaging,” Science 343(6166), 58–61 (2014). [CrossRef]
2. D. Shin, A. Kirmani, V. K. Goyal, and J. H. Shapiro, “Photon-efficient computational 3-D and reflectivity imaging with single-photon detectors,” IEEE Trans. Comput. Imaging 1(2), 112–125 (2015). [CrossRef]
3. Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, and S. McLaughlin, “Lidar waveform-based analysis of depth images constructed using sparse single-photon data,” IEEE Trans. on Image Process. 25(5), 1935–1946 (2016). [CrossRef]
4. J. Rapp and V. K. Goyal, “A few photons among many: Unmixing signal and noise for photon-efficient active imaging,” IEEE Trans. Comput. Imaging 3(3), 445–459 (2017). [CrossRef]
5. D. B. Lindell, M. O’Toole, and G. Wetzstein, “Single-photon 3D imaging with deep sensor fusion,” ACM Trans. Graph. 37(4), 1–12 (2018). [CrossRef]
6. A. M. Pawlikowska, A. Halimi, R. A. Lamb, and G. S. Buller, “Single-photon three-dimensional imaging at up to 10 kilometers range,” Opt. Express 25(10), 11919–11931 (2017). [CrossRef]
7. Z.-P. Li, X. Huang, P.-Y. Jiang, Y. Hong, C. Yu, Y. Cao, J. Zhang, F. Xu, and J.-W. Pan, “Super-resolution single-photon imaging at 8.2 kilometers,” Opt. Express 28(3), 4076 (2020). [CrossRef]
8. Z.-P. Li, X. Huang, Y. Cao, B. Wang, Y.-H. Li, W. Jin, C. Yu, J. Zhang, Q. Zhang, C.-Z. Peng, F. Xu, and J.-W. Pan, “Single-photon computational 3D imaging at 45 km,” Photonics Res. 8(9), 1532–1540 (2020). [CrossRef]
9. D. Shin, F. Xu, F. Wong, J. H. Shapiro, and V. K. Goyal, “Computational multi-depth single-photon imaging,” Opt. Express 24(3), 1873–1888 (2016). [CrossRef]
10. J. Tachella, Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, S. McLaughlin, and J.-Y. Tourneret, “Bayesian 3D reconstruction of complex scenes from single-photon lidar data,” SIAM J. Imaging Sci. 12(1), 521–550 (2019). [CrossRef]
11. R. Tobin, A. Halimi, A. McCarthy, M. Laurenzis, F. Christnacher, and G. S. Buller, “Three-dimensional single-photon imaging through obscurants,” Opt. Express 27(4), 4590–4611 (2019). [CrossRef]
12. A. Halimi, R. Tobin, A. McCarthy, J. M. Bioucas-Dias, S. McLaughlin, and G. S. Buller, “Robust restoration of sparse multidimensional single-photon LiDAR images,” IEEE Trans. Comput. Imaging 6, 138–152 (2020). [CrossRef]
13. A. Maccarone, F. Mattioli Della Rocca, A. McCarthy, R. Henderson, and G. S. Buller, “Three-dimensional imaging of stationary and moving targets in turbid underwater environments using a single-photon detector array,” Opt. Express 27(20), 28437 (2019). [CrossRef]
14. J. Hecht, “Lidar for self-driving cars,” Opt. Photonics News 29(1), 26–33 (2018). [CrossRef]
15. F. Villa, R. Lussana, D. Bronzi, S. Tisa, A. Tosi, F. Zappa, A. Dalla Mora, D. Contini, D. Durini, S. Weyers, and W. Brockherde, “CMOS imager with 1024 SPADs and TDCs for single-photon timing and 3-D time-of-flight,” IEEE J. Sel. Top. Quantum Electron. 20(6), 364–373 (2014). [CrossRef]
16. D. Shin, F. Xu, D. Venkatraman, R. Lussana, F. Villa, F. Zappa, V. K. Goyal, F. N. Wong, and J. H. Shapiro, “Photon-efficient imaging with a single-photon camera,” Nat. Commun. 7(1), 12046 (2016). [CrossRef]
17. S. Burri, H. Homulle, C. Bruschini, and E. Charbon, “LinoSPAD: a time-resolved 256×1 CMOS SPAD line sensor system featuring 64 FPGA-based TDC channels running at up to 8.5 giga-events per second,” Proc. SPIE 9899, 98990D (2016). [CrossRef]
18. D. B. Lindell, M. O’Toole, and G. Wetzstein, “Towards transient imaging at interactive rates with single-photon detectors,” in Proc. IEEE Int. Conf. Comput. Photog., (IEEE, 2018).
19. K. Morimoto, A. Ardelean, M.-L. Wu, A. C. Ulku, I. M. Antolovic, C. Bruschini, and E. Charbon, “Megapixel time-gated SPAD image sensor for 2D and 3D imaging applications,” Optica 7(4), 346–354 (2020). [CrossRef]
20. J. Tachella, Y. Altmann, N. Mellado, A. McCarthy, R. Tobin, G. S. Buller, J. Y. Tourneret, and S. McLaughlin, “Real-time 3D reconstruction from single-photon lidar data using plug-and-play point cloud denoisers,” Nat. Commun. 10(1), 4984 (2019). [CrossRef]
21. J. Rapp, J. Tachella, Y. Altmann, S. McLaughlin, and V. K. Goyal, “Advances in single-photon lidar for autonomous vehicles: Working principles, challenges, and recent advances,” IEEE Signal Process. Mag. 37(4), 62–71 (2020). [CrossRef]
22. J. Richardson, R. Walker, L. Grant, D. Stoppa, F. Borghetti, E. Charbon, M. Gersbach, and R. K. Henderson, “A 32×32 50ps resolution 10 bit time to digital converter array in 130nm CMOS for time correlated imaging,” in Proc. IEEE Custom Integr. Circuits Conf., (2009), pp. 77–80.
23. C. Veerappan, J. Richardson, R. Walker, D.-U. Li, M. W. Fishburn, Y. Aruyama, D. Stoppa, R. K. Henderson, F. Borghetti, M. Gersbach, and E. Charbon, “A 160×128 single-photon image sensor with on-pixel 55ps 10b time-to-digital converter,” in Proc. IEEE Int. Solid-State Circuits Conf., (2011), pp. 312–314.
24. S. Lindner, C. Zhang, I. M. Antolovic, M. Wolf, and E. Charbon, “A 252 × 144 SPAD pixel FLASH LiDAR with 1728 dual-clock 48.8 ps TDCs, integrated histogramming and 14.9-to-1 compression in 180nm CMOS technology,” in IEEE Symp. VLSI Circuits, (2018), pp. 69–70.
25. C. Zhang, S. Lindner, I. M. Antolovic, J. Mata Pavia, M. Wolf, and E. Charbon, “A 30-frames/s, 252 × 144 SPAD flash LiDAR with 1728 dual-clock 48.8-ps TDCs, and pixel-wise integrated histogramming,” IEEE J. Solid-State Circuits 54(4), 1137–1151 (2019). [CrossRef]
26. R. K. Henderson, N. Johnston, F. Mattioli Della Rocca, H. Chen, D. Day-Uei Li, G. Hungerford, R. Hirsch, D. McLoskey, P. Yip, and D. J. Birch, “A 192 × 128 time correlated SPAD image sensor in 40-nm CMOS technology,” IEEE J. Solid-State Circuits 54(7), 1907–1916 (2019). [CrossRef]
27. A. R. Ximenes, P. Padmanabhan, M. J. Lee, Y. Yamashita, D.-N. Yaung, and E. Charbon, “A modular, direct time-of-flight depth sensor in 45/65-nm 3-D-stacked CMOS technology,” IEEE J. Solid-State Circuits 54(11), 3203–3214 (2019). [CrossRef]
28. M. Perenzoni, D. Perenzoni, and D. Stoppa, “A 64 × 64-pixels digital silicon photomultiplier direct TOF sensor with 100-MPhotons/s/pixel background rejection and imaging/altimeter mode with 0.14% precision up to 6 km for spacecraft navigation and landing,” IEEE J. Solid-State Circuits 52(1), 151–160 (2017). [CrossRef]
29. R. K. Henderson, N. Johnston, S. W. Hutchings, I. Gyongy, T. A. Abbas, N. Dutton, M. Tyler, S. Chan, and J. Leach, “A 256×256 40nm/90nm CMOS 3D-stacked 120dB dynamic-range reconfigurable time-resolved SPAD imager,” in Proc. IEEE Int. Solid-State Circuits Conf., (2019), pp. 106–108.
30. I. Gyongy, S. W. Hutchings, A. Halimi, M. Tyler, S. Chan, F. Zhu, S. McLaughlin, R. K. Henderson, and J. Leach, “High-speed 3D sensing via hybrid-mode imaging and guided upsampling,” Optica 7(10), 1253–1260 (2020). [CrossRef]
31. L. R. A. MacColl, Fundamental theory of servomechanisms (D. Van Nostrand Company, Inc., 1945).
32. W. M. Goodall, “Television by pulse code modulation,” Bell Syst. Tech. J. 30(1), 33–49 (1951). [CrossRef]
33. G. G. Furman, “Improving performance of quantizer feedback systems by use of external dither,” M.S. thesis, Massachusetts Institute of Technology (1959).
34. R. C. Jaffe, “Causal and statistical analyses of dithered systems containing three-level quantizers,” M.S. thesis, Massachusetts Institute of Technology (1959).
35. L. Roberts, “Picture coding using pseudo-random noise,” IEEE Trans. Inf. Theory 8(2), 145–154 (1962). [CrossRef]
36. L. Schuchman, “Dither signals and their effect on quantization noise,” IEEE Trans. Commun. 12(4), 162–165 (1964). [CrossRef]
37. N. S. Jayant and L. R. Rabiner, “The application of dither to the quantization of speech signals,” Bell Syst. Tech. J. 51(6), 1293–1304 (1972). [CrossRef]
38. J. Rapp, R. M. A. Dawson, and V. K. Goyal, “Dither-enhanced lidar,” in OSA Imaging and Applied Optics Congress, (Orlando, FL, 2018), p. JW4A.38.
39. J. Rapp, R. M. A. Dawson, and V. K. Goyal, “Improving lidar depth resolution with dither,” in Proc. IEEE Int. Conf. Image Process., (2018), pp. 1553–1557.
40. J. Rapp, R. M. A. Dawson, and V. K. Goyal, “Estimation from quantized Gaussian measurements: When and how to use dither,” IEEE Trans. Signal Process. 67(13), 3424–3438 (2019). [CrossRef]
41. J. Busck and H. Heiselberg, “Gated viewing and high-accuracy three-dimensional laser radar,” Appl. Opt. 43(24), 4705–4710 (2004). [CrossRef]
42. L. F. Gillespie, “Apparent illuminance as a function of range in gated, laser night-viewing systems,” J. Opt. Soc. Am. 56(7), 883–887 (1966). [CrossRef]
43. F. Christnacher, S. Schertzer, N. Metzger, E. Bacher, M. Laurenzis, and R. Habermacher, “Influence of gating and of the gate shape on the penetration capacity of range-gated active imaging in scattering environments,” Opt. Express 23(26), 32897 (2015). [CrossRef]
44. X. Ren, P. W. R. Connolly, A. Halimi, Y. Altmann, S. McLaughlin, I. Gyongy, R. K. Henderson, and G. S. Buller, “High-resolution depth profiling using a range-gated CMOS SPAD quanta image sensor,” Opt. Express 26(5), 5541 (2018). [CrossRef]
45. P. Rehain, Y. M. Sua, S. Zhu, I. Dickson, B. Muthuswamy, J. Ramanathan, A. Shahverdi, and Y.-P. Huang, “Noise-tolerant single photon sensitive three-dimensional imager,” Nat. Commun. 11(1), 921 (2020). [CrossRef]
46. Z. Chen, R. Fan, X. Li, Z. Dong, Z. Zhou, G. Ye, and D. Chen, “Accuracy improvement of imaging lidar based on time-correlated single-photon counting using three laser beams,” Opt. Commun. 429, 175–179 (2018). [CrossRef]
47. A. Raghuram, A. Pediredla, S. G. Narasimhan, I. Gkioulekas, and A. Veeraraghavan, “STORM: Super-resolving transients by OveRsampled measurements,” in Proc. IEEE Int. Conf. Comput. Photog., (2019), pp. 44–54.
48. R. M. Gray and T. G. Stockham, “Dithered quantizers,” IEEE Trans. Inform. Theory 39(3), 805–812 (1993). [CrossRef]
49. S. Mandai, E. Venialgo, and E. Charbon, “Timing optimization utilizing order statistics and multichannel digital silicon photomultipliers,” Opt. Lett. 39(3), 552–554 (2014). [CrossRef]
50. S. Cova, M. Ghioni, M. A. Itzler, J. C. Bienfang, and A. Restelli, Semiconductor-based detectors,” in Single-Photon Generation and Detection, A. Migdall, S. V. Polyakov, J. Fan, and J. C. Bienfang, eds. (Academic Press, 2013), chap. 4, pp. 83–146.
51. W. Becker, Advanced Time-Correlated Single Photon Counting Techniques (Springer, 2005).
52. M. J. Stevens, R. H. Hadfield, R. E. Schwall, S. W. Nam, R. P. Mirin, and J. A. Gupta, “Fast lifetime measurements of infrared emitters using a low-jitter superconducting single-photon detector,” Appl. Phys. Lett. 89(3), 031109 (2006). [CrossRef]
53. Q. Hernandez, D. Gutierrez, and A. Jarabo, “A computational model of a single-photon avalanche diode sensor for transient imaging,” arXiv 1703.02635 (2017).
54. M. Gladney, B. F. Dowden, J. D. Swalen, R. Lusebrink, and C. H. Sederholm, “Computer-assisted gas-liquid chromatography,” Anal. Chem. 41(7), 883–888 (1969). [CrossRef]
55. E. Grushka, “Characterisation of exponentially-modified Gaussian peaks in chromatography,” Anal. Chem. 44(11), 1733–1738 (1972). [CrossRef]
56. M. R. W. Dawson, “Fitting the ex-Gaussian equation to reaction time distributions,” Behav. Res. Methods, Instruments, & Comput. 20(1), 54–57 (1988). [CrossRef]
57. A. P. Greeley, T. A. Neumann, N. T. Kurtz, T. Markus, and A. J. Martino, “Characterizing the system impulse response function from photon-counting LiDAR data,” IEEE Transactions on Geosci. Remote. Sens. 57(9), 6542–6551 (2019). [CrossRef]
58. C. Moret-Tatay, D. Gamermann, E. Navarro-Pardo, and P. F. d. C. Castellá, “ExGUtils: A python package for statistical analysis with the ex-Gaussian probability density,” Front. Psychol. 9, 612 (2018). [CrossRef]
59. B. Zandbelt, “exgauss: a MATLAB toolbox for fitting the ex-Gaussian distribution to response time data,” https://figshare.com/articles/exgauss/971318 (2014).
60. Y. Lacouture and D. Cousineau, “How to use MATLAB to fit the ex-Gaussian and other probability functions to a distribution of response times,” Tutorials Quant. Methods for Psychol. 4(1), 35–45 (2008). [CrossRef]
61. N. C. Beaulieu and Q. Guo, “Novel estimator for the location parameter of the generalized Gaussian distribution,” IEEE Commun. Lett. 16(12), 2064–2067 (2012). [CrossRef]
62. A. Restrepo and A. C. Bovik, “Adaptive trimmed mean filters for image restoration,” IEEE Trans. Acoust. Speech Signal Process. 36(8), 1326–1337 (1988). [CrossRef]
63. P. Daniell, “Observations weighted according to order,” Am. J. Math. 42(4), 222–236 (1920). [CrossRef]
64. P. J. Bickel and E. L. Lehmann, “Descriptive statistics for nonparametric models: II. location,” Ann. Statist. 3(5), 1045–1069 (1975). [CrossRef]
65. E. H. Lloyd, “Least-squares estimation of location and scale parameters using order statistics,” Biometrika 39(1-2), 88–95 (1952). [CrossRef]
66. A. K. Gupta, “Estimation of the mean and standard deviation of a normal population from a censored sample,” Biometrika 39(3-4), 260–273 (1952). [CrossRef]
67. M. F. Wagdy, “Effect of various dither forms on quantization errors of ideal A/D converters,” IEEE Trans. Instrum. Meas. 38(4), 850–855 (1989). [CrossRef]
