Enhancing the spatial resolution of time-of-flight based non-line-of-sight imaging via instrument response function deconvolution

DingJie Wang; DingJie Wang; Wei Hao; Wei Hao; Wei Hao; YuYuan Tian; YuYuan Tian; WeiHao Xu; WeiHao Xu; Yuan Tian; Yuan Tian; HaiHao Cheng; HaiHao Cheng; SongMao Chen; SongMao Chen; SongMao Chen; Ning Zhang; WenHua Zhu; XiuQin Su; XiuQin Su; XiuQin Su

doi:10.1364/OE.518767

1. Introduction

Non-line-of-sight (NLOS) imaging has gained significant attention in recent years due to its potential applications in surveillance, autonomous driving, and rescue operations [1]. Over the past decade, NLOS imaging based on different sensing schemes has been developed, such as thermal imaging [2,3], synthetic wavelength holography (SWH) [4], millimeter-wave synthetic aperture radar (SAR) [5], acoustic [6], speckle coherence [7–9], RGB camera [10,11], electromagnetic waves like doppler [12], polarized infrared [13], multispectral and hyperspectral [14,15], event camera [16], and ordinary camera [17], etc. In particular, time-correlated single photon counting (TCSPC) technology, with its single-photon detection sensitivity and picosecond timing resolution, has greatly promoted the development of NLOS imaging in past few years. Such TCSPC-based NLOS imaging scheme records the time-of-flight (TOF) information of photons reflected from the hidden scene, and its detection mode typically employ picosecond pulsed laser (a few femtosecond laser sources are also used [18,19]) as light source and single photon detector (SPD) as receive terminals, including streak camera [19,20], intensified charge coupled device (ICCD) [21–24], Geiger-mode single-photon avalanche diode (SPAD) [25–27], and superconducting nanowire single-photon detectors (SNSPD) [28,29]. Based on this scheme, various ingenious algorithms have been developed to reconstruct hidden scene from the transient data, such as filtered backprojection (FBP) [19], light-cone transform (LCT) [30], frequency-wavenumber migration (FK) [31], phasor-field (PF) [32], signal–object collaborative regularization (SOCR) [33,34], geometry-based Fermat path [35], etc., which show good results in reconstruction details, robustness to noise and reconstruction efficiency, etc. Some other research combines innovations in hardware has pushed the limits of NLOS imaging, such as long range NLOS imaging and tracking [36–38], high-precision NLOS imaging [18], real-time NLOS imaging [39–41], and keyhole NLOS imaging [42], etc.

Despite the great potential of TCSPC-based NLOS imaging, an important metric of NLOS imaging, i.e., its spatial resolution, is constrained by the TCSPC system’s total time resolution (or total time jitter, TTJ) [18,30,43], which is comprehensive of the pulse width of the laser and the time jitter of each electronic devices (including detectors, counting modules, synchronization signals, etc.) [44]. The most widely used TCSPC-based imaging systems have TTJ of around a few hundred picoseconds [30,36,39,45], limiting their spatial resolution to above $\sim 10\mathrm{cm}$ or slightly below that level. To break through the limitations of the spatial resolution, Wang et al. developed innovations in hardware devices by establishing the imaging system of frequency up-conversion single-photon detector (UCSPD) and femtosecond pulse laser uses pulse pump frequency up-conversion detection technology to achieve a time resolution of 1.4ps in the near-infrared single-photon detector, within the high temporal resolution, high-precision NLOS imaging with lateral spatial resolution of $\sim 2\mathrm{mm}$ and axial spatial resolution of $\sim 180\mathrm{um}$ was successfully achieved [18]. Some other research efforts have focused on the after reconstruction procedure, in 1.43km NLOS imaging [36], Wu et al. modeled temporal and spatial broadening in long-range conditions as 1D and 2D Gaussian model, respectively, and incorporated them into the forward model, which achieved lateral resolution of $\sim 9.4\mathrm{cm}$ after reconstruction at overall resolution of 1.1ns (here the 1.1ns incorporates spatial broadening in addition to the TTJ we have described here). Jin et al. constructed the point spread function (PSF) from the corresponding direct reflection echoes in the confocal scanning scheme and used it as a prior for histogram deconvolution, and improved the spatial resolution to about $1/3$ of the original [46].

The instrument response function (IRF) of the TCSPC system statistically represents the distribution of system’s TTJ. Herein, the IRF deconvolution (IRF-DC) method, originally used in fluorescence lifetime measurement to extract the "true" decay curve from the measured curve [47], was introduced to improve the spatial resolution of NLOS imaging. This is achieved by modeling the mechanisms by which the TTJ degrades the spatial resolution of NLOS imaging as Poisson process of convolution with the normalized IRF as kernel, thus transforming the inverse problem into deconvolution which can be mathematically solved in Bayesian framework.

The primary research of this paper is as follow: in Section 2, we give a statement of the problem, model the degradation mechanism of TTJ on the spatial resolution, and mathematically solve the inverse problem; in Section 3, we customize a numerical simulation model to quantitatively verify the effectiveness of the IRF-DC method in enhancing the spatial resolution of NLOS imaging, in parallel with the advantages for reconstruction at low SNR situation; in Section 4, we implemented an TCSPC-based NLOS imaging system to experimentally validate the effectiveness of the IRF-DC method in real-word scenario; in Section 5, we developed discussion of the methodology; and finally, the conclusion was drawn in Section 6.

