# Robust super-resolution depth imaging via a multi-feature fusion deep
network

Alice Ruget1, Stephen McLaughlin1, Robert K. Henderson2, Istvan Gyongy2,
Abderrahim Halimi1, Jonathan Leach1,∗ 1School of Engineering and Physical
Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK  2School of
Engineering, Institute for Integrated Micro and Nano Systems, The University
of Edinburgh, Edinburgh, EH9 3FF, UK<EMAIL_ADDRESS>
###### Abstract

The number of applications that use depth imaging is increasing rapidly, e.g.
self-driving autonomous vehicles and auto-focus assist on smartphone cameras.
Light detection and ranging (LIDAR) via single-photon sensitive detector
(SPAD) arrays is an emerging technology that enables the acquisition of depth
images at high frame rates. However, the spatial resolution of this technology
is typically low in comparison to the intensity images recorded by
conventional cameras. To increase the native resolution of depth images from a
SPAD camera, we develop a deep network built to take advantage of the multiple
features that can be extracted from a camera’s histogram data. The network is
designed for a SPAD camera operating in a dual-mode such that it captures
alternate low resolution depth and high resolution intensity images at high
frame rates, thus the system does not require any additional sensor to provide
intensity images. The network then uses the intensity images and multiple
features extracted from down-sampled histograms to guide the up-sampling of
the depth. Our network provides significant image resolution enhancement and
image denoising across a wide range of signal-to-noise ratios and photon
levels. Additionally, we show that the network can be applied to other data
types of SPAD data, demonstrating the generality of the algorithm.

##  I Introduction

Figure 1:  Overview of the proposed network. HistNet is designed for a SPAD
camera operating in a dual-mode. The camera provides a histogram of photon
counts and an intensity image. The network takes this as an input and provides
a HR depth map as the output with a resolution that is four times higher in
both spatial dimensions than the initial raw histogram.

Light detection and ranging (lidar), where a pulse of light is used to
illuminate a target and a detector provides time-of-flight information, is one
of the leading technologies for depth imaging. For example, lidar is one of
the key systems for future connected and autonomous vehicles, and it is used
in the latest smartphones and tablet to aid auto-focus assist and enhance
virtual reality. Single-photon avalanche detector (SPAD) arrays are an
emerging technology for depth estimation via lidar. These are devices that are
sensitive to single photons and can provide histograms of times of arrival of
a single photons with respect to a trigger. When used in combination with a
pulsed laser that illuminates a target object, SPAD arrays provide accurate
and fast data that can be converted to depth information.

In the context of lidar, several different SPAD array sensors have been
developed, see Ximenes et al. (2018); Zhang et al. (2018) for recent examples.
They have been used to measure depth in a range of scenarios, including under
water Maccarone et al. (2015); Halimi et al. (2017), long range Chan et al.
(2019); McCarthy et al. (2013); Pawlikowska et al. (2017), at high speed
Hutchings et al. (2019); Henderson et al. (2019); Gyongy et al. (2020); Bronzi
et al. (2014), and providing high-resolution depth information Ren et al.
(2018); Zhang et al. (2018). Recently, Morimoto et al. reported a mega-pixel
SPAD array Morimoto et al. (2020). SPAD array sensors have also been used for
light-in-flight imaging Gariepy et al. (2015) and looking at objects hidden
around corners Gariepy et al. (2016). They have also been used extensively
within the field of biophotonics, see ref Bruschini et al. (2019) for a
review.

Although SPAD arrays are becoming well established in lidar systems, there are
several key challenges to overcome to fully exploit their potential. The
single-photon sensitivity that the SPAD array provides promises depth imaging
at long ranges and in degraded visual environments, but improving the
performance in these scenarios can dramatically increase the range of use of
the detectors. In addition, the native resolution is typically very low in
comparison to conventional image sensors. Ultimately, it is desirable to
operate the SPAD arrays at high frame rates, cover a large field-of-view at
large distances, produce images at high resolutions, and perform well in a
wide range of environmental conditions. Each of these objectives brings
separate challenges that need to be addressed.

Due to the nature of the challenges to single-photon depth imaging,
computational post-processing techniques are known to be a very powerful
method to improve the overall image quality, both in terms of signal-to-noise
and resolution. The latter methods Tachella et al. (2019); Rapp and Goyal
(2017); Halimi et al. (2019); Guo et al. (2019); Sun et al. (2020); Lindell et
al. (2018) take advantage of prior information in one form or another and
attempt to improve the quality of the depth images in the low-photon regime.
There are advantages and disadvantages to each of the methods, often with a
trade off in terms of quality of the reconstruction and for the time taken to
reconstruct.

Several statistical approaches have been implemented to improve the depth maps
of single-photon depth data: Tachella et al. established their images making
use of priors on the depth and intensity estimates, achieving fast
reconstruction and very good performance in the single-photon regime Tachella
et al. (2019); Rapp and Goyal tackled the problem by the creation of super-
pixels that borrow information from relevant neighbours to separate the
contributions of signal and background Rapp and Goyal (2017); and Halimi et
al. implemented an alternating direction method of multipliers (ADMM) to
minimize a cost function using priors on correlations between pixels using
both depth and intensity estimates Halimi et al. (2019). Callenberg et al.
implemented an iterative optimization scheme to increase the spatial
resolution of SPAD array data using a sensor fusion approach Callenberg et al.
(2021). We refer the reader to Wallace et al. (2020); Rapp et al. (2020) for
more details regarding state-of-the-art robust reconstruction algorithms.

Machine learning approaches have also shown good performance when enhancing
the quality of high-resolution depth data. For example, Guo et al developed a
deep neural network to reconstruct a high-resolution depth map from its low
resolution version Guo et al. (2019). In addition, Lindell et al Lindell et
al. (2018) and Sun et al Sun et al. (2020) developed deep networks that
process the whole 3D volume of raw photon counts and output a 2D depth map,
achieving high performance in the low photon regime at the cost of a long
processing time. All references Guo et al. (2019); Lindell et al. (2018); Sun
et al. (2020) make use of an intensity image to guide the reconstruction of
the depth. The networks in references Lindell et al. (2018); Sun et al. (2020)
are specifically designed for SPAD array data, whereas the work in Guo et al.
(2019) processes existing depth maps. Peng et al. Peng et al. (2020)
implemented a network capable of reconstructing depth from SPADs at very long
ranges by exploiting the non local correlations in time and space via the use
of non local blocks and down-sampling operators. Recently Nishimura et al.
Nishimura et al. (2020) demonstrated how a single depth histogram can be used
to resolve the depth ambiguity problem associated with monocular depth
estimation.

This paper proposes and implements a machine learning network to
simultaneously perform up-sampling and de-noising of depth information. We
design and apply the network so that it is suitable for the data provided by
the Quantic 4x4 sensor Hutchings et al. (2019); Henderson et al. (2019);
Gyongy et al. (2020), but the network can be applied to other depth images.
The SPAD camera alternates between two modes at over 1000 frames per second.
It provides high-resolution intensity images at a resolution of 256x128 pixels
followed by low-resolution 64x32x16 histogram of photon counts containing
depth information.

Our approach is to select the essential information from the histogram and
input this directly to the network, without requiring the entire histogram of
counts to be provided to the network. After processing, the final resolution
of our up-sampled depth images from the Quantic 4x4 sensor is increased by a
factor of four to 256x128 pixels. Figure 1 shows an overview of the results of
the proposed network when applied on captured data.

This paper is organized as follows: in Section II, we provide a brief overview
of the SPAD array sensor, the model of photon detection, and we present the
processing done to the SPAD data to extract useful information prior to the
reconstruction via the network. Section III introduces the proposed HistNet in
details. Section IV reports the results on both simulated and real data along
with a comparison to other algorithms, and demonstrates its robustness to
different noise scenarios. Section V presents our conclusions and future work.

##  II SPAD array data

###  II.1 Data acquisition

