# Spline Sketches: An Efficient Approach for Photon Counting Lidar

###### Abstract

Photon counting lidar has become an invaluable tool for 3D depth imaging due
to the fine-precision it can achieve over long ranges. However, high frame
rate, high resolution lidar devices produce an enormous amount of time-of-
flight (ToF) data which can cause a severe data processing bottleneck
hindering the deployment of real-time systems. In this paper, an efficient
photon counting approach is proposed that exploits the simplicity of piecewise
polynomial splines to form a hardware-friendly compressed statistic, or a so-
called spline sketch, of the ToF data without sacrificing the quality of the
recovered image. As each piecewise polynomial spline is a simple function with
limited support over the timing depth window, the spline sketch can be
computed efficiently on-chip with minimal computational overhead. We show that
a piecewise linear or quadratic spline sketch, requiring minimal on-chip
arithmetic computation per photon detection, can reconstruct real-world depth
images with negligible loss of resolution whilst achieving $95\%$ compression
compared to the full ToF data, as well as offering multi-peak detection
performance. These contrast with previously proposed coarse binning histograms
that suffer from a highly nonuniform accuracy across depth and can fail
catastrophically when associated with bright reflectors. Further, by building
range-walk correction into the proposed estimation algorithms, it is
demonstrated that the spline sketches can be made robust to photon pile-up
effects. The computational complexity of both the reconstruction and range
walk correction algorithms scale only with the size of the spline sketch which
is independent to both the photon count and temporal resolution of the lidar
device.

##  1 Introduction

Single-photon light detection and ranging (lidar) is a well-established 3D
depth imaging modality that offers millimeter precision[1] over long
ranges[2]. The technique consists of emitting light pulses and using a single
photon avalanche diode (SPAD) to detect the presence of an incoming photon. By
using a time-correlated single photon counting system (TCSPC), time-of-flight
data can be accumulated over a series of light pulses providing depth and
intensity information for the surfaces in each pixel of the scene. A TCSPC
histogram is most commonly used to represent the ToF data by clustering the
time delay between emitted light pulses and detected photons into time bins
discretized over the whole timing depth window period for each pixel in the
scene (see Figure 1). Rapid hardware and technological advances in recent
years has granted even finer resolution at much quicker frame rates [3, 4].
Modern lidar devices generate massive amounts of data per second that needs to
be transferred off-chip for downstream tasks such as detection, depth
estimation and segmentation. As an example, a high rate, high resolution lidar
device is capable of imaging a scene containing $512\times 512$ pixels at a
frame rate of 50 frames per second (fps). Assuming the timing depth window of
the laser is discretized over $T=1000$ histogram bins of 16 bit precision, the
lidar device would require a data transfer rate of almost 15 GB/s.

To tackle the data-transfer bottleneck of modern lidar devices, practitioners
have often resorted to using coarse timing bin width in the TCSPC histogram.
This technique is often referred to as coarse binning and typically results in
a substantial loss of temporal resolution causing a limiting compression-
resolution trade-off as demonstrated in Figure 1. Gyongy et al. [5] proposed a
practical solution by increasing the width of the impulse response of the
device and using a maximum likelihood estimator to achieve sub-bin resolution
with respect to coarse bins. However, as will be demonstrated in Section 6, we
observe that (1) the worse case error is typically highly dependent on the
position of the signal with respect to the coarse bins, and (2) if there are
highly reflective surfaces (e.g. retroreflectors) in the scene, the measured
IRF can be much narrower than the systems normal impulse response therefore
sub-bin accuracy with respect to the coarse bins cannot be achieved. Recently,
Sheehan et al. [6, 7, 8] proposed a novel solution to tackle the data-transfer
bottleneck whilst circumventing the compression-resolution trade-off. The
technique consists of forming a compact representation, a so-called sketch, of
the ToF data that retains sufficient salient information to accurately
estimate the depth and intensity of multiple surfaces in the scene.
Fundamentally, it was demonstrated that the size of the sketch only needs to
be of the order of the number of surfaces in the scene to achieve negligible
loss of resolution[6]. As there are typically no more than 1 or 2 surfaces
present per pixel, a substantial compression rate can be achieved.
Furthermore, the reconstruction algorithms are based solely on the sketch and
admit both a computational and memory complexity that scales only on the size
of the sketch and is independent of both the resolution of the lidar device as
well as the number of photons detected. However, the method proposed in[6]
requires computing several sinusoidal functions for each photon detection
which may be challenging to execute within the limited dead-time (a period of
insensitivity after a photon detection during which the SPAD cannot register
any incoming photons) of the lidar device.

Fig. 1: True pixelwise TCSPC histogram and its reconstruction (left) compared
with a coarser TCSPC histogram of 50 bins and its reconstruction (right)

In this paper, we introduce new sketch statistics that both captures the
necessary salient information in the data while requiring low on-chip
arithmetic complexity. To this end we propose a spline sketch approach for
photon counting that unites the substantial data transfer compression
capabilities of the Fourier sketched lidar technique with the hardware
friendly computation of a traditional TCSPC histogram. By replacing the
sinusoidal functions of the original sketch in [6] by elementary piecewise
polynomial spline functions, we can construct a compressed representation of
the ToF data that is almost as statistically efficient as the originally
proposed Fourier sketch while having minimal computational overhead on-chip.
The main contributions of this paper are:

  * •

Through the introduction of piecewise polynomial splines, we design low-cost,
hardware-efficient sketches that retain sufficient salient information
required to estimate the depth and intensity parameters of the lidar
observation model and bridge the gap between coarse binning representations
and Fourier sketches whose performance is position independent by definition.

  * •

Both closed-form (for linear splines) and iterative algorithms are proposed
that can accurately estimate the depth and intensity parameters of the scene
while exhibiting a computational and memory complexity that scales solely with
the size of the spline sketch.

  * •

We show that the performance of coarse binning is highly dependent on the
position of the reflector within a given bin. In contrast, linear and
quadratic polynomial spline sketches exhibit only mild variation in such
performance and achieve theoretical error almost equivalent to the Fourier
sketch.

  * •

Finally, we show that the spline sketches carry sufficient information to
enable the correction of range walk associated with bright reflective surfaces
and present a prototype spline sketch range walk correction algorithm.

The paper is organized as follows: In Section 2 the concepts of sketched lidar
and piecewise polynomial splines are introduced. In Section 3, the cost of
computing the spline sketches on-chip is discussed while the spline sketch
based reconstruction algorithms are proposed in Section 4. Both theoretical
Cramér-Rao error bounds and the relative performance of the proposed
algorithms are demonstrated in Section 5 and 6, respectively. Finally, we
conclude the paper in Section 7.

##  2 Background

###  2.1 Lidar Observation Model

