# Compressive sensing for 3D-LiDAR imaging: A pipeline to increase resolution of simulated single-photon camera

Erwan Viala, Paul-Edouard Dupouy, Nicolas Riviere, Laurent Risser

► To cite this version:

Erwan Viala, Paul-Edouard Dupouy, Nicolas Riviere, Laurent Risser. Compressive sensing for 3D-LiDAR imaging: A pipeline to increase resolution of simulated single-photon camera. Journal of Quantitative Spectroscopy and Radiative Transfer, 2024, 314, pp.108866. 10.1016/j.jqsrt.2023.108866. hal-04390950

HAL Id: hal-04390950

https://hal.science/hal-04390950v1

Submitted on 12 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Review

# Compressive sensing for 3D-LiDAR imaging: A pipeline to increase resolution of simulated single-photon camera

Erwan Viala<sup>a,*</sup>, Paul-Edouard Dupouy<sup>a</sup>, Nicolas Riviere<sup>a</sup>, Laurent Risser<sup>b,c</sup>

<figure><img src="image_3.png" /></figure>

ARTICLE INFO

Keywords:

3D-LiDAR

Compressive Sensing

Laser imaging

Photon counting

Signal processing

A B S T R A C T

In this paper, we present a more efficient strategy than existing solutions to enhance the lateral resolution of low photon 3D-LiDAR operating in Geiger mode. Our pipeline makes it possible to reconstruct 3D-images with an unprecedented lateral-resolution, simultaneously at low photon count and Hertz level framerates. It is applied on simulated GmAPD 3D-LiDAR signals. Signals acquired using these kind of sensors are unsuitable for direct applications of Compressive Sensing algorithms. Our contribution focuses on a more efficient strategy for waveform denoising and reconstruction. For each pixel, we reconstruct sub-pixels with a Compressive Sensing approach. After describing our method, we demonstrate its applicability on realistic simulated data.

1. Introduction

Single-photon cameras are increasingly being used in 3D-LiDAR sensing systems, especially in low Size Weight and Power (SWaP) embedded systems. Here we focus on systems based on Geiger mode Avalanche Photodiode (GmAPD) sensors [1]. These sensors are highly sensitive [2–4], capable of statistically detecting fractional photons return levels. This, coupled with centimeter depth-resolution provides a competitive edge against “classical” linear systems. This is especially true in long-range remote sensing (>10 km) and surveillance applications [5,6]. However, GmAPD sensors are limited to a few pixels (up to 128 × 128 pixels) [7,8], which may result in insufficient resolution. Increasing the pixel count at the sensor level might pose its own problems, namely the size of the Focal Plane Array (FPA) and light gathering capabilities. In addition, GmAPD sensors can only detect the first photon [3] without information on photon count. Therefore, the first detection may introduced a loss of later detections, we name it dead-time effect, it is sometime also referred as pile-up [9,10]. Performance in many applications is strongly correlated with the lateral resolution of the imaging system. This paper describes an original method to increase the lateral resolution of GmAPD cameras.

Few methods have been detailed to improve GmAPD lateral-resolution, reduce SWaP and/or acquisition-time for long-range applications. For instance, Gordon et al. [11] performed long-range imaging using a 32 × 32 pixels GmAPD without increasing lateral resolution. Some methods to increase the lateral resolution albeit with different constraints were used successfully in 3D-LiDAR imaging. In

2016, Shin et al. [12] demonstrated photon efficient reconstruction with a visible Single Photon Avalanche Diode (SPAD) by scanning static scenes with a two-axis motorized scanner system. Long-range acquisitions were made at 10 km [13] in 2017 and 45 km in 2020 [14,15] using a scanning system. In 2019, Tachella et al. [16] demonstrated real-time application at medium range (300 m) with improved lateral resolution using a purely numerical method, which can be akin to a 3D interpolation.

Scanning systems are compatible with static ranging but not for dynamic objects observation in an operational context. A FPA delivers higher framerates than scanning systems [17] among 3D-LiDAR. Additionally, Compressive Sensing (CS) demonstrated equivalent results using only 10% of the acquisitions number required for scanning systems [18]. Purely numerical methods such as interpolation are in most cases unable to resolve sub-pixel objects [19]. Moreover, super-resolution tends to be blurrier than CS reconstruction [20]. Thus, CS based methods were shown as less vulnerable to these weaknesses. For these reasons, our study is focused on a CS method [21–23].

Since 2008, CS was used in the 3D-LiDAR community to increase single-pixel GmAPD lateral-resolution [18,24–31]. In most cases, these works use Orthogonal Matching Pursuit [32] or Total Variation [33] algorithms. In 2017, Edgar et al. [28] performed real-time acquisition with CS at a metric range with a single pixel sensor. In this paper, we define a complete pipeline to acquire 3D-images, leveraging CS and InP/InGaAs GmAPD sensors, applied to simulated data. We reconstruct

<footnote><sup>a</sup> ONERA - DOTA, University of Toulouse, Toulouse, FR-31055, France</footnote>

<footnote><sup>b</sup> Institut de Mathématiques de Toulouse (UMR 5219), CNRS, University of Toulouse, Toulouse, FR-31062, France</footnote>

<footnote><sup>c</sup> Artificial and Natural Intelligence Toulouse Institute (ANITI), Toulouse, 31000, France</footnote>

<footnote>* Corresponding author.
E-mail address: erwan.viala@onera.fr (E. Viala).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. 3D-LiDAR setup.</figcaption></figure>

the 3D-images using CS considering each pixel in the FPA as a single-pixel detector. As in previous works on single-pixel imaging [18,24–31], we use a Digital Micro-mirror Device (DMD) [34] to modulate the incoming light onto the FPA in a binary spatial-modulation scheme which we call patterns.

Hereafter, we define a complete pipeline to apply CS to GmAPD sensor acquisition and constraints. This pipeline is fully described in Section 2. Then we numerically demonstrate in Section 3 its feasibility using simulated data and, discuss its possibilities capabilities and limits.

## 2. Methods and materials

### 2.1. General overview

The objective of this work is to numerically reconstruct 3D-images with CS under the constraints imparted by the use of GmAPD sensors. This kind of sensors detects only the first photoelectron without information on the photon count [3]. To apply CS methods, we defined a specific simulated setup in Fig. 1. This setup is a combination of a classical quasi-monostatic active system with a GmAPD array associated with a DMD.

