1.

Introduction

Time-of-flight (ToF) lidar is a powerful technique for collecting accurate three-dimensional (3-D) information about the surrounding area. After the emission of a laser pulse, the signal received at the detector is the aggregated contributions from all reflecting surfaces within the instantaneous field of view (IFOV) of the sensor. By performing time-resolved measurements, the distance to the individual surfaces can be measured and the data represented as a 3-D point cloud, which can be used for higher level tasks, such as object recognition or mapping. The range information improves the chance of separating and locating objects covering only part of the IFOV. The smaller the spatial separation between surfaces, the higher range resolution is needed. Time-correlated single-photon counting (TCSPC) is a lidar technology that can provide very high-range resolution data and hence separate closely spaced objects.¹ TCSPC lidar emits very weak laser pulses and detects the time of arrival of the first photon to be absorbed by the detector with very high temporal accuracy. By performing a large number of measurements, a histogram of arrival times is assembled, which represents the laser radar response of the scene within the IFOV. With a large enough number of measurements, detections will accumulate in the histogram bins corresponding to the range of a reflecting surface, making single surfaces appear as clear peaks above the shot noise of the background detections. The shape of the peak is determined by the instrument response function (IRF), which depends on factors such as the laser pulse width, detector timing jitter, and trigger jitter. Signals from closely spaced surfaces will overlap and may appear as a single, broadened peak in the histogram. In this article, we study signal processing approaches for detecting multiple surfaces at different ranges within one IFOV, with special emphasis on the performance for separating closely spaced surfaces whose contributions to the histogram overlap. The application in mind is detection and segmentation of partly occluded objects, for example, objects behind vegetation. For this application, the TCSPC system must be able to detect and localize multiple surfaces at different ranges within one IFOV. The number of surfaces is unknown a priori and must be estimated in the signal processing.

Several techniques for surface detection in TCSPC lidar data have been proposed in the literature. Most of them try to minimize the number of detected photons, and hence the integration time needed to determine the distance to a single reflecting surface. Examples of algorithms to detect multiple reflecting surfaces have also been published. There are various approaches based on Markov chain Monte Carlo simulations reported in the literature.^2–4 These algorithms were designed for low photon counts, i.e., short integration times or targets with low reflectance. Shin et al.⁵ proposed another peak detection algorithm (SPISTA) for low photon count but based on linear optimization. It was claimed⁵ that it outperforms a mixture of Gaussian (MoG)-based estimator that is using a greedy histogram-data-fitting strategy. There are also approaches based on local maximum detection and cross correlation,⁶ but these approaches are generally not intended for resolving multiple surfaces close to each other in range as the local maxima are distorted when surfaces are so close together that the IRFs overlap.

In our work, we assume longer integration times than what has been previously reported in the literature. The goal of this paper is to make true object surfaces stand out from the background noise. The increased signal-to-noise ratio (SNR) of the histograms is used to detect multiple closely spaced surfaces at different distances. The cost of the increased capability of detecting objects obscured by, e.g., vegetation is longer measurement time. The optimal balance between short integration times corresponding to high area coverage and high SNR allowing detection of multiple surfaces closely spaced in range is obviously application dependent.

In this work, we study a number of peak detection approaches and assess their performance from different perspectives; the possibility of detecting double surfaces within one IFOV, geometrical accuracy of the detections, the number of outliers, and how reduced signal quality affects the performance.

The remainder of this article is organized as follows: in the next section, a number of algorithms for preprocessing and peak detection are described, along with the tested peak detection approaches (mixes of the described algorithms). In Sec. 3, the setup of the tests is described, including a short description of the TCSPC system and the reference measurements performed. Section 4 describes the procedure for measuring range accuracy and in Sec. 5 results are presented. A discussion of the results is included in Sec. 6 and conclusions are given in Sec. 7.

2.

Peak Detection Approaches

A fundamental requirement on a peak detection approach for our intended application is that no a priori assumption of the number of surfaces in the IFOV is made. It should also meet the following quality criteria:

