Photon-efficient imaging with a single-photon camera — TEI-tagged view

Dongeek ShinpersName, Feihu XupersName, Dheera VenkatramanpersName, Rudi LussanapersName, Federica VillapersName, Franco ZappapersName, Vivek K. GoyalpersName, Franco N. C. WongpersName, Jeffrey H. ShapiropersName (2016date). “Photon-efficient imaging with a single-photon camera.” Nature Communicationstitle·j.

Object of study — depth and reflectivity imagingterm·key from very few detected photons with a single-photon cameraterm·key

Research question — whether accurate 3D and reflectivityterm·key images can be recovered from about one photon per pixel with a coarse-timing SPADterm·key array

Contribution — an array-specific reconstruction algorithm exploiting transverse smoothness and longitudinal sparsity of natural scenes

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Photon-efficient imaging with a single-photon cameraterm·key Abstract Reconstructing a scene's 3D structure and reflectivityterm·key accurately with an active imaging system operating in low-light-level conditions has wide-ranging applications, spanning biological imaging to remote sensing. Here we propose and experimentally demonstrate a depth and reflectivity imagingterm·key system with a single-photon cameraterm·key that generates high-quality images from ∼1 detected signal photon per pixel. Previous achievements of similar photon efficiencyterm·key have been with conventional raster-scanning data collection using single-pixel photon counters capable of ∼10-ps time taggingterm·key. In contrast, our camera's detector array requires highly parallelized time-to-digital conversions with photon time-tagging accuracy limited to ∼ns. Thus, we develop an array-specific algorithm that converts coarsely time-binned photon detections to highly accurate scene depth and reflectivityterm·key by exploiting both the transverse smoothness and longitudinal sparsity of natural scenes. By overcoming the coarse time resolution of the array, our framework uniquely achieves high photon efficiencyterm·key in a relatively short acquisition time.

Active optical imaging systemsterm·mth use their own light sources to recover scene information. To suppress photon noise inherent in the optical detection process, they typically require a large number of photon detections. For example, a commercially available flash camera typically collects >10 9 photons (10 3 photons per pixel in a 1 megapixel image) to provide the user with a single photograph 1 . However, in remote sensing of a dynamic scene at a long standoff distance, as well as in microscope imaging of delicate biological samples, limitations on the optical flux and integration time preclude the collection of such a large number of photons 2 3 . A key challenge in such scenarios is to make use of a small number of photon detections to accurately recover the desired scene information. Exacerbating the difficulty is that, for any fixed total acquisition time, serial acquisition through raster scanning reduces the number of photon detections per pixel. Accurate recovery from a small number of photon detections has not previously been achieved in conjunction with parallel acquisition at a large number of pixels, as in a conventional digital camera.

Our interest is in the simultaneous reconstruction of scene three-dimensional (3D) structure and reflectivityterm·key using a small number of photons, something that is important in many real-world imaging scenarios 4 5 6 . Accurately measuring distance and estimating a scene's 3D structure can be done from time-of-flight data collected with a pulsed-source light detection and ranging system 7 8 9 . For applications specific to low-light-level 3D imaging, detectors that can resolve individual photon detections—from either photomultiplication 10 or Geiger-mode avalanche operation 11 —must be used in conjunction with a time correlator. These time-correlated single-photon detectors provide extraordinary sensitivity in time-tagging photon detections, as shown by the authors of ref. 12 , who used a time-correlated single-photon avalanche diode (SPADterm·key) array to track periodic light pulses in flight.

The state-of-the-art in high photon-efficiency depth and reflectivity imagingterm·key was established by the authors of first-photon imaging (FPIterm·mth) 13 , who demonstrated accurate 3D and reflectivityterm·key recovery from the first detected photon at each pixel. Their set-up, which used raster scanning and a time-correlated single SPADterm·key detector, required exactly one photon detection at each pixel, making each pixel's acquisition time a random variable. Consequently, FPIterm·mth is not applicable to operation using a SPAD cameraterm·mth 14 15 16 17 18 19 20 21 —all of whose pixels must have the same acquisition time—thus precluding FPIterm·mth's reaping the marked image-acquisition speedup that the camera's detector array affords 12 22 23 . Although there have been extensions of FPIterm·mth to the fixed acquisition-time operation needed for array detection 24 25 , both their theoretical modeling and experimental validations were still limited to raster scanning, with a single SPADterm·key detector. As a result, they ignored the limitations of currently available SPAD camerasterm·mth—much poorer time-tagging performance and pixel-to-pixel variations of SPADterm·key properties—implying that these initial fixed acquisition-time (pseudo-array) frameworks will yield sub-optimal depth and reflectivityterm·key reconstructions when used with low-light experimental data from an actual SPAD cameraterm·mth.