2. Methodology

2.1 Problem statement

Figure 1(a) illustrates the confocal NLOS (C-NLOS) imaging scenario in which case the sensing point locates same to illumination point. According to the theory proposed by O’Toole [30], the minimum resolvable lateral distance $\Delta x$ and axial distance $\Delta z$ of C-NLOS imaging is:

(1)$$\Delta x\ge \frac{C\sqrt{W^2+Z^2}}{2W}{\gamma} \quad \Delta z\ge \frac{C}{2}{\gamma}$$

where ${\gamma }\approx \sqrt {\sum {{\gamma } _{component}^{2}}}$ represents the TTJ of the TCSPC system (characterized with full width at half maximum, FWHM), and ${\gamma } _{component}$ is time jitter of each component [44]; $C$ is the speed of light; $2W$ is the size of the scan area, and $Z$ is the distance between the hidden scene and relay wall. Eq. (1) reveals that the spatial resolution boundary of C-NLOS is fundamentally determined by the TCSPC system’s TTJ, and the non-confocal (NC-NLOS) configuration also follow the rule, but the lateral resolution is theoretically reduced by a factor of two [43]. Figure 1(b) shows a group of C-NLOS simulation results based on the numerical model in Section 3, the transient data is sampled at $64\times 64$ locations over $1.6\mathrm{m} \times 1.6\mathrm{m}$ region, and the hidden Bunny is $0.8\mathrm{m}$ away from the relay wall, results with the LCT solver indicate that increasing the TTJ from 100ps to 1300ps significantly reduces spatial resolution.

Fig. 1. Schematics of C-NLOS layout and our proposed method. (a) The schematics of C-NLOS imaging. (b) The spatial resolution decreases progressively as the increasing of the system’s TTJ. (c) The pipeline of our IRF-DC method, the LCT reconstructed spatial resolution when TTJ=1300 ps is driven to be close to that of TTJ=100 ps.

Download Full Size | PDF

2.2 Problem formation and mathematically solving

In real NLOS imaging scenario, the optical field of view (FOV) also introduces temporal broadening of the NLOS distribution apart from the TTJ delineated in Section 2.1 [36]. Assuming the NLOS distribution broadened only by the FOV is $\mathbf {g}$, which can be interpreted as the NLOS distribution observed when the IRF is a Dirac function. Then due to the temporal broadening induced by the system’s TTJ, the observed $\mathcal {G}$ can be seen as the convolution of $\mathbf {g}$ with the normalized IRF, similar to the time-resolved fluorescence measurement [44,48], given as:

(2)$$\mathcal{G} (t)= \mathcal{P}\left \{ \mathbf{g}(t)\ast \gamma(t) \right \}$$

where $\mathcal {P}$ represents the Poisson process due to the SPAD in photon starved scenario; "$\ast$" represents convolution operator in time domain, and $t\in \left [ 1,2,\ldots T \right ]$ denotes the $t$th timestamp index of the recorded histogram; $\gamma$ represents the normalized IRF of TCSPC system.

Within the theoretical model in Eq. (2), mathematically we can estimate $\mathbf {g}$ from the recorded $\mathcal {G}$ based on bayesian probabilistic framework, and the inverse problem can be written as:

(3)$$\hat{\mathbf{g}} = \underset{\mathbf{g}}{\text{argmax}} \quad P(\mathbf{g}\mid \mathcal{G})\propto P(\mathcal{G}\mid \mathbf{g})P(\mathbf{g})$$

where $\hat {\mathbf {g}}$ represents the estimation of $\mathbf {g}$. When performing estimation without considering the prior $P(\mathbf {g})$, it leads to a maximum likelihood estimation (MLE) of $g$, and the likelihood function $P(\mathcal {G}\mid \mathbf {g} )$ can be derived based on the Poisson statistical model, given as follow:

(4)$$P(\mathcal{G}\mid \mathbf{g})=\prod_{t=1}^{T} \left ( \frac{(\mathbf{g}\ast\gamma)^{\mathcal{G}}\cdot e ^{-(\mathbf{g}\ast\gamma)} }{\mathcal{G} !} \right )$$

Directly solving the aforementioned Poisson model based joint probability function remains to be challenging. Typically, the problem of maximum the likelihood function $P(\mathcal {G}\mid \mathbf {g})$ equivalents to minimizing its negative log-likelihood, and then the estimation of $\mathbf {g}$ can be transformed into a minimization problem that described as follow:

(5)$$\hat{\mathbf{g}}=\underset{{\mathbf{g}}}{\text{argmin}}\sum_{t=1}^{T} {\left[ -\gamma \cdot \log (\mathbf{g}\ast\gamma) +\left(\mathbf{g}\ast\gamma \right) + \log (\mathcal{G}!)\right]}$$

Since the optimization problem is convex with respect to $\mathbf {g}$, solving the above problem with zero gradient yields an equational constraint on $\mathcal {G}$, $\mathbf {g}$ and $\gamma$ [49]. Then the Richardson-Lucy (RL) algorithm [50] can be employed to solve the problem, RL assumes that the ratio of the optimal solution of $\mathbf {g}$ before and after the iteration tends to 1 as it converges, thus derive a multiplicative iterative solution scheme of $\mathbf {g}$, the iterative formula for $\mathbf {g}$ is given as follows:

(6)$$\hat{\mathbf{g}}_{k+1} (t)=\hat{\mathbf{g}}_{k}(t)\cdot\left [ \bar{\gamma} (t) \ast \frac{\mathcal{G} (t)}{\hat{\mathbf{g}}_{k}(t)\ast \gamma(t)} \right ]$$

where the $\hat {\mathbf {g}}_{k}$ is the estimation of $\mathbf {g}$ after $k$th iterations, and $\bar {\gamma }$ is the conjugate of $\gamma$. Note that the standard RL often lead to a noise-dominated solution as iteration $k\to \infty$, thus an empirical judgment of the iteration interval is needed to prevent noise amplification. Typically the recovered $\hat {\mathbf {g}}$ with IRF-DC remains to be narrower, with higher peaks and more features with respect to measured $\mathcal {G}$, as shown in Fig. 1(c), which enables the reconstruction solvers to enhance the spatial resolution as well as to improve the reconstruction capability at low SNR situation.

2.3 Posterior reconstruction

The IRF-DC is a pre-processing procedure prior to reconstruction by performing deconvolution on the measured transient data, and the only information needed is the IRF of the TCSPC system, which is usually obtained by prior calibration step, thus the IRF-DC is theoretically applicable to both C-NLOS and NC-NLOS schemes, and even SPAD array-based schemes. Futher, typical NLOS reconstruction solvers such as the FBP (for NC-NLOS), LCT and FK(for C-NLOS), and PF (for both NC-NLOS and C-NLOS), etc., can be used for post-reconstruction.

3. Numerical simulation and quantitative assessment

3.1 Numerical model

The numerical model was specifically tailored, which enabled us to precisely control the characteristics of IRF and quantitatively evaluate the effectiveness of the IRF-DC. Note that here we focus on the $3^{rd}$ reflection from the hidden scene, instead of $1^{st}$ reflection from the relay wall.

3.1.1 Energy attenuation of NLOS transfer

The energy attenuation model for three reflections in NLOS scenes is developed based on radiometric theory. Denote the spatial positions of the light source and detector by $R_{l}$ and $R_{c}$, and assume that the reflective volume of hidden target is $\mathcal {S}$ with albedo $\alpha _{T}$, $(r_{l},r_{c})$ denotes a pair of illumination point and sensing point on the relay wall with albedo $\alpha _{W}$. Taking a voxel element $s(s\in \mathcal {S} )$ as example, the echo pulse energy $\mathbf {E}_{r_{l}\to r_{c}}^{s}$ contributed by element $s$ is:

(7)$$\mathbf{E}_{r_{l}\to r_{c}}^{s} =\mathbf{E}\cdot {\frac{\alpha _{W}\cdot \Delta s}{2\pi \cdot \left\| s-r_{l}\right\|^{2}}} \cdot \frac{\alpha _{T} \cdot A_{FOV}}{2\pi \cdot \left\| r_{c}-s\right\|^{2}}\cdot \frac{\alpha _{W}\cdot A_{Rec}}{2\pi\cdot \left\| R_{c}-r_{c}\right\|^{2}}$$

where $\mathbf {E}$ represents the emitted single pulse energy; $\Delta s$, $A_{FOV}$ and $A_{Rec}$ represent the area of element $s$ , the area of field of view (FOV) on the relay wall and the area of optical system’s receiving aperture respectively. And the three fractional terms characterize the $1^{st}$, $2^{nd}$, and $3^{rd}$ reflection decays during transmission respectively.

3.1.2 Observation model

Based on the energy decay model described in Eq. (7) and the Poisson distribution model of the TCSPC system, we can derive a observation model of the $3^{rd}$ reflection for the illumination point to sensing point pair $(r_{l},r_{c})$:

(8)$$\mathcal{G}(r_{l}\to r_{c},t)\sim \mathcal{P}\left \{ \eta\left [ \iiint_{{\mathcal{S} }}^{} \left ( \frac{\mathbf{E}_{r_{l}\to r_{c}}^{s} }{h\nu }\cdot {\gamma} (t-t_{r_{l}\to r_{c}}^{s}) \right )ds +b_{\lambda } \right ]+b_{d} \right \}$$

where $\eta$ represents the quantum efficiency of the detector at working band $\lambda$; $b_{\lambda }$ and $b_{d}$ are background noise and dark counting noise of the detector respectively; $h$ and $\nu$ are Planck’s constant and electromagnetic frequency, and the $\frac {\mathbf {E}_{r_{l}\to r_{c}}^{s}}{h\nu }$ term is the number of echo photons contributed by voxel element $s(s\in \mathcal {S} )$ in single laser pulse period, the integral term indicates that the return signal is contributed by the entire reflective volume $\mathcal {S}$. In particular, $\gamma (t)$ denotes the normalized IRF, and $t_{r_{l}\to r_{c}}^{s}=\frac {\left \|r_{l}-R_{l} \right \|+\left \|s-r_{l} \right \|+\left \|r_{c}-s \right \|+\left \|R_{c}-r_{c} \right \|}{C}$ represents the TOF corresponds to total flight path $R_{l}\to r_{l}\to s\to r_{c}\to R_{c}$. The proposed IRF-DC method has no requirement for the mathematical model of the IRF, here we choose the classical heavy-tailed distribution model [51] to model the IRF in the observation equations, that is:

(9)$$\gamma(t)=\left | (\frac{t}{\tau } )^{2} \cdot e^{-\frac{t}{\tau }} \right | _{Norm},\quad\tau =\frac{\gamma _{{\tiny FWHM} }}{3.5}$$

where $\left | \bullet \right | _{Norm}$ is normalization operator. $\gamma _{{\tiny FWHM} }$ is approximately equal to the FWHM of the IRF curve, and we use this parameter to characterize the TTJ of the TCSPC system.

3.2 Quantitative evaluation

Based on the numerical model, we conducted numerical simulations to thoroughly evaluate the performance of proposed IRF-DC method. First, we validate the proposed method in improving the spatial resolution under sufficiently cumulative conditions, and then, we validate the method to achieve a high quality reconstruction under low SNR conditions. All of the reconstruction was shown with maximum intensity projections (MIP) along the axial direction in the Jet colormap. The structural similarity (SSIM) and peak signal-to-noise ratio (PSNR) are used to quantitative evaluate the quality of the reconstructed MIP images.

3.2.1 Spatial resolution enhancement

We constructed a "E-Chart" target consisting of multiple letter "E" with different size, and all letters were distributed in a plane 0.5m from the relay wall. The transient data is sampled at $64\times 64$ locations on the relay wall over $1\mathrm{m}\times 1\mathrm{m}$ region under C-NLOS settings, and the bin resolution of the histogram is 4ps. Firstly we reconstructed the ground truth with PF solver when the IRF is the Dirac function, i.e., the TOF of the photon event is not affected by jitter, as shown in Fig. 2(a). Then we simulate the TTJs in order of 200ps, 400ps, 600ps, 800ps, 1000ps and 1200ps in the simulation, and the distributions of corresponding IRFs are shown in Fig. 2(b). To eliminate the influence of noise or insufficient accumulation on the reconstruction results, all transient data are fully sampled, i.e., with high SNR.

Fig. 2. The Ground Truth of "E-Chart" and distribution of heavy-tailed IRFs as the TTJ ranges from 200 ps to 1200 ps, the IRFs have been normalized.

Download Full Size | PDF

In Fig. 3, we show the reconstruction results based on the direct LCT, FK and PF solver. From the results in the first, third and fifth row, we can see that at a TTJ of 200 ps, the result of the LCT and FK solver start to show a significant decrease in spatial resolution, and the PF solver is relatively more robust to the increase of TTJ, but when the TTJ reaches to 600 ps, a serious decrease in reconstruction resolution also occurs. The second, fourth and sixth row shows the results of LCT, FK and PF solver with IRF-DC method, showing that the IRF-DC method motivates the three solvers to achieve relatively clearer and more stable reconstructions. Even with TTJ reaching 1200 ps, IRF-DC can drive LCT, FK, and PF solver to reconstruct more reliable results.

Fig. 3. NLOS reconstruction when TTJ ranges from 200 ps to 1200 ps, the first and second row is the result of LCT solver .w/o and .w/ proposed IRF-DC method, the third and fourth row is the result of FK solver .w/o and .w/ IRF-DC method, and the fifth and sixth row is the result of PF solver .w/o and .w/ IRF-DC method.

Download Full Size | PDF

In Fig. 4, we show the quantitative evaluation of reconstruction results in Fig. 3, the dotted line represent the reconstruction using direct LCT, FK and PF solver, and the solid line represent the reconstruction of LCT, FK and PF solver based on the IRF-DC method. We can see that for all three solvers, when TTJ increases from 200ps to 1200ps, the IRF-DC can promote a large improvement in PSNR and SSIM: for the LCT solver, the PSNR increment is greater than 3.2047dB (TTJ= 1200ps), the SSIM increment is greater than 0.2365 (TTJ=200ps); for the FK solver, the PSNR increment is greater than 2.3981dB (TTJ=1200ps), the SSIM increment is greater than 0.261 (TTJ=1200ps); and for the PF solver, the PSNR increment is greater than 2.3326dB (TTJ=400ps), SSIM increment is greater than 0.0182 (TTJ=200ps).

Fig. 4. Quantitative evaluation of the reconstruction results of E-Chart. (a) PSNR of the reconstructed MIPs with LCT, FK, PF solver as TTJ ranges from 200 ps to 1200 ps .w/ and .w/o IRF-DC. (b) SSIM of the reconstructed MIPs with LCT, FK, PF solver as TTJ ranges from 200 ps to 1200 ps .w/ and .w/o IRF-DC.

Download Full Size | PDF

3.2.2 Reconstruction within low SNR

We further verified the effectiveness of the IRF-DC method in low SNR situations. Here we constructed a combination of letters "UCAS" distributed in the same plane 0.5m from the relay wall, and the transient data is also sampled at $64\times 64$ locations on the relay wall over a $1\mathrm{m}\times 1\mathrm{m}$ region. To eliminate the interference of TTJ on reconstruction, the TTJ fixedly to 500ps, and control the SNR by changing the number of noise counts in the observation Eq. (8). Here we define the mean noise count (MNC) of average pixel to characterize the noise intensity, that is, the noise count after summing and averaging the data of all sampling points. Since the accumulation time is fixed, the larger the MNC, the lower the SNR of the data. The ground truth is reconstructed with the PF solver at 0 noise counts and the IRF is a Dirac function, as shown in Fig. 5(a). And then, we sequentially set the values of MNC to 0, 2.7, 5.4, 8.1, 10.8 and 13.5 counts, and Fig. 5(b) shows the photon distribution of average pixel when MNC increased.