In this work, we use the term “approach” for denoting an algorithm for detecting and locating peaks in TCSPC histograms, possibly preceded by an algorithm for preprocessing such histograms.

Two types of preprocessing algorithms are considered: matched filter (MF) and Lucy–Richardson deconvolution (LR).⁷ Both algorithms use the impulse response function (IRF) of the measurement system but for completely different purposes. MF is used for suppressing noise and improving range precision and is obtained by convolving the TCSPC histogram with a time-reversed version of the IRF. Another advantage of MF is that the peak generated by a single surface after filtering is symmetric and hence more easily described by a fitting function. The purpose of the LR deconvolution of TCSPC data is to decompose histograms into modes corresponding to the individual surfaces. LR involves iterative convolutions with both the IRF and its time-reversed version, but still it is computationally fast provided that the number of iterations is kept low. A side effect of LR is that it amplifies noise rather than reducing it, which can be mitigated through imposing constraints on the deconvolved signal or removing spurious points in a postprocessing step.

For the problem of detecting peaks in histograms, the following algorithms are considered:

The local maximum detection (LM) is performed in two steps. First, the range bins corresponding to LM values are located, using thresholds for height and prominence to reduce the number of false alarms. Then, the position of each peak is estimated with subsample accuracy through quadratic interpolation using the LM range bin and its immediate neighbors.

In the EMGM algorithm, data are modeled as a mixture of Gaussian basis functions. This is different compared with the other fitting methods investigated in this paper, which use representations of the IRF. The Gaussian basis functions allow for using computationally efficient standard implementations. The algorithm is based on expectation-maximization, which is an iterative approach for finding the maximum likelihood estimates of parameters in a statistical model. It has a nonconvex cost function, which means that only finding a local minimum can be guaranteed when used for fitting of multimodal distributions. Thus, good initialization is needed.

The SPISTA algorithm is essentially a reformulation of the multiple peak detection problems as an L1 norm problem with a global optimum. It models the input signal and background noise as Poisson distributions. In SPISTA, the single-photon detection statistics from multiple reflectors within a pixel is described as a multidepth profile and expressed as a sparse signal. The multidepth estimation problem from single-photon observations is formulated as a convex optimization problem by combining the statistics of photon detection data with sparsity of multidepth profiles. Linear minimization is known to be less affected by noise and outliers compared with least squares minimization and also computationally faster. SPISTA has a convex cost function, which is different from the EMGM and IRF-LS algorithms.

The IRF-LS algorithm minimizes the difference between the measured data and a linear combination of IRFs in a least square sense, recursively adding more surfaces until no further improvement is found. The approach is good at accurately detecting high-amplitude peaks with high resolution but is relatively computationally demanding and prone to adding low-amplitude surfaces from imperfections of the IRF description. Since, in our implementation, it never removes surfaces as some published MCMC-based methods do,³ it is more prone to overfitting. Furthermore, it has a nonconvex cost function.

The FRI algorithm reformulates the time-of-flight superresolution problem as a parameter estimation problem. Based on the theory of finite rate of innovation (FRI) sampling,¹⁰ a bandlimited approximation of the sampling kernel function is used, which allows for formulating the problem as a linear least square problem in the frequency domain. The surface positions are determined through solving an eigenvalue problem, with a threshold applied on the eigenvalues to determine the number of surfaces in the final result. The method is noniterative and hence computationally attractive. It has shown promising results for both superresolution and resolving overlapping peaks and was reported to provide better estimates than an MCMC-based method and at lower execution time.⁹

The algorithms are combined into the following peak detection approaches:

3.

Tests

3.1.

Lidar System