Single-photon avalanche diode arrays can capture depth and intensity
information of a scene. To achieve this, a short laser pulse is used to
illuminate a target, and the detector records the arrival time of photons
reflected back by the scene with respect to a laser trigger. This data, known
as time tagged data, can be used to generate a temporal histogram of counts
with respect to the time of flights, where the peak in the histogram can be
used to calculate the distance to the target.

In this work, we develop a network suitable for a SPAD array sensor, the
Quantic 4x4 sensor, that generates a histograms of counts on-chip and operates
in a hybrid acquisition mode Hutchings et al. (2019); Henderson et al. (2019);
Gyongy et al. (2020). This hybrid mode alternates between two measurement
modes at a temporal rate exceeding 1000 frames per second. The details of the
modes are: a high-resolution (HR) intensity measurement with a spatial
resolution of 256x128, and a low-resolution (LR) histogram of photon counts
containing depth information at a resolution of 64x32x16 (16 being the number
of time bins of each of the 64x32 histograms). It is the purpose of the
network to increase the resolution of the initial depth data (64x32) to the
same resolution as the intensity data (256x126), while simultaneously
denoising the data.

###  II.2 Pre-processing of data for the network

The SPAD camera provides LR histogram data and HR intensity images in
alternate frames. We will see in the following sections how we select features
from the SPAD array data that maximise the quality of information provided to
the network while minimising the total quantity of data necessary for accurate
super-resolution. There are several different features that we extract from
the data provided by the SPAD array: the first depth map, the second depth
map, the high-resolution intensity image, and the multi-scale depth features
extracted from down-sampled versions of the original histogram. The processing
time to calculate each of these features is minimal, adding very little
computational overhead to our overall procedure. Figure 2 shows SPAD array
data and the different processing steps to compute the arguments of HistNet.

Figure 2:  Representation of the processing of the different arguments of
HistNet. The SPAD array provides alternating LR Histograms of size 64x32x16
and HR intensity images of size 256x128. The 256x128 first and second depth
maps are obtained by computing center of mass around the strongest and second
strongest peaks in the raw histogram and then up-sampling by four in both
spatial dimensions using nearest neighbour interpolation. Depth feature D1 of
size 128x64 is obtained by down-sampling the first depth map by two in both
dimensions. D2 of size 64x32 is obtained by applying center of mass around the
strongest peak in the 64x32x16 raw histogram. D3 and D4 are obtained by down-
sampling the raw 64x32 histogram respectively by two and four, and computing
centre of masses around the strongest peak of subsequent down-sampled
versions.

####  II.2.1 First depth map

The first depth map is calculated directly from the 64x32x16 3D histogram
data. The photon counts data can be assumed to be drawn from the Poisson
distribution $\mathscr{P}\left(.\right)$ as commonly used in Rapp and Goyal
(2017); Halimi et al. (2019). Assuming known background level and a Gaussian
system impulse response, the maximum likelihood estimator of the depth is
obtained as the central mass of the received signal photon time-of-flights
(assuming depths are far from the observation window edges). This estimator is
approximated for each depth pixel $(i,j)$ as

$d_{i,j}=\frac{\sum_{t=\max(1,d^{max}-1)}^{\min(T,d^{max}+1)}t*\max(0,h_{i,j,t}-b_{i,j})}{\sum_{t=\max(1,d^{max}-1)}^{\min(T,d^{max}+1)}\max(0,h_{i,j,t}-b_{i,j})},$
(1)

with $h_{i,j,t}$ being the photon counts acquired in pixel $(i,j)$ for time
bin $t\in[1,T]$, $b_{i,j}$ is the background level of pixel $(i,j)$ estimated
as the median of each pixel, and $d^{max}$ represents the location of the
signal peak estimated as the location of the maximum of the histogram of
counts in photon dense regimes, or using a matched filter in sparse photon
regimes. We integrate over three time bins, between $d^{max}-1$ and
$d^{max}+1$, as this width corresponds approximately to the impulse function
of our system. Before this data is input to the network, we divide the depth
by the total number of time bins to normalize it between 0 and 1. It is then
up-scaled to the size of the desired resolution (four times larger in both
spatial dimensions) using a nearest neighbour interpolation. Nearest neighbour
interpolation is used as it preserves the separate surfaces in the scene. This
is preferable over other interpolation strategies that connect separate
surfaces with new depths that did not exist in the original data.

####  II.2.2 Multi-scale depth features

Multiple resolution scales have been shown to help depth estimation, in
particular in high noise scenarios Gyongy et al. (2020); Rapp and Goyal
(2017); Guo et al. (2019). This information is included to the network by
using four depth features D1, D2, D3 and D4 of different resolution scales.
The dimensions of each feature for our real data (64x32x16 histogram) are
128x64, 64x32, 32x16 and 16x8 respectively. D1 is obtained by down-sampling
the previously obtained 256x128 first depth map by two in both spatial
dimensions using nearest neighbour interpolation, by picking only the top-left
pixel within squares of 2x2 pixels. D2 is obtained by computing center of mass
on the 64x32x16 LR histogram. D3 and D4 are obtained by down-sampling this
histogram by a factor of two and four by summing the neighbouring pixels in
the spatial dimensions, hence obtaining histograms of size 32x16x16 and
16x8x16 respectively. D3 and D4 are then obtained by computing center of mass
on those histograms. Thanks to this process of down-sampling at the level of
histogram, the resultant D3 and D4 have higher signal-to-noise than the first
depth map, albeit with a lower resolution. This helps the network identify
features in images with high levels of noise. All the features are normalized
by dividing with the total number of time bins of the system (i.e., the width
of the range observation window).

####  II.2.3 High resolution intensity

The intensity map of a scene has been used to guide the reconstruction of the
depth in statistical methods Tachella et al. (2019); Gyongy et al. (2020);
Rapp and Goyal (2017); Halimi et al. (2019) and in machine learning methods
Lindell et al. (2018); Guo et al. (2019); Sun et al. (2020). In our case, we
obtain the intensity image directly from the SPAD detector Quantic 4x4. This
intensity image has a spatial resolution of 256x128, which is four times
larger than the 64*32 histogram spatial resolution. We normalize between 0 and
1 this intensity before inputting to the network. The intensity and histogram
data are acquired in alternate frames, so it is possible that objects can move
from frame to frame. However, our system has a high temporal frame rate, so we
assume perfect alignment between the histogram and intensity data.

####  II.2.4 Second depth map

When the Quantic 4x4 sensor operates in the depth mode, each superpixel in the
histogram gathers the photon counts of 4x4 pixels. Therefore, some histograms
might present multiple peaks when observing multiple surfaces located at
different depths. Those histograms contain multiple peaks that correspond to
the different depths involved. While the first depth map is calculated by
identifying the strongest peak, we compute a second depth map based on the
second strongest peak. More precisely, the depth position in the second depth
map is calculated by applying the center of mass of Equation 1 around the
index of the second strongest peak.

We set the following criterion on the minimum number of photon counts a
relevant peak should contain. For each pixel $(i,j)$ we consider a peak at bin
$t$ to be relevant if

$h_{i,j,t}>b_{i,j}+level*\sqrt{b_{i,j}}$ (2)

with $h_{i,j,t}$ being the number photon counts of pixel $(i,j)$ at time bin
$t$; $b_{i,j}$ being the background level at pixel $(i,j)$ estimated by taking
the median value of the histogram; $\sqrt{b_{i,j}}$ represents the standard-
deviation of the Poisson distributed background counts; and $level$ is a
variable adjusted empirically, so that the values in the second depth map do
not come from the noise but mostly represent real signal. In the case of our
captured data, we set $level$ to 12. If the criteria is not met, the second
depth is set to zero.