Here we propose and demonstrate a photon-efficient 3D structure and reflectivity imagingterm·key technique that can deal with the aforementioned constraints that SPAD camerasterm·mth impose. We give the first experimental demonstration of accurate time-correlated SPAD-camera imaging of natural scenes obtained from ∼1 detected signal photon per pixel on average. Unlike prior work, our framework achieves high photon efficiencyterm·key by exploiting the scene's structural information in both the transverse and the longitudinal domains to censor extraneous (background light and dark count) detections from the SPADterm·key array. Earlier works that exploit longitudinal sparsity only in a pixel-by-pixel manner require more detected signal photons to produce accurate estimates 26 27 . Because our new imager achieves highly photon-efficient imaging in a short data-acquisition time, it paves the way for dynamic and noise-tolerant active optical imaging applications in science and technology.

Results Imaging set-up Our experimental set-up is illustrated in Fig. 1 . The illumination source was a pulsed laser diode (PicoQuant LDH series with a 640-nm center wavelength), whose original output-pulse duration was increased to a full-width at half-maximum of ∼2.5 ns (that is, a root mean square (r.m.s.) value of T p ≈1 ns). The laser diode was pulsed at a T r ≈50 ns repetition period set by the SPADterm·key array's trigger output. A diffuser plate spatially spread the laser pulses to flood illuminate the scene of interest. An incandescent lamp injected unwanted background light into the camera. The lamp's power was adjusted so that (averaged over the region that was imaged) each detected photon was equally likely to be due to signal (backreflected laser) light or background light. A standard Canon FL series photographic lens focused the signal plus background light on the SPADterm·key array. Each photon detection from the array was time tagged relative to the time of the most recently transmitted laser pulse and recorded ( Supplementary Methods ).

The SPADterm·key array 18 20 , covering a 4.8 × 4.8-mm footprint, consists of 32 × 32 pixels of fully independent Si SPADsterm·key and complementary metal-oxide-semiconductor (CMOS)-based electronic circuitry that includes a time-to-digital converter for each SPADterm·key detector. The SPADterm·key within each 150 × 150-μm pixel has a 30-μm-diameter circular active region, giving the array a 3.14% fill factor. At the 640-nm operating wavelength, each array element's photon-detection efficiency is ∼20% and its dark count rate is ∼100 Hz at room temperature. To extend the region that could be imaged and increase the number of pixels, we used multiple image scans to form a larger-size composite image. In particular, we mounted the SPADterm·key array on a feedback-controlled, two-axis motorized translation stage to produce images with N x × N y =384 × 384 pixels ( Supplementary Figs 1a,b and 2a–c ).

The SPADterm·key array has a Δ=390 ps time resolution set by its internal clock rate. We set each acquisition frame length to 65 μs, with a gate-on time of 16 μs and a gate-off time of 49 μs for limiting power dissipation of the chip and for data transfer. At the start of each frame, the SPADterm·key array was set to trigger the laser to generate pulses at a ∼20 MHz repetition rate. Hence, in the 16 μs gate-on time of each frame, ∼320 pulses illuminated the scene Supplementary Fig. 3 ).

Observation model We define to be the scene's 3D structure and reflectivityterm·key that we aim to recover, and we let be the average rates of background-light plus dark-count detections. Flood illumination of the scene at time t =0 with a photon-flux pulse s ( t ) then results in the following Poisson-process rate function for ( i , j )-th pixel of the composite image: where η i , j ∈(0,1] is the ( i , j )th detector's photon-detection efficiency and c is the speed of light. Fabrication imperfections of the SPADterm·key array cause some pixels to have inordinately high dark-count rates , making their detection times uninformative in our imaging experimentsterm·mth because they are predominantly from dark counts. Thus we performed camera calibration to determine the set of these ‘hot pixels' (2% of all pixels in our case) so that their outputs could be ignored in the processing of the imaging data.