Fig. 5. The Ground Truth of "UCAS" and average pixel photon distribution as MNC ranges from 0 counts to 13.5 counts.

Download Full Size | PDF

In Fig. 6, we show the reconstruction based on PF solver, and the results were quantitatively evaluated with PSNR. The reconstruction results for MNC=0, 2.7, 5.4, 8.1, 10.8, and 13.5 counts are shown from left to right. The first row is the result of using the direct PF solver, we can see that when MNC = 0 counts, the PF solver produces a low-quality reconstruction, the profile of "UCAS" is contaminated by noise, which is due to insufficient accumulation. As MNC increases to 13.5 counts, we can see that the direct PF solver results in failed reconstructions. The second row is the result of using the PF solver within the IRF-DC method. Although the increase in MNC also caused a slight degradation of the reconstruction, even when MNC = 13.5 counts, a clear "UCAS" profile was reconstructed (PSNR=11.3045dB), and the clarity of the reconstruction was better than that of the reconstruction using the direct PF solver when MNC = 0 counts (PSNR=9.8622dB).

Fig. 6. Reconstruction with PF solver when MNC ranges from 0 to 13.5 counts, the first and second row is the result .w/o and .w/ proposed IRF-DC method respectively.

Download Full Size | PDF

4. Experimental

4.1 Hardware setup

Real-world C-NLOS experiment is conducted to validate the effectiveness of the proposed IRF-DC method. The experimental scenario is shown in Fig. 7(a), the workflow of the system is as follow: the laser source (PicoQuant LDH-D-C-850) emits long pulse (FWHM=500ps) with wavelength of 850nm and average power is 0.1mW@1MHz; the emitted beam is output through free space, adjusted by two mirrors and coupled into the perforated mirror, the direction of laser beam is controlled by 2-axis scanner (GVS-012) to achieve scanning sampling, which is driven by a digital to analog converter (DAC8563), the emitted pulse is reflected on the relay wall, hidden target and relay wall in turn, and finally collected by the optical system along the same path. The photon after three reflections are coupled into the SPAD detector (EXCELITASDTS-SPCMAQRH-16-FC) with quantum efficiency $\sim 40{\% }@850\mathrm{nm}$ and time jitter $\sim 350\mathrm{ps}@\text {Typ}$. Channel 0 of the time-correlated photon counter TCSPC (PicoQuant PicoHarp 300) is the synchronization signal from the pulsed laser, which serves as the reference for the start time, and Channel 1 is the trigger signal from the SPAD. Figure 7(b) shows the details of the confocal light path, to filter out background noise outside the working band, a narrowband filter (Semrock FF01-850/10-2) is added before the lens.

Fig. 7. Experimental setup. (a) Overall perspective of the imaging scenario: the top left corner is the line-of-sight (LOS) perspective, the target "S" is obscured by black curtain, the region encircled by the red dotted line is the scan area. The bottom right corner is the front view of the target "S". Both the scan area and target "S" are made of foam sheet material that exhibits strong diffuse reflective properties. (b) Confocal light path: the coaxiality of the transmitting and receiving optical path is achieved by perforated mirror, the bottom right corner is the close-up of real optical part. (c) The normalized IRF of our TCSPC system: the TTJ of the system is FWHM = 608 ps.

Download Full Size | PDF

4.2 Calibration of system’s IRF

The calibrated IRF (normalized ) of the experimental setup is shown in Fig. 7(c), to reduce the reconstruction time, the invalid bin of the IRF front section is removed to reduce the length of the deconvolution kernel vector. The calibration of the IRF is similar to that of single-photon imaging and includes the following two steps: (i) adjust the optical system so that the laser spot and the FOV coincide to obtain the direct reflected photons. (ii) perform accumulation at safe power condition to reduce the effect of Poisson noise. Since the intensity of the direct reflection is much higher than that of three reflections in NLOS imaging, it may be necessary to reduce the power by increasing the optical path attenuation or adjusting the laser parameters to protect the detector. In our practice, due to the low power of the pulsed laser, we accumulated the IRF for 1 minute under the same conditions as the above experiments (0.1mW), and the calibrated TTJ of our system is about 608ps (FWHM), which is mainly determined by the pulse width of laser, the jitter of SPAD, TCSPC and Synchronization signal from laser source.

4.3 Result and analysis

The hidden target is a letter "S" with dimensions of $0.6\mathrm{m}\times 0.9\mathrm{m}$, and about 0.7m away from the relay wall. We took a $33\times 33$ sweep scanning sample over $1.12\mathrm{m}\times 1.12\mathrm{m}$ region of the relay wall, and the reconstructed region in the depth direction is $z\in [0.6\mathrm{m},1.2\mathrm{m}]$.

Figure 8 illustrates the reconstruction results based on LCT solver, FK solver and PF solver. Figure 8(a) shows direct reconstruction with LCT solver, we can see the general profile of the target "S", but the reconstructed details are lost, for example, the corners of the target "S" are blunted, which is a typical reduction in spatial resolution due to TTJ=608ps. Figure 8(c) is the direct reconstruction result based on the FK solver, shows a worse reconstruction result since only the center part of the target is reconstructed. Figures 8(b) and 8(d) show the reconstruction results of the LCT solver and the FK solver based on IRF-DC respectively, both showing relatively complete reconstruction. Figure 8(e) shows the direct reconstruction based on the PF solver, we can see a relatively complete profile of "S", but from a 3D perspective, we can see that there are a lot of artifacts in the reconstruction area, the reason is that our measurements is not fully accumulated, that is, the SNR of original transient data is very low, but the TTJ=608ps is still within the allowable range of the PF solver, this is the reason why the target profile in Fig. 8(e) is clear but artifacts exist. Figure 8(f) is the reconstruction result of the PF solver based on IRF-DC which shows perfect reconstruction.