We note that in high noise scenarios, discriminating peaks corresponding to
real depths from peaks corresponding to background photons is difficult. In
extremely noisy scenarios, no second depth can be extracted and the second
depth map is set to zero. In the same way as for the first depth map, the
second depth map is 4x4 up-scaled with a nearest neighbour interpolation. The
second depth map can be computed in the case of the Quantic 4x4 camera because
this sensor gathers the photon counts of 4x4 pixels in the depth mode. In the
case of more traditional SPAD sensors, e.g. in section 4.IV.3, the second
depth map is set to zero.

####  II.2.5 Outlier rejection step

For the cases when the objects are contained only within a few bins, e.g. in
section 4.IV.3, we crop the LR histogram in the temporal dimension, so that
the returned signal from the studied objects extends over all time bins. In
that case, all the previously detailed depth features (i.e. first depth map,
second depth map, and the four depth features D1, D2, D3 and D4) are computed
on the cropped LR histogram. This step ensures the maximum contrast of the
depth features and the best performance of the network. We evaluate the range
of time bins in which the objects are present in the following way. For one
scene, the histograms of all pixels are added into one single histogram. On
this single histogram, we use criterion (2) (to estimate the second depth
map), to select the bins where a minimum photon count is received, and then
apply a median filter to erase all isolated peaks. The range of the remaining
time bins is used to crop the LR histogram in the temporal dimension.

##  III Proposed method

The proposed HistNet network increases the spatial resolution of the input
depth map by four in both dimensions (e.g., from 64x32 to 128x64 for our real
data). It is also robust to extremely noisy conditions as highlighted in
following sections. The network structure is independent of the input
dimensions and can therefore reconstruct data of any size. However, for
clarity, we use the dimensions of our real data to present the network
structure.

###  III.1 Network architecture

In the context of guided depth super-resolution, Guo et al. Guo et al. (2019)
developed the U-Net-based Ronneberger et al. (2015) DepthSR-Net algorithm,
which offers state-of-the-art performance. The U-Net architecture is
particularly useful as it requires very few training data and is based on CNNs
which have shown good results for image denoising and super-resolution tasks
Lindell et al. (2018); Sun et al. (2020).

In this paper, we develop a new network to account for the multi-scale
information available from the SPAD array histogram data, making it more
robust to sparse photon regimes and/or high background scenarios, in addition
to exploiting the fine details of the intensity guide. The network makes use
of the different features extracted from the raw histogram (first and second
depth maps, the multi-scale depths) and the intensity image. Our network
performs simultaneous up-sampling and denoising of depth data for a wide range
of scenarios. Details of the number and filter sizes for each layer are
reported in Figure 3. Figure 4 shows a schematic representation of the
network.

With respect to other state-of-the art algorithms, HistNet differs in that it
extracts the key information from a histogram prior to the data being
processed by any neural network. This provides the network with known
informative features, allowing a reduction in processing times.

####  III.1.1 Residual U-Net architecture

The goal of the network is to take the data from the low-resolution histogram
and intensity image and produce a residual map $\mathscr{R}$ that can be added
to an up-scaled version of the low resolution depth map He et al. (2015); Guo
et al. (2019). The sum of the residual map and the low-resolution depth map is
the final high-resolution depth map. The goal of the training is to find the
parameters of the filters, i.e. the weights and biases, that minimises the
$l1$-norm between the residual map and the ground truth.

The network consists of an encoder of five layers, denoted L0 to L4, connected
to a five-layer decoder (L5 to L9) with skip connections. The network includes
a branch that processes the multi-scale depth features (see III.1.2) and a
branch that processes the intensity image (see III.1.3). The main input of the
encoder is the concatenation of the first and second depth maps along a third
dimension. In the case of the real data, this input is therefore of size
128x64x2. Note that each filter of the convolutional and deconvolutional layer
has a height and width of 3.

In the encoder, the main input is passed to layer L0, which consists of two
convolution operations of 64 filters each. Layers L1 to L4 consist of three
steps: first, 2x2 max-pooling that down-samples the data by two in both
spatial dimensions; second, integration of the information of the multi-scale
depth features by concatenation with the layer of the depth guidance that has
the same shape (see III.1.2); and finally, two convolution operations with a
number of filters of 128, 256, 512 and 1024 for L1, L2, L3 and L4,
respectively.

In the decoder, the layers L5 to L8 consist of four steps: first,
deconvolution operations filter and up-sample the image by two in both spatial
dimensions; second, skip connections between encoder and decoder are computed
by concatenating the decoder layer with the layer in the encoder that has the
same shape; third, guidance with the intensity is provided by concatenation
with the right layer in the intensity guidance branch (see III.1.3); and
finally, two convolution operations are performed with a number of filters of
512, 256, 128 and 64 for L5 to L8, respectively. L9 is a convolutional layer
with one filter that provides the output. The predicted HR depth map is
obtained by adding the output of the network $\mathscr{R}$, which is known as
the residual map, to the LR depth input, i.e. the first depth map.

####  III.1.2 Guidance with the noise-tolerant multi-scale depth features

The multi-scale depth features are of size 64x32, 32x16, 16x8 and 8x4. Each
feature passes through a convolutional layer of 64, 128, 256 and 512 filters
of size 3x3 respectively, creating cubes of size 64x32x64, 32x16x128, 16x8x256
and 8x4x512. These cubes are integrated in the encoder by concatenation along
the filter dimension with the layer of corresponding size (see III.1.1).

####  III.1.3 Guidance with the intensity map

We use intensity guidance in the same manner as Guo et al. Guo et al. (2019).
The guidance branch consists of convolutional operations followed by 2x2 max-
pooling. The number of filters of the convolution operation of each layer is
64, 128, 256 and 512. For our real data, the outputs of the convolutional
operations of each layer in the guidance branch are of size 128x64x64,
64x32x128, 32x16x256 and 16x8x512. These outputs are integrated along the
decoder part of the network by concatenation along the filter dimension with
the layer of corresponding size (see III.1.1).

Figure 3:  Details of network architecture and associated parameters. The
”output shapes” of the layers specified in the fourth column correspond to the
case of processing our real data (i.e. histograms of spatial resolution 64x32
and 128x64 intensity images). The filters of each layer are described with
four parameters: N stands for the number of filters; w the width; h the height
and c the number of inner channels. cv stands for convolutional layer; mp for
max-pooling layer; cat() for concatenation with the layers specified in the
brackets and dcv for deconvolutional layer. Figure 4:  Representation of the
HistNet. Input of the network consists of the first and the second depth maps.
Multi-resolution depths features are integrated along the contracting path of
the U-Net. The intensity image is processed at multiple resolution and
integrated along the expansive part of the U-Net. Skip connections between the
contracting and the expansive paths are displayed as red arrows.

###  III.2 Loss

We minimize the $l1$ loss to help reconstruct separate surfaces in depth. The
$l1$-norm is known to promote sparsity Tibshirani (1996); Donoho and Logan
(1992). During the training, a batch-mode learning method with a batch size of
$M$ = 64 was used, and the loss was defined as

$\mathscr{L}(\theta)=\dfrac{1}{MN}\sum_{m=1}^{M}\sum_{n=1}^{N}\left|{\mathscr{R}_{m,n}(\theta)+d_{m,n}-d^{ref}_{m,n}}\right|,$
(3)

with $M$ the number of images within one batch, $N$ the number of pixels of
each image, $\theta$ the trainable parameters of the network, the
$\mathscr{R}$ the residual map predicted by HistNet, $d$ the first depth map,
and $d^{ref}$ the ground truth depth.

###  III.3 Simulated datasets for training, validation and testing

We simulate realistic SPAD array measurements (LR histogram and HR intensity)
from 23 scenes of the MPI Sintel Depth Dataset Butler et al. (2012); Wulff et
al. (2012) for the training and validation dataset, and from six scenes of the
Middlebury dataset Scharstein and Pal (2007); Hirschmüller and Scharstein
(2007) for the test dataset.