In most cases, signals acquired using a GmAPD array are reconstructed by using histograms of return-times to increase the dynamic range. While this approach bears good results [35], its use with CS is precluded by a fixed upper bound on the number of acquisitions (less than $10^3$ per pattern [2]) required to maintain a sufficient final framerate. In our case, the signal-to-noise ratio can be low and time-correlated due to the GmAPD nature. As shown in Section 3.1, under these constraints, applying iteratively the noise reduction and then the waveform estimation (resp. waveform estimation then denoising) leads to poor waveform estimation (resp. dramatic increase in noise over time). So, to reconstruct the waveform, denoising and waveform estimation are made at the same time.

The state-of-the-art of ruggedized fiber laser-sources compatible with low SWaP systems using sufficient energy and nanosecond pulse length is around 20 kHz [36] at 1.55 µm. Commercially available GmAPD cameras can achieve framerates of approximately 150 kHz [37].

This difference will be exploited (Section 2.2.2) to acquire frames without laser illumination during the source's dead time. These additional acquisitions can be leveraged to partially compensate the GmAPD constraints. The idea is to perform a statistical test on each time-bin for each pixel in order to separate the signal from the noise and directly estimate its support.

As illustrated in Fig. 2, the whole process can be split in several steps by (i) doing a compressed acquisition, (ii) building the histogram of return-times for each pattern and, (iii) simultaneously retrieving the signal support and correcting waveform before (iv) sequentially applying a classical CS solver in each bin to reconstruct the final 3D-image.

We need to address the temporal correlation induced by a dead-time effect. This can occur if the backscattered light generates a signal with multiple echoes. The pattern-coefficient is defined as the acquired signal for a bin and a given pattern. Let us consider two acquisitions with two different patterns but identical cross-section. In that case, their pattern-coefficients must be equal independently of detections in previous bins. Not meeting this condition may lead to an erroneous estimate of pattern-coefficients and thus makes 3D-image reconstruction by CS incorrect.

Fig. 3 presents the temporal correlation effect. Two patterns are selected. The $\Phi_1$ pattern hides only one object. Both objects are visible with the $\Phi_2$ pattern. In the first case, the signal peaks at twice the level of the second case. Correctly, recovered signal should exhibit identical peak level for the “triangle” object as they have an identical cross-section. This illustrates the effect of the pattern selection on the temporal signal while dead-time occurs. The same effect will occur for a slope or any target combination producing multiple echoes in the retrieved signal.

### 2.2. Preprocessing

#### 2.2.1. Full waveform reconstruction

To address the temporal correlation, there are two ways commonly used. First, decreasing the signal level to reach the domain where GmAPD have a linear behavior, typically at 0.01 average photon
<figure><img src="image_5.png" /><figcaption>Fig. 2. Illustration of the pipeline and its associated notations by considering a reconstruction example.</figcaption></figure>

counts [38]. Second, define a process to correct the dead-time effect inducing the temporal correlation. In recent advances for dead-time effect correction, one imply inverse method base on bayesian framework to revert the dead-time effect and demonstrate sub-picosecond- accurate 3D-image using it [39]. Another recent approach [40], defines a probabilistic model correction using Markov-chain under high-flux conditions for a non-paralyzable and asynchronous sensor enabling higher dynamic range GmAPD measurement.

We define a simpler and faster signal estimator that works when the ratio of pulse width to bin length is small (typically less than or equal to one), as demonstrated in Section 3.1. When using a larger pulse, consider using a slower but much more accurate algorithm in this case, such as [39,40]. Our method is inspired by Coates pile-up correction process [9,10,41] extended for an asynchronous GmAPD sensor in a CS setup. GmAPD sensors are designed to detect events at the single-photon level. Their dynamic ranges are effectively of one since they can only detect an event with no information on the photon count. O'Brien and Fouche [3] defined and validated the following model: Eq. (1).

$$P(p; t_1; t_2) = \frac{1}{p!} [Y(t_1, t_2)]^p e^{-Y(t_1, t_2)} \quad (1)$$

where $P(p; t_1; t_2)$ is the probability for a given pattern that an amount of p events occurs between times $t_1$ and $t_2$ and $Y(t_1, t_2)$ is the rate function of the process [3] which varies with noise, scene (including optical and geometrical properties) and setup. Based on the state-of-the-art of the cameras, we only consider equally sized bins. As GmAPD sensors can only measure whether at least one photon was detected, we define two

probabilities for each bin in Eq. (2).

$$ \begin{cases} P(p_k = 0) = P(0; t_k; t_{k+1}) = e^{-Y(t_k, t_{k+1})} \\ P(p_k > 0) = 1 - P(0; t_k; t_{k+1}) = 1 - e^{-Y(t_k, t_{k+1})} \end{cases} \quad (2) $$

where $p_k$ is the number of events that occurs between $t_k$ and $t_{k+1}$. This probability assumes that all bins are independent and only the first return is detected. Thus, we define the probability $D_k$ that a detection is made at bin k in Eqs. (3) and (4) taking the time dependency into account.

$$D_k = P(p_k > 0) \prod_{j=0}^{k-1} P(p_j = 0) \quad (3)$$

$$\Rightarrow D_k = (1 - e^{-Y_k}) \prod_{j=0}^{k-1} e^{-Y_j} \quad (4)$$

where $Y_k = Y(t_k, t_{k+1})$ is the rate function of the process between times $t_k$ and $t_{k+1}$. Thus, to retrieve the correct waveform, we must approximate the rate function $Y_k$ for each bin k.

The signal is retrieved from the return-times histogram. Let $H_k$ be the normalized histogram at bin k:

$$H_k = \frac{1}{N} \sum_{j=0}^{N-1} \mathbb{1}_{z_j=k} \quad (5)$$

where N is the number of acquisitions for each used pattern, and $z_j$ is the bin where the signal is detected during acquisition j. From the law of large numbers, $H_k$ converges to the probability of a detection at bin k, which is $D_k$ in Eq. (6).

$$H_k \underset{N \to \infty}{\longrightarrow} P[\mathbb{1}_{z=k}] = P(z=k) = D_k \quad (6)$$

Thus, $H_k$ approximates $D_k$. We use Eq. (6) to infer an approximation of all $(\hat{Y}_k)_{k \in \{0, k-1\}}$, which is the purpose of this step. If $D_k \in [0, 1]$, we deduce from Eq. (4) that:

$$e^{-Y_k} = 1 - D_k \times \left( \prod_{j=0}^{k-1} Y_j \right)^{-1} \in [0, 1] \quad (7)$$

We have $H_k \in [0, 1]$ by definition. Using the logarithm of Eq. (7) and by approximating $D_k$ by $H_k$, we define $\hat{Y}_k$ an estimator of $Y_k$ by recurrence in Eq. (8).

$$\hat{Y}_k = -\ln \left( 1 - e^{\sum_{j=0}^{k-1} \hat{Y}_j} \times H_k \right) \quad (8)$$

Eq. (8) allows us to estimate the signal using histograms. The noise needs to be removed to accurately reconstruct the signal. As mentioned before, this step must be concomitant to the reconstruction.

### 2.2.2. Function support estimation

We introduce a new method to denoise a signal using the constraint of a state-of-the-art setup (Section 2.1) at our advantage. Noise rate is nearly time invariant for each pixel, but variable between GmAPD array's pixels [42]. As mentioned before, it is possible to acquire a few frames (at least one) of noise between each useful frame with concomitant laser illumination. By statistically comparing bin-sets (all pattern coefficients for a given bin k) with $H_k^1$ and without $H_k^0$ concomitant laser illumination, we determine the support of the signal directly from the waveform, with set probability $\alpha$. For all patterns, we calculate the support of the signal and the corrected waveform. Then, we multiply both as illustrated in Fig. 2, step (iii).

We statistically test each bin-set of potential signal $H_k^1$ against pure-noise bin-set $H_k^0$ to define which population is the larger. The function's support is built from the bin-sets where $H_k^1$ is stochastically larger.

We define the null hypothesis as: “the distributions of both population are identical”. Usually a Student's t-test should be used to test this statistical hypothesis. In our case, the distributions of $H_k^1$ and $H_k^0$ are not normal. So, the use of non-parametric tests, without additional assumption on the noise or data distribution, is preferable. We choose
<figure><img src="image_6.png" /><figcaption>Fig. 3. Illustration of case where the temporal correlation is especially strong leading to an inaccurate 3D-image.</figcaption></figure>

the Mann-Whitney U rank test [43] among them, it is defined in Eq. (10). Here, we test the whole waveform for each pixel, considering each pattern-coefficient as an element of the population. In practice, we concatenate $H_k^1$ and $H_k^0$ to form an overall population. Then the resulting list is sorted by photon count. The U statistic is defined as the number times $H_k^0$ precedes $H_k^1$ as in Eq. (9). Finally U is used to compare the rank of the elements of both populations.

$$U = \sum_{m=1}^{M} \sum_{l_1=1}^{N_{pulse}} \sum_{l_2=1}^{N_{noise}} \mathbb{1}_{(H_k^1)_{m,l_1} > (H_k^0)_{m,l_2}} \quad (9)$$

where M is the number of patterns, $\langle \cdot \rangle_{m,l}$ the lth coefficient acquisitions of pattern m, $N_{pulse}$ and $N_{noise}$ are respectively the number of acquisitions with and without a concomitant laser illumination.

$$P(U \le \bar{U}) = \alpha, \text{ under the null hypothesis} \Rightarrow \text{the hypothesis of identical distribution is rejected} \quad (10)$$

Thus, we consider the test as significant with level $\alpha$ if the probability that U is inferior or equal to the $\bar{U}$ under the null hypothesis, i.e. if the distributions were identical. Under the null hypothesis, both distributions are equal so the ranks are uniformly random. $\alpha$ determines the tolerance level at which we can affirm that a bin is in the support by rejecting the null hypothesis. The practical implementation used is defined in [44].

## 2.3. Compressive sensing

After reconstruction through the first two steps (Sections 2.2.1 and 2.2.2), each bin is time-independent from the others and, noise is only present where signal is detected. Practically, we simulate a set of small DMD pixels imaged by a GmAPD sensor. Reconstruction is performed by considering each GmAPD focal-plane-array pixel with a set of DMD pixels as a single-pixel camera configuration. For each camera pixel, we solve a CS optimization problem for each bin to reconstruct local 2D-images that we concatenate to obtain a global 3D-image.

Let $X_k$ be the vectorized reconstructed image at the pixel $(i, j)$ and the bin $k$) and $Y_k$ the corresponding measurement. Let $\Phi$ be the measurement basis. Due to the binary nature of the DMD, a linear transformation to convert this basis into binary $\{0, 1\}$ must exist. Let $\Psi$ be the basis where $X_k$ is K-sparse (i.e. only K components of $X_k$ are non-zero). Let N the number of DMD pixels in one group corresponding to the dimension of $X_k$. We can express the CS optimization problem as:

$$(P_{CS}) \left\{ \begin{array}{ll} \min_{s \in \mathbb{R}^N} & \|s\|_0 \\ \text{s.t.} & \|Y - \Phi \Psi s\|_2^2 \le \epsilon, \quad \epsilon \ge 0 \end{array} \right. \quad (11)$$

where $s_k = \Psi X_k$ and $\epsilon$ is an error term. To recover accurately $s_k$ from $Y_k$ and then $X_k$, it is established [21,23] that $\Phi$ must follow

the 2K-RIP (Restricted Isometry Property). This means that $s_k$ can be inferred from $Y_k$, i.e. $s_k$ must be K-sparse and $\Psi$ must be highly incoherent with $\Phi$. In practice, solving exactly the problem with $l_0$ norm is known to be NP-hard [45]. Therefore, the use of a different approach is required. The optimization problem can be relaxed by using the $l_1$ norm. Another approach is to approximate the solution with a greedy algorithm [32]. Here, we choose the later option with Orthogonal Matching Pursuit (OMP) since it provides the best precision / performance ratio [46]. In this paper, the measurement process $\Phi$ is defined as an undersampled Hadamard transform and the sparsifying basis $\Psi$ as the 2D-Haar transform. Thus, as a measurement process, each row of $\Phi$ represents one vectorized pattern.

### 2.3.1. Orthogonal matching pursuit using cholesky decomposition

The following algorithm is an efficient way to implement OMP as defined in [32].

We can solve the $l_2$ minimization problem evaluated at each step of the 'while' loop of Algorithm 1 by the orthogonal projection $(\Phi_{|\Lambda}^T \Phi_{|\Lambda})^{-1} \Phi_{|\Lambda}^T$. The projection's computation is accelerated by approximating the inverse matrix using Cholesky or QR decomposition. This strategy was proposed in [47,48]. This implementation is efficient when the size of the matrix is low enough ($<10^2$). Otherwise, an approximation of the projection can be made with a similar method, like the Batch-OMP algorithm [48]. We achieved a reduction of the execution time by two orders of magnitude with this implementation.

#### Algorithm 1 Orthogonal Matching Pursuit

1: **procedure** OMP($N \times d$ measurement matrix $\Phi$, $N - dim$ measurement vector $Y, N \times N$ Sparsifying matrix $\Psi$)

2: $r \leftarrow Y$

3: $\Lambda \leftarrow \emptyset$

4: **while** halting criterion not True **do**

5: $\Lambda \leftarrow \Lambda \cup \left\{ \underset{j=1,...,d}{\operatorname{argmax}} \left| \langle r, (\Phi \Psi)_j \rangle \right| \right\}$

6: $S \leftarrow \operatorname{argmin} \|Y - (\Phi \Psi)_{|\Lambda} s\|_2$

7: $r \leftarrow Y - (\Phi \Psi)_{|\Lambda} S_{|\Lambda}$

8: **end while**

9: **return** $\Psi S$

10: **end procedure**

### 2.3.2. Graphics Processing Unit (GPU) acceleration

The problem is well suited for parallel implementation since we execute our algorithm on each pixel independently. Commercially available sensors are 32 × 32 or 32 × 128 pixels array which is compatible with the number of blocks of current GPU generation. With the Cholesky version of OMP, we also reduce the necessary memory.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Selected scenes for simulation tests. (a) A synthetic scene at 13 km from the sensor with a field of view of 0.8 mrad in both horizontal and vertical directions equivalent to a resolution of 4 cm at 13 km. (b) The French motto, a bust of the French Marianne and a 1 m wide "ONERA" letters at 10 km from the 3D-LiDAR system with a square field of view of 0.4 mrad equivalent to a resolution of 2 cm.</figcaption></figure>

This enable the computation of the complete pipeline for all pixels simultaneously. One can notice that this version of OMP is almost only vector multiplication that could more optimized to work in parallel using the different threads inside each block.

When pipeline ran only on GPU, we gain another two orders of magnitude in terms of execution time. A short processing time allows low latency reconstruction that is an important advantage, eg. for surveillance applications (close to real-time).

## 2.4. Simulation of 3D-LiDAR signal

Simulated data are used to design, test and validate our method. The MATLAB software — an ONERA software [49] allows us to simulate the signal from a GmAPD sensors paired with a laser source. This end-to-end tool models 3D laser imagers, using physically realistic sensor models, validated through laboratory and field experiments. It relies on a description of the optronic scene, including the geometry and the optical properties of the materials for each facet of the target [50,51]. Moreover, simulation provides full control of input parameters and information unavailable on field data to design, test and validate our approach:

* The same mathematical model [3] can be used to both simulate detection and reconstruction of the correct waveform;

* Realistic signal features are artificially accessible eg. photon count or signal/noise classes.

The pipeline is sensitive to the geometry of the scene. To test its robustness against different geometrical features, various 3D objects were created. The following Fig. 4 illustrates the 3D targets used to generate the simulated data. The (a) scene provides multiple-types of planar surfaces including slopes. From these, we evaluate spatial frequency, lateral and temporal resolutions and, responses of our pipeline.

The (b) scene provides finer details, with broader spatial frequency content distribution. The French motto and the bust of Marianne are used to illustrate behavior of our process in different context.

The pipeline is also sensitive to signal level to a lower extend. We define the surface of the object covered by a 10% Lambertian hemispherical reflectance material to get a uniform response from all sensor's pixels. In this work, we simulate the camera high enough to neglect the atmospheric turbulence effect. We introduce specifications (Table 1) compatible with a single-photon 3D-LiDAR using commercially available components. We consider two sources of noise in the simulation, solar detection and Dark Count Rate. We only give the combine contribution of noise for each pixel as Noise Count Rate.

The simulations are carried out at 13 km range. In a real scene, the start of the gate could be a real burden, it is not alleviated in the work.

<table><caption>**Table 1**
Specification of the simulated 3D-LiDAR system using commercially available components.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Laser</strong></td><td></td></tr><tr><td>Pulse width</td><td>0.25 ns</td></tr><tr><td>Pulse energy</td><td>100–150 μJ</td></tr><tr><td>Pulse repetition frequency</td><td>20 kHz</td></tr><tr><td>Per pattern pulse acquisition</td><td>1000</td></tr><tr><td>Wavelength</td><td>1550 nm</td></tr><tr><td>Divergence</td><td>Matching the field of view</td></tr><tr><td>Aperture</td><td>5 mm</td></tr><tr><td><strong>Sensor</strong></td><td></td></tr><tr><td>GmAPD acquisition frequency</td><td>186 kHz</td></tr><tr><td>Number of GmAPD array pixels</td><td>32 × 32</td></tr><tr><td>DMD pattern frequency</td><td>5 kHz</td></tr><tr><td>DMD maximum/used pixels</td><td>1024 × 768/256 × 256</td></tr><tr><td>Bin length</td><td>0.25 ns</td></tr><tr><td>Noise Count Rate</td><td>1 MHz</td></tr><tr><td>Quantum efficiency</td><td>0.4</td></tr></tbody></table>

A lot of applications will need to consider complementary approach, to this method to find a correct gate starts (e.g., another sensor for detection purpose only).

All simulations made in this paper, aims a 64-folds resolution enhancement. The physical factors which could limits the lateral resolution, will be linked to the instrument size due to diffraction. This is not restrictive in a simulation, but limiting the resolution enhancement to 64-folds leaves enough margin to consider a realistic instrument size at low SWaP.

## 3. Results and discussion

The complete pipeline is applied to multiple simulated scenes by the MATLAB software to validate our approach. To cover most types of real objects that could be imaged, we select four representative categories: (i) objects spread over the 2-lateral dimensions (flat surface with large cross-section) (ii) object(s) producing multiple echoes in a single pixel, (iii) objects spread over the temporal dimension (slanted surface), (iv) sub-pixel objects.

### 3.1. Waveform reconstruction

Results obtained after the reweighting step (see Section 2.2.1) are plotted in Figs. 5 and 6. Outputs of the support estimation and their comparisons with ground truth (Section 2.2.2) are plotted in Fig. 7 and summarized in Table 2. In Figs. 5 to 7, $Y_k$ and its support are obtained from artificially available data thanks to the simulation. These
<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of the mean histogram and the mean reconstructed waveform (Section 2.2.1) of 100 simulated signals. Each sub-picture corresponds to a different pixel of the GmAPD camera. Blue is the signals' ground-truth, red the histograms' mean and its 95% confidence interval and, green the estimated signals' mean and 95% confidence interval. The gate starts at 13 km from the camera.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Noise-effect illustration on a large-gate signal. The blue, red and green curves represent the true waveform $Y_k$, the estimated signal $\hat{D}_k$ as defined in Section 2.2.1 and the histogram $H_k$, respectively. The gate starts at 10 km from the camera.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. From top to bottom: the 'ground truth' $Y_k$, the raw histogram $H_k$ and $\hat{Y}_k$ the estimated support using the test defined in Section 2.2.2 applied to the estimated waveform $\hat{D}_k$. Filled bins are in the function's support. Each column is the signal processed for different GmAPD pixel. From left to right: case (B), (C) and (D). Min-Max normalization, by column, is applied to the three signals. The gate starts at 13 km from the camera.</figcaption></figure>

data are not accessible to the algorithm, as these variables are just considered for direct simulations in MATLIS software. Nevertheless,

we store this values for future exploitation and comparisons between ground truth and reconstructed point clouds. $\hat{D}_k$ is the reconstructed,
<figure><img src="image_11.png" /><figcaption>Fig. 8. From top to bottom: ‘ground truth’ $Y_k$, estimated waveform $\hat{D}_k$ with thresholding and then $\hat{Y}_k$ estimated support applied to the estimated waveform $\hat{D}_k$. Filled bins are in the function’s support. Signal is normalized with Min–Max normalization. The gate starts at 10 km from the camera.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Results obtained after: (a) and (c) 3D-LIDAR acquisition plotted with only the point labeled as signal; (b) and (d) a CS acquisition with all patterns, processed with the complete pipeline, without any supplementary post-processing. The grayscale represents the normalized intensity from the pipeline.</figcaption></figure>

time decorrelated waveform, as defined in Section 2.2.1. $H_k$ is the return-times histogram.

We use the Peak Signal-to-Noise Ratio (PSNR, Eq. (12)) to evaluate the quality of the waveform-estimation (Fig. 6).

$$PSNR_{Y_k}(\hat{Y}_k) = 20\log_{10}\left(\frac{\max_{Y}}{\sqrt{\frac{1}{N}\sum_{i=0}^{N}\|Y_{k,i} - \hat{Y}_{k,i}\|_{2}^{2}}}\right) \quad (12)$$

We process one hundred signals to better represent the trends in the estimation. Then we calculate the mean of the signal and an estimated 95% confidence interval for each bin.

We select four pixels of interest to illustrate the results on different specific situations:

(A) Plane surface perpendicular to the camera optical axis, with a size superior to a camera pixel

(B) Plane surface perpendicular to the camera optical axis, with a size inferior to a projected DMD pixel
<caption>**Table 2**</caption>

Comparison of confusion matrices, (a) is the confusion matrix of support($H_k$) versus support($Y_k$); (b) is the confusion matrix of support($H_k$) using a threshold versus support($Y_k$) and (c) is the confusion matrix of support($\hat{Y}_k$) versus support($Y_k$). Green hue corresponds to the true positive/negative and red to the false positive/negative, and saturation lowers with the quality of the estimation.

<table><thead><tr><th>Predicted</th><th>out</th><th>in</th></tr></thead><tbody><tr><th>True</th><td></td><td></td></tr><tr><th>out</th><td>216464</td><td>53181</td></tr><tr><th>in</th><td>393</td><td>1322</td></tr></tbody></table>

<table><thead><tr><th>Predicted</th><th>out</th><th>in</th></tr></thead><tbody><tr><th>True</th><td></td><td></td></tr><tr><th>out</th><td>263652</td><td>5992</td></tr><tr><th>in</th><td>565</td><td>1151</td></tr></tbody></table>

<table><thead><tr><th>Predicted</th><th>out</th><th>in</th></tr></thead><tbody><tr><th>True</th><td></td><td></td></tr><tr><th>out</th><td>269274</td><td>371</td></tr><tr><th>in</th><td>164</td><td>1551</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** 16 pattern From top to bottom, result obtained after: removing the waveform reconstruction step, replacing the support estimation step by a threshold chosen to avoid more losses for signal recovery and, using the whole pipeline. From left to right, result using: the synthetic scene with 16, 32, 48 pattern measurement. The grayscale represents the intensity normalized from the pipeline for a given scene.</figcaption></figure>

(C) Two planes perpendicular to the camera optical axis, at two different ranges, each filling approximately half of a sensor pixel field-of-view

(D) A 75° slope, with respect to the camera optical axis, filling all the pixels

While the histogram's mean PSNR is 64.6 with a variance of 180.4, that of the estimation's is 71.3 with a variance of 2.6 both calculated on the synthetic scene (Fig. 4(a)). The estimated waveform error is less sensitive to the spatial distribution of the objects and signal level (see Fig. 5). High ($6 \times 10^{-1}$) and low ($6 \times 10^{-3}$) photon-counts estimated waveform errors are more stable to low photon count ($6 \times 10^{-3}$)
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** From top to bottom, result obtained after: (a) removing the waveform reconstruction step, (b) replacing the support estimation step by a threshold chosen to avoid more losses for signal recovery and, (c) using the whole pipeline. From left to right, result using the ONERA scene with 16, 32, 48 pattern measurement. The grayscale represents the intensity normalized from the pipeline for a given scene.</figcaption></figure>

estimation error than for the histogram. Coates method [9,10] gives nearly the same result as the method define in Section 2.2.1. These methods are sensitive to saturation: if the sensor saturates (average photon count is close to one), the estimation tends to be unstable due to the exponential term in the reconstruction formula. If a measurement of the average photon number close to one is required, more complex and computational intensive method as [39,40] should be used to replace this step.

The estimated waveform is less time-correlated, which is the main goal of this step. It is clear in most cases except for the second one, where signal levels are too low to observe dead-time effect. Indeed, the error increases with time while the signal decreases. This effect is illustrated in case (D). If signal level is low enough, the following bins will not be greatly affected, as in (B). Temporal response is correlated to pulse length. Long pulses with respect to the bin-length will cause temporal distortion. In this case, adding a deconvolution step using available knowledge on the pulse should be considered. Here, we only consider short pulses, with full width at half maximum shorter than half the temporal bin size, where this effect can be neglected.

A new case with multiple targets along the line of sight and increased gate length is considered to illustrate the effect of noise on

the support estimation. Fig. 6 highlights a drawback of our method: the noise is increasing along the gate. Here, discriminating signal from noise is especially challenging. For example, the signal around 23 m is at the same level as noise at 40 m. Using a simple threshold is insufficient to insure both low false alarm and good detection rates.

## 3.2. Support estimation

Estimated function support is compared against ground truth and histogram based method in Fig. 7 for cases (B), (C) and (D).

Each point of the estimated waveform is considered as a statistical observation, labeled in or out with respect to ground truth. Thus, we construct the confusion matrix (Table 2) to evaluate the estimation's quality. In most cases, support is correctly retrieved with a tendency to slightly broaden it around peaks.

Estimating the support using this method increases true positives/negatives and decreases false positives/negatives numbers. No absolute threshold is necessary. Only a significance level for the test must be set. Arbitrarily, we set it at 0, 1%. Our method provides better support estimation than the classical histogram method, with a higher true positive and lower false negative. Statistics in Table 2 are made on
all pixels, bins and patterns. Our method provides higher statistical
confidence than evaluating support for each pattern separately.

We apply the current step to the case presented in Fig. 6. In Fig. 8, we illustrate that our approach works well even with lower signal to noise ratio. For example, support at 23 m and 40 m is correctly evaluated.

3.3. Complete process

As shown in Fig. 9, the pipeline allows the recovery of many details in the scene in comparison with the raw data thanks to the increased number of points.

One can notice reconstruction artefacts (e.g. holes, noise), stemming
from the resolution of the CS optimization problem. Highest spatial
frequency (spoke target) and plane surface with a cross-section smaller
than one projected DMD pixel-size are not recovered (synthetic scene).

Our method has clear advantages for recovering details by recon-
structing better-defined edges like in cases (A), (C) and (D) (Figs. 9(a)
and 9(b)) or the ‘ONERA’ letters (Figs. 9(c) and 9(d)). Another advan-
tage is about high spatial frequency objects (French Motto) and with
depth variations (Marianne’s bust). These objects are identifiable in
Figs. 9(b) and 9(d) where they were only recognizable in Figs. 9(a) and
9(c). The lateral resolution is increased by a factor of 8 in each lateral
direction (Fig. 9), depth resolution remains unchanged.

In Figs. 10 and 11, we investigate the variation of quality in the
reconstruction process while removing different parts of the pipeline
on the different scenes and different numbers of pattern-measurements.
Without the waveform reconstruction the overall intensity decreases
and additional artefacts show up. These artefacts depend on sparsifying
and measurement basis. In the second row, an optimal threshold was
define to avoid losses for CS optimization recovery and lets a lot
of noisy acquisition. As classical CS algorithm, it allows a 3D-image
acquisition with increased lateral resolution using only few compressed
measurement. One can notice that the quality decreases faster relatively
to the compression rate than others classical CS algorithms [18]. The
decrease in quality with the compression ratio does not decrease lin-
early with the image size [52]. However, to highlight the interest of our
method, it only needs 16 full-waveform pattern-measurements with a
32 × 32 pixels GmAPD array to reconstruct an 256 × 256 × 256 vox-
els 3D-image. Moreover, using the complete pipeline limits unwanted
artefacts emerging from CS recovery while efficiently denoising.

4. Conclusion

This paper describes a complete pipeline to increase lateral reso-
lution of GmAPD 3D-LiDAR systems using Compressive Sensing (CS).
Our pipeline considers the nature of the signal using a statistical
approach to denoise and estimate a corrected waveform before solv-
ing the CS optimization problem. The feasibility of our method is
numerically demonstrated on simulated data. In our case, the lateral
resolution is enhanced 64-folds while preserving the depth resolution
and minimizing artefacts. This pipeline is suitable under daylight at
very long-range. We validated on static targets. Future work will deal
with the reconstruction of moving targets by evaluating the use of
pattern registration.

Funding

This work is supported by the Occitanie Region through a doctoral
fellowship.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Data availability

Data will be made available on request.

References

[1] Gatt P, Johnson S, Nichols T. Geiger-mode avalanche photodiode ladar receiver performance characteristics and detection statistics. Appl Opt 2009;48(17):3261-76. http://dx.doi.org/10.1364/AO.48.003261, URL http://ao.osa.org/abstract.cfm?URI=ao-48-17-3261.

[2] Hiskett PA, Gordon KJ, Copley JW, Lamb RA. Long range 3D imaging with a 32 × 32 Geiger mode InGaAs/InP camera. In: Advanced photon counting techniques VIII. 9114, SPIE; 2014, p. 67–79. http://dx.doi.org/10.1117/12.2050540.

[3] O'Brien ME, Fouche DG. Simulation of 3D laser radar systems. Linc Lab J 2005;15(1):37–60, doi: None, URL https://www.ll.mit.edu/sites/default/files/page/doc/2019-01/15_1simulation.pdf.

[4] Entwistle M, Itzler MA, Chen J, Owens M, Patel K, Jiang X, et al. Geiger-mode APD camera system for single-photon 3D LADAR imaging. In: Advanced photon counting techniques VI, vol. 8375. SPIE; 2012, p. 78–89. http://dx.doi.org/10.1117/12.921004.

[5] Hespel L, Riviere N, Fraces M, Dupouy PE, Coyac A, Barillot P, et al. 2D and 3D flash laser imaging for long-range surveillance in maritime border security: detection and identification for counter UAS applications. In: Laser radar technology and applications XXII, vol. 10191. SPIE; 2017, p. 52–62. http://dx.doi.org/10.1117/12.2269141.

[6] Felix T. Maritime protection of critical infrastructure assets in the Campeche sound [Ph.D. thesis], Naval Postgraduate School; 2005.

[7] Aull B. Geiger-mode avalanche photodiode arrays integrated to all-digital CMOS circuits. Sensors 2016;16(4):495. http://dx.doi.org/10.3390/s16040495, URL https://www.mdpi.com/1424-8220/16/4/495.

[8] Bradley CP, Mukherjee SS, Reinhardt AD, McManamon PF, Lee AO, Dhulla V. 3D imaging with 128 × 128 eye safe InGaAs P-I-N LiDAR camera. In: Laser radar technology and applications XXIV, vol. 11005. SPIE; 2019, p. 273–9. http://dx.doi.org/10.1117/12.2521981.

[9] Coates PB. The correction for photon 'pile-up' in the measurement of radiative lifetimes. J Phys E 1968;1(8):878. http://dx.doi.org/10.1088/0022-3735/1/8/437.

[10] Coates P. Pile-up corrections in the measurement of lifetimes. J Phys E 1972;5(2):148. http://dx.doi.org/10.1088/0022-3735/5/2/018.

[11] Gordon KJ, Hiskett PA, Lamb RA. Advanced 3D imaging LiDAR concepts for long range sensing. In: Advanced photon counting techniques VIII, vol. 9114. SPIE; 2014, p. 49–55. http://dx.doi.org/10.1117/12.2051013.

[12] Shin D, Xu F, Venkatraman D, Lussana R, Villa F, Zappa F, et al. Photon-efficient imaging with a single-photon camera. Nature Commun 2016;7(1):12046. http://dx.doi.org/10.1038/ncomms12046.

[13] Pawlikowska AM, Halimi A, Lamb RA, Buller GS. Single-photon three-dimensional imaging at up to 10 kilometers range. Opt Express 2017;25(10):11919–31. http://dx.doi.org/10.1364/OE.25.011919, URL http://www.opticsexpress.org/abstract.cfm?URI=oe-25-10-11919.

[14] Liu X, Shi J, Sun L, Li Y, Fan J, Zeng G. Photon-limited single-pixel imaging. Opt Express 2020;28(6):8132–44. http://dx.doi.org/10.1364/OE.381785, URL http://www.opticsexpress.org/abstract.cfm?URI=oe-28-6-8132.

[15] Li Z-P, Huang X, Cao Y, Wang B, Li Y-H, Jin W, et al. Single-photon computational 3D imaging at 45km. Photon Res 2020;8(9):1532–40. http://dx.doi.org/10.1364/PRJ.390091, URL http://www.osapublishing.org/prj/abstract.cfm?URI=prj-8-9-1532.

[16] Tachella J, Altmann Y, Mellado N, McCarthy A, Tobin R, Buller GS, et al. Real-time 3D reconstruction from single-photon LiDAR data using plug-and-play point cloud denoisers. Nature Commun 2019;10(1):4989. http://dx.doi.org/10.1038/s41467-019-12943-7.

[17] Villa F, Severini F, Madonini F, Zappa F. SPADs and SiPMs arrays for long-range high-speed light detection and ranging (LiDAR). Sensors 2021;21(11):3839. http://dx.doi.org/10.3390/s21113839, URL https://www.mdpi.com/1424-8220/21/11/3839.

[18] Duarte MF, Davenport MA, Takhar D, Laska JN, Sun T, Kelly KF, et al. Single-pixel imaging via compressive sampling. IEEE Signal Process Mag 2008;25(2):83–91. http://dx.doi.org/10.1109/MSP.2007.914730, URL https://ieeexplore.ieee.org/document/4472247.

[19] Willett RM, Marcia RF, Nichols JM. Compressed sensing for practical optical imaging systems: A tutorial. Opt Eng 2011;50(7):072601. http://dx.doi.org/10.1117/1.3596602.

[20] Sun Y, Gu G, Sui X, Liu Y, Yang C. Single image super-resolution using compressive sensing with a redundant dictionary. IEEE Photonics J 2015;7(2):1-11. http://dx.doi.org/10.1109/JPHOT.2015.2409063, URL https://ieeexplore.ieee.org/document/7055220.

[21] Candès EJ, Romberg JK, Tao T. Stable signal recovery from incomplete and inaccurate measurements. Comm Pure Appl Math 2006;59(8):1207–23. http://dx.doi.org/10.1002/cpa.20124, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.20124.
[22] Candes EJ, Tao T. Near-optimal signal recovery from random projections: Universal encoding strategies? IEEE Trans Inform Theory 2006;52(12):5406-25. http://dx.doi.org/10.1109/TIT.2006.885507, URL https://ieeexplore.ieee.org/document/4016283.

[23] Candes E, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inform Theory 2006;52(2):489-509. http://dx.doi.org/10.1109/TIT.2005.862083, URL https://ieeexplore.ieee.org/document/1580791.

[24] Romberg J. Imaging via compressive sampling. IEEE Signal Process Mag 2008;25(2):14-20. http://dx.doi.org/10.1109/MSP.2007.914729, URL https://ieeexplore.ieee.org/document/4472239.

[25] Gibson GM, Johnson SD, Padgett MJ. Single-pixel imaging 12 years on: A review. Opt Express 2020;28(19):28190-208. http://dx.doi.org/10.1364/OE.403195, URL http://www.opticsexpress.org/abstract.cfm?URI=oe-28-19-28190.

[26] Howland GA, Dixon PB, Howell JC. Photon-counting compressive sensing laser radar for 3D imaging. Appl Opt 2011;50(31):5917-20. http://dx.doi.org/10.1364/AO.50.005917, URL http://ao.osa.org/abstract.cfm?URI=ao-50-31-5917.

[27] Howland GA, Lum DJ, Ware MR, Howell JC. Photon counting compressive depth mapping. Opt Express 2013;21(20):23822-37. http://dx.doi.org/10.1364/OE.21.023822, URL http://www.opticsexpress.org/abstract.cfm?URI=oe-21-20-23822.

[28] Edgar MP, Johnson S, Phillips DB, Padgett MJ. Real-time computational photon-counting LiDAR. Opt Eng 2017;57(3):1-7. http://dx.doi.org/10.1117/1.OE.57.3.031304.

[29] Sun M-J, Tong M, Edgar M, Padgett M, Radwell N. A Russian dolls ordering of the Hadamard basis for compressive single-pixel imaging. Sci Rep 2017;7(1):3464. http://dx.doi.org/10.1038/s41598-017-03725-6, URL https://www.nature.com/articles/s41598-017-03725-6.

[30] Higham C, Murray-Smith R, Padgett M, Edgar M. Deep learning for real-time single-pixel video. Sci Rep 2018;8(1):2369. http://dx.doi.org/10.1038/s41598-018-20521-y, URL https://www.nature.com/articles/s41598-018-20521-y.

[31] Edgar MP, Sun M-J, Gibson GM, Spalding GC, Phillips DB, Padgett MJ. Real-time 3D video utilizing a compressed sensing time-of-flight single-pixel camera. In: Optical trapping and optical micromanipulation XIII, vol. 9922. SPIE; 2016, p. 171-8. http://dx.doi.org/10.1117/12.2239113.

[32] Tropp JA, Gilbert AC. Signal recovery from random measurements via orthogonal matching pursuit. IEEE Trans Inform Theory 2007;53(12):4655-66. http://dx.doi.org/10.1109/TIT.2007.909108, URL https://ieeexplore.ieee.org/document/4385788.

[33] Chengbo L. Compressive sensing for 3D data processing tasks applications, models and algorithms [Ph.D. thesis], Rice University; 2011.

[34] Lee B. Introduction to ±12 degree orthogonal digital micromirror devices (dmds). Tech. rep., Texas Instruments; 2008.

[35] McCarthy A, Krichel NJ, Gemmell NR, Ren X, Tanner MG, Dorenbos SN, et al. Kilometer-range, high resolution depth imaging via 1560 nm wavelength single-photon detection. Opt Express 2013;21(7):8904-15. http://dx.doi.org/10.1364/OE.21.008904, URL http://opg.optica.org/oe/abstract.cfm?URI=oe-21-7-8904.

[36] by Lumibird K. PEFL-KULT pulsed erbium fiber laser brochure V1.0. 2019, URL https://www.keopsys.com/wp-content/uploads/PDF/PEFL-KULT_brochure-v1-0.pdf.

[37] Itzler MA, Entwistle M, Owens M, Patel K, Jiang X, Slomkowski K, et al. Geiger-mode avalanche photodiode focal plane arrays for three-dimensional imaging LADAR. In: Infrared remote sensing and instrumentation XVIII, vol. 7808. SPIE; 2010, p. 75-88. http://dx.doi.org/10.1117/12.861600.

[38] Walker JG. Iterative correction for 'pile-up' in single-photon lifetime measurement. Opt Commun 2002;201(4):271-7. http://dx.doi.org/10.1016/S0030-4018(01)01663-7, URL https://www.sciencedirect.com/science/article/pii/S0030401801016637.

[39] Heide F, Diamond S, Lindell DB, Wetzstein G. Sub-picosecond photon-efficient 3D imaging using single-photon sensors. Sci Rep 2018;8(1):17726. http://dx.doi.org/10.1038/s41598-018-35212-x.

[40] Rapp J, Ma Y, Dawson RMA, Goyal VK. High-flux single-photon LiDAR. Optica 2021;8(1):30-9. http://dx.doi.org/10.1364/OPTICA.403190, URL https://opg.optica.org/optica/abstract.cfm?URI=optica-8-1-30.

[41] Liu X, Ma Y, Li S, Yang J, Zhang Z, Tian X. Photon counting correction method to improve the quality of reconstructed images in single photon compressive imaging systems. Opt Express 2021;29(23):37945-61. http://dx.doi.org/10.1364/OE.443084, URL https://opg.optica.org/oe/abstract.cfm?URI=oe-29-23-37945.

[42] Itzler MA, Krishnamachari U, Chau Q, Jiang X, Entwistle M, Owens M, et al. Statistical analysis of dark count rate in geiger-mode APD FPAs. In: Electro-optical remote sensing, photonic technologies, and applications VIII; and military applications in hyperspectral imaging and high spatial resolution sensing II, vol. 9250. SPIE; 2014, p. 1-12. http://dx.doi.org/10.1117/12.2068744.

[43] Mann HB, Whitney DR. On a test of whether one of two random variables is stochastically larger than the other. Ann Math Stat 1947;18(1):50-60. http://dx.doi.org/10.1214/aoms/1177730491.

[44] Bucchianico A. Combinatorics, computer algebra and the wilcoxon-mann-whitney test. J Statist Plann Inference 1999;79(2):349-64. http://dx.doi.org/10.1016/S0378-3758(98)00261-4, URL https://www.sciencedirect.com/science/article/pii/S0378375898002614.

[45] Natarajan BK. Sparse approximate solutions to linear systems. SIAM J Comput 1995;24(2):227-34. http://dx.doi.org/10.1137/S0097539792240406.

[46] Arjoune Y, Kaabouch N, El Ghazi H, Tamtaoui A. Compressive sensing: Performance comparison of sparse recovery algorithms. In: 2017 IEEE 7th annual computing and communication workshop and conference, vol. 1. IEEE; 2017, p. 1-7. http://dx.doi.org/10.1109/CCWC.2017.7868430, URL https://ieeexplore.ieee.org/document/7868430.

[47] Davis G, Mallat M. Adaptive greedy approximations. Constr Approx 1997;13(1):57-98. http://dx.doi.org/10.1007/BF02678430.

[48] Rubinstein R, Zibulevsky M, Elad M. Efficient implementation of the K-SVD algorithm and the Batch-OMP method. Tech. rep., Technion - Computer Science Departement; 2008.

[49] Anna G, Hespel L, Riviere N, Hamoir D, Tanguy B. Physical modelling of point-cloud (3D) and full-wave-form (4D) laser imaging. In: Electro-Optical Remote Sensing, Photonic Technologies, and Applications IV. 7835, SPIE; 2010, p. 95-109. http://dx.doi.org/10.1117/12.864688.

[50] Coyac A, Riviere N, Hespel L, Briottet X. Comparison of simulated and experimental 3D laser images using a gmapd array: application to long range detection. In: Laser radar technology and applications XXI, vol. 9832. SPIE; 2016, p. 46-53. http://dx.doi.org/10.1117/12.2223603.

[51] Coyac A, Hespel L, Riviere N, Briottet X. Performance assessment of simulated 3D laser images using geiger-mode avalanche photo-diode: tests on simple synthetic scenarios. In: Electro-optical and infrared systems: technology and applications XII; and quantum information science and technology, vol. 9648. SPIE; 2015, p. 22-31. http://dx.doi.org/10.1117/12.2194304.

[52] Roman B, Hansen A, Adcock B. On asymptotic structure in compressed sensing. 2014, http://dx.doi.org/10.48550/ARXIV.1406.4178, URL https://arxiv.org/abs/1406.4178.