During the course of $N$ laser pulses, the SPAD registers $n\leq N$ photon
detections for a given pixel in the scene, where each photon has a detection
time denoted by $x_{j}$ for $1\leq j\leq n$. Assuming there are $K$ distinct
surfaces in the field of view each at a timing depth of $t_{k}\in[0,T-1]$,
then a simple model for the detection time $x_{j}$ of a single photon for a
given pixel is given by a mixture distribution [9]111This model does not
capture all aspects of the SPAD device. However it serves as a useful
idealisation. The sketching framework is not reliant on this specific model
and indeed, a more complicated one will be considered in Section 4.3 when
discussing range walk effects.

$\pi(x_{j}|{\boldsymbol{\theta}})=\sum^{K}_{k=1}\alpha_{k}\pi_{s}(x_{j}|t_{k})+\alpha_{0}\pi_{b}(x_{j}),$
(1)

where $\alpha_{1},\alpha_{2},\dots,\alpha_{K}$ denote the probability of
detection from the $k$th surface, $\alpha_{0}$ denotes the probability of
detection from an ambient background source (e.g. sunlight) and
$\sum_{k=0}^{K}\alpha_{k}=1$. Given that the lidar system has a known impulse
response function (IRF) $h(t)$ then the distribution of photons originating
from a surface is defined as $\pi_{s}(x|t)=h(x-t)/H$ where
$H=\sum_{t=0}^{T-1}h(t)$. The distribution of background photon detections can
be modelled uniformally over the whole timing depth window $\pi_{b}(x)=1/T$.
As a result, single photon counting lidar reduces to estimating the parameters
${\boldsymbol{\theta}}=(\alpha_{0},\alpha_{1},\dots,\alpha_{K},t_{1},\dots,t_{K})$
given a collection of ToF data (e.g. TCSPC histogram) for each pixel in the
scene. As discussed in Section 1, the complexities of transferring ToF data
are often dependent on the depth resolution of the lidar device which is
determined by the discrete number of time bins $T$ used over the whole timing
depth window. The memory requirement for a single pixel in a frame is
typically $\mathcal{O}(Tb)$ where $b$ denotes the integer bit precision which
can be memory expensive in modern day devices with fine depth resolution.

###  2.2 Sketched Lidar

Recently, Sheehan et al. [6, 8] developed a novel solution to capture
sufficient salient information from the ToF data $\\{x_{j}\\}^{n}_{j=1}$ to
estimate the parameters of the observation model in (1). The TCSPC histogram
(see Figure 1) can be seen as a statistic that captures enough information
from the ToF data required for estimation. However it is a very inefficient
representation with memory size scaling proportional to required resolution.
The authors in [6] designed a much more compact representation of the ToF data
that has a size independent of the resolution $T$ of the lidar device and the
number of photon detections $n$.

The basic concept of sketching is to introduce a nonlinear mapping $\Phi(x)$
that defines a sketch statistic as its empirical average,
$\mathbf{z}=\frac{1}{n}\sum_{j=1}^{n}\Phi(x_{j})$ with respect to data $x_{j}$
drawn from a probability density $\pi(x|\theta)$. Then we can estimate the
parameters $\theta$ by comparing the empirical sketch, $\mathbf{z}$ with the
hypothesised expectation: $\mathbb{E}_{x\sim\pi(\cdot|\theta)}\Phi(x)$. In the
original sketched lidar work[6] a particular Fourier sketch proposed that was
inspired by compressive learning theory [10]. Let
$\Phi_{\mathcal{F}}(x)=[e^{{\rm i}\omega_{\ell}x}]_{l=1}^{m}$ denote $m$
complex Fourier features, then the Fourier sketch of size $m$ can be defined
as

$\mathbf{z}_{\mathcal{F}}=\frac{1}{n}\sum_{j=1}^{n}\Phi_{\mathcal{F}}(x_{j}).$
(2)

Due to the summation in (2), the sketch can be updated throughout the pulse
cycle for each photon detection. Once the sketch has been computed, the
parameters of the observation model in (1) can be estimated, e.g. through a
generalized method of moments [11] scheme

$\operatorname*{arg\,min}_{{\boldsymbol{\theta}}}\lVert\mathbf{z}_{\mathcal{F}}-\mathbb{E}_{x\sim\pi}\Phi_{\mathcal{F}}(x)\rVert_{\mathbf{W}}^{2},$
(3)

for some positive definite weighting matrix $\mathbf{W}\in\mathbb{C}^{m\times
m}$. Advantageously for the Fourier sketch, the sketch in expectation is
simply the characteristic function (CF) of the observation model in (1)
sampled at the $m$ frequencies $\omega_{1},\dots,\omega_{m}$. The CF exists
for all distributions and, in the case of the observation model in (1), has a
simple closed form solution which can be easily computed. Fundamentally, the
size of the sketch (i.e. the number of CF samples) required to accurately
estimate ${\boldsymbol{\theta}}$ needs only to be of the order of the number
of surfaces in the pixel, for instance $m=\mathcal{O}(K)$. Moreover, in[6, 8],
the authors propose several sketch-only reconstruction algorithms, including
exploiting powerful spatial denoisers, that have both a memory and
computational complexity of $\mathcal{O}(m)$.

One of the challenges with the original sketched lidar approach is that the
feature function $\Phi_{\mathcal{F}}$ requires nontrivial on-chip processing.
Firstly, the feature function in (2) consists of $2m$ sinusoidal
functions222The original Fourier sketch is a complex number, however it is
possible to exploit the complex structure to simultaneously calculate sin and
cos within FPGA implementations. which is resource intensive, despite the
existence of various FPGA implementations, e.g. [12]. Secondly, the Fourier
sketch is a global sketch as each sinusoidal function has a support equal to
the whole pulse cycle and therefore all $2m$ components of the sketch must be
updated for each photon detection. This is in contrast to the TCSPC histogram
where only a single timing bin needs to be updated for each photon detection.
In this paper we propose to replace the Fourier features $\Phi_{\mathcal{F}}$
with a low-cost, semi-local spline sketch consisting of $M$ piecewise
polynomial spline feature functions that have a limited support over the
timing depth window. While one can intuitively think of the spline functions
as approximating the Fourier representation, it is important to stress that
they define valid sketches in their own right (we are not simply approximating
the Fourier sketch) and as we will see in Section 4.1 have additional pleasing
properties that are unrelated to the previous Fourier sketch.

###  2.3 Cardinal Basis Splines

Cardinal splines are piecewise polynomial functions of degree $p$ that are
used extensively in approximation theory [13]. Let $\phi_{0}$ denote the
spline of degree $p=0$ defined as

$\phi_{0}(x)=\begin{cases}1&x\in[0,1)\\\ 0&\text{otherwise}.\end{cases}$ (4)