From the HR depth and intensity images provided by both datasets, we create
histograms of counts of 16 time bins. Here we simulate the case where the
return signal from objects extends over all time bins of the histogram, and
therefore no outlier rejection step is necessary (see II.2.5). The photon
counts of the histogram are assumed to be taken from a Poisson distribution
and the impulse response of our SPAD camera is approximated by a Gaussian
function $\mathscr{G}(m,\sigma)$ with average m and standard deviation
$\sigma=0.5714$ histogram bins Henderson et al. (2019); Gyongy et al. (2020).
For each pixel $(i,j)$, the photon counts $h_{i,j,t}$ of the acquired
histogram at time bin $t$ can be expressed as function of the intensity
$r_{i,j}$ and the depth $d_{i,j}$ as follows

$h_{i,j,t}=\mathscr{P}\left(r_{i,j}*\mathscr{G}(d_{i,j},\sigma)+b_{i,j,t}\right),$
(4)

with $b_{i,j,t}$ the background level which is assumed constant for all time
bins of a given pixel. To reconstruct the real measurements, we estimate the
background $b_{i,j,t}$ of the simulated training dataset as the median of the
histograms of the real measurements. For the simulated measurements, the
background $b_{i,j,t}$ of the simulated training dataset is estimated from the
SBR value. The LR histograms are simulated by down-sampling the HR histogram
by integrating over four by four pixels in the spatial dimensions. We consider
two metrics to assess the level of noise: the signal to background ratio (SBR)
and the number photon counts per pixel ($ppp$). We define the pixel
$SBR_{i,j}$ as

$SBR_{i,j}=\frac{\sum_{t=\max(1,d^{max}-1)}^{\min(T,d^{max}+1)}(h_{i,j,t}-b_{i,j})}{b_{i,j}*(\min(T,d^{max}+1)-\max(1,d^{max}-1))}.$
(5)

The image SBR is then the average value over all pixels, i.e.
$SBR=\frac{1}{N}\sum_{i,j}{SBR_{i,j}}$. In a similar fashion, we define
$ppp_{i,j}$, the number of photon counts reflected from the target in a pixel
(i,j), as

$ppp_{i,j}=\sum_{t=\max(1,d^{max}-1)}^{\min(T,d^{max}+1)}h_{i,j,t}-b_{i,j}$
(6)

The average $ppp$ for the image is then the average value over all the pixels,
i.e. $ppp=\frac{1}{N}\sum_{i,j}{ppp_{i,j}}$.

It should be noted that the 436x1024 MPI Sintel Depth Dataset are pre-
processed to remove all NaN and Inf values using median filtering. We use 21
images for training and two are saved for the validation. We increase the
number of training images by a factor of eight by using all possible
combinations of 90° image rotations and flips. Furthermore, the images are
split into overlapping patches of size 96x96 with a stride of 48.

###  III.4 Implementation details

We implemented HistNet within the Tensorflow framework and use the
ProximalAdagradOptimizer optimizer Singer and Duchi (2009), as this enables
the minimization of the $l1$ loss function. The learning rate was set to 1e-1.
The training was performed on a NVIDIA RTX 6000 GPU. We trained for 2000
epochs, which took about 10 hours.

##  IV Results

###  IV.1 Evaluation on simulated measurements (4-times up-sampling)

####  IV.1.1 Noise scenarios

Three different noise scenarios were considered: a scenario mimicking the
lighting conditions of Gyongy et al. (2020) with ppp = 1200 counts and SBR = 2
denoted ”high signal-to-noise scenario”, a scenario corresponding to a lower
photon count and lower signal to noise with ppp = 4 counts and SBR = 0.02
denoted ”medium signal-to-noise scenario”, and a scenario corresponding to a
lower photon count and much lower signal to noise with ppp = 4 counts and SBR
= 0.006 denoted ”low signal-to-noise scenario” . We trained a separate network
for the three different noise scenarios.

####  IV.1.2 Evaluation metrics

We quantify the performance of the network on simulated data by comparing the
reconstruction with the ground truth depth using the root mean squared error
metric RMSE $=\sqrt{\frac{1}{N}\|{\mathscr{R}+d-d^{ref}}\|^{2}}$ and the
absolute depth error ADE $=\left|{\mathscr{R}+d-d^{ref}}\right|$,
$\mathscr{R}$ being the residual map predicted by HistNet after the training,
$d$ the up-scaled version of the low resolution depth map, i.e the first depth
map, and $d^{ref}$ the ground truth.

####  IV.1.3 Comparison algorithms

We compare the results of HistNet with the following methods:

  * •

Nearest-neighbour interpolation Depth is estimated with center of mass on the
276x344x16 histogram of counts and up-sampled with nearest neighbour
interpolation to a 1104x1376 image. Note that this is the method used to
produce the first depth argument of HistNet. We choose nearest-neighbour
interpolation to avoid joining spatially separated surfaces.

  * •

Guided Image Filtering of He et al. 2013 He et al. (2013) We perform further
processing to the estimated depth from the nearest-neighbour interpolation by
applying the guided filtering algorithm with the HR intensity image as a
guide.

  * •

DepthSR-Net of Guo et al. 2019 Guo et al. (2019) We retrained this network
using the same training datasets for our network. This network outputs a 4x
upsampled depth map from a LR depth map using a HR intensity map to guide the
reconstruction. In Guo et al. (2019), the LR depth map is first up-sampled to
the desired size with a bicubic interpolation. However, we want to reconstruct
surfaces that are well separated from one another. Therefore, nearest
neighbour interpolation instead of bicubic interpolation is used to upsample
the input LR depth map.

  * •

Algorithm of Gyongy et al. 2020 Gyongy et al. (2020) This algorithm is
designed to process the Quantic 4 x4 SPAD array. It consists of various steps
of guided filtering and up-sampling with low computational cost. One part of
the algorithm is designed to compensate the inherent misalignement between the
depth and the intensity information that the SPAD provides. Since our
synthetic data consists of perfectly aligned intensity and depth, we do not
use this part of the algorithm.

####  IV.1.4 Results for the different signal-to-noise scenarios

Figure 5 shows the reconstruction for the different reconstruction methods for
high signal-to-noise scenario: (ppp = 1200 counts and SBR = 2). We see that
HistNet is able to produce sharp and clean boundaries. Guided Filtering and
DepthSR-Net introduce blurred details around the edges and nearest
interpolation leads to pixelated images. We report the root mean squared error
and the absolute error in Table 1 for the two Middlebury scenes reconstructed
with the different methods. This table indicates that HistNet outperforms the
other methods in both performance categories. The processing time of the
different methods is reported in Table 1. Guided Filtering and nearest
interpolation have a very low computational cost and process the image the
fastest, in a few milliseconds for a 1104x1376 input. The algorithm of Gyongy
et al. Gyongy et al. (2020) reconstructs the image in about 4 seconds. The
reconstructions of HistNet and DepthSR-Net were performed on a NVIDIA RTX 6000
GPU. Each 1104x1376 Middlebury scene took about 7 seconds to reconstruct with
the two networks.

Figure 6 and 7 show the results for measurements simulated with an medium
signal-to-noise (ppp=4 counts, SBR =0.02) and low signal-to-noise (ppp=4
counts, SBR =0.006), respectively. Visually, we see that our method performs
the best in high noise scenarios. Quantitative comparison can be found in
Table 1. For both scenes, HistNet performs better in terms of the RMSE and
ADE.