We define to be the total number of time bins in which the photon detections can be found, and let C i , j , k be the observed number of photon counts in the k th time bin for pixel ( i , j ) after n s pulsed-illumination trials. By the theory of photon counting 28 , we have that C i , j , k 's statistical distribution is for k =1,2,…, N z , where we have assumed that the pulse repetition period is long enough to preclude pulse aliasing artifacts. Also, we operate in a low-flux condition such that ∑ k =1 N z C i , j , k , the total number of detections at a pixel, is much less than n s , the total number of illumination pulses to avoid pulse-pileup distortions. Our imaging problem is then to construct accurate estimates, and , of the scene's reflectivityterm·key A and 3D structure Z , using the sparse photon-detection data .

3D structure and reflectivityterm·key reconstruction In the low-flux regime, wherein there are very few detections and many of them are extraneous, an algorithm that relies solely on the aforementioned pixelwise photodetection statistics has very limited robustness. We aim to achieve high photon efficiencyterm·key by combining those photodetection statistics with prior information about natural scenes.

Most natural and man-made scenes have strong spatial correlations among neighbouring pixels in both transverse and longitudinal measurements, punctuated by sharp boundaries 29 . While conventional works normally treat each pixel independently, our imaging framework exploits these correlations to censor/remove extraneous (and randomly distributed) photon-detection events due to background light and detector dark counts. It should be noted that unlike noise mitigation via spatial filteringterm·mth and averaging, in which fine spatial features are washed out due to oversmoothing, our technique retains the spatial resolution set by the SPADterm·key array. Our reconstruction algorithm optimizes between two constraints for a given set of censored measurements: that the 3D and reflectivityterm·key image estimates come from a scene that is correlated in both the transverse and longitudinal domains, and that the estimates employ the Poisson statistics of the raw single-photon measurements.

The implementation of our reconstruction algorithm can be divided into the following three steps ( Fig. 2 ; Supplementary Note 1 ).

Step 1: natural scenes have reflectivities that are spatially correlated—the reflectivityterm·key at a given pixel tends to be similar to the values at its nearest neighbours—with abrupt transitions at the boundaries between objects. We exploit these correlations by imposing a transverse-smoothness constraint using the total-variation (TV) norm 30 on our reflectivityterm·key image. In this process, we ignore data from the hot-pixel set . The final reflectivityterm·key image is thus obtained by solving a regularized optimization problem.

Step 2: natural scenes have a finite number of reflectors that are clustered in depth. It follows that in an acquisition without background-light or dark-count detections, the set of detection times collected over the entire scene would have a histogram with N z bins that possesses non-zero entries in only a small number of small subintervals. This longitudinal sparsity constraint is enforced in our algorithm by solving a sparse deconvolution problem from the coarsely time-binned photon-detection data, which is specific to the array imaging set-up, to obtain a small number of representative scene depths. Raw photon-detection events at times corresponding to depths differing by more than cT p /2 from the representative scene depths are censored. As step 2 has identified coarse depth clusters of the scene objects, the next step of the algorithm uses the filtered set of photon detections to determine a high-resolution depth image within all identified clusters.

Step 3: similar to what was done in step 1 for reflectivityterm·key estimation, we impose a TV-norm spatial smoothness constraint on our depth image, where data from the hot-pixel set and censored detections at the remaining pixels are ignored. Thus, we obtain by solving a regularized optimization problem.

Reconstruction results Figure 3 shows experimental results of 3D structure and reflectivityterm·key reconstructions for a scene comprised of a mannequin and sunflower when, averaged over the scene, there was ∼1 signal photon detected per pixel and ∼1 extraneous (background light plus dark count) detection per pixel. In our experimentsterm·mth, the per-pixel average photon-count rates of backreflected waveform and background-light plus dark-count response were 1,089 counts/s and 995 counts/s, respectively. The image resolution was 384 × 384 for this experimentterm·mth. We compare our proposed method with the baseline pixelwise imaging method that uses filtered histograms 22 and the state-of-the-art pseudo-array imaging method 25 .