Then the degree $p$ spline $\phi_{p}$ is defined as the $p+1$ convolution of
$\phi_{0}$:

$\phi_{p}(x)=(\underbrace{\phi_{0}*\phi_{0}*\dots*\phi_{0}}_{p+1\text{
times}})(x)$ (5)

where $(f*g)(t)=\int_{\mathbb{R}}f(t-s)g(s)\,ds$ denotes the convolution
operator. In this paper, we will only consider the first three degree splines
due to their overall computational simplicity. It can be easily deduced from
(4) that the linear ($p=1$) and quadratic ($p=2$) cardinal splines are given
as

$\phi_{1}(x)=\begin{cases}x&x\in[0,1)\\\ 2-x&x\in[1,2)\\\
0&\text{otherwise}\end{cases}$ (6)

and

$\phi_{2}(x)=\begin{cases}\frac{1}{2}x^{2}&x\in[0,1)\\\
\frac{1}{2}\;+\;(x-1)\;-\;(x-1)^{2}&x\in[1,2)\\\
\frac{1}{2}\;-\;(x-2)\;+\;\frac{1}{2}(x-2)^{2}&x\in[2,3)\\\
0&\text{otherwise}\end{cases}$ (7)

respectively. Cardinal splines of degree $p$ have the appealing properties
that (1) they are $p$ times differentiable and (2) are compactly supported in
$\phi_{p}\in[0,p+1]$. Figure 2 depicts the first three degree cardinal
splines.

Fig. 2: The $p=0,1$ and 2 degree cardinal B splines.

##  3 Spline Sketch

In this paper, we are interested in constructing a sketch consisting of
several splines that together cover the whole support of the timing depth
window $T$, which we will treat as having periodic boundary conditions,
exactly as in the Fourier setting. To do so, let
$\xi_{0},\xi_{1},\dots,\xi_{M}$ be a sequence of equispaced points, called
knot points, over the whole timing depth window $T$ such that $\xi_{0}=0$ and
$\xi_{M}=T$. In addition, let $\Delta\coloneqq T/M$ denote the knot interval
distance between two consecutive knot points. Then the scaled and translated
spline of degree $p$ is defined as

$\phi_{i,p}(x)\;=\;\phi_{p}\left(\frac{x}{\Delta}-i\right),$ (8)

for $i=0,1,\dots,M-1$. The spline $\phi_{i,p}$ is supported on
$\left[i\Delta,(i+p+1)\Delta\right]$. As a result, a spline sketch of size $M$
can be expressed as

$\mathbf{z}_{p}=\frac{1}{n}\sum^{n}_{i=1}\Phi_{p}(x_{i})$ (9)

where $\Phi_{p}(x)=\left[\phi_{i,p}(x)\right]_{i=0}^{M-1}$ is the set of
spline feature functions of degree $p$. Figure 3 illustrates the spline
feature function $\Phi_{p}$ of size $M=4$ for $p=0,1,2$ as well as the real
and imaginary components of the original Fourier sketch.

Fig. 3: The $p=0,1,2$ spline and Fourier sketches of size $M=4$. The dotted
black lines illustrate the partition of the timing depth window into 4 knot
intervals defined by the knot points $\xi_{i}$.

One can see straightaway that for $p=0$ the spline feature function is
equivalent to forming $M$ coarse bins over the whole timing depth window. As
the degree $p$ increases, we observe two important points (1) the support of
each spline feature increases linearly with respect to $p$, (2) each knot
interval $\Delta$ contains $p+1$ spline features. These two observations will
become important when we consider the on-chip implementations of the spline
sketches in Section 3.1. Each of the Fourier features
($\phi_{\omega_{j}}(x)=e^{{\rm i}\omega_{j}x}$) are depicted in the bottom
right of Figure 3. In contrast to the spline features, each individual Fourier
feature has a support equal to the whole timing depth window $T$.

###  3.1 Hardware Implementation

Algorithm 1 details a prototype on-chip algorithm for computing a spline
sketch of degree $p$ in an online manner. Once the appropriate knot interval
of the photon detection has been established, only $p+1$ of the total $M$
sketch features need to be updated. We detail further the total number of
arithmetic operations required per photon detection for each of the $p=0,1,2$
spline sketches.

Initialisation: Degree $p$ , $\mathbf{z}_{p}=0,n=0$

while Acquisition Window do

if New Photon Arrival with Time Stamp $x$ then

Detect Knot interval $x\in\left[i\Delta,(i+1)\Delta\right)$

for $q\xleftarrow{}0\text{ to }p$ do

$\mathbf{z}_{p}=\mathbf{z}_{p}+\phi_{i-q,\Delta,p}(x)$

end for

$n\xleftarrow{}n+1$

end if

end while

$\mathbf{z}_{p}\xleftarrow{}\mathbf{z}_{p}/n$

The sketch $\mathbf{z}_{p}$ is transferred off-chip for post-processing.

Algorithm 1 Sketch Online Processing

For efficiency it makes sense to set $T$ and $M$ to integer powers of 2 so
that identifying the knot interval can be read from the $\log(M)$ leading bits
of the $\log(T)$-bit digital representation of the detection time. The
remaining $b$-bits then represent the relative position within the knot
interval which we denote by $r$. We further note that identifying the knot
interval is the same for $p=0,1$ and $2$. Therefore in examining the different
arithmetic operations we will only focus on the calculation of accumulated
spline function values. For the spline function calculation we assume that the
spline functions are suitably scaled such that exact integer arithmetic can be
used throughout. For the $p=0$ degree spline function $\Phi_{0}$ (i.e. coarse
binning), once we have located the appropriate bin all that is required is to
update the bin count which involves a single one bit addition (e.g.
implemented via a ripple counter). In the case of the $p=1$, two neighboring
sketch values have to be calculated. From (6) we can see the first value is
simply $r$ and requires no further calculation. The second value is
$2^{b+1}-r$ and requires a single addition. These values then have to be added
to the current sketch values requiring a further two additions. Therefore the
linear spline needs to execute a total of 3 addition/subtractions per photon
detection.

For $p=2$, there are three spline values that must be calculated per photon
detection. In terms of the relative knot interval position, $r$, the three
values that require calculation are: $r^{2}$, $2^{2b}-2^{b}r-r^{2}$, and
$2^{2b}-2^{b+1}r+r^{2}$. Notice here that all scaling involved is by powers of
2 and therefore only involves a binary shift (whose computation we ignore).
Once calculated, these values will need to be added to the associated three
sketch components. Therefore the quadratic spline needs to execute a total of
1 multiplication and 7 addition/subtractions per photon detection. In
contrast, while there are various ways to compute the original Fourier sketch
in (2), e.g. CORDIC or local polynomial approximation333i.e.
$\sin(\omega_{j}x)\approx\omega_{j}x-\frac{\omega_{j}^{3}x^{3}}{6}$ and
$\cos(\omega_{j}x)\approx
1-\frac{\omega_{j}^{2}x^{2}}{2}+\frac{\omega_{j}^{4}x^{4}}{24}$ for
$j=1,2,\dots,M/2$.[12], depending on the latency and/or the computation
desired, the computational cost is typically much higher. Furthermore, in each
case all $M$ sketch values need to be updated per photon detection. Specific
details of the arithmetic operations required for each spline sketch per
photon detection is summarized in Table 1.