| Proposed HistNet |  Gyongy et al. Gyongy et al. (2020) |  Guo et al. Guo et al. (2019) | NNI | Guided Image Filtering
---|---|---|---|---|---
Rec time per scene | 7s (on GPU) | 4s | 7s (on GPU) | 1ms | 0.4s
Training on high signal-to-noise data with second depth; ppp=1200 counts and
SBR=2
Scene |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE
Art |  0.023 |  0.0027 |  0.043 |  0.0076 |  0.026 |  0.0080 |  0.053 |  0.038 |  0.046 |  0.039
Reindeer |  0.012 |  0.0018 |  0.023 |  0.0040 |  0.015 |  0.0051 |  0.040 |  0.035 |  0.037 |  0.035
Training on medium signal-to-noise data without second depth; ppp=4 counts and
SBR=0.02
Scene |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE
Art |  0.053 |  0.019 |  0.11 |  0.050 |  0.054 |  0.023 |  0.32 |  0.22 |  0.22 |  0.17
Reindeer |  0.040 |  0.019 |  0.12 |  0.060 |  0.047 |  0.024 |  0.31 |  0.21 |  0.21 |  0.16
Training on low signal-to-noise data without second depth; ppp=4 counts and
SBR=0.006
Scene |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE
Art |  0.082 |  0.055 |  0.248 |  0.187 |  0.102 |  0.064 |  0.363 |  0.276 |  0.27 |  0.22
Reindeer |  0.075 |  0.050 |  0.234 |  0.168 |  0.083 |  0.053 |  0.357 |  0.272 |  0.259 |  0.206
Table 1:  Quantitative comparison of the different reconstruction methods for
4x up-sampling on simulated measurements with a high, medium and low signal-
to-noise. RMSE is the root-mean-square error; ADE is the absolute depth error.
Figure 5:  High signal-to-noise. Comparison of reconstruction techniques for
measurements simulated with an average of 1200 signal photons per pixel and a
signal-to-background ratio of 2. ADE is the absolute depth error calculated
with normalized data between 0 and 1. (a) displays the ground truth depth
image; (b) displays the ground truth of the intensity; (c) is the
reconstruction of the depth with nearest interpolation; (d) is the up-sampled
version of the first depth data with Guided Image Filtering He et al. (2013);
(e) is the reconstruction of the first depth data with the algorithm DepthSR-
Net of Guo et al. 2019 Guo et al. (2019); (f) is the reconstruction of the
first depth data with the algorithm of Gyongy et al. 2020 Gyongy et al.
(2020); (g) is the reconstruction via our proposed method HistNet. (a*), (b*),
(c*), (d*), (e*), (f*) and (g*) are the closeup views of (a)-(g). (h)-(n) are
the corresponding images from another simulated measurements with (h*)-(n*)
the closeup views.  Figure 6:  Medium signal-to-noise. Comparison of
reconstruction techniques for measurements simulated with an average of 4
signal photons per pixel and a signal-to-background ratio of 0.02. ADE is the
absolute depth error calculated with normalized data between 0 and 1. (a)
displays the ground truth depth image; (b) displays the ground truth of the
intensity; (c) is the reconstruction of the depth with nearest interpolation ;
(d) is the up-sampled version of the first depth data with Guided Image
Filtering He et al. (2013); (e) is the reconstruction of the first depth data
with the algorithm DepthSR-Net of Guo et al. 2019 Guo et al. (2019); (f) is
the reconstruction of the first depth data with the algorithm of Gyongy et al.
2020 Gyongy et al. (2020); (g) is the reconstruction via our proposed method
HistNet. (a*), (b*), (c*), (d*), (e*), (f*) and (g*) are the closeup views of
(a)-(g). (h)-(n) are the corresponding images from another simulated
measurements with (h*)-(n*) the closeup views. Figure 7:  Low signal-to-noise.
Comparison of reconstruction techniques for measurements simulated with an
average of 4 signal photons per pixel and a signal-to-background ratio of
0.006. ADE is the absolute depth error calculated with normalized data between
0 and 1. (a) displays the ground truth depth image; (b) displays the ground
truth of the intensity; (c) is the reconstruction of the depth with nearest
interpolation; (d) is the up-sampled version of the first depth data with
Guided Image Filtering He et al. (2013); (e) is the reconstruction of the
first depth data with the algorithm DepthSR-Net of Guo et al. 2019 Guo et al.
(2019); (f) is the reconstruction of the first depth data with the algorithm
of Gyongy et al. 2020 Gyongy et al. (2020); (g) is the reconstruction via our
proposed method HistNet. (a*), (b*), (c*), (d*), (e*), (f*) and (g*) are the
closeup views of (a)-(g). (h)-(n) are the corresponding images from another
simulated measurements with (h*)-(n*) the closeup views.

####  IV.1.5 Robustness to noise

We study how well a network trained on data with a specific SBR and ppp levels
can reconstruct data with other noise levels. HistNet trained according to the
high and medium signal-to-noise scenarios was tested on data with ppp levels
ranging from $1*10^{-1}$ to $7*10^{5}$ and SBR ranging from $1*10^{-5}$ to
$70$. Figure 8 shows the RMSE value between the ground truth and the
reconstruction of HistNet with respect to the SBR and the $ppp$ of the testing
data. We see that HistNet shows good performance on a variety of SBR and $ppp$
levels. HistNet trained with the medium signal-to-noise scenario is able to
reconstruct data that presents higher noise than when trained with the high
signal-to-noise scenario. However, the best performance is always achieved
when the $ppp$ and SBR of the testing data approximately matches the one
implemented in the training dataset. Figure 9 shows the reconstruction of data
of low signal-to-noise scenario by HistNet trained on data with median signal-
to-noise scenario.

Figure 8:  RMSE between predictions of network and ground truth for different
SBR and ppp. The noise level implemented in the training dataset is
represented with a red marker. In (a), medium signal-to-noise, the network was
trained on data with an average signal count ($ppp$) of 4 per pixel and a
signal to background ratio (SBR) of 0.02. In (b), high-signal-to-noise, the
network was trained on data with an average signal count of 1200 and a signal
to background ratio (SBR) of 2. Both networks were tested on simulated data
with $ppp$ levels ranging from $1*10^{-1}$ to $7*10^{5}$ and SBR ranging from
$1*10^{5}$ to $70$. The network performs best when the training and testing
noise matches. It is also robust to data that have a lower noise. However, the
performance drops when the testing data presents higher noise than the
training data.  Figure 9: Reconstruction of data with low signal-to-noise
(ppp=4 counts and SBR=0.006) with HistNet and DepthSR-Net of Guo et al. Guo et
al. (2019) trained on data with medium signal-to-noise (ppp=4 and SBR=0.02).
(a) is the ground truth depth; (b) shows the first depth map input of the
network ; (c) is the 4x-upsampled reconstruction with DepthSR-Net of Guo et
al. ; (d) s the 4x-upsampled reconstruction with HistNet. RMSE is the root-
mean-square error.

####  IV.1.6 Ablation study

In this section, we investigate the relative benefits of the different input
features of the network, i.e. the intensity image, the multi-scale depth
features, and the second depth map. For each of these different feature
inputs, we used three versions of HistNet: one that does not use the intensity
(w/o Intensity); one that does not use the multi-scale depth features (w/0
Depth Features); and one that does not use the second depth maps (w/0 Second
Depth). We trained each of these versions of HistNet for 4x up-sampling for
the medium and high noise scenarios.

Quantitative results of different simulated measurements from the Middlebury
dataset are displayed in Tables 2 and 3 for the high and medium signal-to-
noise scenarios, respectively. In the high signal-to-noise scenario, the
absence of the second depth map results in a noticeable reduction in the
quality of the reconstruction. In the medium signal-to-noise scenario, the
absence of the multi-scale depth features and the intensity image
significantly reduces the performance of the network. These results
demonstrate the value that each of the features brings to the network.