From the visualization of reflectivityterm·key overlaid on depth, we observe that the baseline pixelwise imaging method ( Fig. 3a,e ) generates noisy depth and reflectivityterm·key images without useful scene features, owing to the combination of low-flux operation and high-background detections plus detector dark counts. In contrast, the existing pseudo-array method—which exploits transverse spatial correlations, but presumes constant B i , j —gives a reflectivityterm·key image that captures overall object features, but is oversmoothed to mitigate hot-pixel contributions ( Fig. 3b ). Furthermore, because the pseudo-array method presumes the 10-ps-class time taggingterm·key of a single-element SPADterm·key that is used in raster-scanning set-ups, its depth image fails to reproduce the 3D structure of the mannequin's face from the ns-class time taggingterm·key afforded by our SPAD cameraterm·mth's detector array. In particular, it overestimates the head's dimensions and oversmooths the facial features ( Fig. 3f ), whereas our array-specific method accurately captures the scene's 3D structure and reflectivityterm·key ( Fig. 3c,g ). This accuracy can be seen by comparing our framework's result with the high-flux pixelwise depth and reflectivityterm·key images ( Fig. 3d,h )—obtained by detecting 550 signal photons per pixel and performing time-gated pixelwise processing—that serve as ground-truth proxies for the scene's actual depth and reflectivityterm·key. For the fairest comparisons in Fig. 3 , each algorithm—baseline pixelwise processing, pseudo-array processing and our new framework—had its parameters tuned to minimize the mean-squared degradation from the ground-truth proxies.

The depth error maps in Fig. 3i–k quantify the resolution improvements from our imager over the existing ones for this low-flux imaging experimentterm·mth. Although the mean number of signal photon detections is ∼1 per pixel, in the high-reflectivity facial regions of the mannequin the average is ∼8 signal photon detections per pixel, while the number of signal photons detected at almost every pixel in the background portion of the scene is 0. Despite this fact, the pseudo-array imaging technique leads to a face estimate with high depth error due to oversmoothing incurred in its effort to mitigate background noise. It particularly suffers at reconstructing the depth boundaries with low-photon counts as well. Compared with conventional methods, our framework gives a much better estimate of the full 3D structure. We also study how the depth error of our framework depends on the number of photon counts at a given pixel. For example, our framework gives errors of 4.4 and 0.9 cm at pixels (260, 121) and (107, 187), which correspond to 1-photon-count depth boundary region and 8-photon-count mannequin face, respectively. Overall, we observe a negative correlation between the depth error and the number of photons detected at a pixel for our method ( Supplementary Fig. 4 ). Recall that the time bin duration of each pixel of the SPAD cameraterm·mth is Δ=390 ps, corresponding to c Δ/2≈6-cm depth resolution. Overall, our imager successfully recovers depth with mean absolute error of 2 cm and thus with sub-bin-duration resolution, while existing methods fail to do so.

In terms of measuring the improvements in gray-scale imaging, we can compute the peak signal-to-noise ratio (PSNR) between the reference ground-truth reflectivityterm·key and the estimated reflectivityterm·key. While the PSNR values of the conventional pixelwise estimation and pseudo-array imaging are 19.3 and 24.9 dB, respectively, that of our method is 29.1 dB; it improves over both existing methods by at least 4 dB. We emphasize the difficulty of single-photon imaging in our set-up by computing that the SNR of the time-of-flight of a single photon ranges from −2.2 to 7.8 dB ( Supplementary Note 2 ). We also note that varying the regularization parameters also affects the imaging performance ( Supplementary Fig. 5 ). Lastly, the robustness of our reconstruction algorithm is evaluated by imaging an entirely different scene consisting of watering can and basketball, using regularization parameters pre-trained on the mannequin scene ( Supplementary Fig. 6 ).

