Few-frame depth image reconstruction algorithm for GM-APD LiDAR based on three-dimensional spatiotemporal kernel density estimation

1 Introduction

Geiger-mode Avalanche Photodiode (GM-APD) LiDAR features single-photon detection sensitivity and picosecond-level temporal resolution, and has been widely applied to long-range target detection, autonomous driving, high-precision 3D topographic mapping and other fields [1, 2]. In practical outdoor scenarios, solar background light and atmospheric scattering introduce substantial environmental noise. When the Signal-to-Background Ratio (SBR) drops to a low level (e.g., SBR <0.1) [3], weak target echoes are easily overwhelmed by background noise [4]. Furthermore, for tasks such as fast beam scanning and airborne high-dynamic detection, the dwell time of the laser beam on each individual pixel is limited. Accordingly, the system can only acquire echo data within an extremely short time window consisting of dozens to hundreds of frames [5, 6]. Under the combined constraints of low SBR and few frames, the detected echo photons from targets are highly sparse and discrete, posing great challenges to accurate depth extraction.

One-dimensional temporal statistical algorithms based on single pixels are commonly adopted for depth extraction from single-photon data. For instance, the peak-picking method determines target depth by locating the maximum count in the photon time-of-flight histogram [7]. The matched filter method takes the laser emission waveform as a template and conducts cross-correlation with received signals [8]. The center-of-mass method identifies targets by calculating the weighted average of photon counts within the time window [9]. These one-dimensional algorithms achieve satisfactory depth extraction performance when sufficient frames are accumulated and the signal quality is high. Nevertheless, GM-APD detectors operating in the single-stop mode suffer from the inherent dead-time effect. Specifically, the first arriving photon within the time gate triggers an avalanche and deactivates the detector, which fails to respond to subsequent photons in the same detection cycle. This physical response mechanism gives rise to an asymmetric distribution of the photon count histogram along the time axis, characterized by a sharp rising edge and a slowly decaying tail [10, 11]. Under low SBR conditions, abundant background photons occupy most trigger opportunities, so that faint target peaks are severely masked by accumulated front-end noise. Meanwhile, the few-frame setting drastically limits the accumulation of valid photons, making it difficult to establish a reliable statistical distribution for individual pixels. In such cases, the peak-picking method is prone to false alarms caused by random background noise; the cross-correlation peaks of the matched filter degrade significantly; and the center-of-mass method suffers from range walk due to the asymmetric distribution. Therefore, pure single-pixel one-dimensional processing cannot simultaneously guarantee high target detection accuracy and ranging precision under low-SBR and few-frame conditions.

To address the problem that individual pixels fail to form valid peaks due to insufficient accumulated photons, spatial features of adjacent pixels are introduced for spatiotemporal joint filtering. The Sparse Poisson Intensity Reconstruction Algorithm (SPIRAL) proposed by Shin et al. [12] adopts Total Variation (TV) regularization and achieves promising contour recovery under extremely low signal quality. However, its iterative optimization based on high-dimensional spatiotemporal matrices incurs enormous computational overhead. In recent years, deep learning has been widely used in single-photon imaging. Lindell et al. developed a U-Net based convolutional neural network to process single-photon 3D data [13]. Several researchers adopted Residual Network (ResNet) and spatiotemporal Transformer models to extract spatial and temporal features from point clouds, achieving excellent denoising performance in specific scenarios [14–16]. Generally, deep learning models rely on large-scale training datasets with a distribution consistent with practical application scenarios. However, considering the high cost and operational difficulty of collecting real-world field single-photon data, large-scale public open datasets are still lacking. This greatly restricts the generalization ability of deep learning models when dealing with unknown complex outdoor environments and unstructured background interference. To reduce the dependence of algorithms on prior datasets and global iterative optimization, Kernel Density Estimation (KDE) based on pixel neighborhoods has been introduced into point cloud denoising [17, 18]. This method aggregates photons in local spatial regions to extract temporal peaks. Nevertheless, under low-SBR and few-frame conditions, conventional neighborhood aggregation is prone to misclassify intense background noise as target peaks. Moreover, traditional local maximum extraction strategies fail to computationally correct asymmetric ranging errors induced by the dead-time effect, resulting in substantial ranging deviations in the reconstructed depth images.