Fig. 8. Experimental results. (a),(b) Reconstruction based on LCT solver .w/o and .w/ IRF-DC method. (c),(d) Reconstruction based on FK solver .w/o and .w/ IRF-DC method. (e),(f) Reconstruction based on PF solver .w/o and .w/ IRF-DC method.

Download Full Size | PDF

In our experiments, for the LCT and FK solver, the IRF-DC offers more benefits in enhancing spatial resolution, which can be observed by comparing Fig. 8(a) with Fig. 8(b), or Fig. 8(c) with Fig. 8(d); while for the PF solver, the IRF-DC has obvious effects in improving the reconstructed SNR, seen from the comparison between Fig. 8(e) and Fig. 8(f).

5. Discussion

Here we develop discussion about the IRF-DC method, including two points: computing time cost and applicability, it is specified as follows:

• The IRF-DC is a preprocessing procedure that requires additional computational time, the inverse problem is solved using the RL algorithm, which gradually approaches the optimal solution through iterations. In our study, we found that the larger the TTJ is, the more iterations are required, in the simulation in Section 3.2, the number of iterations are $k$=50 (TTJ=200ps, 400ps) and $k$=70 (TTJ=600ps, 800ps, 1000ps, 1200ps), and the computation time is also related to the size of the data volume. In the real-world experiments in Section 4, we have a data volume of $33\times 33\times 512$, with iteration number of $k$=20, and the runtime is 2.74s with Matlab platform on an Intel- i7-10875H CPU @ 2.30GHz with 16GB.
• For TCSPC-based NLOS imaging schemes, the manufacturing defects of the hardware will make time-resolved detection modes suffer from TOF measurement errors, i.e., the problem of time jitter will always exist. In contrast to method proposed in [46], which relies on the $1^{st}$ reflection echo in C-NLOS schemes to construct the PSF, here we can obtain the IRF of the TCSPC system through a prior calibration step, which makes the IRF-DC method applicable to both single-pixel scanning schemes (including C-NLOS and NC-NLOS) and detector array-based NLOS imaging schemes.

6. Conclusion

In summary, we propose an instrumental response function deconvolution (IRF-DC) method to overcome the bottleneck of spatial resolution of NLOS imaging caused by TCSPC system’s TTJ. This is achieved by modeling the degradation mechanism of spatial resolution as as Poisson convolution process with the normalized IRF as convolution kernel, and then solve the inverse problem with iterative deconvolution method. To verify the effectiveness of the method, we first establish a numerical model, the simulation results show that at TTJ=1200ps, IRF-DC can promote the LCT and FK solvers to reconstruct better than the direct reconstruction at TTJ=200ps, as well as can promote the PF solver to still reconstruct better than the direct reconstruction at TTJ=600ps; meanwhile, the reconstruction of the measured data at low SNR by the IRF-DC method also shows obvious advantages. Besides, we carried out indoor C-NLOS imaging experiments to further validate the performance of the method, the experimental results demonstrate that at TTJ = 608 ps, the IRF-DC motivates the LCT, FK, and PF solvers to achieve sharper reconstruction for low SNR measurements, whereas the direct reconstruction all have difficulty in recognizing the target contour. The proposed IRF-DC can effectively break through the spatial resolution limitation of NLOS imaging due to hardware defects and promote the development of high-resolution NLOS imaging. In future, we will investigate how to find the optimal solution quickly rather than relying on experience and iterative testing of the optimal number of iterations.

Funding

National Natural Science Foundation of China (62305375); Strategic High Technology Innovation Fund of The Chinese Academy of Sciences (GQRC-19-19); China Postdoctoral Science Foundation (2020M683600); Open Research Fund for Development of High-end Scientific Instruments and Core Components of the Center for Shared Technologies and Facilities (E32931Q101).

Acknowledgments

The authors would like to acknowledge Dr. Wenhua Zhu, Dr. Haihao Cheng and Dr. Ningbo Xie for their technical assistance in optical system design.

Disclosures

The authors declare no conflicts of interest.

Data availability

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

References

1. Y. Altmann, S. McLaughlin, and M. J. Padgett, “Quantum-inspired computational imaging,” Science 361(6403), eaat2298 (2018). [CrossRef]

2. M. Kaga, T. Kushida, and T. Takatani, “Thermal non-line-of-sight imaging from specular and diffuse reflections,” IPSJ Transactions on Computer Vision and Applications 11(1), 8 (2019). [CrossRef]

3. T. Sasaki, C. Hashemi, and J. R. Leger, “Passive 3d location estimation of non-line-of-sight objects from a scattered thermal infrared light field,” Opt. Express 29(26), 43642–43661 (2021). [CrossRef]

4. F. Willomitzer, P. V. Rangarajan, and F. Li, “Fast non-line-of-sight imaging with high-resolution and wide field of view using synthetic wavelength holography,” Nat. Commun. 12(1), 6647 (2021). [CrossRef]

