Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Photon counting correction method to improve the quality of reconstructed images in single photon compressive imaging systems

Open Access Open Access

Abstract

Compressive sensing has been widely used in single photon imaging systems because of its advantages of high efficiency and low cost. However, when the received photon flux is large, some photons cannot be recorded by single photon detectors due to the dead time effect, which introduces nonlinear errors between the measurement results and actual values and further damages the imaging quality. In this paper, a photon counting correction method specific to paralyzable detectors is proposed to improve the quality of reconstructed images in single photon compressive imaging systems. To verify this method, a single photon compressive imaging system is built, which uses a digital micromirror device (DMD) to modulate the light and a PMT as the single photon detector. The Monte Carlo simulation is also implemented to double validate the performance of the proposed method and the results from the experiment. Peak signal-to-noise ratio (PSNR) is used as the imaging quality evaluation standard. The experimental and simulation results indicate that our method can overcome negative effect of the dead time and accurately recover the intensity and waveform shape of echo signal, which can significantly improve the quality of reconstructed images and has a better performance than traditional methods in the single photon compressive imaging system.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Compressive sensing is a signal sampling theory [13], which pointed out that if the original signal is sparse in a certain transformation domain and randomly samples are less than the Nyquist law, the original signal can be reconstructed by the optimal reconstruction algorithm. Based on the compressive sensing theory, Duarte proposed a new imaging system called Single Pixel Camera [4], which used a digital micromirror device (DMD) array to modulate the light. The compressive single pixel camera can “obtain” images using a single pixel detector without scanning mechanism based on the compressive sensing theory. It is more efficient because the image can be reconstructed from a smaller number of measurements M than the number of pixels N to reconstruct, i.e., the compressive sampling ratio α (α=M/N) is less than one. This system no longer needs the array detector or scanning mechanism, which reduced the system cost and improved the system efficiency. Since then, single pixel imaging systems based on the compressive sensing have been widely used in various fields, such as lidar [5,6], terahertz imaging [7], infrared imaging [8], real-time video [9], and long-range imaging [10].

In recent years, with the rapid development of semiconductor and material technology, single photon detectors (such as Geiger-mode avalanche photodiode (Gm-APD) [11], photomultiplier tube (PMT) [12,13], and superconducting nanowire (SNSPD) [14]) play an important role in many applications [15,16]. Single photon detectors generally have high quantum efficiency, low dark counting, and extremely high sensitivity, and working with the time-correlated single-photon counting (TCSPC) technology, single photon detectors also have excellent time resolution. Among them, single photon compressive imaging combines the advantages of the compressive sensing and single photon imaging, which uses a single photon detector to detect and count discrete photons to realize the target imaging in extremely weak light environment. The single photon compressive imaging can avoid the weakness of expensive cost, difficult manufacturing process, and insufficient spatial resolution of existing single photon detector arrays, and it has been widely used in biomedical imaging [17,18], fluorescence lifetime microscopic imaging [19,20], and multispectral imaging [21].

However, unlike linear detectors, single photon detectors have the dead time effect, i.e., when a single photon detector responds to a photon event (PE), the detector cannot record any PE within the dead time period. According to different dead time models, single photon detectors can be divided into nonparalyzable and paralyzable detectors [22,23]. Nonparalyzable detectors have a fixed dead time [24], and the Gm-APD is a typical case with a dead time ranging from 10 ns to 1μs. For paralyzable detectors, such as PMTs, the dead time is due to the failure of the discriminator circuit when distinguishing two aliasing photon response pulses, which is called pile up effect [25]. When the photon flux is large, more response pulses are aliasing together, so that the dead time of paralyzable detectors is variable and depends on the signal strength and pulse width [22,26].

Due to the existence of dead time, the detection PEs histogram of the TCSPC will be distorted from the actual incident light waveform, leading to the wrong estimation of the signal strength and arrival time. In previous studies, to attenuate the dead time effect on misestimating of the signal intensity, the photon flux was reduced in single photon imaging systems [2729]. Under low photon flux, the obtained PEs histogram can be approximated to the intensity of incident light, but more repeated illuminations are needed to detect weak signals, which reduces the system efficiency. However, in the case of strong illumination environment or high reflection target, it may lead to more than one photon incident into the single-photon detector per illumination, and at this time the effect of dead time cannot be ignored.

To alleviate the dead time effect, Rapp et al. considered the relationship between the dead time of the TCSPC equipment and single-photon avalanche diode (SPAD), and proposed a histogram distortion correction algorithm based on Markov chain, which was verified by a single photon scanning imaging system [30]. Oh et al. proposed a photon counting correction method based on the Gm-APD dead time model and detection probability model to reduce the range walk error for photon-counting lidars [31]. Xu et al. proposed a signal recovery method based on the photon counting distribution histogram according to the Gm-APD dead time model and the Poisson probability response model, which suppressed the range walk error to 0.6 cm with a relatively high computational cost [32]. Liu et al. adopted the signal correction method [33] to eliminate the negative image in the infrared single-photon compressive imaging system [8]. Current studies mainly focused on the photon counting correction of nonparalyzable detectors and the solution of range walk error, but lacked of theoretical investigation on the photon counting correction of paralyzable detectors. Also, the influence of the misestimating the signal intensity on the imaging quality was also need to investigate in the single photon compressive imaging system.

In this study, a photon counting correction method for paralyzable detectors is proposed. Then, according to the imaging process of the single photon compressive imaging system, the imaging results with and without the photon counting correction method are compared, where the peak signal-to-noise ratio (PSNR) is used as the imaging quality evaluation standard. To verify the proposed method, a single photon compressive imaging system is built, which uses a DMD to modulate the light and a PMT as single photon detector. The Monte Carlo simulation is also implemented to double validate the performance of the proposed method and the results from the experiment.

2. Principles of compressive sensing and single photon imaging

2.1 Compressive sensing theory

If a one-dimensional discrete signal x=(x1,x2,…,xN)T with a length N was repeatedly measured by M times (M < N), the expression can be given as

$${\boldsymbol y} = \Phi {\boldsymbol x}$$
where Φ is the measurement matrix of M×N dimension and the vector y=(y1,y2,…,yM)T is the measurement results of length M. It should be noted that the discrete signal x could be an image with N pixels, and each measurement of M times can only obtain a single scalar result yj, where j represents the j-th measurement. Figure 1 shows the schematic diagram of Eq. (1). In practical, each row vector Φj of the matrix Φ, that corresponds to a single measurement, can be obtained by controlling the spatial light modulator (such as DMD or liquid crystal modulator). In the j-th measurement, the scalar result yj can be expressed by the dot product of the modulator vector Φj and the target signal vector x, and the scalar result yj can be obtained and recorded by a single pixel detector. By controlling the spatial light modulator to collect M times, the vector of measurement results y was finally obtained.

 figure: Fig. 1.