Aiming at the issues of severe background interference and large ranging errors in depth extraction under few-frame and low-SBR conditions, this paper proposes a depth image reconstruction algorithm based on three-dimensional spatiotemporal kernel density estimation (3D-STKDE). First, a spatiotemporal correlation weighted preprocessing strategy is adopted. Continuous weights are assigned according to the clustering characteristics of photons in local spatiotemporal neighborhoods to preliminarily suppress randomly distributed background noise. Second, a 3D Gaussian kernel with separable spatial and temporal components is constructed to perform convolution smoothing on discrete photons. The algorithm compensates for insufficient sampling along the single temporal axis under the few-frame constraint by exploiting local spatial correlation, and further derives a continuous 3D probability density distribution. In the depth calculation stage, an adaptive density threshold is used to eliminate residual noise and extract valid peaks. Afterwards, a density-weighted centroid calculation model is applied to effectively mitigate range walk caused by asymmetric distributions. Finally, the Inverse Distance Weighted (IDW) spatial interpolation algorithm is employed to fill local holes in depth images. Comparative field experiments demonstrate that in a field scenario with an extremely low SBR of approximately 0.0311 and only 100 accumulated frames, the proposed algorithm achieves a Target Recovery (TR) of 82.92% and a Mean Absolute Error (MAE) of 1.24 m. The experimental results fully verify the reliable reconstruction capability and high ranging precision of the proposed algorithm under data-limited conditions.

2 Proposed algorithm

This work extends the conventional processing scheme relying on one-dimensional temporal accumulation at individual pixels to a photon density clustering and calculation framework in the three-dimensional spatiotemporal domain. The overall workflow of the proposed depth image reconstruction algorithm is illustrated in Figure 1. Taking the raw sparse point cloud as input, the algorithm first suppresses background noise via spatiotemporal correlation weighted preprocessing. It then extracts spatially continuous density features using spatially and temporally separable 3D-STKDE. Density-weighted centroid computation is subsequently performed to correct ranging errors induced by the detector dead-time effect. Finally, high-fidelity two-dimensional depth images are generated through spatial interpolation.

2.1 GM-APD detection mechanism and peak degradation of few-frame data

GM-APD single-photon LiDAR achieves range measurement based on the time-of-flight (ToF) principle [19]. The system mainly consists of a Main Control and Signal Processing Unit, a 1064 nm Pulsed Laser, Transmitting and Receiving Optical systems, a Time Delay and Gate Control Unit, a 64 × 64 Gm-APD Detector, and a high-precision time-to-digital converter (TDC) [20]. A schematic diagram of its operating principle is presented in Figure 2.

During each detection cycle, the Main Control and Signal Processing Unit acts as the synchronization center. It sends a Trigger Signal to the 1064 nm Pulsed Laser to emit pulses toward the target through the Transmitting Optical system. Concurrently, it activates the Time Delay and Gate Control Unit. To suppress near-field backscattering interference, this control unit generates a Range Gating signal to activate the detector only within a specific temporal window, while simultaneously sending a Start Signal to initialize the TDC for timing.

When echo photons reflected by the target pass through the Receiving Optical system and arrive at the active Gm-APD Detector, pixels operating in Geiger mode generate avalanche currents upon absorbing a single photon. The resulting discrete logic pulses are immediately output to the TDC as a Stop Signal. The target range R can be calculated from the round-trip flight time of photons:where c denotes the speed of light in vacuum, and t is the time interval recorded by the TDC.

R = c · t / 2 , (1)

According to the above ranging principle, target depth can be accurately retrieved by accumulating abundant trigger timestamps and constructing a one-dimensional probability distribution histogram under ideal conditions. Nevertheless, for practical field detection, photon arrival events statistically obey a Poisson distribution. Meanwhile, environmental factors such as solar background light introduce massive random noise photons into the system. Under the combined influence of the detector’s intrinsic dead-time effect and single-stop operation mode, noise photons arriving prior to signal photons will trigger pixels preferentially. Consequently, the detector remains in its recovery phase when genuine target echoes arrive and cannot generate a response, which further gives rise to the photon pile-up effect [21].