The TCSPC lidar data subject to peak detection was recorded using the system described in detail by Henriksson et al.¹ The lidar system uses a $1.55\text{-}\mu \mathrm{m}$ fiber laser source with 22 ps pulse length and pulse repetition frequency of 8 MHz. The optical transceiver is a 500-mm focal length Cassegrain mirror telescope, where the laser light is emitted via a mirror in front of the central obscuration to provide a monostatic configuration. The collected photons are focused on a graded index fiber and, after bandpass filtering, focused on an InGaAs single-photon avalanche detector. The times of detections, measured relative to the emissions of laser pulses, were recorded with 32-ps precision using a Picoquant Picoharp 300. The detector field of view was continuously swept row by row at $0.04\mathrm{deg}/\mathrm{s}$ during the measurement. The signal is subsequently divided into pixels for small angular intervals during the postprocessing. Due to performing the measurements while sweeping, there will be a spatial broadening of objects, but, with small angular intervals for each pixel, this decrease in spatial resolution can be kept at an acceptable level.¹ The measurements are aimed at accurately finding targets in vegetation or behind camouflage and hence involve longer integration times than in Refs. 23.4.5.–6. The effective integration time per pixel was $126\mathrm{ms}/\text{pixel}$, which means that $\sim {10}^6$ measurements were performed per pixel. The laser power, however, was low, meaning that on average 1120 detections were recorded per pixel. The background rate, caused by detector dark counts and solar illumination, was estimated to be $0.22\text{detections}/(\text{range}\mathrm{b}\mathrm{i}\mathrm{n}\times \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l})$ or in total $870\text{detections}/\text{pixel}$. Approximately, 78% of the detections are thus uncorrelated to the laser and only contribute noise to the histograms.

The measured IRF was 392-ps full width at half maximum (FWHM). The fitting function described in Ref. 1 consisting of a Gaussian peak with exponential tails was used to describe the IRF in the signal processing. As seen in Fig. 1, the general shape is well described, but some high-frequency features near the peak leave a residual error.

Fig. 1

The measured IRF and the fitting function used in the signal processing. The FWHM of the IRF is 392 ps.

3.2.

Target Setup and Measurements

Measurements were performed against a reference target with holes behind which metal plates were placed at different distances: 33, 53, 73, and 93 mm (see Fig. 2). In the further analysis the front surface is designated as $S_1$, and the back plates are designated as $S_2$, $S_3$, $S_4$, and $S_5$, in order of increasing distance from the front surface. The incidence angle deviated from the normal of the target surfaces with about 7-deg horizontally and 4-deg vertically, respectively.

Fig. 2

The reference target.

Scanning across the target means that the histogram of every pixel corresponds to an IFOV in a different position on the target. We are, especially, interested in the cases where the IFOV covers the edge of a hole in the front plate so that the histogram is composed of photons reflected both from the front surface and plate behind the hole. The four different distances 33 mm ($S_1$ to $S_2$), 53 mm ($S_1$ to $S_3$), 73 mm ($S_1$ to $S_4$), and 93 mm ($S_1$ to $S_5$) allow us to investigate what distance between reflecting surfaces within the IFOV can be resolved. Figure 3 shows examples of collected histograms for a point at the edge of each hole, in addition to a model consisting of a linear combination of two IRFs fitted to the measured data. The histograms shown in Fig. 3 are created using a different dataset with higher SNR than the data used for peak detection analysis and only serve to illustrate the phenomenon. At different positions at the edges of the holes varying portions of the IFOV will be on the front and back plates, leading to varying proportions of the signal belonging to the first and second peaks in the histogram. The effects of different obscuration coverages or different reflectivities on the obscuration and target will have a similar effect. If the obscuration is caused by foliage or other vegetation, it will probably be in the form of several small surfaces with a distribution in range to the target.

Fig. 3

Measurements at four different points at edges of each hole in the characterization target. The black lines are measured data and the red dotted lines are two fitted IRF peaks. Point number 1 includes surfaces $S_1$ and $S_2$, point 2 $S_1$ and $S_3$, point 3 $S_1$ and $S_4$, and point 4 $S_1$ and $S_5$.

4.

Measuring Relative Range Accuracy