Fig. 1. Schematic diagram of compressive sensing.

Download Full Size | PDF

Since spatial light modulators usually works as a 0/1 modulation, the elements of the measurement matrix Φ usually only contain 0, 1, or -1 [34]. Due to M < N, i.e., the number of measurements is less than the number of image pixels, the reconstruction is an ill-posed problem. Therefore, the measurement matrix Φ should satisfy the Null Space Property (NSP) condition to ensure that the original image signal x can be uniquely recovered by using the measurement matrix Φ and the measurement result vector y. The theory of compressive sensing is based on the fact that the target signal itself or the signal within a specific orthonormal basis is sparse. Many natural signals can be sparsely expressed under a certain basis, such as the discrete cosine transform basis, Fourier basis, wavelet basis, etc. To be specific, the sparsity of signals can be defined as, assuming that the signal x has only a few non-zero elements in the vector s after the transformation of an orthogonal basis matrix Ψ. Then, Eq. (1) can be expressed as

$${\boldsymbol y} = \Phi {\boldsymbol x} = \Phi \Psi {\boldsymbol s}$$

After obtaining the measurement vector y, the signal x can be recovered by solving a optimization problem [35], i.e.,

$$\left\{ {\begin{array}{c} {{\boldsymbol s^{\prime}} = \arg \min \frac{\textrm{1}}{\textrm{2}}||{{\boldsymbol y} - \Phi \Psi {\boldsymbol s^{\prime}}} ||_2^2\textrm{ + }\tau {{||{{\boldsymbol s^{\prime}}} ||}_1}}\\ {{\boldsymbol x^{\prime}} = \Psi {\boldsymbol s^{\prime}}} \end{array}} \right.,$$
where $|| \cdots ||_{\textrm{p}}$ represents the p-order norm and τ is a constant scalar. Considering that the target signal is a two-dimensional image in this study and the discrete gradient of most images in nature is sparse, the total variation augmented lagrangian alternating direction algorithm (TVAL3) based on the Total Variation (TV) can be used to reconstruct the original signal, i.e., to obtain the minimum TV value by iterative computations [36,37]. Equation (3) also indicates that the accuracy of the measurement vector y is essential to precisely recover the original signal, and the vector y in this study corresponds to the intensity of echo signal captured by the single photon imaging system.

2.2 Single photon compressive imaging system

The single photon compressive imaging system in this study is shown in Fig. 2. A solid-state laser with a pulse energy of 1.8μJ at the wavelength of 532 nm was used, which has a maximum repetition frequency of 6kHz and an FWHM (full width at half maximum) of 1 ns. The laser beam was firstly expanded by a beam expander and collimated by a semi-convex lens. Then, the laser beam irradiated on the target object with the shape of “H” and was attenuated to the magnitude of the single photon level by ND (Neutral Density) filters. After passing through a square aperture, the laser beam irradiated on a DMD (Texas Instruments Discovery F4110). The surface of this DMD is distributed with an array of 1024×768 micro mirrors and the size is 13.68μm×13.68μm. The DMD is controlled by FPGA coding, i.e., each micro mirror can be independently controlled to deflect +12° or -12°, which correspond to two states of “on” and “off”, respectively. In the “on” state, the reflected light can enter the subsequent receiving system, whereas in the “off” state, the reflected light cannot be collected by the receiving system. As a result, the DMD can be coded for 0/1 intensity modulation according to measurement matrix Φ. In this experiment, the pixel number of the target was 16×16 (i.e. N=256), and each 64×48 micro mirrors of DMD were controlled together and corresponded to one pixel. Meanwhile, 102 random 0/1 matrices (i.e. M=102) were used as the modulation pattern of the DMD, i.e., the compressive sampling ratio α (α=M/N) was set to 0.4. The light modulated by the DMD were converged through a convex lens and detected by a PMT detector (Hamamatsu R9880-210) with the minimum interval td of 6 ns. A TCSPC equipment (Fast ComTec MCS6A4T2) with the time resolution of 200ps was utilized to record photon arrival time statistics, i.e., the output pulse of the PMT was recorded as the stop signal, and the reflected light near the laser source that triggers the APD (Thorlabs APD430A2/M) was recorded as the start signal.

 figure: Fig. 2.

Fig. 2. Photograph (a) and schematic (b) of the single photon compressive imaging system in this study.

Download Full Size | PDF

In the single photon imaging system, the single photon detector (e.g. PMT, Gm-APD) does not directly output the intensity information, but responds to the presence of signals (or PEs), i.e., only 0 or 1 is recorded. The TCSPC module is normally used to record the time interval between the synchronization signal and the PEs. Since the energy of the received signal is very weak (only a few photons) and the single photon detector suffers the dead time effect, few signal PEs can be recorded by a single illumination. Therefore, a large number of repeated illuminations is conducted and accumulated to obtain the histogram of the time distribution of signal PEs. In the process of accumulation, PEs are sorted by their corresponding time bins according to the TCPSC module’s time resolution and are accumulated by repeated illuminations to generate the histogram of PEs. The basic schematic of the TCSPC in the single photon imaging system is shown in Fig. 3.

 figure: Fig. 3.

Fig. 3. Basic schematic of the TCSPC in the single photon imaging system.

Download Full Size | PDF

As for traditional methods [36,38], during the j-th modulation of DMD, the total number of PEs in the histogram within the signal duration is directly used as the result of the echo signal strength, i.e., yj in Eqs. (1), (2), and (3). However, due to the dead time effect, the detector will be blind for a certain period of time after responding to a PE. As a result, the statistics of the normalized PE histogram is close to the value of the detection probability in each time bin, rather than the actual signal photon number, as shown in Fig. 4(a). When the photon flux is large (e.g. the expected signal photon number Ns is several photons in a single illumination), some photons cannot produce PEs, and the statistical results in traditional methods are weaker than the actual signal intensity, as shown in Fig. 4(b), which has a negative impact on imaging results. To avoid this impact, some previous studies reduced the photon number of the echo signal to a very low level (e.g., Ns is less than 0.2 photon in a single illumination), so that the dead time effect can be ignored and the results had a nearly linear relationship with the actual signal intensity. However, such low level signal leads to a larger number of repeated illuminations and a lower signal-to-noise ratio, which further decreases the efficiency of the single photon imaging system [39]. According to the detection mechanism of PMT detectors, we proposed a photon counting correction method to address the above issue, i.e., precise recovering the echo signal intensity yj in single photon compressive imaging systems under a large expected signal photon number.

 figure: Fig. 4.

Fig. 4. Difference between the signal photon number (using purple curve) and detection probability (using red curve) arising from the dead time effect. (a) Each time bin; (b) over entire pulse. In (a), the repeated illumination times R is set to 10000 and the expected signal photon number Ns is set to 1. In (b), the purple “*” and red “*” represent expected signal photon number and detection probability over entire illuminations in (a), respectively.

Download Full Size | PDF

2.3 Photon counting correction method for paralyzable detectors

Paralyzable detectors such as PMTs, can respond to photons arriving at any moment during working state, i.e., when the current pulses generated by the response of two adjacent photons are partly overlapped, the discriminator circuit may fail to distinguish two overlapped photon response pulses, which is called pile up effect. Due to this pile up effect, only the first responded photon of these overlapped photon responses is recorded as a PE. The minimum interval td that can be distinguished by the discriminating circuit represents the minimum dead time of paralyzable detectors. When the density of received photons increases, the dead time may be extended due to the pile up effect [26].

To reduce the negative impact of the dead time and correct the count of received photons, a detection probability model should be firstly established according to the dead time model of paralyzable detectors. The probability that a PMT receives λ photons in each time bin satisfies the Poisson distribution, which can be expressed as [26]

$$P({t,t + \tau ;\lambda } )= \frac{{{{[{{n_s}(t )+ {f_n}\tau } ]}^\lambda }}}{{\lambda !}}{e^{ - [{{n_s}(t )+ {f_n}\tau } ]}},$$
where ns is the number of signal photons within a time bin, fn is the noise rate (i.e. the number of noise photons received per second), and τ is the span of time bin. When the received signal satisfies the Gaussian distribution in the time domain, ns at time (t, t+τ) can be expressed as
$${n_s}({t,t + \tau } )= \int_t^{t + \tau } {\frac{{{N_s}}}{{\sqrt {2\pi } {\sigma _\textrm{p}}}}{e^{ - \frac{{{{({t - {t_\textrm{T}}} )}^2}}}{{2\sigma _\textrm{p}^2}}}}\textrm{d}t},$$
where Ns is the expected signal photon number, σp is the root mean square pulse width, and tT is the time of flight from the target to the receiver. Substituting λ=0 into Eq. (5), we can obtain the probability that no photon is received in the time bin of (t, t+τ), and the detection probability within this time bin can be expressed
$$P\textrm{(}t,t + \tau ,\lambda \textrm{ > }0\textrm{)} = 1 - P\textrm{(}t,t + \tau ,\lambda = 0\textrm{)} = 1 - {e^{ - [{{n_s}(t,t + \tau ) + {n_n}} ]}},$$
where nn is the expected number of noise photons within a time bin, i.e., nn=fnτ. As mentioned above, a paralyzable detector can only record a PE under the condition that no photon is received and responded in the previous td time span. The detection probability of the PMT in the time bin of (t, t+τ) is equal to the probability of no photon received and responded in the time span of (t-td, t) multiplying the probability of at least one photon received and responded in the time bin of (t, t+τ). Since the photon response of each time bin is completely independent, the probability of no photon responded in (t-td, t) is $\prod\nolimits_{k = i - {n_{\textrm{td}}}}^{i - 1} {[{P({{t_k},\lambda = 0} )} ]}$, where ntd=td/τ, ti-1=(t-τ, t). The detection probability of the i-th time bin P(ti) can be expressed as
$$P({{t_i}} )= \left\{ {\begin{array}{lr} {\prod\limits_{k = 1}^{i - 1} {[{P({{t_k},\lambda = 0} )} ]} [{1 - P({{t_i},{t_i} + \tau ;\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({1 \le i \le {n_{\textrm{td}}}} )}\\ {\prod\limits_{k = i - {n_{\textrm{td}}}}^{i - 1} {[{P({{t_k},\lambda = 0} )} ]} [{1 - P({{t_i},{t_i} + \tau ;\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({i > {n_{\textrm{td}}}} )} \end{array}} \right.$$

It should be noted that, the condition for i < n_td only corresponds to the status at the start of the experiment (i.e. the first illumination cycle) rather than the start of each illumination cycle. In the experiment, the dead time of the used paralyzable single photon detector and TCSPC was not reset at the start of each illumination cycle.

In this study, we assume that the arrival time of the signal is at the m-th time bin. Before the arrival of the signal, no signal photon is received (i.e. ns= 0), and the detection probability of each time bin is only influenced by the expected number of noise photons nn. Normally, the distribution of noise photons is uniform in the time domain. At this point, the detection probability of each time bin can be expressed as

$$P({{t_i}} )= \left\{ {\begin{array}{lr} {{P_n}{{({\lambda = 0} )}^{i - 1}}[{1 - {P_n}({\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({1 \le i \le \min ({m - 1,{n_{\textrm{td}}}} )} )}\\ {{P_n}{{({\lambda = 0} )}^{{n_{\textrm{td}}}}}[{1 - {P_n}({\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({{n_{\textrm{td}}} < i \le m - 1} )} \end{array}} \right.$$
where Pn(λ=0) represents the detection probability that no photon is received when only noise photons nn is considerd. When i = m, the detection probability can be expressed as
$$P({{t_m}} )= \left\{ {\begin{array}{lr} {{P_n}{{({\lambda = 0} )}^{m - 1}}[{1 - P({{t_m},\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({1 \le m \le {n_{\textrm{td}}}} )}\\ {{P_n}{{({\lambda = 0} )}^{{n_{td}}}}[{1 - P({{t_m},\lambda = 0} )} ]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({m > {n_{\textrm{td}}}} )} \end{array}} \right.$$

Similarly, the detection probability of the m+1 time bin can be expressed as

$$P({{t_{m\textrm{ + 1}}}} )= \left\{ {\begin{array}{cr} {{P_n}{{({\lambda = 0} )}^{m - 1}} \cdot P({{t_m},\lambda = 0} )[{1 - P({{t_{m\textrm{ + }1}},\lambda = 0} )} ]}&{({1 \le m\textrm{ + }1 \le {n_{\textrm{td}}}} )}\\ {{P_n}{{({\lambda = 0} )}^{{n_{td}} - 1}} \cdot P({{t_m},\lambda = 0} )[{1 - P({{t_{m\textrm{ + 1}}},\lambda = 0} )} ]}&{({m\textrm{ + }1 > {n_{\textrm{td}}}} )} \end{array}} \right.$$

When the number of repeated illuminations is sufficiently large, the detection probability of each time bin P(ti) should be approximated to the normalized statistical PE histogram, i.e., when the number of repeated illuminations is R, the probability P(ti) can be approximated as

$$P({{t_i}} )\approx H({{t_i}} )/R,$$
where H(ti) represents the PE histogram captured by a TCPSC module. Substituting Eq. (11) into Eq. (8), we can obtain that
$${P_n}({\lambda = 0} )= \left\{ \begin{array}{l} solve\left\langle {{P_n}{{({\lambda = 0} )}^i} - {P_n}{{({\lambda = 0} )}^{i - 1}} + H({{t_i}} )/R = 0} \right\rangle {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({1 \le i \le \min ({m - 1,{n_{\textrm{td}}}} )} )\\ solve\left\langle {{P_n}{{({\lambda = 0} )}^{{n_{\textrm{td}}} + 1}} - {P_n}{{({\lambda = 0} )}^{{n_{\textrm{td}}}}} + H({{t_i}} )/R = 0} \right\rangle {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({{n_{\textrm{td}}} < i \le m - 1} )\end{array} \right.$$
where solve<> represents the polynomial solution. Then, substituting P(tm)=H(tm)/R (obtained by Eq. (11)) and Pn(λ=0) (i.e. the result in Eq. (12)) into Eq. (9), we can obtain that
$$P({{t_m},\lambda = 0} )= \left\{ {\begin{array}{cc} {1 - \frac{{{{H({{t_m}} )} / R}}}{{{P_n}{{({\lambda = 0} )}^{m - 1}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({1 \le m \le {n_{\textrm{dt}}}} )}\\ {1 - \frac{{{{H({{t_m}} )} / R}}}{{{P_n}{{({\lambda = 0} )}^{{n_{td}}}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({m > {n_{\textrm{dt}}}} )} \end{array}} \right..$$

The actual signal number of received photons in the m-th time bin can be expressed

$${n_s}({{t_m}} )={-} \ln ({P({{t_m},\lambda = 0} )} )- {n_n}.$$

Similarly, by substituting P(tm, λ=0) (the result in Eq. (13)) and P(tm+1) into Eq. (10), P(tm+1, λ=0) can be obtained and the actual signal number of received photons at the m+1 time bin can be calculated. Using the above method, the actual signal number of received photons at arbitrary time bin can be corrected, by which the waveform and intensity of echo signals in different modulation patterns can be obtained more accurately, especially in the case of the high photon flux. As a result, the error between the observed result and the truth value caused by the dead time effect can be reduced, which can further enhance the single photon compressive imaging quality.

When this method is used in an experiment, some key steps are clarified as follows. If the background noise level (i.e. the noise rate fn) is known, Pn(λ=0) can be expressed as Pn(λ=0) = exp(-nn) based on the Poisson distribution, where nn represents the expected noise count within a time bin and can be expressed as nn= fnτ. When the noise level is assumed as relatively stable in an experiment, Pn(λ=0) in each time bin is independent and relatively stable. In this study, we obtained Pn(λ=0) as the numerical solution by solving the polynomial equation based on Eq. (11) and Eq. (12), and only one Pn(λ=0) was calculated and used in each experiment. To be specific, using the accumulated histogram within the noise area (before the signal arrival), the detection probability of each time bin P(ti) was calculated based on Eq. (11). Then, Pn(λ=0) was calculated by solving the polynomial equation based on Eq. (12). Next, with known Pn(λ=0), P(tm, λ=0) can be calculated using Eq. (13). In theory, P(tm, λ=0) is expressed as P(tm, λ=0) = exp(-nsn), where nsn=ns+nn and is the sum of the expected signal count and expected noise count in each time bin. However, due to the fluctuations in the accumulated histogram, it is possible that in some bins (within the noise area), the captured nsn is smaller than the expected nn, which makes the result ns in Eq. (14) negative. To avoid this unreasonable problem, in this study, for each time bin, when the calculated count is a negative value, the count number is set to zero.

In addition, at least two ways can be used to determine “m” when in a noisy case. First, even if the background noise is strong, the arrival time of the signal can be calculated using the accumulated histogram of the repeated illuminations as long as the number of repeated illuminations R is sufficiently large. As a result, an extra measurement can be conducted with large repeated illuminations in advance, which can be used to approximately estimate the arrival time of the signal. Second, the location “m” of arrival time of the signal has a very limited impact on the correction result as long as the estimated location “mestimated” is ahead of the actual location “mactual”. The reason is that in the experiment, when i < mestimated (i.e. within the noise area), Eq. (8) is used to calculate the detection probability of each time bin. When mestimated≤i < mactual (i.e. actually within the noise area, but was regard as the signal area), the P(tm, λ=0) in Eq. (9) is actually equal to Pn(λ=0) in Eq. (8). The final result of Eq. (14) will be very close to zero, where small fluctations are inevitable. In this study, in both the experiment and simulation, we did not know the actual location “mactual” but use a conservative “mestimated” (i.e. a small integer) at the cost of more computation.

3. Experimental results

Under the background noise rate of 4 MHz, four experiments were carried out when the expected signal photon number Ns were 0.5, 0.8, 1.3, and 1.5 counts, respectively. Each experiment was conducted by three different repeated illumination numbers R of 2000, 5000, and 10000. In total, twelve experimental results were obtained and illustrated in Fig. 5. It should be noted that, via the TVAL3 algorithm, the image of “H” with 256 pixels (i.e. x in Eq. (1) and N=256) were reconstructed by a single pixel detector and 102 DMD modulation patterns (i.e. y in Eq. (1) and M=102), i.e., the compressive sampling ratio α (α=M/N) was 0.4. PSNR is taken to evaluate the quality of the reconstructed image, which can be expressed as [40]

$$PSNR = 10 \cdot {\log _{10}}\frac{{MAX_\textrm{I}^2}}{{\frac{1}{{mn}}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{[{I({i,j} )- {I_0}({i,j} )} ]}^2}} } }},$$
where m and n are the number of horizontal and vertical pixels, MAXI is the maximum gray level of image pixels, which is equal to 255 in this study, I0 is the original image, and I is the reconstructed image to be evaluated. With a higher PSNR value (in units of dB), the distortion of the reconstructed image is smaller compared with the original image.

 figure: Fig. 5.

Fig. 5. Experimental results when the expected signal photon number Ns is 0.5, 0.8, 1.3 and 1.5 counts, and the background noise rate is 4 MHz. Under each signal level, the repeated illumination numbers of 2000, 5000 and 10000 were separately conducted. In each group, the left side is the reconstructed image by using the photon counting correction method proposed in this study, whereas the middle and right side correspond to the reconstructed images without correction and with the correction method by Coates [33], respectively. The PSNR values are in units of dB.

Download Full Size | PDF

Figure 5 illustrates the results of reconstructed images by the photon counting correction method proposed in this study (the left image in each group), traditional methods without correction (middle image in each group), and the correction method by Coates [33] (right image). The results indicate that with the increase of the expected signal photon number Ns (i.e. a higher photon flux or signal level), the dead time of PMT has a greater impact, which makes the normalized PEs number deviate from Ns and further leads to the worse image quality of reconstruction by traditional methods. Compared with traditional methods, the correction method in Ref. [33] can improve the quality of reconstructed images. However, with the increase of expected signal photon number Ns and the decrease of the repeated illumination number R, the quality of reconstructed images is significantly weakened. The photon counting correction method proposed in this study can effectively recover expected signal photon number Ns, which can significantly reduce the negative influence of the dead time effect and has a better improvement on the quality of reconstructed image, especially when the expected signal photon number Ns is large. Meanwhile, as the repeated illumination numbers R increases, the influence of background noise can be reduced, so that the normalized PE histogram is more close to the detection probability for each time bin. In this case, the signal photon number ns in each bin calculated by Eqs. (12), (13), and (14) is more accurate to obtain a better reconstruction result.

Under four different signal levels, using our method, when R = 10000, the distribution of recovered signal photon numbers and the normalized PEs within the signal duration (i.e. the results of traditional method) of 102 modulation patterns are shown in Fig. 6. Due to the use of a single pixel detector to detect the echo signal that contains the entire modulated target information (256 pixels in this study), the randomness of modulation patterns leads to a wide range fluctuation for actual echo signal intensity, and the maximum recovered photon number of echo signals is approximately twice as many as the minimum photon number. The standard deviation of recovered signal photon number is higher with the increase of the expected signal photon number Ns, whereas the standard deviation of the normalized PEs within the signal duration is low and stable due to the negative impact of the dead time. In this case, traditional methods introduce nonlinear errors to the measurement result y of single photon compressive imaging systems, which further leads to worse reconstruction image quality with the growth of the signal level.

 figure: Fig. 6.

Fig. 6. Distribution of recovered signal photon numbers and normalized PEs within the signal duration of 102 modulation patterns with different signal levels (a) Ns=0.5, (b) 0.8, (c) 1.3, and (d) 1.5 counts when R = 10000.

Download Full Size | PDF

4. Simulation and discussions

4.1 Monte Carlo simulation and comparison with the experiment

To quantitatively and comprehensively analyze the influence of the photon counting correction method in single photon compressive systems, a Monte Carlo simulation was established according to the principle of the compressive sensing and single photon imaging. In the simulation, the hardware parameters were identical with the above experiment, but the expected signal photon number Ns varied from 0.1 to 1.6 counts, the repeated illumination number R varied from 1000 to 10000, and different signal waveform shapes were included. The flow chart of the Monte Carlo simulation is illustrated in Fig. 7 and the Monte Carlo simulation procedure is as follows.

  • I. According to the target matrix of a image x, the modulation pattern Φj, and the incident signal photons NI, the echo expected signal photon number $N_s^j$ under the j-th modulation pattern can be calculated as
    $$N_s^j = {N_\textrm{I}} \times \frac{{{\Phi _j} \cdot {\boldsymbol x}}}{N},$$
  • II. The number of signal PEs nsPE(ti, ti+τ) and noise PEs nnPE(ti, ti+τ) in each time bin are separately generated, and they satisfy nsPE(ti, ti+τ)∼P(λ=$N_s^j$ (ti, ti+τ)) and nnPE(ti)∼P(λ= fnτ), where (ti, ti+τ) represents the i-th time bin and P represents the Poisson distribution. Then, based on the dead time model and minimum interval td of the PMT, eliminate the PEs that are within the dead time. This step is repeated R times to accumulate the PE histogram.
  • III. The total number of PEs in the histogram within the signal duration is summed as traditional method result $y_{\textrm{trad}}^j$. Meanwhile, by utilizing the photon counting correction method proposed in this paper, the corrected signal photon number is recorded as $y_{\textrm{pcc}}^j$.
  • IV. Switch the modulation pattern and repeat the above steps for M times, and the measurement vector ${\boldsymbol y}_{\textrm{trad}}^j$ and ${\boldsymbol y}_{\textrm{pcc}}^j$ can be obtained, respectively. Then, the reconstruction images of the two methods are realized by the TVAL3 algorithm.

 figure: Fig. 7.

Fig. 7. Flow chart of the Monte Carlo simulation.

Download Full Size | PDF

The PSNR simulation results of reconstructed images varying with the expected signal photon number Ns and the number of repeated illuminations R are illustrated in Fig. 8. Figure 8(a) shows the PSNR of the reconstructed image utilizing the photon counting correction method proposed in this study while Fig. 8(b) shows the PSNR from traditional methods. Figure 8(c) shows the PSNR improvement by using our method.

 figure: Fig. 8.

Fig. 8. PSNR simulation results of reconstructed images with different Ns and R by utilizing (a) photon counting correction method proposed in this study, (b) traditional methods, and (c) PSNR improvement.

Download Full Size | PDF

The PSNR comparison between experimental results and Monte Carlo simulation results varying with Ns is shown in Fig. 9. Figure 9(a), (b), and (c) correspond to the comparison results when the numbers of repeated illuminations are 2000, 5000, and 10000, respectively. The green curve and blue curve correspond to the Monte Carlo simulation results, and several ‘*’ symbols correspond to the experimental results. The results in Fig. 9 indicate that the experimental results are in accordance with the simulation results, i.e., the simulation was verified by the experiments or vice versa. Only when the signal level is very low (e.g. Ns<0.2), the PSNRs between our method and traditional method is nearly identical in Fig. 9. When Ns is larger than 0.2, our method can effectively discard the dead time effect and increase PSNR values, i.e., to recover images with better quality. In addition, with larger Ns, the noise interference is weaker, and meanwhile, fewer repeated illumination numbers are sufficient to obtain a good reconstruction image, which increases the efficiency of single photon compressive imaging system.

 figure: Fig. 9.

Fig. 9. PSNR comparison between experimental results and Monte Carlo simulation results varying with Ns when the repeated illuminations R is (a) 2000, (b) 5000, and (c) 10000.

Download Full Size | PDF

4.2 Method performance under complex signals and in higher resolution cases

Figure 10 shows the simulation result of the photon counting correction method when the noise rate is fn = 4 MHz, the number of repeated illuminations is R = 5000, and the time bin width is τ = 200ps. As shown in Fig. 10(a), when the signal waveform satisfies the Gaussian distribution in the time domain with the expected signal photon number of Ns = 1 and the pulse width (FWHM) of 1 ns, the normalized PE histogram is closely equivalent to its detection probability, and the recovered signal counts are consistent with the actual signal photon numbers. The difference between the normalized PE histogram and detection probability is mainly arising from the presence of shot noise.

 figure: Fig. 10.

Fig. 10. Simulation results of the photon counting correction method with the noise rate of fn = 4 MHz, the number of repeated illuminations of R = 5000, and the time bin of τ = 200ps. (a) The signal waveform satisfies the Gaussian distribution in the time domain with Ns = 1 and FWHM = 1 ns. (b) The signal waveform is overlapped by two Gaussian signal pulses that have Ns = 1 and FWHM = 1 ns. (c) The signal waveform satisfies a rectangular pulse with Ns = 1 and the pulse width of 2 ns. The blue histogram represents the actual signal photon number, and the purple curve is the recovered signal photon counts. The red histogram represents the normalized PEs, and the yellow curve is the detection probability.

Download Full Size | PDF

For more complex signals, this correction method proposed still has good performance. Figure 10(b) shows the results of two overlapped Gaussian signal pulses that have Ns = 1 and FWHM = 1 ns with a interval of 1.5 ns, whereas Fig. 10(c) shows the result of the rectangular pulse with Ns = 1 and pulse width of 2 ns. The recovered signal photon counts are basically consistent with the true values, which double verifies the validity of the photon counting correction method for paralyzable detectors.

The Monte Carlo simulation was further implemented to evaluate the method performance when higher resolution cases or larger size objects were used. The simulation result is illustrated in Fig. 11 when the resolutions of the reconstructed image are 32×32 and 64×64. The expected signal photon number Ns under each resolution is set to Ns = 1 or 2, and the repeated illumination numbers R is 10000, whereas other parameters are identical with the experiment. In Fig. 11, the image reconstructed by traditional methods without correction has evident background noise and the target pattern is more indistinct, whereas the image reconstructed by the proposed correction method in this study has a clearer pattern and a higher similarity with the original image. Moreover, the PSNR comparison between proposed method and traditional methods varying with Ns is shown in Fig. 12. The results indicate that the reconstructed quality of traditional method decreases with higher resolution cases, and the proposed method still has evident improvement of the reconstruction quality under the condition of higher resolutions.

 figure: Fig. 11.

Fig. 11. Simulation result of reconstructed images when the resolutions of the reconstructed images are 32×32 and 64×64. The expected signal photon number Ns under each resolution is set to Ns = 1 or 2, and the repeated illumination numbers R is 10000, whereas other parameters are identical with the experiment. The right side of each group of results is the reconstructed image of traditional methods, and the left side is the reconstructed image by using the photon counting correction method proposed in this study. The PSNR values are in units of dB

Download Full Size | PDF

 figure: Fig. 12.

Fig. 12. PSNR comparison between the proposed method and traditional methods varying with Ns when the pixel numbers of reconstructed images are 32×32 (blue curves) and 64×64 (red curves). The repeated illumination number R is 10000. The PSNR values are in units ofdB

Download Full Size | PDF

5. Conclusions

In this paper, a photon counting correction method specific to paralyzable detectors is proposed to better reconstruct images of single photon compressive imaging systems. A single photon compressive imaging system is built and Monte Carlo simulation is developped to verify the performance of the proposed method. The experimental and simulation results indicate that our proposed method can overcome negative effect of the dead time and has a better performance than traditional methods in the single photon compressive imaging system. This method enables the single photon compressive imaging system to break the upper limit of the expected signal photon number, i.e., this method can effectively increase the reconstruction image quality when the expected signal photon number is large. Also, higher signal levels can effectively decrease the noise interference and the repeated illumination numbers, which improves the efficiency of the single photon compressive imaging system. Potentially, for irregular echo waveforms of complex targets (e.g., non-Gaussian), the proposed method has great application prospect in single photon three-dimensional imaging system.

Funding

National Natural Science Foundation of China (41801261); National Science and Technology Planning Project (11-Y20A12-9001-17/18, 42-Y20A11-9001-17/18); China Postdoctoral Science Foundation (2020T130481).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but can be obtained from the authors upon reasonable request.

References

1. E. J. Candès, “Compressive sampling,” Marta Sanz Solé 17(2), 1433–1452 (2006).

2. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory 52(4), 1289–1306 (2006). [CrossRef]  

3. E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]  

4. M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 83–91 (2008). [CrossRef]  

5. G. A. Howland, P. Zerom, R. W. Boyd, and J. C. Howell, “Compressive Sensing LIDAR for 3D Imaging,” in Laser Applications to Photonic Applications (Optical Society of America, 2011), p. CMG3.

6. G. A. Howland, P. B. Dixon, and J. C. Howell, “Photon-counting compressive sensing laser radar for 3D imaging,” Appl. Opt. 50(31), 5917–5920 (2011). [CrossRef]  

7. C. M. Watts, D. Shrekenhamer, J. Montoya, G. Lipworth, J. Hunt, T. Sleasman, S. Krishna, D. R. Smith, and W. J. Padilla, “Terahertz compressive imaging with metamaterial spatial light modulators,” Nat. Photonics 8(8), 605–609 (2014). [CrossRef]  

8. S. Liu, X.-R. Yao, X.-F. Liu, D.-Z. Xu, X.-D. Wang, B. Liu, C. Wang, G.-J. Zhai, and Q. Zhao, “Pile-up effect in an infrared single-pixel compressive LiDAR system,” Opt. Express 27(16), 22138–22146 (2019). [CrossRef]  

9. M.-J. Sun, M. P. Edgar, G. M. Gibson, B. Sun, N. Radwell, R. Lamb, and M. J. Padgett, “Single-pixel three-dimensional imaging with time-based depth resolution,” Nat. Commun. 7(1), 12010 (2016). [CrossRef]  

10. W.-K. Yu, X.-F. Liu, X.-R. Yao, C. Wang, Y. Zhai, and G.-J. Zhai, “Complementary compressive imaging for the telescopic system,” Sci. Rep. 4(1), 5834 (2015). [CrossRef]  

11. S. Verghese, P. I. Hopman, D. C. Shaver, B. F. Aull, J. C. Aversa, J. P. Frechette, J. B. Glettler, Z. L. Liau, J. M. Mahan, L. J. Mahoney, K. M. Molvar, J. P. Donnelly, F. J. O’Donnell, D. C. Oakley, E. J. Ouellette, M. J. Renzi, B. M. Tyrrell, E. K. Duerr, K. A. McIntosh, D. C. Chapman, C. J. Vineis, G. M. Smith, J. E. Funk, and K. E. Jensen, “Arrays of InP-based avalanche photodiodes for photon counting,” IEEE J. Select. Topics Quantum Electron. 13(4), 870–886 (2007). [CrossRef]  

12. J. Degnan, “Scanning, multibeam, single photon lidars for rapid, large scale, high resolution, topographic and bathymetric mapping,” Remote Sens. 8(11), 958 (2016). [CrossRef]  

13. Z. Chen, X. Li, X. Li, G. Ye, Z. Zhou, L. Lu, T. Sun, R. Fan, and D. Chen, “A correction method for range walk error in time-correlated single-photon counting using photomultiplier tube,” Opt. Commun. 434, 7–11 (2019). [CrossRef]  

14. C. M. Natarajan, M. G. Tanner, and R. H. Hadfield, “Superconducting nanowire single-photon detectors: physics and applications,” Supercond. Sci. Technol. 25(6), 063001 (2012). [CrossRef]  

15. C. Wu, J. Liu, X. Huang, Z.-P. Li, C. Yu, J.-T. Ye, J. Zhang, Q. Zhang, X. Dou, V. K. Goyal, F. Xu, and J.-W. Pan, “Non–line-of-sight imaging over 1.43 km,” P. Natl. Acad. Sci. USA 118(10), e2024468118 (2021). [CrossRef]  

16. Z. Zhang, X. Liu, Y. Ma, N. Xu, W. Zhang, and S. Li, “Signal photon extraction method for weak beam data of ICESat-2 using information provided by strong beam data in mountainous areas,” Remote Sens. 13(5), 863 (2021). [CrossRef]  

17. K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting x-ray detectors in medical imaging: Vision 20/20: Photon counting detectors,” Med. Phys. 40(10), 100901 (2013). [CrossRef]  

18. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proceedings of the National Academy of Sciences 109(26), E1679–E1687 (2012). [CrossRef]  

19. W. Becker, A. Bergmann, M. A. Hink, K. König, K. Benndorf, and C. Biskup, “Fluorescence lifetime imaging by time-correlated single-photon counting: Fluorescence Lifetime Imaging by TCPSC,” Microsc. Res. Techniq. 63(1), 58–66 (2004). [CrossRef]  

20. Y. Chen and A. Periasamy, “Time-correlated single-photon counting fluorescence lifetime imaging–FRET microscopy for protein localization,” in Molecular Imaging (Elsevier, 2005), pp. 239–259.

21. X.-F. Liu, W.-K. Yu, X.-R. Yao, B. Dai, L.-Z. Li, C. Wang, and G.-J. Zhai, “Measurement dimensions compressed spectral imaging with a single point detector,” Opt. Commun. 365, 173–179 (2016). [CrossRef]  

22. J. W. Müller, “Dead-time problems,” Nuclear Instruments and Methods 112(1-2), 47–57 (1973). [CrossRef]  

23. J. Rapp, Y. Ma, R. M. A. Dawson, and V. K. Goyal, “Dead Time Compensation for High-Flux Ranging,” IEEE Trans. Signal Process. 67(13), 3471–3486 (2019). [CrossRef]  

24. F. Ichino, “Efficient model of photon-counting laser radar for distance error calibration,” Opt. Commun. 427, 278–287 (2018). [CrossRef]  

25. C. W. Helstrom, “Output distributions of electrons in a photomultiplier,” J. Appl. Phys. 55(7), 2786–2792 (1984). [CrossRef]  

26. Z. Zhang, Y. Ma, S. Li, P. Zhao, Y. Xiang, X. Liu, and W. Zhang, “Ranging performance model considering the pulse pileup effect for PMT-based photon-counting lidars,” Opt. Express 28(9), 13586–13600 (2020). [CrossRef]  

27. Y. Liu, J. Shi, and G. Zeng, “Single-photon-counting polarization ghost imaging,” Appl. Opt. 55(36), 10347–10351 (2016). [CrossRef]  

28. Y. Yang, J. Shi, F. Cao, J. Peng, and G. Zeng, “Computational imaging based on time-correlated single-photon-counting technique at low light level,” Appl. Opt. 54(31), 9277–9283 (2015). [CrossRef]  

29. W.-K. Yu, X.-F. Liu, X.-R. Yao, C. Wang, S.-Q. Gao, G.-J. Zhai, Q. Zhao, and M.-L. Ge, “Single photon counting imaging system via compressive sensing,” arXiv:1202.5866 [physics.optics] (2012).

30. J. Rapp, Y. Ma, R. M. A. Dawson, and V. K. Goyal, “High-flux single-photon lidar,” Optica 8(1), 30–39 (2021). [CrossRef]  

31. M. S. Oh, H. J. Kong, T. H. Kim, K. H. Hong, and B. W. Kim, “Reduction of range walk error in direct detection laser radar using a Geiger mode avalanche photodiode,” Opt. Commun. 283(2), 304–308 (2010). [CrossRef]  

32. L. Xu, Y. Zhang, Y. Zhang, L. Wu, C. Yang, X. Yang, Z. Zhang, and Y. Zhao, “Signal restoration method for restraining the range walk error of Geiger-mode avalanche photodiode lidar in acquiring a merged three-dimensional image,” Appl. Opt. 56(11), 3059–3063 (2017). [CrossRef]  

33. P. B. Coates, “The correction for photon `pile-up’ in the measurement of radiative lifetimes,” J. Phys. E: Sci. Instrum. 1(8), 878–879 (1968). [CrossRef]  

34. X.-F. Liu, X.-R. Yao, C. Wang, X.-Y. Guo, and G.-J. Zhai, “Quantum limit of photon-counting imaging based on compressed sensing,” Opt. Express 25(4), 3286–3296 (2017). [CrossRef]  

35. C. Yuan, Q. Yan, Y. Wu, Y. Wang, and Y. Wang, “Single Photon Compressive Imaging Based on Digital Grayscale Modulation Method,” Photonic Sens. 11(3), 350–361 (2021). [CrossRef]  

36. H. Wang, Q. Yan, B. Li, C. Yuan, and Y. Wang, “Sampling time adaptive single-photon compressive imaging,” IEEE Photonics J. 11(3), 1–10 (2019). [CrossRef]  

37. C. Li, W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented Lagrangian method with applications to total variation minimization,” Comput Optim Appl 56(3), 507–530 (2013). [CrossRef]  

38. A. Gupta, A. Ingle, A. Velten, and M. Gupta, “Photon-flooded single-photon 3D cameras,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (IEEE, 2019), pp. 6763–6772.

39. M. Patting, M. Wahl, P. Kapusta, and R. Erdmann, “Dead-time effects in TCSPC data analysis,” in International Congress on Optics and Optoelectronics (SPIE, 2007), p. 658307.

40. A. Hore and D. Ziou, “Image quality metrics: PSNR vs. SSIM,” in20th International Conference on Pattern Recognition (IEEE, 2010), pp. 2366–2369.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but can be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1.
Fig. 1. Schematic diagram of compressive sensing.
Fig. 2.
Fig. 2. Photograph (a) and schematic (b) of the single photon compressive imaging system in this study.
Fig. 3.
Fig. 3. Basic schematic of the TCSPC in the single photon imaging system.
Fig. 4.
Fig. 4. Difference between the signal photon number (using purple curve) and detection probability (using red curve) arising from the dead time effect. (a) Each time bin; (b) over entire pulse. In (a), the repeated illumination times R is set to 10000 and the expected signal photon number Ns is set to 1. In (b), the purple “*” and red “*” represent expected signal photon number and detection probability over entire illuminations in (a), respectively.
Fig. 5.
Fig. 5. Experimental results when the expected signal photon number Ns is 0.5, 0.8, 1.3 and 1.5 counts, and the background noise rate is 4 MHz. Under each signal level, the repeated illumination numbers of 2000, 5000 and 10000 were separately conducted. In each group, the left side is the reconstructed image by using the photon counting correction method proposed in this study, whereas the middle and right side correspond to the reconstructed images without correction and with the correction method by Coates [33], respectively. The PSNR values are in units of dB.
Fig. 6.
Fig. 6. Distribution of recovered signal photon numbers and normalized PEs within the signal duration of 102 modulation patterns with different signal levels (a) Ns=0.5, (b) 0.8, (c) 1.3, and (d) 1.5 counts when R = 10000.
Fig. 7.
Fig. 7. Flow chart of the Monte Carlo simulation.
Fig. 8.
Fig. 8. PSNR simulation results of reconstructed images with different Ns and R by utilizing (a) photon counting correction method proposed in this study, (b) traditional methods, and (c) PSNR improvement.
Fig. 9.
Fig. 9. PSNR comparison between experimental results and Monte Carlo simulation results varying with Ns when the repeated illuminations R is (a) 2000, (b) 5000, and (c) 10000.
Fig. 10.
Fig. 10. Simulation results of the photon counting correction method with the noise rate of fn = 4 MHz, the number of repeated illuminations of R = 5000, and the time bin of τ = 200ps. (a) The signal waveform satisfies the Gaussian distribution in the time domain with Ns = 1 and FWHM = 1 ns. (b) The signal waveform is overlapped by two Gaussian signal pulses that have Ns = 1 and FWHM = 1 ns. (c) The signal waveform satisfies a rectangular pulse with Ns = 1 and the pulse width of 2 ns. The blue histogram represents the actual signal photon number, and the purple curve is the recovered signal photon counts. The red histogram represents the normalized PEs, and the yellow curve is the detection probability.
Fig. 11.
Fig. 11. Simulation result of reconstructed images when the resolutions of the reconstructed images are 32×32 and 64×64. The expected signal photon number Ns under each resolution is set to Ns = 1 or 2, and the repeated illumination numbers R is 10000, whereas other parameters are identical with the experiment. The right side of each group of results is the reconstructed image of traditional methods, and the left side is the reconstructed image by using the photon counting correction method proposed in this study. The PSNR values are in units of dB
Fig. 12.
Fig. 12. PSNR comparison between the proposed method and traditional methods varying with Ns when the pixel numbers of reconstructed images are 32×32 (blue curves) and 64×64 (red curves). The repeated illumination number R is 10000. The PSNR values are in units ofdB

Equations (16)

Equations on this page are rendered with MathJax. Learn more.

y = Φ x
y = Φ x = Φ Ψ s
{ s = arg min 1 2 | | y Φ Ψ s | | 2 2  +  τ | | s | | 1 x = Ψ s ,
P ( t , t + τ ; λ ) = [ n s ( t ) + f n τ ] λ λ ! e [ n s ( t ) + f n τ ] ,
n s ( t , t + τ ) = t t + τ N s 2 π σ p e ( t t T ) 2 2 σ p 2 d t ,
P ( t , t + τ , λ  >  0 ) = 1 P ( t , t + τ , λ = 0 ) = 1 e [ n s ( t , t + τ ) + n n ] ,
P ( t i ) = { k = 1 i 1 [ P ( t k , λ = 0 ) ] [ 1 P ( t i , t i + τ ; λ = 0 ) ] ( 1 i n td ) k = i n td i 1 [ P ( t k , λ = 0 ) ] [ 1 P ( t i , t i + τ ; λ = 0 ) ] ( i > n td )
P ( t i ) = { P n ( λ = 0 ) i 1 [ 1 P n ( λ = 0 ) ] ( 1 i min ( m 1 , n td ) ) P n ( λ = 0 ) n td [ 1 P n ( λ = 0 ) ] ( n td < i m 1 )
P ( t m ) = { P n ( λ = 0 ) m 1 [ 1 P ( t m , λ = 0 ) ] ( 1 m n td ) P n ( λ = 0 ) n t d [ 1 P ( t m , λ = 0 ) ] ( m > n td )
P ( t m  + 1 ) = { P n ( λ = 0 ) m 1 P ( t m , λ = 0 ) [ 1 P ( t m  +  1 , λ = 0 ) ] ( 1 m  +  1 n td ) P n ( λ = 0 ) n t d 1 P ( t m , λ = 0 ) [ 1 P ( t m  + 1 , λ = 0 ) ] ( m  +  1 > n td )
P ( t i ) H ( t i ) / R ,
P n ( λ = 0 ) = { s o l v e P n ( λ = 0 ) i P n ( λ = 0 ) i 1 + H ( t i ) / R = 0 ( 1 i min ( m 1 , n td ) ) s o l v e P n ( λ = 0 ) n td + 1 P n ( λ = 0 ) n td + H ( t i ) / R = 0 ( n td < i m 1 )
P ( t m , λ = 0 ) = { 1 H ( t m ) / R P n ( λ = 0 ) m 1 ( 1 m n dt ) 1 H ( t m ) / R P n ( λ = 0 ) n t d ( m > n dt ) .
n s ( t m ) = ln ( P ( t m , λ = 0 ) ) n n .
P S N R = 10 log 10 M A X I 2 1 m n i = 1 m j = 1 n [ I ( i , j ) I 0 ( i , j ) ] 2 ,
N s j = N I × Φ j x N ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.