To quantitatively characterize the aforementioned physical degradation process, the detection time gate is discretized into N time bins. Let λ i denote the average photon arrival rate (including both signal and background noise) within the i − t h bin. The theoretical trigger probability of receiving at least one photon in this bin is:

P D i = 1 − e − λ i . (2)

Due to the dead-time mechanism, the actual trigger probability P d i of this bin is strictly constrained by the cumulative survival probability of the preceding i − 1 bins that have not been triggered [22]:

P d i = P D i ∏ j = 1 i − 1 1 − P D j = 1 − e − λ i exp − ∑ j = 1 i − 1 λ j . (3)

As indicated by Equation 3, since GM-APD detectors typically operate in the single-stop mode, background noise photons randomly distributed in the front part of the time gate tend to trigger the detector in advance and force it into the dead-time state. This physical constraint leads to an exponential decline in the detector’s response probability to valid target signals in the middle and rear sections of the time gate. SBR is introduced as a fundamental metric to quantify the relative intensity between background noise and target echo. It is defined as the ratio of valid signal photon counts reflected by the target to background noise photon counts within the target-containing range gate.

Based on the above definition and probability model, Monte Carlo simulations are conducted to analyze the actual triggering probability distribution of a single pixel under different frame accumulation numbers. The simulation parameters are set as follows: the time gate contains 1000 bins, the target is located at the 600th bin, and the average SBR is 0.02. Figure 3 compares the actual triggering probability and theoretical distribution with accumulated frames of 10,000, 1000, 100 and 50, respectively.

For 10,000 accumulated frames (Figure 3a), the actual probability distribution agrees well with the theoretical model, and the target peak is clearly distinguishable. When the frame number decreases to 100 (Figure 3c) and 50 (Figure 3d), the combined effects of front-end noise preemption and insufficient sampling frames result in highly sparse and discrete one-dimensional histograms for individual pixels. Such severe degradation of statistical characteristics demonstrates that the true target depth cannot be accurately reconstructed merely relying on one-dimensional temporal statistics of single pixels under low SBR and few-frame conditions. For this reason, the proposed algorithm extends one-dimensional temporal processing to the three-dimensional spatiotemporal domain. By exploiting spatial correlation among adjacent pixels, the shortage of accumulated samples along the single temporal axis is compensated, enabling high-fidelity depth reconstruction from sparse data.

2.2 Spatiotemporal correlation weighting and three-dimensional kernel density estimation

Under low SBR and few-frame conditions, numerous background noise points induced by ambient light are randomly distributed across the full field-of-view in the 3D discrete dataset. Direct implementation of global 3D spatiotemporal clustering is highly susceptible to random noise and will lead to substantial errors in depth estimation. Physically, background noise photons follow a uniform Poisson distribution across the spatiotemporal domain and exhibit strong randomness and discreteness. In contrast, echo photons reflected from real targets are jointly governed by the spatial continuity of the target itself, the point spread function (PSF) of the optical receiving system [23], and the temporal profile of the laser pulse. These factors collectively yield prominent clustering features of photons within local spatiotemporal neighborhoods [24].

According to the above spatiotemporal distribution properties of signals, a spatiotemporal correlation weighting preprocessing scheme is introduced prior to the core clustering algorithm to denoise the raw point cloud data. For any candidate photon event in the dataset, the number of clustered photons within its local spatiotemporal sliding window Ω x , y , t is first counted N. A nonlinear activation function is then adopted to assign a confidence weight [0,1] ranging within to the corresponding photon W c :

W c = 1 − exp − max 0 , N − λ b g β . (4)

In the formula, λ b g denotes the expected baseline that quantifies the local background noise intensity, while β is a noise-dependent smoothing parameter governing the weight decay rate. Consequently, a larger β is applied under lower SBR regimes to enforce stronger smoothing for effective noise outlier suppression. The value of β is determined by the background noise level via a pre-calibrated mapping. This mathematical model indicates that isolated or sparsely distributed photons, where the local photon count N ≤ λ bg is low, predominantly originate from background noise. Their weights are suppressed and approach zero accordingly. For photons with evident local clustering features, their weights rise in a nonlinear manner. Different from hard threshold truncation, this scheme adopts continuous weight mapping. It effectively suppresses discrete noise while preserving weak target echo photons that are submerged by noise to the greatest extent.