Choice of laser pulse root mean square time duration For a transform-limited laser pulse, such as the Gaussian s ( t ) that our imaging framework presumes, the r.m.s. time durationterm·mth T p is a direct measure of system bandwidth. As such, it has an impact on the depth-imaging accuracy in low-flux operation. This impact is borne out by the simulation results in Fig. 4 , where we see that the pulse waveform with the shortest r.m.s. duration does not provide the best depth recovery. Thus, in our experimentsterm·mth, we broadened the laser's output pulse to T p ≈1 ns. The full-width at half-maximum is then 2.4 ns. This pulse duration allowed our framework to have a mean absolute depth error of 2 cm and resolve depth features well below the c Δ/2≈6 cm value set by the SPADterm·key array's 390-ps-duration time bins (see Fig. 4 for details on depth-recovery accuracy versus r.m.s. pulse duration).

For application of our framework to different array technology, the optimal pulsewidth should scale with the SPAD cameraterm·mth's time-bin duration. For example, if our SPADterm·key hardware were replaced to improve the timing resolution to the 50-ps range 17 , we would want to make sure the r.m.s. pulsewidth remains approximately three times longer than the time bin (or the full-width at half-maximum is approximately six times longer) based on our method for choosing optimal pulsewidth from Fig. 4 . Thus, for accurate single-photon imaging with a 50-ps time-binning SPADterm·key array, we would shorten our pulse from 1.1 ns to 140 ps.

Discussion We have proposed and demonstrated a SPAD-camera imaging framework that generates highly accurate images of a scene's 3D structure and reflectivityterm·key from ∼1 detected signal photon per pixel, despite the presence of extraneous detections at roughly the same rate from background light and dark counts. By explicitly modeling the limited single-photon time-tagging resolution of SPAD-array imagers, our framework markedly improves reconstruction accuracy in this low-flux regime as compared with what is achieved with existing methods. The photon efficiencyterm·key of our proposed framework is quantified in Fig. 5 , where we have plotted the sub-bin-duration depth error it incurs in imaging the mannequin and sunflower scene versus the average number of detected signal photons per pixel. For this task, our algorithm realizes centimetre-class depth resolution down to <1 detected signal photon per pixel, while the baseline pixelwise imager's depth resolution is more than an order of magnitude worse because of its inability to cope with extraneous detections.

Because our framework employs a SPAD cameraterm·mth for highly photon-efficient imaging, it opens up new ways to image 3D structure and reflectivityterm·key on very short time scales, while requiring very few photon detections. Hence, it could find widespread use in applications that require fast and accurate imaging using extremely small amounts of light, such as remote terrestrial mapping 31 , seismic imaging 32 , fluorescence profiling 2 and astronomy 33 . We emphasize, in this regard, that our framework affords automatic rejection of ambient-light and dark-count noise effects without requiring sophisticated time-gating hardware. It follows that our imager could also enable rapid and noise-tolerant 3D vision for self-navigating advanced robotic systems, such as unmanned aerial vehicles and exploration rovers 34 .

Methods Reconstruction algorithm Before initiating our three-step imaging algorithm, we first performed calibration measurements to: identify , the SPADterm·key array's set of hot pixels; obtain the average background-light plus dark-count rates for the remaining pixels; and determine the laser pulse's r.m.s. time durationterm·mth ( Supplementary Fig. 7a–c ). For hot-pixel identification, we placed the SPAD cameraterm·mth in a dark room, and identified pixels with dark-count rate >150 counts/s, since the standard pixel of our camera should have a dark-count rate of 100 counts/s. The background-light plus dark-count rate at each pixel was identified by simply measuring the count rate, when the background-light source was on but the laser was off. Finally, the laser-pulse shape was calibrated by measuring the time histogram of 3,344 photon detections from a white calibration surface target that was placed ∼1 m away from the imaging set-up. It turned out that: ∼2% of our camera's 1,024 pixels were placed in ; the background-light plus dark-count rates were indeed spatially varying across the remaining pixels; and the laser pulse's time duration was T p ≈1 ns and reasonably approximated as a Gaussian.