5. Y. Wen, S. Wei, J. Wei, et al., “Non-line-of-sight imaging of hidden moving target using millimeter-wave inverse synthetic aperture radar,” in International Geoscience and Remote Sensing Symposium (IEEE, 2022), pp. 555–558.

6. D. B. Lindell, G. Wetzstein, and V. Koltun, “Acoustic non-line-of-sight imaging,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition (2019), pp. 6780–6789.

7. S. Divitt, D. F. Gardner, and A. T. Watnik, “Imaging around corners in the mid-infrared using speckle correlations,” Opt. Express 28(8), 11051–11064 (2020). [CrossRef]

8. O. Katz, P. Heidmann, M. Fink, et al., “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nat. Photonics 8(10), 784–790 (2014). [CrossRef]

9. M. Zhou, Y. Zhang, and P. Wang, “Speckle-correlation-based non-line-of-sight imaging under white-light illumination,” Optics & Laser Technology 170, 110231 (2024). [CrossRef]

10. C.-Y. Tsai, A. C. Sankaranarayanan, and I. Gkioulekas, “Beyond volumetric albedo–a surface optimization framework for non-line-of-sight imaging,” in IEEE/CVF conference on computer vision and pattern recognition (2019), pp. 1545–1555.

11. K. L. Bouman, V. Ye, A. B. Yedidia, et al., “Turning corners into cameras: Principles and methods,” in International Conference on Computer Vision (IEEE, 2017), pp. 2270–2278.

12. N. Scheiner, F. Kraus, F. Wei, et al., “Seeing around street corners: Non-line-of-sight detection and tracking in-the-wild using doppler radar,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition (2020), pp. 2068–2077.

13. H. Liu, P. Wang, and X. He, “Pi-nlos: polarized infrared non-line-of-sight imaging,” Opt. Express 31(26), 44113–44126 (2023). [CrossRef]

14. C. Hashemi, R. Avelar, and J. Leger, “Clutter rejection in passive non-line-of-sight imaging via blind multispectral unmixing,” Opt. Express 32(2), 2132–2146 (2024). [CrossRef]

15. C. Hashemi, R. Avelar, and J. Leger, “Isolating signals in passive non-line-of-sight imaging using spectral content,” IEEE Transactions on Pattern Analysis and Machine Intelligence (2023).

16. C. Wang, Y. He, and X. Wang, “Passive non-line-of-sight imaging for moving targets with an event camera,” Chin. Opt. Lett. 21(6), 061103 (2023). [CrossRef]

17. R. Czajkowski and J. Murray-Bruce, “Two-edge-resolved three-dimensional non-line-of-sight imaging with an ordinary camera,” Nat. Commun. 15(1), 1162 (2024). [CrossRef]

18. B. Wang, M.-Y. Zheng, and J.-J. Han, “Non-line-of-sight imaging with picosecond temporal resolution,” Phys. Rev. Lett. 127(5), 053602 (2021). [CrossRef]

19. A. Velten, T. Willwacher, and O. Gupta, “Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging,” Nat. Commun. 3(1), 745 (2012). [CrossRef]

20. O. Gupta, T. Willwacher, and A. Velten, “Reconstruction of hidden 3d shapes using diffuse reflections,” Opt. Express 20(17), 19096–19108 (2012). [CrossRef]

21. F. Heide, L. Xiao, W. Heidrich, et al., “Diffuse mirrors: 3d reconstruction from diffuse indirect illumination using inexpensive time-of-flight sensors,” in Conference on Computer Vision and Pattern Recognition (IEEE, 2014), pp. 3222–3229.

22. M. Laurenzis and A. Velten, “Non-line-of-sight active imaging of scattered photons,” in Electro-Optical Remote Sensing, Photonic Technologies, and Applications VII; and Military Applications in Hyperspectral Imaging and High Spatial Resolution Sensing, vol. 8897 (SPIE, 2013), pp. 47–53.

23. M. Laurenzis and A. Velten, “Feature selection and back-projection algorithms for nonline-of-sight laser–gated viewing,” Journal of Electronic Imaging 23(6), 063003 (2014). [CrossRef]

24. M. Laurenzis and A. Velten, “Nonline-of-sight laser gated viewing of scattered photons,” Opt. Eng. 53(2), 023102 (2014). [CrossRef]

25. M. Buttafava, J. Zeman, and A. Tosi, “Non-line-of-sight imaging using a time-gated single photon avalanche diode,” Opt. Express 23(16), 20997–21011 (2015). [CrossRef]

26. J. Zhao, F. Gramuglia, P. Keshavarzian, et al., “A gradient-gated spad array for non-line-of-sight imaging,” IEEE J. Selected Topics in Quantum Electronics 30(1), 8000110 (2023). [CrossRef]

27. Y. Zheng, W. Wang, and C. Zhang, “Converting non-confocal measurements into semi-confocal ones with timing-accuracy improving for non-line-of-sight imaging,” Optics and Lasers in Engineering 176, 108067 (2024). [CrossRef]

28. T. Li, Y. Luo, and K. Zhao, “Weak non-line-of-sight target echoes extraction without accumulation,” Opt. Express 31(22), 36209–36218 (2023). [CrossRef]