After weight assignment for discrete photon events, continuous density distribution in the three-dimensional spatiotemporal domain is further calculated to extract core depth features. For single-photon LiDAR systems, the spatial dispersion of target echoes is mainly determined by the point spread effect of the receiving optics, while the temporal span is governed by the inherent laser pulse width. Given the essential physical differences between spatial and temporal dimensions, a spatially and temporally separable 3D spatiotemporal Gaussian kernel is constructed to perform convolution smoothing on the weighted discrete photons:

K x , y , t = exp − x 2 + y 2 2 σ s 2 − t 2 2 σ t 2 . (5)

The spatial kernel bandwidth σ s and temporal kernel bandwidth σ t are determined according to the optical PSF of the receiver optics and the physical laser pulse duration of the system, respectively. Under few-frame and low SBR conditions, data-driven adaptive bandwidth estimation can become unreliable due to photon sparsity. By adopting deterministic physical priors as fixed kernel bandwidths, the algorithm maintains clustering stability and reduces density estimation bias caused by random noise. Driven by real system parameters, this smoothing procedure converts discrete photon timestamps with continuous weights into a continuous three-dimensional probability density distribution. As illustrated in Figure 4, the proposed clustering method fully exploits spatial correlation of local neighborhoods. It effectively compensates for insufficient sampling along the single temporal axis under few-frame conditions, enabling accurate extraction of true depth peaks from sparse discrete point clouds.

2.3 Density-weighted centroid calculation and spatial interpolation for depth images figures

After obtaining the continuous three-dimensional probability density distribution, peak coordinates are further calculated to determine the precise depth of targets. To suppress interference from residual ambient background noise, an adaptive density threshold is formulated based on the statistical property of Poisson-distributed noise, whereby the standard deviation of noise fluctuations grows proportionally to the square root of accumulated frames N [25]:where ρ ¯ bg represents the estimated average density of background noise, and k is a threshold proportional coefficient dynamically configured according to the actual background noise level to strike a balance between noise removal and signal preservation. During peak extraction, subsequent depth calculation is only performed on the set of candidate pixels whose kernel density exceeds this threshold.

ρ t h = k · N · ρ ¯ b g , (6)

In single-photon detection, the detector dead-time effect generally leads to an asymmetric photon count histogram along the temporal axis, featuring a sharp rising edge and a slowly decaying tail. When dealing with such asymmetric distributions, the conventional arithmetic mean method tends to shift the estimated center toward the trailing edge. This centroid shift introduces severe ranging errors, directly contributing to the range walk phenomenon [26]. As explicitly illustrated in Figure 5, the raw photon counts exhibit this asymmetric behavior, drawing the arithmetic mean away from the ground truth depth. To mitigate this positioning deviation, we adopt the continuous kernel density values derived from 3D-STKDE as confidence weights and implement density-weighted centroid calculation within the peak neighborhood. For the valid time-of-flight candidate set t i and corresponding kernel density values ρ i at a given pixel x , y, the depth value D x , y is calculated as:where γ denotes the nonlinear amplification factor constrained by the dead-time characteristics of the detector. By assigning larger nonlinear weights to the high-density core region of peaks, the proposed model effectively compensates for asymmetric range walk deviations induced by the dead-time effect.

D x , y = ∑ i t i · ρ i γ ∑ i ρ i γ , (7)

After depth calculation, local empty pixels frequently appear in the reconstructed two-dimensional depth images due to the lack of signal photons, which is inevitable under the combined constraints of low SBR and few accumulated frames. To acquire intact depth images, the IDW spatial interpolation method is adopted for post-processing and hole filling [27].

Let the coordinate of a depth hole to be filled in the 2D depth image be x i , y i. The algorithm first searches its local neighborhood to collect all valid pixels with solved depth values. For an arbitrary valid pixel x 0 , y 0 within this neighborhood with a reconstructed depth of D i, the physical Euclidean distance d i between this pixel and the center of the depth hole is defined as:

d i = x 0 − x i 2 + y 0 − y i 2 . (8)