Table 1: Table highlighting the arithmetic operations required for each sketch per photon detection. Feature Function | Active Features | Add/Sub | Mult.
---|---|---|---
$\Phi_{0}$ | 1 | 1 (1-bit) | 0
$\Phi_{1}$ | 2 | 3 | 0
$\Phi_{2}$ | 3 | 7 | 1

##  4 Spline Sketch Reconstruction

###  4.1 Closed Form Solution

We begin by proposing a simple closed form solution for the piecewise linear
spline sketch444In principle this could also be done for p=2 but we have not
pursued this here that accurately estimates the parameters of the observation
model in (1). The method assumes an IRF with local support and utilises
adjacent linear splines located in the region of the reflected pulse to form a
local mean estimator. Due to the semi-local nature that is inherent to splines
(see Section 3) one can easily establish a subset of splines that are needed
to construct the local mean estimator. First we locate the individual spline
sketch component with the maximum magnitude,

$l\;=\;\operatorname*{arg\,max}_{1\leq j\leq M}\,\mathbf{z}^{j}_{1}$ (10)

where $\mathbf{z}^{j}_{1}$ is the $j$th entry of the sketch $\mathbf{z}_{1}$.
Given that the largest individual sketch $\mathbf{z}_{1}^{l}$ has been located
and assuming that the support of the IRF is less than one knot interval, then
there are only 3 possible scenarios that can exist. Figure 4 depicts each of
the 3 possible scenarios.

Fig. 4: Left column: the 3 (out of 8) spline sketches in the region of the
signal. Middle column: a bar chart of the spline sketches. Left column: local
means estimator of each scenario.

In scenario 1, the IRF is contained within 1 knot interval on the left side of
the $l$th spline. This is reflected in the bar chart (middle) of the spline
sketch values. Similarly, in scenario 2, the IRF is contained within 1 knot
interval on the right side of the $l$th spline. This is also reflected in the
bar chart of the spline sketch values where the second largest sketch value is
to the right of $\mathbf{z}_{1}^{l}$. Finally in scenario 3, the IRF straddles
a knot point and is contained within both knot intervals, as captured by the
bar chart of the spline sketch values. In each scenario, we can build a simple
local mean estimator by taking the difference of adjacent sketch values about
$\mathbf{z}_{1}^{l}$ as is captured in the third column in Figure 4. Notice
that when the IRF straddles a knot point, the local mean estimator is
calculated by the difference between $\mathbf{z}_{1}^{l-1}$ and
$\mathbf{z}_{1}^{l+1}$ hence the estimator is linear over 2 knot intervals.
Specifically, the 3 estimators are calculated as follows

Scenario 1:
$\displaystyle\quad\hat{t}\;=\;\xi_{l-1}+\frac{\Delta}{2}\;+\;\frac{T(\mathbf{z}^{l}_{1}-\mathbf{z}^{l-1}_{1})}{2M\alpha}$
(11) Scenario 2:
$\displaystyle\quad\hat{t}\;=\;\xi_{l}+\frac{\Delta}{2}\;+\;\frac{T(\mathbf{z}^{l+1}_{1}-\mathbf{z}^{l}_{1})}{2M\alpha}$
(12) Scenario 3:
$\displaystyle\quad\hat{t}\;=\;\xi_{l}\;+\;\frac{T(\mathbf{z}^{l+1}_{1}-\mathbf{z}^{l-1}_{1})}{M\alpha}$
(13)

In each case, the depth estimator $\hat{t}$ is dependent on knowing or
accurately estimating $\alpha$. However,this can easily be estimated from
neighbouring spline values that are assumed to only contain background
photons. For the features located in positions that consist solely of photon
detections from ambient sources it can be easily seen (see (1)) that

$\mathbb{E}_{x\sim\pi}\,\phi_{i,p}(x)\;=\;\frac{1-\alpha}{M}.$ (14)

One can therefore use the sketch values $\mathbf{z}^{i}_{1}$ that are known to
only contain background counts to estimate $\alpha$ as follows:

$\hat{\alpha}\;=\;1-M\mathbf{z}^{i}_{1}.$ (15)

If one has access to multiple sketch values that contain photon detections
solely from ambient sources then taking the average over those estimates can
lead to a more accurate $\hat{\alpha}$.

Given the 3 estimators above, one can choose the scenario which produces the
lowest spline sketch loss:

$\lVert\mathbf{z}_{1}-\mathbb{E}_{x\sim\pi}\Phi_{1}(x)\rVert_{2}.$ (16)

As one has to only compute the loss in (16) a total of 3 times, the
computational complexity of this closed form solution is minimal. As the
method returns the (local) mean of the IRF distribution, it is independent of
the IRF shape (incorporating the appropriate leading-edge to mean correction
as appropriate). For additional accuracy one could further fine-tune the
estimates by running a few iterations of a spline version of the SMLE
algorithm proposed in [6] or additionally incorporate powerful spatial
denoisers as proposed in [8]. We evaluate the performance of this method in
Section 6 and compare the error it produces to the optimal Cramér-Rao lower
bound.

###  4.2 Matching Pursuit Algorithm

The closed form estimator proposed in Section 4.1 produces a simple and
computationally cheap solution. However, the solution is dependent on an
accurate estimate of the reflectivity $\alpha$ which can become problematic if
there are multiple reflectors in a scene. In the case of multiple reflectors
per pixel, we propose a robust matching pursuit (MP) inspired algorithm,
similar to the one proposed by Keriven et al. in [14], that iteratively
selects a reflector that is most correlated to the current residual error. In
addition, it can comfortably handle spline sketches of any degree $p$. Details
of the proposed algorithm are given in Algorithm 2.

Initialisation: degree $p$, number of reflectors $K$, $\theta=\emptyset$,
residual $\mathbf{r}=\mathbf{z}_{p}$

for $k\xleftarrow{}1\text{ to }K$ do

Find a reflector highly correlated with the residual
$t^{*}_{k}\xleftarrow{}\,\max_{t}\Big{\langle}\frac{\mathbb{E}_{x\sim\pi_{s}(t)}\,\Phi_{p}(x)}{\lVert\mathbb{E}_{x\sim\pi_{s}(t)}\,\Phi_{p}(x)\rVert_{2}},\mathbf{r}\Big{\rangle}$