| HistNet | HistNet w/o Intensity | HistNet w/o Depth Features | HistNet w/o Second Depth
---|---|---|---|---
Scene |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE
Art |  0.0230 |  0.0027 |  0.0210 |  0.0027 |  0.0230 |  0.0029 |  0.0260 |  0.0041
Books |  0.0057 |  0.0014 |  0.0061 |  0.0014 |  0.0060 |  0.0014 |  0.0082 |  0.0017
Dolls |  0.0050 |  0.0015 |  0.0055 |  0.0017 |  0.0051 |  0.0016 |  0.0061 |  0.0020
Laundry |  0.0128 |  0.0033 |  0.0124 |  0.0034 |  0.0128 |  0.0034 |  0.0150 |  0.0043
Moebius |  0.0129 |  0.005 |  0.0153 |  0.0064 |  0.0140 |  0.0056 |  0.0150 |  0.0063
Reindeer |  0.0126 |  0.0019 |  0.0125 |  0.0019 |  0.0123 |  0.0019 |  0.0140 |  0.0023
AVG |  0.012 $\pm$ 0.0065 |  0.026 $\pm$ 0.014 |  0.012 $\pm$0.0058 |  0.0029 $\pm$ 0.0019 |  0.012 $\pm$ 0.0065 |  0.0028 $\pm$ 0.0015 |  0.014 $\pm$ 0.0069 |  0.0095 $\pm$ 0.0153
Table 2:  High signal-to-noise. Quantitative comparison on different simulated measurements with ppp = 1200 counts and SBR = 2 between versions of HistNet with and without using the intensity maps, the down-sampled features or the second depth map. The version of HistNet that integrates the intensity, the down-sampled features and the second depth map has the best performance in terms of RMSE and ADE for most of the simulated measurements. RMSE is the root-mean-square error; ADE is the absolute depth error.  | HistNet | HistNet w/o Intensity | HistNet w/o Depth Features
---|---|---|---
Scene | RMSE | ADE | RMSE | ADE | RMSE | ADE
Art | 0.053 | 0.019 | 0.057 | 0.020 | 0.064 | 0.023
Books | 0.017 | 0.009 | 0.021 | 0.010 | 0.019 | 0.010
Dolls | 0.019 | 0.012 | 0.029 | 0.013 | 0.061 | 0.018
Laundry | 0.024 | 0.013 | 0.026 | 0.012 | 0.026 | 0.013
Moebius | 0.025 | 0.015 | 0.024 | 0.013 | 0.045 | 0.020
Reindeer | 0.040 | 0.019 | 0.050 | 0.021 | 0.055 | 0.023
AVG | 0.030 $\pm$ 0.013 | 0.015 $\pm$ 0.004 |  0.034 $\pm$ 0.014 | 0.015 $\pm$ 0.004 |  0.045 $\pm$ 0.017 |  0.018 $\pm$ 0.005
Table 3:  Medium signal-to-noise. Quantitative comparison on different
simulated measurements with average ppp = 4 counts and SBR = 0.02 between
versions of HistNet with and without using the intensity maps, the down-
sampled features or the second depth maps. The version of HistNet that
integrates both the intensity and the down-sampled features performs better
than the one without it in terms of RMSE and ADE for most of the simulated
measurements. RMSE is the root-mean-square error; ADE is the absolute depth
error.

###  IV.2 Results for the Quantic 4x4 camera (4-times up-sampling)

####  IV.2.1 Quantic 4x4 camera

We test the performance of HistNet on real measurements captured by the
Quantic 4x4 camera Gyongy et al. (2020). The spatial resolution of the
histogram data is of 32x64 and the number of time bins is of 16. The
resolution of the intensity image is of 128x256. The data is first
interpolated to the size of the intensity image with nearest interpolation and
is calibrated using a compensation frame. We estimated a number of photon
counts per pixel of 1200 and a signal-to-background ratio of 2 in this data.
Therefore, we use HistNet trained with this scenario to reconstruct the depth
maps. Figure 10 show the reconstruction of Quantic 4x4 data via HistNet
together with a comparison with different reconstruction algorithms. We see
that HistNet leads to more accurate image with sharper edges.

Figure 10:  Reconstruction of Quantic 4x4 data. (a) displays the reflectivity
image from the SPAD Gyongy et al. (2020); (b) is the reconstruction of the
depth using nearest interpolation ; (c) is the reconstruction of the depth
using guided image filtering et al. 2013 He et al. (2013); (d) is the
reconstruction via the algorithm DepthSR-Net of Chunle Guo et al. 2019 Guo et
al. (2019); (e) displays the reconstruction of the algorithm of Gyongy et al.
2020 Gyongy et al. (2020) ; (f) is the reconstruction of our proposed HistNet.
The second row (g)-(l) displays the equivalent for another captured
measurement and (g*)-(l*) displays closeup of (g)-(l).

###  IV.3 Generalisation to a new dataset (8-times up-sampling)

In this section, we consider the case of 8-fold up-sampling of depth from data
acquired by another single-photon detector array. Here we apply our algorithm
for up-sampling the dataset collected and presented in Lindell et al. Lindell
et al. (2018). This shows that our network can be applied to other data
formats, not just data from the Quantic 4x4 camera.

####  IV.3.1 Description of the data

The SPAD array sensor in ref Lindell et al. (2018) acquires histograms of
photon counts with a spatial resolution of 72x88 with 1024 time bins. A
conventional high-resolution camera acquires a corresponding HR intensity
image at a resolution of 576x704 pixels (8x higher spatial resolution than the
original data). The histograms are very sparse, i.e. the return signal from
objects is only contained within a range of 60 time bins, and the remaining
924 time bins contain only background photons.

####  IV.3.2 Pre-processing of data for the network

The LR histograms (of size 72x88x1024) are cropped in the temporal dimension
to include only the return signal from objects as described in II.2.5. We note
though that the rejection step removes data in temporal bins where there is no
signal, so we do not anticipate that it has any impact on the signal-to-
background ratio (5) where any signal is present.

The first depth map and the multi-scale depth features are computed from this
cropped version of the LR histogram. The first depth map is computed using a
matched filter on this cropped histogram of photon counts, which is then up-
scaled with the nearest neighbour interpolation to the desired 576x704
resolution (8x larger in both spatial dimensions in this section). The first
depth map is divided by the number of bins of the cropped histogram to be
normalize between zero and one. Unlike the Quantic 4x4, the sensor considered
in this section does not integrate the histograms of 4x4 pixels. Therefore,
the second depth map is set to zero.

The multi-scale depth features D1 and D2 are obtained by down-sampling by two
and four the previously obtained first depth map. D3 is obtained by using a
matched filter on the cropped histogram of photon counts. D4 is obtained by
down-sampling the cropped histogram of photon counts by two by summing the
neighbouring pixels in the spatial dimension and by using matched filter on
the down-sampled histogram. All D1-4 are divided by the number of bins of the
cropped histogram. The dimensions of the features are 288x352, 144x276, 72x88
and 36x44 for D1, D2, D3 and D4 respectively. The HR intensity acquired by the
camera is normalized between zero and one. Note that the algorithm of Lindell
et al. does not require any pre-processing as it operates directly on
histogram data.

####  IV.3.3 Training of HistNet for 8-times up-sampling

We simulate the SPAD measurements (LR Histogram and HR intensity) from 23
depth images extracted from the MPI dataset. The HR depth and intensity images
are first down-sampled by a factor of eight with bicubic interpolation.
Histograms with 60 time bins are computed from these LR depth and LR intensity
images (as detailed in III.3), so that the depth is spread over the entire
range of bins. The noise level of our training dataset is chosen to match the
noise level of the Lindell dataset as close as possible. We use a ppp of 2
counts and a SBR of 40.

####  IV.3.4 Comparison algorithms

  * •

Nearest-neighbour interpolation Depth is estimated with matched filter on the
cropped 276x344x60 histogram of counts and up-sampled with nearest neighbour
interpolation to a 576x704 image.

  * •

Guided Image Filtering of He et al. 2013 He et al. (2013) We perform further
processing to the estimated depth from the nearest-neighbour interpolation by
applying the guided filtering algorithm with the HR intensity image as a
guide.

  * •

DepthSR-Net of Guo et al. 2019 Guo et al. (2019) We trained this network using
the same training datasets for our network. This network outputs an 8x up-
sampled depth map from a LR depth map using a HR intensity map to guide the
reconstruction.

  * •