Based on the general prior that depth correlation between adjacent pixels decays with increasing physical distance, spatial weight coefficients w i are assigned to all valid neighboring pixels:where p is the distance power exponent. The final filled depth value D fill of the hole pixel is obtained via normalized weighted summation of valid depths in the neighborhood:

w i = 1 d i p , (9)

D f i l l = ∑ w i · D i ∑ w i . (10)

This method fully exploits the local spatial continuity of target surfaces. It effectively fills depth holes without global iterative optimization, and further improves the overall integrity of reconstructed depth images.

3 Experiments and analysis

3.1 Definition of evaluation metrics

To objectively quantify the depth image reconstruction performance of the proposed algorithm under few-frame and low SBR conditions, Target Recovery and Mean Absolute Error are selected as the core evaluation metrics in this work.

TR is used to evaluate the completeness of the target geometric morphology [28] recovered by the algorithm, defined as the ratio of successfully reconstructed valid pixels to the total number of true target pixels in the reference. Its formula is given by:where N r e c is the number of successfully reconstructed valid pixels, and N t o t a l is the total number of true target pixels in the ground-truth depth image. A pixel at spatial coordinates x , y is considered to be effectively reconstructed if the deviation between its reconstructed depth D r e c x , y and the corresponding ground-truth depth D g t x , y satisfies:where L p denotes the laser pulse width. A higher TR value indicates more complete preservation of the target morphology.

T R = N r e c N t o t a l × 100 % , (11)

∣ D r e c x , y − D g t x , y ∣ ≤ L p , (12)

MAE is adopted to quantify the actual ranging error of correctly reconstructed pixels, which directly reflects the ability of the algorithm to correct the range walk phenomenon caused by asymmetric distributions [29]. Its formula is defined as:where V is the set of coordinates of all valid reconstructed pixels, and Δ R is the physical distance resolution corresponding to a single bin, given by Δ R = c · Δ t / 2 (where Δ t is the resolution of TDC). In this work, MAE is converted to the actual physical distance (in meters). A smaller MAE value indicates higher reconstruction accuracy.

M A E = 1 N r e c ∑ x , y ∈ V D r e c x , y − D g t x , y · Δ R , (13)

3.2 Simulation experiments and result analysis

To evaluate the depth extraction performance of the algorithm under controlled conditions with absolute distance references, this paper adopts a Monte Carlo model to generate simulation data for a 64 × 64 pixel array GM-APD LiDAR system. Two targets at different depths are set in the scene, corresponding to the 600th bin and the 1200th bin respectively. The overall SBR is regulated by adjusting the generation probabilities of signal photons and background noise points. Three groups of test datasets are finally obtained with SBR values of 0.08, 0.05 and 0.02. Figure 6 illustrates the ground-truth depth image of the simulated scene, where pseudo-colors are used to represent specific depth bins, as well as the one-dimensional photon triggering histogram acquired by global pixel accumulation under the conditions of SBR = 0.02 and 200 accumulated frames.

In the simulation model, the temporal profile of the laser echo is characterized by a Gaussian pulse with a full width at half maximum (FWHM t) of 11.56 ns, corresponding to 5.78 bins under the system’s temporal resolution of 2.0 ns/bin. To ensure a fair and physically consistent comparison, the parameter configurations and underlying mathematical models of all evaluated algorithms are explicitly defined. For the proposed 3D-STKDE method, the spatial kernel bandwidth σ s is configured as 0.8 pixels to match the optical PSF (FWHM s ≈ 1.88 pixels, yielding σ s = 1.88/2.355 ≈ 0.8 pixels). The temporal kernel bandwidth σ t is set to 2.45 bins, which is explicitly matched to the theoretical signal standard deviation (σ t = FWHM t/2.355 ≈ 2.45 bins) to perform matched-filtering temporal noise suppression. Meanwhile, the spatiotemporal sliding window, which covers the effective support domain of the Gaussian kernel, is set to 3 × 3, and the neighborhood size for spatial interpolation in the hole-filling stage is specified as 5 × 5. The noise-smoothing parameter β is selected according to the background noise level via a pre-calibrated mapping, with values of 1.5, 2.0, and 2.5 corresponding to SBR levels of 0.08, 0.05, and 0.02, respectively.