$\mathbf{z}_{\theta_{k}}\,\xleftarrow{}\,\mathbb{E}_{x\sim\pi_{s}(t^{*}_{k})}\,\Phi_{p}(x)$

$\alpha_{k}\,\xleftarrow{}\,\frac{\langle\mathbf{z}_{p},\mathbf{z}_{\theta_{k}}\rangle}{\lVert\mathbf{z}_{\theta_{k}}\rVert_{2}}$

$\theta\xleftarrow{}\theta\cup\left\\{\alpha_{k},t_{k}\right\\}$

$\mathbf{r}\xleftarrow{}\mathbf{z}_{p}-\alpha_{k}\mathbf{z}_{\theta_{k}}$

end for

Normalise $\alpha$ such that $\sum_{k=0}^{K}\alpha_{k}=1$.

Estimates for the depth and intensity information of $K$ reflectors in a given
pixel.

Algorithm 2 Spline sketch matching pursuit reconstruction algorithm.

The computational complexity of the MP spline sketch reconstruction algorithm
is $\mathcal{O}(TMK)$ and only takes $K$ iterations to complete. Due to the
semi-local nature of the spline sketches, one can easily initialise $t_{k}$
with a coarse grid in the region of the maximum spline $\mathbf{z}_{p}^{l}$
(see the middle column of Figure 4) and in the specific case of linear
splines, one can use the simple closed form solution discussed in Section 4.1
for a good initialisation. In doing so, the computational complexity reduces
to $\mathcal{O}(T^{\prime}MK)$ where $T^{\prime}\leq T$ is the size of the
local coarse grid. The MP algorithm can handle empirical IRF functions by
approximating $\mathbb{E}_{x\sim\pi}\,\Phi_{p}(x)$ using simple numerical
integration schemes, for instance the trapezium rule [15], which can be stored
within a look-up table. One can further fine-tune the estimates by running a
few iterations of a spline version of the SMLE algorithm proposed in [6] or
additionally incorporate powerful spatial denoisers as proposed in [8]. In
Section 6, we evaluate the MP algorithm and compare the error it produces to
the optimal Cramér-Rao lower bounds.

###  4.3 Range Walk Correction

Lidar systems are subject to systematic range errors that arise due to a
difference between the estimated results and the ground truth distance of the
object in the scene [16, 17]. This can happen due to many aspects of the
detector’s response, but it is particularly linked to imaging highly
reflective objects such as a retroreflector. Retroreflectors are ultra
reflective surfaces that lead to a high probability of photon detection per
laser pulse. Due to the intrinsic deadtimes of the lidar device, the SPAD
detection deviates from the classical Poisson statistical model and more
sophisticated models need to be considered. Here we will use a model similar
to the one presented in[17] to examine the associated range walk effects and
possible sketch corrections. The photon flux $\lambda_{i}$ on the SPAD during
the $i$th time bin of the timing depth window, denoted by the interval
$I_{i}$, can be modelled by an inhomogeneous Poisson process with rate
function

$\lambda_{i}\;=\;\int_{I_{i}}\mu\left(\beta\frac{h(x-t)}{H}+s\right)\;dx$ (17)

where $\mu$ is the SPAD photon detection probability ($\sim 1\%$),
$\beta\in[0,1]$ is the target reflectivity and $s$ is the level of detections
from background sources. In the medium and large flux regimes, the SPAD tends
to detect photon arrivals at earlier bins along the timing depth window,
resulting in an observed IRF that is distorted towards the leading edge of the
original pulse. This phenomena is commonly referred to as photon pile-up and
can cause large range errors for highly reflective surfaces.

Fig. 5: Distorted IRFs as a function of the target reflectivity.

Figure 5 depicts the distortion of a real life IRF from the face dataset in
[18] simulated using the multinomial distribution in [17] for different target
reflectivity $\beta$ and with a SPAD photon detection probability of
$\mu=0.01$. It can be seen that as the target reflectivity increases to 1,
photon detections accumulate more frequently towards the rising edge of the
original IRF with larger intensity. As the location of the reflecting surface
can be identified with the leading edge of the observed IRF, any estimator
based on the peak or the mean of the distribution will need a reflectivity
dependent correction, often called range walk correction, in order to
accurately locate the leading edge. Intensity based correction models are a
well-established technique in range walk correction [16, 19] and can be easily
integrated into a spline sketch estimation algorithm as any degree spline
sketch contains intensity information. The left of Figure 6 shows the range
walk correction as a function of the normalized intensity calculated using a
spline sketch of size $M=50$. As with all intensity based range walk models,
the range walk error can be uniquely corrected up until a point of saturation.
This saturation occurs when the probability of detection is essentially 1
(i.e. $n\approx N$) and happens between the medium and high flux regimes.
After the saturation point has been hit, it is impossible to distinguish
between different range-walk errors using only intensity which can lead to
large estimation errors (as it will be demonstrated in Figure 13).

Fig. 6: Left: A look-up table evaluating intensity to range walk error. Right:
A look-up table evaluating the standard deviation of the observed IRF to range
walk error.

In order to accurately correct for range walk within the saturated regime, it
is necessary to capture higher order statistics related to the shape of the
observed IRF. This can be easily done using a quadratic spline sketch as they
also capture local second order statistics of the ToF data allowing one to
distinguish between the variances of the IRFs. It can be seen in Figure 5 that
as the target reflectivity increases, the observed IRF becomes ever more
concentrated about the leading edge. The right of Figure 6 shows the range
walk error as a function of the local estimate of the standard deviation which
can be calculated using the quadratic spline sketch ($M=50$) following the
same principle as the local mean estimators presented in Section 4.1. In
contrast to the intensity based range walk model, the shape based range walk
model defines a robust look-up table correction over the full range of
reflectivities. We therefore have a unique and well-defined correction scheme.

In a similar fashion to the closed form estimator in Section 4.1, one can
construct a local standard deviation estimator in closed form that is both
simple and computationally cheap. First we locate the individual spline sketch
with the maximum magnitude

$l\;=\;\operatorname*{arg\,max}_{1\leq j\leq M}\,\mathbf{z}^{j}_{2}$ (18)

where $\mathbf{z}^{j}_{2}$ is the $j$th entry of the quadratic spline sketch
$\mathbf{z}_{2}$. Next, we remove the noise contribution from the individual
sketches by defining the noise corrected sketches as

$\tilde{\mathbf{z}}^{j}_{2}\;=\;\frac{1}{\alpha}\left(\mathbf{z}^{j}_{2}-\frac{1-\alpha}{M}\right).$
(19)

By letting

$\tilde{\mathbf{z}}_{\text{sub}}\;=\;\left[\tilde{\mathbf{z}}^{l-2}_{2},\tilde{\mathbf{z}}^{l-1}_{2},\tilde{\mathbf{z}}^{l}_{2},\tilde{\mathbf{z}}^{l+1}_{2},\tilde{\mathbf{z}}^{l+2}_{2}\right]$
(20)