To measure the range accuracy, the detected peaks were converted into 3-D points in Euclidean space using the estimated range values in combination with the angle obtained through interpolation, assuming constant angular velocity between the start and stop angle values. The algorithms may introduce a slightly different bias to the absolute range estimates. However, we are only interested in the relative accuracy within the point cloud obtained by the respective approach. As a consequence, the relative position of the reference target in each of the respective datasets has to be determined first. This is achieved through fitting of a plane to the extracted 3-D data points. First, the orientation of the front plane is estimated using RANSAC¹¹-based plane fitting. For this purpose, a point cloud was selected after empiric testing with different parameters and visual inspection of both data and the resulting plane to ensure a good fit. The fitted plane was then used as an input to the assessment of performance for all other approaches, yet shifted in range to compensate for the differences in range bias.

The position of the front plane implicitly gives the position and orientation of the back plates, which are parallel to the front plate and at known distances (33, 53, 73, and 93 mm). Using the equation of the plane, the normal distance from each point to the front surface is then computed. Based on that, the points are segmented into five clusters corresponding to the front and the four back plates (Fig. 5).

A point is regarded as belonging to the nearest surface provided that the distance is small enough. This is to avoid having outliers far from the surfaces affecting the range accuracy estimation. We introduce a binary weight parameter $w_{p,S_i}$

5.

Results

In this section, we present the results obtained from the tests, studying the performance of the approaches regarding the criteria listed in Sec. 2 such as detection of double surfaces, range precision, and the number of outliers. Furthermore, we study how the performance of the respective approaches varies when reducing the signal strength.

Figure 6 shows the points along a line transect across the reference target detected with the respective approach. Ideally, all points should lie on the solid lines representing the position of the planar plates of the reference target relative to the data. The color of the dots shows the estimated amplitude of the peak. All approaches except EMGM-MF produce the same overall, expected pattern with stronger peaks for single surfaces and weaker peaks in the double surface regions. Though, the LR-LM approach shows somewhat larger intensity variations between surfaces and is not as reliable in terms of drawing conclusions about the reflectivity of the surfaces. This is an expected effect of the LR preprocessing step as it amplifies noise.

Fig. 4

Range histograms from a pixel within a double surface region computed with 100%, 50%, and 25% of the collected data, respectively.

Fig. 5

Extracted points segmented according to their distance to the surfaces. $S_1$: large front surface (dark blue), $S_2$: surface nearest background (light blue), $S_5$: surface farthest in back (yellow). Axes in meters.

Fig. 6

The figure shows 2-D plots of points along a line transect across the reference target. The line segments show the position of the surfaces of the reference target relative to the data. Each dot corresponds to one detected peak, with dots on the same lateral position belonging to the same histogram. The color of the dots corresponds to the estimated peak amplitude, with the color scale going from blue (weak) via green to yellow (strong).

The plot shows that none of the approaches accurately resolve the step to the left-most surface ($S_2$) situated 33 mm behind the front plate but instead they all produce spurious points in between. We can also see that LR-LM and SPISTA cannot resolve the step to the second closest surface $S_3$ (at 53-mm distance) and that EMGM-MF still has difficulties at 73 mm. We also note that the range noise is higher for EMGM-MF and FRI-MF than for the other methods.

6.

Discussion

Parameter values and thresholds in the peak detection approaches affect the number of correctly detected points, the number of outliers, and RMSE values. In this work, the values were chosen empirically rather than with an exhaustive analysis of the relation between parameter settings, system parameters, detection probability, false alarm rate, and accuracy. For some approaches, the number of outliers can easily be decreased by increasing the minimum signal strength threshold but that would impair the ability to find weak surfaces. Analogously, the thresholds could be adjusted to find more multiple surfaces at the expense of adding noise points. The FRI-MF approach may suffer in resolution performance because of too hard truncation of the Fourier series describing the IRF. Including more components has the cost of more false alarms and longer execution time. Finding the optimal parameters to maximize performance under different circumstances and with different requirements will be part of the next phase of our work. Nevertheless, we notice that for the IRF-LS, IRF-LS-MF, and SPISTA approaches, it was possible to find parameters that give a good trade-off between detection, false alarm rate, accuracy, and signal quality, which we consider a strength in itself. In cases with a high contrast between the highest peak and the threshold for the lowest peak, which is counted as a surface, the IRF-LS algorithm has been found to add extra peaks to compensate for imperfections in the description of the IRF. Because the threshold was set higher compared with the peak value than the error in the IRF description ($\sim 3\%$) this did not generate any false alarms here.