Network of Lindell et al. 2018 Lindell et al. (2018) We use the end-to-end
trained network for depth estimation and guided 8x up-sampling presented in
Lindell et al. (2018). This network takes a LR histogram of photon counts and
a HR intensity image as the input, and it outputs an eight-fold up-sampled
depth. The network consists of two parts: first, a denoising branch that
estimates a clear depth image from the histogram; and second, a guided up-
sampling branch that up-samples the depth image to the desired resolution.
Both branches consist of a series of 3D convolutions and make use of the
intensity image as guidance. The denoising branch processes the histogram of
photon counts at multiple scales.

####  IV.3.5 Results

Figure 11 shows the reconstruction for the different reconstruction methods,
and we report quantitative results for eight Middlebury scenes in Table 4.
This table shows that HistNet compares well with the state-of-the-art
algorithms DepthSR-Net Guo et al. (2019) and that of Lindell et al. Lindell et
al. (2018). The RMSE and ADE metrics are all within one standard deviation of
each other, thus all algorithms are extremely close in terms of image quality
in this scenario. This similarity in performance is to be expected due to the
high signal-to-noise level of this dataset. In terms of processing time,
HistNet and the network of Guo et al. Guo et al. (2019) take 3.6s to
reconstruct the image, the network of Lindell et al. Lindell et al. (2018)
takes 11.7s and guided filtering and nearest interpolation processes the image
in a few milliseconds. The networks were tested using a NVIDIA RTX 6000 GPU.
We note that our proposed HistNet and DepthSR-Net of Guo et al. Guo et al.
(2019) act on pre-processed histograms by extracting informative depth
features. Hence, the processing time of those learning-based methods depends
on the spatial desired resolution of the depth maps and does not depend on the
number of time bins in the histogram.

| Proposed HistNet |  Lindell et al. Lindell et al. (2018) |  Guo et al. Guo et al. (2019) | NNI | Guided Image Filtering
---|---|---|---|---|---
Rec time per scene | 3.6s (on GPU) | 11.7s (on GPU) | 3.6s (on GPU) | 0.004s | 0.05s
Scene |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE |  RMSE |  ADE
Art |  0.092 |  0.043 |  0.111 |  0.060 |  0.085 |  0.040 |  0.313 |  0.182 |  0.264 |  0.169
Books |  0.042 |  0.024 |  0.039 |  0.022 |  0.042 |  0.023 |  0.315 |  0.177 |  0.265 |  0.165
Dolls |  0.046 |  0.024 |  0.037 |  0.024 |  0.039 |  0.023 |  0.280 |  0.155 |  0.227 |  0.138
Laundry |  0.056 |  0.024 |  0.052 |  0.027 |  0.062 |  0.028 |  0.232 |  0.128 |  0.185 |  0.112
Moebius |  0.042 |  0.024 |  0.043 |  0.027 |  0.039 |  0.024 |  0.273 |  0.152 |  0.221 |  0.139
Reindeer |  0.110 |  0.067 |  0.085 |  0.043 |  0.117 |  0.067 |  0.383 |  0.230 |  0.334 |  0.213
Bowling |  0.065 |  0.028 |  0.070 |  0.030 |  0.074 |  0.035 |  0.274 |  0.157 |  0.222 |  0.139
Plastic |  0.046 |  0.023 |  0.037 |  0.019 |  0.052 |  0.026 |  0.252 |  0.140 |  0.197 |  0.121
AVG |  0.062 $\pm$ 0.023 |  0.032 $\pm$ 0.015 |  0.059 $\pm$ 0.027 |  0.031 $\pm$ 0.013 |  0.064 $\pm$ 0.026 |  0.033 $\pm$ 0.014 |  0.290 $\pm$ 0.047 |  0.165 $\pm$ 0.032 |  0.240 $\pm$ 0.047 |  0.150 $\pm$ 0.032
Table 4:  Quantitative comparison of the different reconstruction methods for
8x up-sampling on simulated measurements with a ppp of 2 counts and an SBR of
40. RMSE is the root-mean-square error; ADE is the absolute depth error. The
error metrics were calculated on normalized data such that range of the ground
truth depth images was between zero and one.  Figure 11:  Comparaison of
reconstruction techniques for 8-fold up-sampling ADE is the absolute depth
error calculated with normalized data between 0 and 1. (a) displays the ground
truth depth image; (b) displays the ground truth of the intensity; (c) is the
reconstruction of the depth with nearest interpolation ; (d) is the up-sampled
version of the first depth data with Guided Image Filtering He et al. (2013);
(e) is the reconstruction estimated by the CNN network of Lindell et al. 2018
Lindell et al. (2018); (f) is the reconstruction of the first depth data with
the algorithm DepthSR-Net of Guo et al. 2019 Guo et al. (2019); (g) is the
reconstruction via our proposed method HistNet. (a*), (b*), (c*), (d*), (e*),
(f*) and (g*) are the closeup views of (a)-(g). (h)-(n) are the corresponding
images from another simulated measurements with (h*)-(n*) the closeup views.

##  V Discussion and conclusion

In this work we present a deep network for up-scaling and denoising depth
images. The network is designed for measurements provided by SPAD array
sensors and provides significant improvements for image quality and resolution
over a wide range of noise scenarios. Our method exploits the SPAD array data
in a simple and efficient manner, i.e. we extract multi-scale features and
multiple depths from the histogram data, and these are provided directly to
the network. Additionally we made use of the $l1$ loss to promote sparsity in
the reconstructed images.

The combination of the additional features and the loss function result in a
network with state-of-the-art performance in terms of image metrics and
processing times. The network performs well with respect to other up-sampling
algorithms, especially when applied to data with low signal-to-noise ratios
and low photon levels. Moreover, our method shows robustness to a wide range
of different noise scenarios, so that the noise statistics of the training
dataset do not need to closely match the input data’s. Future work will focus
on the high frame rate of the SPAD array sensor and use information in the
temporal domain to achieve better spatial resolutions for depth images. We
also propose to tackle the misalignment between the histogram and the
intensity image, which is inherent to the operating mode of our SPAD detector
that acquires them alternately.

## Funding

Engineering and Physical Science Research Council (EP/T00097X/1, EP/S001638/1
and EP/L016753/1); UK Royal Academy of Engineering through the Research
Fellowship Scheme under Grant RF/201718/17128; DSTL Dasa project
DSTLX1000147844

## Acknowledgements

The authors are grateful to STMicroelectronics and the ENIAC-POLIS project for
chip fabrication for the Quantic 4x4 sensor. We are grateful to the authors of
the DepthSRNet paper for sharing their code Guo et al. (2019). The code for
this work can be found at https://github.com/HWQuantum/HistNet. We thank
Frédéric Ruget for helpful discussions relating to the simulation of the
datasets.

## Disclosures

The authors declare no conflicts of interest.