define the collection of noise corrected sketches in the neighbourhood of
$\mathbf{z}_{2}^{l}$, then it can be easily shown that the depth and shape
estimates have the closed form solution

$\hat{t}\;=\;\xi_{l}+\mathbf{c}_{1}^{T}\tilde{\mathbf{z}}_{\text{sub}},$ (21)

and

$\hat{\sigma}^{2}\;=\;\mathbf{c}_{2}^{T}\tilde{\mathbf{z}}_{\text{sub}}\;-\;\left(\mathbf{c}_{1}^{T}\tilde{\mathbf{z}}_{\text{sub}}\right)^{2},$
(22)

respectively, where

$\displaystyle\mathbf{c}_{1}$
$\displaystyle=[-2\Delta,-\Delta,0,\Delta,2\Delta],$
$\displaystyle\mathbf{c}_{2}$
$\displaystyle=\left[\frac{15}{4}\Delta^{2},\frac{3}{4}\Delta^{2},-\frac{1}{4}\Delta^{2},\frac{3}{4}\Delta^{2},\frac{15}{4}\Delta^{2}\right]$

and recalling that $\Delta=\frac{T}{M}$.

Here we define the standard deviation estimate in terms of 5 neighbouring
spline sketches due to the slightly larger support of the original IRF for the
face dataset, however one can easily construct equivalent closed form
estimates for larger or smaller neighbourhoods555See Section 8 for details on
MATLAB examples. . As with the local mean estimates we need to use a
background level estimate. Further details can be found in the code.

##  5 Cramér-Rao Error Bounds

In this section, the Cramér-Rao error bounds of the spline sketches for
$p=0,1,2$ are compared with the original Fourier sketch as well as the full
TCSPC histogram (i.e. no compression). The Cramér-Rao bound (CRB) gives a
lower bound for the root mean squared error (RMSE) of an estimator
$\hat{\theta}$ and can therefore provide the best case performance one could
achieve through a reconstruction algorithm. Given the observation model $\pi$
in (1) and the corresponding Fisher information matrix (FIM), defined as

$\mathcal{I}_{\text{data}}(\theta)\;=\;n\mathbb{E}_{x\sim\pi}\Bigg{[}\Bigg{(}\frac{\partial\log\pi(x\mid\theta)}{\partial\theta}\Bigg{)}^{2}\Bigg{]},$
(23)

then the optimal Cramér-Rao RMSE, in terms of the full data, is defined as

$\text{RMSE}_{n}\;=\;\sqrt{\sum^{2K}_{k=1}[\mathcal{I}_{\text{data}}(\theta)^{-1}]_{\\{kk\\}}}.$
(24)

Equivalently, for the sketched case the FIM is defined as (see for instance
[6])

$\big{(}\mathcal{I}_{\text{sketch}}(\theta)\big{)}_{ij}\;=\;n\dfrac{\partial\mathbf{z}_{\theta}}{\partial\theta_{i}}\Sigma^{-1}_{\theta}\dfrac{\partial\mathbf{z}_{\theta}}{\partial\theta_{j}},$
(25)

where for shorthand we denote
$\mathbf{z}_{\theta}=\mathbb{E}_{x\sim\pi}\Phi_{p}(x)$ and where
$\Sigma_{\theta}$ denotes the covariance matrix of the spline sketch.
Similarly, we define the optimal sketched Cramér-Rao RMSE as

$\text{RMSE}_{M}\;=\;\sqrt{\sum^{2K}_{k=1}[\mathcal{I}_{\text{sketch}}(\theta)^{-1}]_{\\{kk\\}}}.$
(26)

In this section, the IRF is a Gaussian function with pulse width $\sigma$.
This ensures that the formula in (26) has a closed form solution, however, one
can expect similar results to hold for more general IRFs. In addition, as many
of the environmental factors need to be fixed, we choose ones that reflect
typical lidar scenes. As the CRB is dependent on many factors, the following
results vary one factor at a time so that one can understand the effect.

###  5.1 Spatially Averaged CRB Performance

Fig. 7: CRB for varying SBR for a Gaussian pulse width fixed at 64cm. Sketch
size fixed at $M=8$.

We begin exploring the spatially averaged CRB performance of the different
sketches. Each result is averaged along 1000 uniform depths across the timing
depth window $[0,T-1]$. In our first experiment, we vary the SBR from a low
($10^{-1}$) to a high ($10$) SBR regime. Here the acquisition time corresponds
to a depth range $24$m where each of the $T$ timing bins equates to a
precision of $4$cm. The sketch size for each of the spline and Fourier
sketches is fixed at $M=8$ and the number of photon detections is $n=1000$.
Figure 7 shows the CRB for the different sketches for a Gaussian pulse width
of $\sigma=64$cm for spline sketches of degree $p=0,1,2$, a Fourier sketch and
a full data TCSPC histogram (i.e. no compression). First, while the degree
$p=0$ spline sketch (i.e. coarse binning) is able to achieve good performance
at high SBR this reduces as the SBR gets less. Furthermore the $p=0$ RMSE is
the highest of all the sketches. There is a significant reduction of RMSE as
the degree of the spline is increased to $p=1$. Notably, the improvement of
using a degree $p=2$ spline sketch is negligible in comparison to its linear
counterpart. This is potentially due to the fact that locally first order
statistics are almost sufficient in capturing the depth information and
therefore the RMSE saturates for $p\geq 1$. There is a slight loss of
information in comparison to using the original Fourier sketch of the same
size, however one gains in the computational simplicity of constructing the
spline sketches. Nonetheless, for a moderate SBR of 1, the RMSE of all the
estimators apart from coarse binning are approximately the same at 5cm.

Figure 8 shows the results for the same set up as before but now with a
narrower Gaussian pulse width of $\sigma=24$cm (as might occur in medium-high
flux scenarios). In this case, one can clearly see that the coarse binning
method fails catastrophically even at high SBR levels (the reason for this
will become clear in the next subsection). In contrast, the other methods
achieve similar performance to before and in fact are slightly more accurate.
In both Figure 7 and Figure 8, we see that each of the $p=1,2$ degree splines
and the Fourier sketches converge to the devices resolution (4cm in this case)
as the SBR becomes high.

Fig. 8: CRB for varying SBR for a Gaussian pulse width fixed at 24cm. Sketch
size fixed at $M=8$.