The low performance for the EMGM-MF approach is due to the fact that an MoG is not a good model for the wide peaks resulting after the MF preprocessing step. Shin⁵ stated that SPISTA would have better performance than an MoG approach, which is the case also in our study.

We are ultimately aiming for peak detection that runs in real time to allow for human-in-the-loop experiments in realistic scenarios. The current implementations of the algorithms are made in Matlab with varying degrees of optimization, and it is not relevant to report execution time at this stage. However, we notice that the approaches involving least square optimization with recursive addition of components (IRF-LS and IRF-LS-MF) are noticeably slower than the other approaches.

In this work, we have focused on the range accuracy of the detected points. When creating a point cloud, the 3-D coordinates of the detected points are defined by the viewing angle of the sensor and not the actual lateral position of the surfaces within the IFOV. A surface may be detected as soon as it is within the footprint of the laser beam and as a consequence, objects typically appear wider than they really are. This effect cannot be mitigated by analyzing each range histogram individually. Instead, it requires taking the data in the vicinity of the point into account, either through postprocessing of the extracted 3-D point cloud or by including neighboring histograms in the peak detection process. Future work will include efforts to address this problem. Similar effects would appear from cross talk between pixels in an array detector.

The experiments used for this investigation were performed with a maximum of two discrete surfaces within the IFOV corresponding to one pixel, to have control over the true distances between peaks. One interesting application is targets behind foliage, where the obscuration consists of several small objects at varying distances to the target surface. The results shown here should still be valid for the minimum distance between different surfaces, which is resolvable, regardless of whether there are two or more surfaces. The important factors to accurately measure the distance to the target surface are thus the number of photons detected reflected off the target surface and the distance from the last obscuring object to the target surface. The distance between different obscuring objects will impact whether they can be resolved or not, but this is often of little or no importance.

In this work, we assumed a static scene and a stationary sensor, which allows us to create and analyze range data as one-dimensional (1-D) histograms. Adding platform motion requires a different approach to handle the collected data. Future work will address efficient data structures for accumulation and analysis of TCSPC data as well as advanced positioning techniques for accurate 3-D registration of TCSPC data acquired with a moving sensor. This will probably need significantly higher data collection rate than the results shown here. Higher laser power that gives a detection probability of up to 50% for every laser pulse, which has been shown to be the optimum power to detect an obscured surface,¹² could be one solution, which may also require a faster scanning mechanism. Array detectors perform multiple measurements for every laser pulse may be another solution. Higher data collection rate will also allow for small field of regard 3-D imaging of moving objects.

7.

Conclusions

Of the approaches discussed herein, IRF-LS is our top candidate for peak detection in TCSPC data. It has the best performance regarding detection of multiple surfaces, which gives good range of accuracy and a low false alarm rate. However, it is computationally demanding in its current implementation and depending on what degree of optimization is achievable, other approaches may be considered in practical applications. Among the faster approaches, SPISTA gives the best results. Hence, a combination of several approaches is worthwhile investigating, where an advanced approach may be applied when multiple surfaces are present or when high accuracy is needed, and a less computationally demanding technique may be used to process the signal received from single surfaces. The width of a peak in a histogram might serve as an indicator of which approach to apply.

The processing of the data from the TCSPC lidar system with 392 ps FWHM IRF, indicating resolution of 59-mm, could with several methods resolve contributions from two surfaces at 53-mm range difference but not separate the contributions at 33-mm range difference.