For the comparative methods, the specific models and parameter implementations are as follows:

• Single-pixel peak-picking method [7]: This method identifies the target depth by directly locating the time bin with the maximum photon count in the raw one-dimensional histogram. To reflect the original data statistics, no additional thresholding or temporal smoothing is applied.
• Matched filter method [8]: This approach utilizes a normalized reference template S t to perform cross-correlation with the received histogram t :

D M F = arg ⁡ max t ∑ τ H τ S τ − t . (14)

To ensure a consistent physical constraint, the normalized reference template S t is modeled as a Gaussian function matching the echo pulse shape with a FWHM of 11.56 ns, corresponding to 5.78 bins under 2.0 ns/bin resolution, i.e., σ ≈ 2.45 bins. where L Y | z represents the data fidelity term utilizing the observed photons Y. The T V z penalty weight η is set to 0.15 to balance noise suppression and structural fidelity, and the maximum number of iterations is limited to 50. Figure 7 visualizes the depth reconstruction results of the single-pixel peak-picking method, matched filter method, high-order SPIRAL and the proposed algorithm under various SBR levels.

• SPIRAL method [12]: This algorithm recovers the depth map z by solving a penalized negative Poisson log-likelihood objective function with Total Variation (TV) regularization:

z ^ = arg min z − L Y | z + η T V z } . (15)

From a visual perspective, the peak-picking method fails to reconstruct valid target structures across all SBR conditions, as it cannot distinguish target peaks from background noise using sparse counting data. The matched filter can roughly identify the main building area as the frame number increases, yet substantial shot noise remains in the background and around building edges. Although the SPIRAL algorithm achieves moderate smoothing for sparse data, obvious structural fractures and pixel holes appear on the reconstructed building surfaces. In contrast, benefiting from the effective fusion of spatiotemporal correlation via 3D-STKDE, the proposed algorithm effectively suppresses background noise and accurately extracts target contours even in the high-noise environment with SBR as low as 0.02. It exhibits superior noise immunity and depth extraction accuracy.

Table 1 summarizes the quantitative metrics of all competing algorithms at different SBR values. The results reveal that the reconstruction performance of all methods degrades to varying degrees as SBR decreases. Nevertheless, the proposed algorithm outperforms all comparative approaches in terms of TR and MAE. At SBR = 0.05, the proposed algorithm achieves a TR of 91.85% and an MAE of 0.59 m. When SBR drops further to 0.02, its TR still remains at 81.79%, which is 21.46% higher than that of the SPIRAL algorithm, while the MAE is maintained at 0.66 m. The quantitative results demonstrate that the proposed algorithm can effectively suppress ranging errors and realize high-precision depth image reconstruction under few-frame and low-SBR conditions.

TABLE 1. Quantitative metrics of different algorithms under various SBR levels in simulation.

3.3 Field experiments and result analysis

To further validate the effectiveness of the proposed algorithm on practical hardware platforms and in complex optical environments, field experiments are conducted with a single-photon LiDAR system. All tests are performed under natural daylight with clear weather, and the atmospheric visibility is greater than 10 km. In such conditions, Intense stray solar light constitutes the primary source of background noise for the system, which provides a realistic test environment to evaluate the algorithm performance under low SBR.

A focal-plane array GM-APD LiDAR system is employed in the experiments, and a photograph of the actual experimental setup is shown in Figure 8. Its key hardware specifications are presented as follows. At the transmitting end, a pulsed laser with a central wavelength of 1064 nm is used, featuring a pulse width of 11.56 ns and a repetition rate of 25.0 kHz. At the receiving end, the system integrates a 64 × 64 GM-APD detector array equipped with a high-precision TDC (time resolution Δ t = 2.0 ns/bin), together with a receiving optical lens with a focal length of 300 mm and an F-number of 5.6.

Two typical scenes with different SBR values and spatial distributions are selected for testing. As shown in Figure 9, Scene one targets residential buildings approximately 440 m away from the system, which contain multiple depth layers due to different object positions. The average SBR of this scene is calculated as 0.0635. Scene two comprises a residential building 500 m away and a quadrotor drone hovering in front of the building. Affected by strong solar radiation and atmospheric scattering, this scene suffers from more severe background noise, with an average SBR of 0.0311.