To further understand the relationship between the pulse width and the RMSE of
the estimators, we fix the SBR at 1 and vary only the pulse width within the
range of 2 and 100 cm. Figure 9 demonstrates that the coarse binning RMSE is
extremely dependent on the pulse width of the IRF, whilst, in contrast, the
other methods are agnostic to the pulse shape. Some practitioners, for
instance [5], have proposed increasing the pulse width of the device’s laser
to achieve a coarse binning RMSE that is only slightly worse than using either
the sketches or full data. However, this is only practical in certain
circumstances, and if the device experiences photon pile-up (see Section 4.3)
then the width of the detected signal can collapse towards a Dirac function.
In such an event, coarse binning would fail catastrophically as illustrated in
Figure 9.

Fig. 9: CRB for varying Gaussian pulse width with the SBR fixed at 1 and a
sketch size fixed at $M=8$.

###  5.2 CRB Spatial Dependence

By analysing the poor performance of coarse binning with narrow IRFs reveals a
more troubling phenomenon. The results presented in the previous subsection
averaged the CRBs over the whole timing depth window. However, the CRBs, with
the exception of the Fourier sketch, do vary as a function of target depth.
Therefore the average CRBs may be hiding a substantially worse performance at
certain depths. We therefore next investigate the spatial dependency of the
CRBs for the different sketches. To do so, the SBR and pulse width are fixed
to 1 and 64cm, respectively, and the depth of the reflector is varied along
the range of the lidar device. Figure 10 shows the results. As expected, the
CRBs for the full data and Fourier sketch are independent of the depth of the
reflector. In contrast, the coarse binning estimator ($p=0$) varies from a
small RMSE at the knot point to a significantly larger maximum error ($\sim
28$cm) at the centre of the knot interval. This can be intuitively understood
as follows. If the return pulse straddles a knot point then its precise
position can be estimated accurately from the relative magnitudes of the
neighboring histogram bins. However, if the pulse lies solely in one knot
interval then the estimator is blind to where the pulse is within the knot
interval. This dependency also occurs to a certain extent for $p=1,2$ degree
spline sketches but at a much dampened state. Interestingly, the $p=1$ spline
sketch achieves the largest RMSE at the knot points and the smallest RMSE at
the centre of the knot intervals which is in contrast to the even degree
splines sketches of $p=0,2$. This result importantly highlights the worst case
performance one can expect to achieve from each spline sketch which, as
demonstrated, can be extreme for coarse binning yet is controllable for both
the linear and quadratic spline sketches.

Fig. 10: CRB for varying depth positions with the SBR fixed at 1, the Gaussian
pulse width fixed at 64cm and a sketch size fixed at $M=8$.

##  6 Simulations

In the section we evaluate the proposed algorithms on both synthetic and real
datasets and compare their performance to the CRB discussed in Section 5.

###  6.1 Algorithm Dependent Performance

While the CRBs provide us with performance bounds for any estimator using the
associated statistics, it is also necessary to understand how a practical
algorithm performs with respect to the associated CRB. We therefore conclude
this section by showing how the different algorithms presented in section 4
compare against the CRBs. The set up for this experiment is identical to the
one presented in Figure 7. We vary the SBR from a $10^{-1}$ to $10$. The depth
range is $24$m with $T$ timing bins equating to a precision of $4$cm each. The
sketch size is $M=8$, the number of photon detections is $n=1000$, and the IRF
is modelled as a Gaussian pulse with width $\sigma=64$cm. Figure 11 shows the
that all algorithms perform equivalently at high SBR. While for low SBR all
algorithms exhibit a small drop in performance compared with their associated
CRB. The MP algorithm appears to perform similarly for both $p=1$ and $p=2$
while the closed form estimate is slightly worse. This is most probably due to
the assumption that the support of the IRF does not exceed a single know
interval. This could be remedied by additional fine-tuning as discussed in
(11).

Fig. 11: A comparison of the proposed reconstruction algorithms and their
respective CRBs with the SBR fixed at 1, the Gaussian pulse width fixed at
64cm and a sketch size fixed at $M=8$.

###  6.2 Real Dataset

We evaluate the proposed algorithms in Section 4 on a real dataset consisting
of a polystyrene head placed at a distance of 40 metres[18]. We compare the
proposed methods with the cross correlation algorithm[1] that works directly
on the full TCSPC (i.e. no compression), as well as the SMLE algorithm in [6]
that uses Fourier features to construct the sketch. The polystyrene head
dataset has size of $141\times 141$ pixels with $T=4613$. Most of the pixels
in this scene contain exactly $K=1$ surface and the average number of photon
detections per pixel is 337. Figure 12 and table 2 show the reconstruction of
the scene and the corresponding RMSE metrics, respectively, using a Fourier
sketch and the SMLE algorithm in [6], coarse binning ($p=0$), the closed form
solution in Section 4.1, and the $p=1,2$ spline sketches using the MP
algorithm in Section 4. These compression algorithms are compared to the
matched filtering algorithm [1] that acts on the original uncompressed TCSPC
histogram. For the compression techniques, we consider a sketch size of
$M=10,20,30$ and 40. Initially, we see that the coarse binning method fails to
accurately recover the scene for the smaller sketch sizes. In comparison, both
$p=1,2$ degree splines manage to retain enough information from the ToF data
to allow for accurate reconstruction. One can observe a slight drop in
performance between the closed form solution and the matching pursuit
algorithm especially around the border of the head. Interestingly, a linear
and quadratic spline sketch of just $M=20$ allows for very accurate
reconstruction in comparison to both the Fourier sketch and the full data
reconstruction. For this specific dataset, the linear and quadratic spline
sketches correspond to a compression of approximately $95\%$ in comparison to
storing the raw photon detection times without sacrificing the overall quality
of the reconstructed image.

Fig. 12: Reconstruction of the face dataset using a Fourier sketch, coarse binning ($p=0$), the local means estimator in Section 4.1, and the $p=1,2$ spline sketches using the MP algorithm. Sketch Size $M$ | Cross Corr. | Fourier | Coarse Bin. | LME ($p=1$) | Spline MP ($p=1$) | Spline MP ($p=2$)
---|---|---|---|---|---|---
10 |  | 8.2 | 74.5 | 15.3 | 12.1 | 11.7
20 | 4.4 | 6.2 | 22.8 | 11.4 | 8.4 | 8.5
30 |  | 4.8 | 18.1 | 8.6 | 6.2 | 6.4
40 |  | 4.6 | 15.1 | 7.0 | 5.7 | 5.9
Table 2: RMSE (bins) of the different compression techniques for a sketch of
size $M=6,12$ and 24, respectively.

###  6.3 Range Walk Correction Simulations