We then proceed to step 1 of the reconstruction algorithm: we estimate reflectivityterm·key by combining the Poisson statistics of photon counts ( Supplementary Fig. 8a ) with a TV-norm smoothness constraint on the estimated reflectivityterm·key—while censoring the set of hot pixels—to write the optimization as a TV-regularized, Poisson image inpainting problem. This optimization problem is convex in the reflectivityterm·key image variable A , which allows us to solve it in a computationally efficient manner with simple projected gradient methods 35 . This step of the algorithm inputs a parameter τ A that controls the degree of spatial smoothness of the final reflectivityterm·key image estimate. For a 384 × 384 image, the processing time of step 1 was ∼6 s on a standard laptop computer ( Supplementary Note 1 ).

For step 2 of the reconstruction algorithm, we filtered the photon-detection data set to impose the longitudinal constraint that the scene has a sparse set of reflectors. This is because the scaled detection-time histogram hist( c T /2) that has been corrected for the average background-light plus dark-count detections per bin is a proxy solution for hist( Z Δ ), where hist( Z Δ ) is a size– N z histogram that bins the scene's 3D structure at the camera's c Δ/2 native range resolution. We used orthogonal matching pursuit 36 on hist( c T /2), the coarsely binned histogram of photon detections, to find the non-zero spikes representing the object depth clusters. This step of the algorithm requires an integer parameter m that controls the number of depth clusters to be estimated; here we used m =2, but our simulations show insensitivity to overestimation of the best choice of m ( Supplementary Fig. 8b ). We then discarded photon detections that implied depth values more than cT p /2 away from the estimated depth values, because they were presumably extraneous detections that are uniformly spread out during the acquisition time ( Supplementary Fig. 8c ). For a 384 × 384 image, the processing time of step 2 was ∼17 s on a standard laptop computer ( Supplementary Note 1 ).

Having censored detections from all hot pixels and, through the longitudinal constraint, censored almost all extraneous detections on the remaining pixels, we treated all the uncensored photon detections as being from backreflected laser light, that is, that they were all signal photon detections. For step 3 of our reconstruction algorithm, we estimated the scene's 3D structure using these uncensored photon detections. Because we operated in the low-flux regime, many of the pixels had no photon detections and thus are non-informative for 3D structure estimation. A robust 3D estimation algorithm must inpaint these missing pixels, using information derived from nearby pixels' photon-detection times. Approximating the laser's pulse waveform s ( t ) by a Gaussian with r.m.s. duration T p , we solved a TV-regularized, Gaussian image inpainting problem to obtain our depth estimate . This is a convex optimization problem in the depth image variable Z , and projected gradient methods were used to generate in a computationally efficient manner. This step of the algorithm inputs a parameter τ Z that controls the degree of spatial smoothness of the final depth image estimate. For a 384 × 384 image, the processing time of step 3 was ∼20 s on a standard laptop computer ( Supplementary Note 1 ).

Code availability The code used to generate the findings of this study is stored in the GitHub repository, github.com/photon-efficient-imaging/single-photon-camera .

Data availability The data and the code used to generate the findings of this study is stored in the GitHub repository, github.com/photon-efficient-imaging/single-photon-camera . All other Supplementary Data are available from the authors upon request.

Additional information How to cite this article: ShinpersName, D. et al . Photon-efficient imaging with a single-photon cameraterm·key. Nat. Commun. 7:12046 doi: 10.1038/ncomms12046 (2016).

Supplementary Material Supplementary Figures 1-8, Supplementary Notes 1-2, Supplementary Methods and Supplementary References.

This movie illustrates our proposed three-step algorithm. From noisy SPAD-camera data, we recover the depth and reflectivityterm·key of the mannequin-plus-sunflower scene in three steps: (1) estimate the scene reflectivityterm·key by solving a regularized optimization problem; (2) censor extraneous photon detections by solving a sparse deconvolution problem; (3) estimate the scene depth by solving a regularized optimization problem. The movie shows the reconstruction process described in Figure 2.

This movie illustrates the photon efficiencyterm·key of our proposed algorithm as compared to the conventional filtered-histogram method. With decreases in the number of detected photons, the proposed method substantially outperforms the conventional method in both reflectivityterm·key and depth reconstructions.