To quantitatively evaluate the reconstruction quality in subsequent analyses using predefined metrics, such as TR and MAE, reliable reference ground-truth depth images were generated for the two scenes, as illustrated in Figures 9b,d. Specifically, we accumulated 20,000 consecutive frames of raw data acquired by the system to guarantee the statistical dominance of the true signal peak [30], allowing the conventional peak-picking method to accurately extract the raw ToF depth. Subsequently, a standard morphological closing operation with a fixed 3 × 3 template was applied to restore the 2D spatial connectivity of the targets and repair minor holes. This morphological step operates exclusively on the 2D spatial plane without altering the original ToF depth values. Consequently, the reference depths remain strictly anchored to the raw photon statistics, preserving the objectivity of the evaluation metrics and completely avoiding potential bias from spatial smoothing.

To evaluate algorithm performance under few-frame conditions, datasets with 40, 70, 100, and 120 accumulated frames are extracted. Comparative methods (Peak, Matched Filter, and SPIRAL) retain the same model configurations as in the simulation section. For the proposed 3D-STKDE, the spatial and temporal kernel bandwidths (σ s = 0.8 pixel, σ t = 2.45 bins) as well as the 3 × 3 spatiotemporal window are strictly maintained to align with the system physics (FWHM t = 11.56 ns at 2.0 ns/bin resolution). Based on the measured ambient SBR, β is set to 1.8 for Scene one and 2.2 for Scene 2. In the hole-filling stage, a 5 × 5 spatial neighborhood is applied, effectively bridging the 2-to-3-pixel hole spans typical of few-frame data without smoothing out subtle target boundaries.

The reconstruction results are illustrated in Figure 10. Limited by the sparsity of photon trigger events, the single-pixel peak-picking method generates numerous discrete noise points at all frame numbers and fails to reconstruct clear depth contours of buildings. The matched filter can roughly distinguish the main building areas as the frame number increases, while discrete noise still remains in the background and around building edges. Although the SPIRAL algorithm achieves certain smoothing effects on sparse data, obvious structural fractures and pixel holes exist on the reconstructed building surfaces.

In comparison, the proposed algorithm maintains stable reconstruction performance across different frame numbers. It can effectively recover layered contours of various buildings even with only 40 frames. When the accumulated frame number increases to 120, local holes on target surfaces are fully filled, and the overall geometric structures and boundaries remain complete. The proposed method achieves higher reconstruction completeness than all comparison algorithms.

Figure 11 presents the reconstruction results of Scene 2 with a much lower SBR of 0.0311. This scene includes a drone with elaborate structures and building backgrounds. The results of the peak-picking method and matched filter degrade drastically, and the drone is almost completely obscured by background noise. The SPIRAL algorithm exhibits large-scale pixel missing in target regions, particularly on wall surfaces and the drone body. In contrast, the proposed algorithm well preserves the morphological integrity of both the drone and background buildings. With the increase of accumulated frames, background noise is thoroughly removed. The hovering position of the drone and planar features of buildings are accurately restored, demonstrating superior performance over the other three methods.

Table 2 lists the quantitative metrics, including TR and MAE, of all algorithms under different frame numbers for the two field scenes. For Scene 1 (SBR = 0.0635), the proposed algorithm achieves a TR of 71.68% and an MAE of 0.56 m using 40 frames. When the frame number rises to 120, the TR increases to 90.47% and the MAE decreases to 0.49 m. By contrast, the SPIRAL algorithm only reaches a TR of 75.66% at the same frame number.

TABLE 2. Quantitative metrics of different algorithms with various accumulated frames for Field Scene one and Scene 2.

The strong robustness of the proposed algorithm is further verified in Scene 2 with a lower SBR of 0.0311. With 100 frames, the proposed algorithm obtains a TR of 82.92% and an MAE of 1.24 m. Meanwhile, the SPIRAL algorithm only achieves a TR of 60.16% due to extensive pixel missing in target areas. Restricted by severe background noise, the matched filter yields a TR of merely 43.50%. Both methods have an MAE above 1.62 m. The quantitative analysis verifies that the proposed algorithm possesses excellent target contour recovery capability and high spatial positioning accuracy under few-frame and low-SBR conditions.