To evaluate our shaped based range walk model discussed in Section 4.3, we use
the IRF in[18] where $T=4613$ (see Section 6.2 for more details) and build a
look-up table as in Figure 6 for a sketch size of $M=25,50,\dots,125$ using
the model in [17]. Once an intensity or standard deviation value has been
estimated from the spline sketches, we consult the aforementioned look-up
tables using the estimates to correct for the range walk. For 50 different
target reflectivities (on a log-scale between $10^{-4}$ and 1) and SBR values
of 100,10,1 and 0.1, we compute the root mean squared error (RMSE) calculated
over 250 random target distance initialisations. Figure 13 shows the RMSE for
each of the SBR values as a function of target reflectivity. One can initially
see that as the point of intensity saturation is hit, the RMSE grows rapidly
as the intensity based range walk model is unable to accurately distinguish
the correct range walk. In contrast, the shape based range walk model stays at
quite a constant error level throughout despite the increase in target
reflectivity. For the lowest SBR of 0.1, the smaller sketch sizes of $M=25$
and $M=50$ have a larger RMSE for the higher target reflectivity. This is
because in challenging imaging scenes, a larger sketch size is required to
accurately estimate and remove the contribution from background sources. These
results demonstrate additional statistical information captured by the
quadratic spline sketch provide further benefits in complex environments.

Fig. 13: The RMSE of both the intensity and shape based range walk correction
models for different SBR levels.

##  7 Conclusion

In this paper, we proposed a semi-local spline sketch that can be efficiently
implemented in hardware and that achieves the compression capabilities of the
original Fourier sketch whilst enjoying minimal overhead computations on-chip.
Simulations demonstrated that both the linear and quadratic splines sketches
achieve almost the same performance as the original sketch with only a slight
drop in performance for very low SBR scenes. Increasing the degree of the
spline sketch to $p=2$ also enabled us to capture local 2nd order statistics
that allow us to accurately correct for range walk error where intensity based
range walk models fail to do so.

##  8 Code Availability

A MATLAB implementation of all the algorithms discussed are available at the
repository https://gitlab.com/mpsheehan1995/spline-sketch-lidar.

## References

  * [1] A. McCarthy, R. J. Collins, N. J. Krichel, V. Fernández, A. M. Wallace, and G. S. Buller, “Long-range time-of-flight scanning sensor based on high-speed time-correlated single-photon counting,” _Appl. Opt._ , vol. 48, no. 32, pp. 6241–6251, Nov 2009.
  * [2] A. M. Pawlikowska, A. Halimi, R. A. Lamb, and G. S. Buller, “Single-photon three-dimensional imaging at up to 10 kilometers range,” _Opt. Express_ , vol. 25, no. 10, pp. 11 919–11 931, May 2017.
  * [3] R. K. Henderson, N. Johnston, H. Chen, D. D. Li, G. Hungerford, R. Hirsch, D. McLoskey, P. Yip, and D. J. S. Birch, “A 192×128 time correlated single photon counting imager in 40nm CMOS technology,” in _ESSCIRC 2018 - IEEE 44th European Solid State Circuits Conference (ESSCIRC)_ , 2018, pp. 54–57.
  * [4] S. W. Hutchings, N. Johnston, I. Gyongy, T. Al Abbas, N. A. W. Dutton, M. Tyler, S. Chan, J. Leach, and R. K. Henderson, “A reconfigurable 3-D-stacked SPAD imager with in-pixel histogramming for flash lidar or high-speed time-of-flight imaging,” _IEEE Journal of Solid-State Circuits_ , vol. 54, no. 11, pp. 2947–2956, 2019.
  * [5] I. Gyongy, S. W. Hutchings, A. Halimi, M. Tyler, S. Chan, F. Zhu, S. McLaughlin, R. K. Henderson, and J. Leach, “High-speed 3D sensing via hybrid-mode imaging and guided upsampling,” _Optica_ , vol. 7, no. 10, pp. 1253–1260, 2020.
  * [6] M. P. Sheehan, J. Tachella, and M. E. Davies, “A sketching framework for reduced data transfer in photon counting lidar,” _IEEE Transactions on Computational Imaging_ , vol. 7, pp. 989–1004, 2021.
  * [7] ——, “Surface detection for sketched single photon lidar,” in _2021 29th European Signal Processing Conference (EUSIPCO)_ , 2021, pp. 621–625.
  * [8] J. Tachella, M. P. Sheehan, and M. E. Davies, “Sketched RT3D: How to reconstruct billions of photons per second,” in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 1566–1570.
  * [9] Y. Altmann and S. McLaughlin, “Range estimation from single-photon lidar data using a stochastic em approach,” in _2018 26th European Signal Processing Conference (EUSIPCO)_ , 2018, pp. 1112–1116.
  * [10] R. Gribonval, A. Chatalic, N. Keriven, V. Schellekens, L. Jacques, and P. Schniter, “Sketching data sets for large-scale learning: Keeping only what you need,” _IEEE Signal Processing Magazine_ , vol. 38, no. 5, pp. 12–36, 2021.
  * [11] A. Hall, _Generalized Method of Moments_. Oxford University Press, 11 2007, pp. 230 – 255.
  * [12] F. de Dinechin, M. Istoan, and G. Sergent, “Fixed-point trigonometric functions on FPGAs,” vol. 41, no. 5, p. 83–88, jun 2014. [Online]. Available: https://doi.org/10.1145/2641361.2641375
  * [13] M. Unser, “Splines: a perfect fit for signal and image processing,” _IEEE Signal Processing Magazine_ , vol. 16, no. 6, pp. 22–38, 1999.
  * [14] N. Keriven, A. Bourrier, R. Gribonval, and P. Pérez, “Sketching for large-scale learning of mixture models,” _Information and Inference: A Journal of the IMA_ , vol. 7, no. 3, pp. 447–508, 2018.
  * [15] A. H. Stroud, “Methods of numerical integration,” _SIAM Review_ , vol. 28, no. 1, pp. 128–128, 1986. [Online]. Available: https://doi.org/10.1137/1028042
  * [16] A. Tontini, L. Gasparini, and M. Perenzoni, “Numerical model of spad-based direct time-of-flight flash lidar cmos image sensors,” _Sensors_ , vol. 20, no. 18, 2020. [Online]. Available: https://www.mdpi.com/1424-8220/20/18/5203
  * [17] F. Heide, S. Diamond, D. B. Lindell, and G. Wetzstein, “Sub-picosecond photon-efficient 3D imaging using single-photon sensors,” _Scientific reports_ , vol. 8, no. 1, pp. 1–8, 2018.
  * [18] Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, and S. McLaughlin, “Lidar waveform-based analysis of depth images constructed using sparse single-photon data,” _IEEE Trans. Image Process._ , vol. 25, no. 5, pp. 1935–1946, 2016.
  * [19] M. S. Oh, H. J. Kong, T. H. Kim, K. H. Hong, and B. W. Kim, “Reduction of range walk error in direct detection laser radar using a geiger mode avalanche photodiode,” _Optics Communications_ , vol. 283, no. 2, pp. 304–308, 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0030401809009663