## References

  * Ximenes et al. (2018) A. R. Ximenes, P. Padmanabhan, M. Lee, Y. Yamashita, D. N. Yaung, and E. Charbon, in _2018 IEEE International Solid - State Circuits Conference - (ISSCC)_ (2018), pp. 96–98.
  * Zhang et al. (2018) C. Zhang, S. Lindner, I. M. Antolovic, J. Pavia, M. Wolf, and E. Charbon, IEEE Journal of Solid-State Circuits PP, 1 (2018).
  * Maccarone et al. (2015) A. Maccarone, A. McCarthy, X. Ren, R. E. Warburton, A. M. Wallace, J. Moffat, Y. Petillot, and G. S. Buller, Opt. Express 23, 33911 (2015).
  * Halimi et al. (2017) A. Halimi, A. Maccarone, A. McCarthy, S. McLaughlin, and G. S. Buller, IEEE Transactions on Computational Imaging 3, 472 (2017).
  * Chan et al. (2019) S. Chan, A. Halimi, F. Zhu, I. Gyongy, R. K. Henderson, R. Bowman, S. McLaughlin, G. S. Buller, and J. Leach, Scientific Reports 9, 8075 (2019), URL https://doi.org/10.1038/s41598-019-44316-x.
  * McCarthy et al. (2013) A. McCarthy, X. Ren, A. D. Frera, N. R. Gemmell, N. J. Krichel, C. Scarcella, A. Ruggeri, A. Tosi, and G. S. Buller, Opt. Express 21, 22098 (2013), URL http://www.opticsexpress.org/abstract.cfm?URI=oe-21-19-22098.
  * Pawlikowska et al. (2017) A. M. Pawlikowska, A. Halimi, R. A. Lamb, and G. S. Buller, Opt. Express 25, 11919 (2017), URL http://www.opticsexpress.org/abstract.cfm?URI=oe-25-10-11919.
  * Hutchings et al. (2019) S. W. Hutchings, N. Johnston, I. Gyongy, T. Al Abbas, N. A. W. Dutton, M. Tyler, S. Chan, J. Leach, and R. K. Henderson, IEEE Journal of Solid-State Circuits 54, 2947 (2019).
  * Henderson et al. (2019) R. K. Henderson, N. Johnston, S. W. Hutchings, I. Gyongy, T. A. Abbas, N. Dutton, M. Tyler, S. Chan, and J. Leach, in _2019 IEEE International Solid- State Circuits Conference - (ISSCC)_ (2019), pp. 106–108.
  * Gyongy et al. (2020) I. Gyongy, S. W. Hutchings, A. Halimi, M. Tyler, S. Chan, F. Zhu, S. McLaughlin, R. K. Henderson, and J. Leach, Optica (2020).
  * Bronzi et al. (2014) D. Bronzi, F. Villa, S. Tisa, A. Tosi, F. Zappa, D. Durini, S. Weyers, and W. Brockherde, IEEE Journal of Selected Topics in Quantum Electronics 20, 354 (2014).
  * Ren et al. (2018) X. Ren, P. W. R. Connolly, A. Halimi, Y. Altmann, S. McLaughlin, I. Gyongy, R. K. Henderson, and G. S. Buller, Opt. Express 26, 5541 (2018), URL http://www.opticsexpress.org/abstract.cfm?URI=oe-26-5-5541.
  * Morimoto et al. (2020) K. Morimoto, A. Ardelean, M.-L. Wu, A. C. Ulku, I. M. Antolovic, C. Bruschini, and E. Charbon, Optica 7, 346 (2020), URL http://www.osapublishing.org/optica/abstract.cfm?URI=optica-7-4-346.
  * Gariepy et al. (2015) G. Gariepy, N. Krstajić, R. Henderson, C. Li, R. R. Thomson, G. S. Buller, B. Heshmat, R. Raskar, J. Leach, and D. Faccio, Nature Communications 6, 6021 (2015), URL https://doi.org/10.1038/ncomms7021.
  * Gariepy et al. (2016) G. Gariepy, F. Tonolini, R. Henderson, J. Leach, and D. Faccio, Nature Photonics 10, 23 (2016), URL https://doi.org/10.1038/nphoton.2015.234.
  * Bruschini et al. (2019) C. Bruschini, H. Homulle, I. M. Antolovic, S. Burri, and E. Charbon, Light: Science & Applications 8, 87 (2019), URL https://doi.org/10.1038/s41377-019-0191-5.
  * Tachella et al. (2019) J. Tachella, Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, S. McLaughlin, and J.-Y. Tourneret, SIAM Journal on Imaging Sciences 12, 521–550 (2019), ISSN 1936-4954.
  * Rapp and Goyal (2017) J. Rapp and V. K. Goyal, IEEE Transactions on Computational Imaging 3, 445 (2017).
  * Halimi et al. (2019) A. Halimi, R. Tobin, A. McCarthy, J. Bioucas-Dias, S. McLaughlin, and G. Buller, IEEE Transactions on Computational Imaging 6, 138 (2019), ISSN 2333-9403.
  * Guo et al. (2019) C. Guo, C. Li, H. F. Jichang Guo, and P. Han, IEEE Transactions on Image Processing 28, 2545 (2019).
  * Sun et al. (2020) Z. Sun, D. B. Lindell, O. Solgaard, and G. Wetzstein, Opt. Express 28, 14948 (2020), URL http://www.opticsexpress.org/abstract.cfm?URI=oe-28-10-14948.
  * Lindell et al. (2018) D. B. Lindell, M. O’Toole, and G. Wetzstein, ACM Trans. Graph. (2018).
  * Callenberg et al. (2021) C. Callenberg, A. Lyons, D. d. Brok, A. Fatima, A. Turpin, V. Zickus, L. Machesky, J. Whitelaw, D. Faccio, and M. B. Hullin, Scientific Reports 11, 1689 (2021), URL https://doi.org/10.1038/s41598-021-81159-x.
  * Wallace et al. (2020) A. Wallace, A. Halimi, and G. Buller, IEEE Transactions on Vehicular Technology PP, 1 (2020).
  * Rapp et al. (2020) J. Rapp, J. Tachella, Y. Altmann, S. McLaughlin, and V. K. Goyal, IEEE Signal Processing Magazine 37, 62 (2020).
  * Peng et al. (2020) J. Peng, Z. Xiong, X. Huang, Z.-P. Li, D. Liu, and F. Xu, in _Computer Vision – ECCV 2020_ , edited by A. Vedaldi, H. Bischof, T. Brox, and J.-M. Frahm (Springer International Publishing, Cham, 2020), pp. 225–241.
  * Nishimura et al. (2020) M. Nishimura, D. B. Lindell, C. Metzler, and G. Wetzstein, in _Computer Vision – ECCV 2020_ , edited by A. Vedaldi, H. Bischof, T. Brox, and J.-M. Frahm (Springer International Publishing, Cham, 2020), pp. 139–155.
  * Ronneberger et al. (2015) O. Ronneberger, P. Fischer, and T. Brox, CoRR abs/1505.04597 (2015), eprint 1505.04597, URL http://arxiv.org/abs/1505.04597.
  * He et al. (2015) K. He, X. Zhang, S. Ren, and J. Sun, _Deep residual learning for image recognition_ (2015), eprint 1512.03385.
  * Tibshirani (1996) R. Tibshirani, Journal of the Royal Statistical Society: Series B (Methodological) 58, 267 (1996), URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.2517-6161.1996.tb02080.x.
  * Donoho and Logan (1992) D. L. Donoho and B. F. Logan, SIAM Journal on Applied Mathematics 52, 577 (1992), URL http://www.jstor.org/stable/2102426.
  * Butler et al. (2012) D. J. Butler, J. Wulff, G. B. Stanley, and M. J. Black, in _European Conf. on Computer Vision (ECCV)_ , edited by A. Fitzgibbon et al. (Eds.) (Springer-Verlag, 2012), Part IV, LNCS 7577, pp. 611–625.
  * Wulff et al. (2012) J. Wulff, D. J. Butler, G. B. Stanley, and M. J. Black, in _ECCV Workshop on Unsolved Problems in Optical Flow and Stereo Estimation_ , edited by A. Fusiello et al. (Eds.) (Springer-Verlag, 2012), Part II, LNCS 7584, pp. 168–177.
  * Scharstein and Pal (2007) D. Scharstein and C. Pal, IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2007) (2007).
  * Hirschmüller and Scharstein (2007) H. Hirschmüller and D. Scharstein, IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2007) (2007).
  * Singer and Duchi (2009) Y. Singer and J. C. Duchi, in _Advances in Neural Information Processing Systems 22_ , edited by Y. Bengio, D. Schuurmans, J. D. Lafferty, C. K. I. Williams, and A. Culotta (Curran Associates, Inc., 2009), pp. 495–503, URL http://papers.nips.cc/paper/3793-efficient-learning-using-forward-backward-splitting.pdf.
  * He et al. (2013) K. He, J. Sun, and X. Tang, IEEE Transactions on Pattern Analysis and Machine Intelligence 35, 1397 (2013).