In addition to reconstruction accuracy, computational efficiency is also a considered metric for practical depth detection systems. For fair comparison, all evaluated algorithms were executed on a computing platform equipped with an Intel Core i5-13400F CPU and 32 GB of RAM. Table 2 presents the execution time of each algorithm evaluated at 100 accumulated frames. While the foundational methods (Peak and Matched Filter) require minimal execution time (e.g., 0.003 s and 0.045 s in Scene 1), their quantitative performance degrades in low-SBR scenarios. Among the compared methods, SPIRAL relies on global regularized iterative minimization, which imposes a computational complexity of O I · V ⁡ log ⁡ V, where I represents the number of required iterations. Its execution time increases from 2.521 s in Scene one–3.228 s in Scene 2, as lower SBR conditions typically require more iterations to converge. In contrast, the proposed 3D-STKDE algorithm is non-iterative. Its core operations consist of photon mapping with a complexity of O N and local 3D spatial convolution with a complexity O V · M of, where N, V, and M represent the number of detected photons, total spatiotemporal voxels, and the fixed-size local kernel volume, respectively. Consequently, the proposed method maintains an execution time of 0.258 s in Scene one and 0.282 s in Scene 2, demonstrating relatively stable computational efficiency across different noise levels while ensuring contour recovery and positioning accuracy.

4 Discussion

This paper proposes a depth image reconstruction algorithm based on 3D spatiotemporal kernel density estimation with separated spatial and temporal domains to extract target depth for Geiger-mode single-photon LiDAR under few-frame and low-SBR conditions. The method suppresses background noise using spatiotemporal correlation weighted preprocessing, and applies a separated 3D Gaussian kernel to obtain a continuous probability density distribution. Then, an adaptive density threshold removes residual noise, and a density-weighted centroid calculation mitigates the asymmetric ranging error caused by the detector dead-time effect. Finally, local spatial interpolation repairs depth voids.

Both simulation and field experiments validate the proposed algorithm. In Monte Carlo simulations (SBR = 0.02), the algorithm achieves a TR of 91.85% and an MAE of 0.66 m. Its TR is 21.46% higher than that of the SPIRAL algorithm. In field tests with severe solar background interference (average SBR ≈0.0311, 100 effective frames), conventional single-pixel one-dimensional statistical methods produce extensive discrete noise, and existing spatial reconstruction algorithms cause structural fractures. In contrast, the proposed method preserves the geometric shapes and boundary features of the targets, yielding a TR of 82.92% and an MAE of 1.24 m.

Compared with previous studies, the proposed method shows clear differences in its mechanism. Conventional single-pixel one-dimensional methods ignore spatial correlation among adjacent pixels and are easily affected by background noise. Global iterative optimization algorithms (such as SPIRAL) use TV regularization, which tends to cause structural fractures and severe pixel loss when effective photons are extremely limited. The proposed 3D-STKDE method sets the spatial and temporal kernel bandwidths based on the point spread function and laser pulse width. This design preserves the physical clustering of signal photons without global iterative optimization. Additionally, the density-weighted centroid calculation compensates for the backward shift of depth estimation caused by the dead-time effect, effectively reducing the MAE.

In summary, the proposed algorithm provides a practical solution for 3D depth image reconstruction under low SBR and few-frame conditions. However, the current implementation employs fixed spatial-temporal bandwidths and a constant spatial interpolation window. Under elevated background illumination or for targets with small spatial footprints (e.g., spanning only a few pixels), the limited availability of neighboring signal photons can weaken spatial correlation gains, which may lead to a localized decrease in TR. Similarly, in regions with complex depth discontinuities or sharp density variations, this static mechanism may introduce slight edge blurring and a minor increase in MAE. Future work will investigate adaptive bandwidths and interpolation window sizes based on local photon density to improve reconstruction accuracy across dynamic scenarios. Given the non-iterative nature and computational efficiency of the algorithm, hardware acceleration on parallel platforms (such as FPGA) will also be explored for real-time applications.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Qiyin Fang, McMaster University, Canada

Reviewed by: Siyuan Wang, Chinese Academy of Sciences (CAS), China

Tamer Elkhatib, Cairo University, Egypt