29. Y. Feng, X. Cui, and Y. Meng, “Non-line-of-sight imaging at infrared wavelengths using a superconducting nanowire single-photon detector,” Opt. Express 31(25), 42240–42254 (2023). [CrossRef]

30. M. O’Toole, D. B. Lindell, and G. Wetzstein, “Confocal non-line-of-sight imaging based on the light-cone transform,” Nature 555(7696), 338–341 (2018). [CrossRef]

31. D. B. Lindell, G. Wetzstein, and M. O’Toole, “Wave-based non-line-of-sight imaging using fast fk migration,” ACM Transactions on Graphics (TOG) 38(4), 1–13 (2019). [CrossRef]

32. X. Liu, I. Guillén, and M. La Manna, “Non-line-of-sight imaging using phasor-field virtual wave optics,” Nature 572(7771), 620–623 (2019). [CrossRef]

33. X. Liu, J. Wang, and Z. Li, “Non-line-of-sight reconstruction with signal–object collaborative regularization,” Light: Science & Applications 10(1), 198 (2021). [CrossRef]

34. X. Liu, J. Wang, and L. Xiao, “Non-line-of-sight imaging with arbitrary illumination and detection pattern,” Nat. Commun. 14(1), 3230 (2023). [CrossRef]

35. S. Xin, S. Nousias, K. N. Kutulakos, et al., “A theory of fermat paths for non-line-of-sight shape reconstruction,” in IEEE/CVF conference on computer vision and pattern recognition (2019), pp. 6800–6809.

36. C. Wu, J. Liu, X. Huang, et al., “Non–line-of-sight imaging over 1.43 km,” Proc. Natl. Acad. Sci. 118, e2024468118 (2021). [CrossRef]

37. G. Gariepy, F. Tonolini, and R. Henderson, “Detection and tracking of moving objects hidden from view,” Nat. Photonics 10(1), 23–26 (2016). [CrossRef]

38. S. Chan, R. E. Warburton, and G. Gariepy, “Non-line-of-sight tracking of people at long range,” Opt. Express 25(9), 10109–10117 (2017). [CrossRef]

39. C. Pei, A. Zhang, and Y. Deng, “Dynamic non-line-of-sight imaging system based on the optimization of point spread functions,” Opt. Express 29(20), 32349–32364 (2021). [CrossRef]

40. J. H. Nam, E. Brandt, and S. Bauer, “Low-latency time-of-flight non-line-of-sight imaging at 5 frames per second,” Nat. Commun. 12(1), 6526 (2021). [CrossRef]

41. M. O’Toole, D. B. Lindell, and G. Wetzstein, “Real-time non-line-of-sight imaging,” in SIGGRAPH Emerging Technologies (Association for Computing Machinery, 2018).

42. C. A. Metzler, D. B. Lindell, and G. Wetzstein, “Keyhole imaging: non-line-of-sight imaging and tracking of moving objects along a single optical path,” IEEE Transactions on Computational Imaging 7, 1–12 (2020). [CrossRef]

43. D. Faccio, A. Velten, and G. Wetzstein, “Non-line-of-sight imaging,” Nat. Rev. Phys. 2(6), 318–327 (2020). [CrossRef]

44. W. Becker, Advanced Time-correlated Single Photon Counting Applications, vol. 111 (Springer, 2015).

45. Z. Li, X. Liu, and J. Wang, “Fast non-line-of-sight imaging based on first photon event stamping,” Opt. Lett. 47(8), 1928–1931 (2022). [CrossRef]

46. C. Jin, J. Xie, and S. Zhang, “Richardson–lucy deconvolution of time histograms for high-resolution non-line-of-sight imaging based on a back-projection method,” Opt. Lett. 43(23), 5885–5888 (2018). [CrossRef]

47. D. O’Connor, W. Ware, and J. Andre, “Deconvolution of fluorescence decay curves. a critical comparison of techniques,” J. Phys. Chem. 83(10), 1333–1343 (1979). [CrossRef]

48. U. P. Wild, A. R. Holzwarth, and H. P. Good, “Measurement and analysis of fluorescence decay curves,” Rev. Sci. Instrum. 48(12), 1621–1627 (1977). [CrossRef]

49. N. Dey, L. Blanc-Féraud, C. Zimmer, et al., “3d microscopy deconvolution using richardson-lucy algorithm with total variation regularization,” Ph.D. thesis, INRIA (2004).

50. W. H. Richardson, “Bayesian-based iterative method of image restoration,” JoSA 62(1), 55–59 (1972). [CrossRef]

51. T. Carlsson, O. Steinvall, and D. Letalick, “Signature simulation and signal analysis for 3-d laser radar,” Sweden: Swedish Defence Research Agency, pp. 7–8 (2001).

Enhancing the spatial resolution of time-of-flight based non-line-of-sight imaging via instrument response function deconvolution

Abstract

1. Introduction

2. Methodology

2.1 Problem statement

2.2 Problem formation and mathematically solving

2.3 Posterior reconstruction

3. Numerical simulation and quantitative assessment

3.1 Numerical model

3.1.1 Energy attenuation of NLOS transfer

3.1.2 Observation model

3.2 Quantitative evaluation

3.2.1 Spatial resolution enhancement

3.2.2 Reconstruction within low SNR

4. Experimental

4.1 Hardware setup

4.2 Calibration of system’s IRF

4.3 Result and analysis

5. Discussion

6. Conclusion

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (8)

Equations (9)

Optics Express