Bayesian 3D Reconstruction of Subsampled Multispectral Single-photon Lidar Signals

TL;DR

This paper introduces a Bayesian 3D reconstruction method for multispectral single-photon Lidar that efficiently reconstructs multiple surfaces with fewer photons and measurements, reducing acquisition time and computational load.

Contribution

The novel Bayesian algorithm jointly processes spectral bands and models spatial correlations, enabling accurate multi-surface reconstruction from subsampled data with less photon detection.

Findings

Better reconstruction accuracy than existing methods.

Effective subsampling schemes reduce measurements without performance loss.

Demonstrated on synthetic and real data.

Abstract

Light detection and ranging (Lidar) single-photon devices capture range and intensity information from a 3D scene. This modality enables long range 3D reconstruction with high range precision and low laser power. A multispectral single-photon Lidar system provides additional spectral diversity, allowing the discrimination of different materials. However, the main drawback of such systems can be the long acquisition time needed to collect enough photons in each spectral band. In this work, we tackle this problem in two ways: first, we propose a Bayesian 3D reconstruction algorithm that is able to find multiple surfaces per pixel, using few photons, i.e., shorter acquisitions. In contrast to previous algorithms, the novel method processes the jointly all the spectral bands, obtaining better reconstructions using less photon detections. The proposed model promotes spatial correlation between neighbouring points within a given surface using spatial point processes. Secondly, we account for different spatial and spectral subsampling schemes, which reduce the total number of measurements, without significant degradation of the reconstruction performance. In this way, the total acquisition time, memory requirements and computational time can be significantly reduced. The experiments performed using both synthetic and real single-photon Lidar data demonstrate the advantages of tailored sampling schemes over random alternatives. Furthermore, the proposed algorithm yields better estimates than other existing methods for multi-surface reconstruction using multispectral Lidar data.

Tables

Table 1. TABLE I: Hyperparameters values

Hyperparameter	Coarse scale	Fine scale
${\gamma }_a$	$e^2$	$e^3$
${\lambda }_a$	${\left(N_r{N}_r/N_{\text{bin}}^2\right)}^{1.5}$	${\left(N_r{N}_r\right)}^{1.5}$
$N_b$	$8N_{\text{bin}}{\Delta }_p/{\Delta }_{\text{bin}}$	$8{\Delta }_p/{\Delta }_{\text{bin}}$
$d_{\min }$	$2N_b+1$	$2N_b+1$
${\sigma }^2$	${0.6}^2$	${0.6}^2/N_{\text{bin}}$
$\beta$	${\sigma }^2/100$	${\sigma }^2/100$

Table 2. TABLE II: IAE, background NMSE and execution time of Depth TV and the proposed method for the blocks and leaves dataset with and without the backplane.

Backplane		Yes			No
Acq. time [ms]		1	0.1	0.01	1	0.1	0.01
IAE [photons]	Depth TV	36	3.7	0.5	45	4.9	0.8
	MuSaPoP	14	1.7	0.3	19	2.4	0.4
Bkg. NMSE	Depth TV	0.12	>1	>1	0.12	>1	>1
	MuSaPoP	0.04	0.11	0.36	0.04	0.11	0.26
Execution time [h]	Depth TV	19.8	17.9	17.7	19.3	17.9	17.7
	MuSaPoP	2.8	1.4	1.0	1.5	1	0.8

Equations

$$\mathbf{\Phi }=\{({\mathbf{c}}_n,{\mathbf{r}}_n),n=1,\dots ,N_{\Phi }\}$$

$$z_{i,j,\ell ,t}|\mathbf{\Phi },b_{i,j,\ell }\sim \mathcal{P}\left(g_{i,j,\ell }\left(\sum _{{\mathcal{N}}_{i,j}}{r}_{n,\ell }{h}_{\ell }(t-t_n)+b_{i,j,\ell }\right)\right)$$

$$p(\mathbf{Z}|\mathbf{\Phi },\mathbf{B})=\prod _{i=1}^{N_c}\prod _{j=1}^{N_r}\prod _{\ell =1}^L\prod _{t=1}^{T}p(z_{i,j,\ell ,t}|\mathbf{\Phi },b_{i,j,\ell }).$$

$$f({\mathbf{\Phi }}_c)\propto f_1({\mathbf{\Phi }}_c)f_2({\mathbf{\Phi }}_c|{\gamma }_a,{\lambda }_a)$$

$$f_{1}(\boldsymbol{\Phi}_{c})\propto\left\{\begin{array}[]{ll}0&\text{if }\exists~{}n\neq n^{\prime}:x_{n}=x_{n^{\prime}},y_{n}=y_{n^{\prime}}\\ &\text{and }|t_{n}-t_{n^{\prime}}|<d_{\min}\\ 1&\mbox{otherwise}\end{array}\right.$$

$$f_2({\mathbf{\Phi }}_c|{\gamma }_a,{\lambda }_a)\propto {\lambda }_a^{N_{\Phi }}{\gamma }_a^{-m\left({\bigcup }_{n=1}^{N_{\Phi }}S({\mathbf{c}}_n)\right)}$$

$$m_{n,\ell }=\log (r_{n,\ell })$$

$$p({\mathbf{\Phi }}_r|{\sigma }^2,\beta )=\prod _{\ell =1}^{L}p({\mathbf{m}}_{\ell }|{\mathbf{\Phi }}_c,{\sigma }^2,\beta ).$$

$${\mathbf{m}}_{\ell }|{\sigma }^2,\beta ,{\mathbf{\Phi }}_c\sim \mathcal{N}(0,{\sigma }^2{\mathbf{P}}^{-1})$$

$$[\mathbf{P}{]}_{n,n^{\prime }}=\{\begin{array}{ll}\beta +{\sum }_{\tilde{n}\in {\mathcal{M}}_{pp}({\mathbf{c}}_n)}\frac{1}{d({\mathbf{c}}_n;{\mathbf{c}}_{\tilde{n}})}& \text{if n=n′}\\ -\frac{1}{d({\mathbf{c}}_n;{\mathbf{c}}_{n^{\prime }})}& \text{if cn∈Mpp(cn′)}\\ 0& \text{otherwise}\end{array}$$

$$\{\begin{array}{rl}p(\mathbf{B}|\mathbf{K},\mathbf{\Theta })& =\prod _{i=1}^{N_r}\prod _{j=1}^{N_c}\prod _{\ell =1}^{L}p(b_{i,j,\ell }|k_{i,j,\ell },{\theta }_{i,j,\ell })\\ b_{i,j,\ell }|k_{i,j,\ell },{\theta }_{i,j\ell }& \sim \mathcal{G}(k_{i,j,\ell },{\theta }_{i,j\ell })\end{array}$$

$$(\mathbf{K},\mathbf{\Theta })=\underset{\mathbf{K},\mathbf{\Theta }}{\mathrm{argmin}}\text{KL}(q||p).$$

$$\{\begin{array}{rl}{\bar{z}}_{i,j,\ell ,b}|b_{i,j,\ell }& \sim \mathcal{P}(g_{i,j,\ell }{v}_{i,j,\ell }{b}_{i,j,\ell })\\ b_{i,j,\ell }& =\exp ({\tilde{b}}_{i,j,\ell }+{\mu }_{\ell })\\ {\tilde{\mathbf{b}}}_{\ell }|{\alpha }_B& \sim \mathcal{N}(0,{\alpha }_B{\mathbf{D}}^{-1})\end{array}$$

$$\begin{gathered} p({\mathbf{\Phi }}_c,{\mathbf{\Phi }}_r,\mathbf{B}|\mathbf{Z},\Psi )\propto p(\mathbf{Z}|{\mathbf{\Phi }}_c,{\mathbf{\Phi }}_r,\mathbf{B})p({\mathbf{\Phi }}_r|{\mathbf{\Phi }}_c,{\sigma }^2,\beta )\times \\ {f}_1({\mathbf{\Phi }}_c|{\gamma }_a,{\lambda }_a)f_2({\mathbf{\Phi }}_c|{\gamma }_{st})\pi ({\mathbf{\Phi }}_c)p(\mathbf{B}|\mathbf{K},\mathbf{\Theta }) \end{gathered}$$

$$\hat{\mathbf{\Phi }}=\underset{\mathbf{\Phi }}{\mathrm{argmax}}p(\mathbf{\Phi },\mathbf{B}|\mathbf{Z},\Psi )$$

$$\hat{\mathbf{B}}=\mathbb{E}\{\mathbf{B}|\mathbf{Z},\Psi \}.$$

$$\{{\mathbf{\Phi }}^{(s)},{\mathbf{B}}^{(s)}\forall s=0,1,\dots ,N_m-1\}.$$

$$\hat{\mathbf{\Phi }}$$

$$\hat{\mathbf{B}}$$

$$r\left(\mathbf{\theta },{\mathbf{\theta }}^{\prime }\right)=\frac{p({\mathbf{\theta }}^{\prime }|\mathbf{Z},\Psi )K(\mathbf{\theta }|{\mathbf{\theta }}^{\prime })p({\mathbf{u}}^{\prime })}{p(\mathbf{\theta }|\mathbf{Z},\Psi )K({\mathbf{\theta }}^{\prime }|\mathbf{\theta })p(\mathbf{u})}\left|\frac{\partial g(\mathbf{\theta },\mathbf{u})}{\partial (\mathbf{\theta },\mathbf{u})}\right|$$

$$\left\{\begin{aligned} {\color[rgb]{0,0,0}u}_{\ell}&\sim\mathcal{U}(0,1),b_{i,j,\ell}^{\prime}={\color[rgb]{0,0,0}u}_{\ell}b_{i,j,\ell}\\ e^{m_{N_{\Phi}+1,\ell}}&=w_{\ell}(1-{\color[rgb]{0,0,0}u}_{\ell})b_{i,j,\ell}\frac{T}{\sum_{t=1}^{T}h_{\ell}(t)}\end{aligned}\right.,$$

$$b_{i,j,\ell }^{\prime }=b_{i,j,\ell }+e^{m_{N_{\Phi }+1,\ell }}\frac{{\sum }_{t=1}^T{h}_{\ell }(t)}{w_{\ell }T}\forall \ell =1,\dots ,L.$$

$$m_{n,\ell }^{\prime }\sim \mathcal{N}\left(m_{n,\ell },{\delta }_m\right)\forall \ell =1,\dots ,L.$$

$$t_n^{\prime }\sim \mathcal{N}\left(t_n,{\delta }_t\right)$$

$$\{\begin{array}{rl}{u}_{\ell }& \sim \mathcal{B}(\eta ,\eta )\\ m_{k_1,\ell }^{\prime }& =m_{n,\ell }+\log (u_{\ell })\\ m_{k_2,\ell }^{\prime }& =m_{n,\ell }+\log (1-u_{\ell })\end{array}$$

$$\{\begin{array}{rl}{s}_{\ell }& \sim \mathcal{B}\text{e}(0.5)\\ \Delta & \sim \mathcal{U}(d_{\min },L_h)\\ t_{k_1}^{\prime }& =t_n+(-1{)}^{s_{\ell }}\Delta \frac{{\sum }_{\ell =1}^L(1-u_{\ell })e^{m_{n,\ell }}}{{\sum }_{\ell =1}^L{e}^{m_{n,\ell }}}\\ t_{k_2}^{\prime }& =t_n-(-1{)}^{s_{\ell }}\Delta \frac{{\sum }_{\ell =1}^L{u}_{\ell }{e}^{m_{n,\ell }}}{{\sum }_{\ell =1}^L{e}^{m_{n,\ell }}}\end{array}$$

$$d_{\min }<\left|t_{k_1}-t_{k_2}\right|\le L_h\forall \ell =1,\dots ,L$$

$$\{\begin{array}{rl}& e^{m_{n,\ell }^{\prime }}=e^{m_{k_1,\ell }}+e^{m_{k_2,\ell }}\forall \ell =1,\dots ,L\\ & t_n^{\prime }=t_{k_1}\frac{{\sum }_{\ell =1}^L{e}^{m_{k_1,\ell }}}{{\sum }_{\ell =1}^L{e}^{m_{k_1,\ell }}+e^{m_{k_2,\ell }}}+t_{k_2}\frac{{\sum }_{\ell =1}^L{e}^{m_{k_2,\ell }}}{{\sum }_{\ell =1}^L{e}^{m_{k_1,\ell }}+e^{m_{k_2,\ell }}}\end{array}$$

$$z_{i,j,t,\ell }=$$

$${\tilde{z}}_{i,j,t,\ell ,b}\sim$$

$${\tilde{z}}_{i,j,t,\ell ,n}\sim$$

$$\{\begin{array}{rl}{\tilde{z}}_{i,j,t,\ell ,b}& \sim \mathcal{B}\left(z_{i,j,l,t},\frac{b_{i,j,\ell }}{{\sum }_{\begin{array}{c}n:(x_n,y_n)=(i,j)\end{array}}{r}_{n,\ell }h(t-t_n)+b_{i,j,\ell }}\right)\\ b_{i,j,\ell }& \sim \mathcal{G}\left(r_{i,j,\ell }+\sum _{t=1}^T{\tilde{z}}_{i,j,t,\ell ,b},\frac{{\theta }_{i,j,\ell }}{g_{i,j,\ell }T{\theta }_{i,j,\ell }+1}\right)\end{array}$$

Full text

I Introduction

Single-photon Lidar devices provide accurate range information by constructing, for each pixel, a histogram of time delays between emitted light pulses and photon detections. Using time correlated single-photon counting (TCSPC) technology, Lidar systems are able to provide accurate depth information (of the order of centimetres) over very long ranges, while allowing the use of laser sources of low power [1]. The acquired range information (3D structure) has many important applications, such as self-driving cars [2], the study of structures behind dense forests [3] and environmental monitoring [4]. Multispectral Lidar (MSL) systems gather measurements at many spectral bands, making it possible to distinguish distinct materials, as illustrated in Figure 1. For example, spectral diversity was used in [5] to differentiate leaves from trunks and in [6] to estimate plant area indices and abundance profiles. The MSL modality consists of constructing one histogram of time delays per wavelength, as shown in Figure 2. The spectral diversity can be obtained either using a supercontinuum laser source [6, 7] or multiple lasers [8]. The detector generally consists of a spatial form of wavelength routing to demultiplex the channels [6, 7, 8] or wavelength-to-time codification [9]. Recovering spatial and spectral information from MSL data is a challenging task, specially in scenes with strong ambient illumination (i.e., multiple spurious detections) or when the acquisition time is very low (i.e., very few photon detections per histogram). Moreover, in a general setting, it is possible to find more than one object per pixel. This phenomenon occurs in scenes where the light goes through semi-transparent materials (e.g., glass), or when the laser footprint is such that multiple surfaces appear in the field of view. Thus, several signal processing algorithms have been proposed to address these challenges: while many algorithms are available for single-wavelength Lidar, either assuming a single surface per pixel [10, 11, 12] or multiple surfaces per pixel [13, 14, 15], to the best of our knowledge, only the single-surface-per-pixel case was studied in the multispectral case [7, 16, 17].

This single-depth assumption greatly simplifies the reconstruction problem, as it significantly reduces memory requirements and overall complexity. Datasets containing dozens of wavelengths can be prohibitively large for practical 3D reconstruction algorithms, both in terms of memory and computing requirements. For example, a typical MSL hypercube with 32 wavelengths has more than $10^{9}$ data voxels. To alleviate this problem, some compressive acquisition strategies have been proposed. While TCSPC technology hinders compressive techniques along the depth axis111A coarse time-of-flight gating is applied to the photon detections, hindering measurements of an arbitrary subset of histogram bins or linear combination of them. However, other alternatives such as gated cameras [18] can provide such measurements., reducing the number of measurements can be achieved by integrating multiple wavelengths in a single histogram [9] or measuring fewer histograms (i.e., subsampling) [16]. The wavelength-to-time approach proposed by Ren et al. [9] is not well-suited in the presence of multiple surfaces per pixel. Indeed, this method compresses $L$ histograms (associated with $L$ wavelengths) into a single waveform by shifting in time the photon detections according to each measured wavelength. While significantly reducing the data size, the resulting likelihood becomes highly multimodal and extremely difficult to handle in the presence of multiple surfaces. Different random subsampling schemes were studied in [16] without obtaining any significant differences in terms of reconstruction quality in the low-photon count regime.

In this work, we investigate a new pseudo-random subsampling scheme for low-photon count MSL data based on ideas from coded aperture design [19, 20]. By choosing the pixels measured for each wavelength in a more principled way, we achieve better results than the completely random schemes of Altmann et al. [16]. Furthermore, the proposed subsampling strategy can be easily implemented in many Lidar systems, reducing the total number of measurements, i.e., the time to acquire an MSL frame. Raster-scan systems using a laser supercontinuum source [6, 7] can be easily modified to measure only a subset of pixels per wavelength. More interestingly, single-wavelength array technology [21] can be combined with coded apertures [19], which acquire different wavelengths at each pixel.

Furthermore, we propose a new method to perform 3D reconstruction from subsampled MSL data, which is capable of finding multiple surfaces per pixel. The novel method draws ideas from a recently published algorithm named ManiPoP [15], provided state-of-the-art reconstructions in the multiple surface, single-wavelength case. However, due to the significantly larger dimensionality of multispectral data, we propose to modify the Bayesian model and estimation strategy of ManiPoP. Adopting a Bayesian framework similar to [15], we assign a spatial point process prior to promote spatial correlation and a Gaussian Markov random field prior to regularize the spectral reflectivity within surfaces/objects. Inference using the resulting posterior distribution is performed using a reversible jump Markov chain Monte Carlo algorithm (RJ-MCMC), coupled with a multiresolution approach that improves the convergence speed and reduces the total computing time of the algorithm. We introduce new RJ-MCMC proposals, which take into account the additional spectral dimension and improve the acceptance ratio, compared to the ones proposed in [15]. Moreover, we propose an empirical Bayes approach to build the prior distribution associated with the background detections, which further improves the convergence of the RJ-MCMC sampler. Contrary to multi-depth methods that require storage of dense volumetric estimates (e.g., [13, 14] in the single-wavelength case and [22] in the multi-temporal case), the memory requirements of the proposed method are minimal (just the 3D points and spectral signatures are stored in memory), enabling the acquisition and processing of very large datasets (dozens of wavelengths and hundreds of pixels).

The main contributions of this work are

• •

a new Bayesian 3D reconstruction algorithm for subsampled MSL data with multiple surfaces per pixel

• •

the analysis of an appropriate subsampling scheme based on coded aperture design, which provides better results than completely random alternatives.

The new algorithm is referred to as MuSaPoP, as it models MultiSpectrAl Lidar signals using POint Processes. The remainder of the paper is organized as follows. Section II recalls the classical observation model for single-photon MSL data. Sections III and IV present the Bayesian 3D reconstruction algorithm and the associated RJ-MCMC sampler. Section V details the principled subsampling strategy. Experiments performed with synthetic and real Lidar data are introduced and discussed in Section VI. Conclusions and future work are finally reported in Section VII.

II Single-photon multispectral Lidar

A full multispectral Lidar data hypercube $\boldsymbol{Z}\in\mathbb{Z}_{+}^{N_{r}\times N_{c}\times L\times T}$ consists of discrete photon count measurements, where $\mathbb{Z}_{+}=\{0,1,\dots\}$ is the set of positive integers, $N_{r}$ and $N_{c}$ are the numbers of vertical and horizontal pixels respectively, $L$ is the number of acquired wavelengths and $T$ is the histogram length (i.e., the number of time bins). The 3D reconstruction algorithm estimates a point cloud, referred to as $\boldsymbol{\Phi}$, from the data hypercube $\boldsymbol{Z}$. The 3D point cloud is represented by an unordered set of points, that is

[TABLE]

where $N_{\Phi}$ is the total number of points, and $\boldsymbol{c}_{n}=[x_{n},y_{n},t_{n}]^{T}\in[1,N_{r}]\times[1,N_{c}]\times[1,T]\subset\mathbb{Z}_{+}^{3}$ and $\boldsymbol{r}_{n}=[r_{n,1},\dots,r_{n,L}]^{T}\in\mathbb{R}_{+}^{L}$ are the coordinate vector and the spectral response of the $n$th point, respectively. The observed photon count at pixel $(i,j)$, bin $t$ and spectral band $\ell$ follows a Poisson distribution [16], whose intensity is a mixture of the background level, denoted by $b_{i,j,\ell}$, and the responses of the surfaces present in that pixel, i.e.,

[TABLE]

where $\mathcal{P}(\cdot)$ denotes the Poisson distribution, $g_{i,j,\ell}\in\{0,1\}$ is a known binary variable that indicates whether wavelength $\ell$ at pixel $(i,j)$ has been acquired/observed or not, $\mathcal{N}_{i,j}=\{n:(x_{n},y_{n})=(i,j)\}$ is the set of points corresponding to pixel $(i,j)$ and wavelength $\ell$, and $h_{\ell}(t)$ is the impulse response of the Lidar device at wavelength $\ell$, which can be measured using a spectralon panel during a calibration step. Note that to lighten notations, we assume that the impulse responses are only wavelength-dependent but the algorithm can be easily adapted when the responses vary among pixels. Assuming mutual independence between noise realizations in different bins, pixels and spectral bands [1], the likelihood of the proposed model can be written as

[TABLE]

For clarity in the notation, we will also denote the set of point coordinates as $\boldsymbol{\Phi}_{c}=\{\boldsymbol{c}_{n},n=1,\dots,N_{\Phi}\}$ and the set of spectral responses as $\boldsymbol{\Phi}_{r}=\{\boldsymbol{r}_{n},n=1,\dots,N_{\Phi}\}$. The set of all background levels is denoted by $\boldsymbol{B}=[\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{L}]\in\mathbb{R}_{+}^{N_{r}\times N_{c}\times L}$, which is the concatenation of $L$ images $\boldsymbol{b}_{\ell}$, one for each wavelength. The cube of binary measurements is designated by $\boldsymbol{G}\in\{0,1\}^{N_{r}\times N_{c}\times L}$, where $[\boldsymbol{G}]_{i,j,\ell}=g_{i,j,\ell}$. Note that the model used in the ManiPoP algorithm [15] can be obtained from (2) by setting all the binary variables to 1, and considering only one band, i.e., $L=1$.

III Multiple-return multiple-wavelength 3D reconstruction

Recovering the position of the objects $\boldsymbol{c}_{n}$, their spectral signatures $\boldsymbol{r}_{n}$ and background levels $\boldsymbol{B}$ from the subsampled MSL data $\boldsymbol{Z}$ is an ill-posed inverse problem, as many solutions can explain the observed photon counts. Thus, prior regularization is necessary to promote reconstructions in a set of feasible 3D point clouds. In this work, we adopt a Bayesian framework, which allows us to include prior knowledge about the scene through tailored prior distributions assigned to the parameters of interest.

III-A Prior distributions

The proposed model considers prior regularization for the point positions and reflectivity. As explained in Section III-A3, an empirical Bayes prior [23] is assigned to the background levels.

III-A1 Spatial configuration

We adopt the spatial prior distribution of 3D points developed in the ManiPoP algorithm. This distribution is designed to promote spatial correlation between points within one surface and repulsion between points belonging to different surfaces. While only a brief summary is included below, we refer the reader to [15] for a detailed discussion about this prior model. The spatial point process prior for the position of the points is modelled by a density defined with respect to a Poisson point process reference measure [24, Chapter 9], i.e.,

[TABLE]

where $f_{1}(\boldsymbol{\Phi}_{c})$ and $f_{2}(\boldsymbol{\Phi}_{c}|\gamma_{a},\lambda_{a})$ are the Strauss and area interaction [25] processes respectively. The repulsive Strauss process is written as

[TABLE]

where $d_{\min}$ is the minimum separation between two surfaces in the same pixel. Attraction between points within the same surface is promoted by the area interaction process, that is

[TABLE]

where $m(\cdot)$ denotes the standard Lebesgue measure, $S(\boldsymbol{c}_{n})$ defines a convex set around the point $\boldsymbol{c}_{n}$, and $\gamma_{a}$ and $\lambda_{a}$ are two hyperparameters, accounting for the amount of attraction and total number of points, respectively. Both densities define Markovian interactions between points, only correlating points in a local neighbourhood. Moreover, the combination of both processes implicitly defines a connected-surface structure, which is used to model 2D manifolds in a 3D space.

III-A2 Reflectivity

The spectral signatures are related to the materials of the surfaces [17]. Neighbouring points corresponding to a surface composed of a specific material show similar spectral signatures. This prior information is added to the model using a Gaussian Markov random field distribution, where the neighbours are defined by the connected-surface structure of the point process prior. First, in order to avoid the positivity constraint on the intensity $r_{n,\ell}$, we use the canonical form [26, 27, 28],

[TABLE]

where $m_{n,\ell}\in\mathbb{R}$ denotes the log-intensity of the $n$th point at band $\ell$. As multispectral devices only acquire dozens of well-separated wavelengths, the spectral measurements within a pixel do not show significant correlation. Hence, although potential correlations between wavelengths could be modelled, we choose here to neglect this correlation to keep the estimation strategy tractable. As a consequence, we consider the following reflectivity prior model

[TABLE]

Spatial correlations between log-intensity values in neighbouring pixels are defined according to the distribution

[TABLE]

where $\sigma^{2}$ and $\beta$ are hyperparameters controlling the level of smoothness. The precision matrix $\boldsymbol{P}$ is very sparse due to the Markovian structure, being defined by

[TABLE]

where $\mathcal{M}_{pp}(\boldsymbol{c}_{n})$ is the set of neighbours of point $\boldsymbol{c}_{n}$, which is obtained using the connected-surface structure illustrated in Figure 3, and $d(\boldsymbol{c}_{n};\boldsymbol{c}_{\tilde{n}})$ is the Euclidean distance between two points, normalized according to the camera parameters of the scene to have a physical meaning [29].

III-A3 Background levels

Background detections are due to detector dark counts and ambient illumination, as explained in [13, 30]. If the transceiver system is mono-static222In mono-static Lidar systems, the laser and detector are coaxial, whereas in bi-static systems, the source and detector do not share the same axis. [1], the set of background levels can be interpreted as a multispectral passive image of the scene, as background detections generally come from ambient illumination reflected onto the target and collected by the single-photon detector. In this case, the background levels are strongly spatially correlated within each wavelength. However, in bi-static systems [13], the transmit and receive channels do not share the same objective lens aperture, yielding weakly or uncorrelated background detections. Although not showing a strong spatial correlation in this second case, all the background levels have similar values, which also serves as prior information.

Independent prior distributions

In order to simultaneously model potential spatial correlation and ensure the positivity of the background levels, ManiPoP uses a 2D gamma Markov random field, which was introduced by Dikmen et al. in [31]. However, this prior is not well suited for MSL data as it introduces an undesired penalization for large background levels, whose negative effects are amplified when considering multiple bands (see Appendix A for details). Other alternatives such as Gaussian Markov random fields [26] cannot be sampled directly in closed form, requiring proposals with a rejection step, whose mixing and convergence scale badly with the dimension of the spectral cube, as shown in [32]. To alleviate these problems, we assign independent gamma priors, i.e.,

[TABLE]

where $[\boldsymbol{\Theta}]_{i,j,\ell}=\theta_{i,j,\ell}$ and $[\boldsymbol{K}]_{i,j,\ell}=k_{i,j,\ell}$ are the shape and scale hyperparameters of the gamma distributions. Despite using independent priors, we can capture the spatial correlation by setting the hyperparameters $(\boldsymbol{K},\boldsymbol{\Theta})$ appropriately. More precisely, in a similar fashion to variational Bayes [33] or expectation propagation [34] methods, in order to simplify the estimation of $\boldsymbol{B}$, we specify (8) such that $p(\boldsymbol{B}|\boldsymbol{K},\boldsymbol{\Theta})$ is similar to another distribution $q(\boldsymbol{B})=\prod_{\ell=1}^{L}q_{\ell}(\boldsymbol{b}_{\ell})$ which explicitly correlates the background levels in neighbouring pixels and assumes mutual independence between spectral bands. Here, we use as a similarity criterion the Kullback-Leibler divergence

[TABLE]

As discussed in Appendix B, solving (9) can be achieved by computing expectations with respect to $q(\boldsymbol{B})$.

Empirical Bayes approach

To ensure that the prior model (8) is informative, a suitable distribution $q(\boldsymbol{B})$ should be chosen. Assuming that we have a coarse estimate of the point cloud (this information will be obtained using the multiresolution approach detailed in Section IV-B), we build the distribution $q(\boldsymbol{B})$ following an empirical Bayes approach, as illustrated in Figure 4. First, one can discard almost all the signal photons in the dataset by removing the photons detected in the compact support of $h_{\ell}(t)$ around each point (see Figure 4, central subplot). The number of bins that is not excluded in each pixel is referred to as $v_{i,j,\ell}$. Secondly, we integrate the remaining photons of each pixel, obtaining a coarse estimate of the per-pixel background photon levels, denoted by $\bar{z}_{i,j,\ell,b}$. We then define $q_{\ell}(\boldsymbol{b}_{\ell})\propto p(\boldsymbol{\bar{z}}_{\ell}|\boldsymbol{b}_{\ell})p(\boldsymbol{b}_{\ell}|\alpha_{B})$ with

[TABLE]

where $\boldsymbol{b}_{\ell}$ a vectorized image of background levels at wavelength $\ell$, $\boldsymbol{D}\in\mathbb{R}^{NrNc\times NrNc}$ a positive semidefinite matrix, $\alpha_{B}$ is a fixed hyperparameter controlling the degree of smoothness. In mono-static systems, a two-dimensional Laplacian filter is chosen for $\boldsymbol{D}$ to promote spatial correlation [26], whereas in bi-static systems, $\boldsymbol{D}$ is replaced by the identity matrix, only penalizing large background levels. $\mu_{\ell}$ is a translation parameter centring $\tilde{b}_{i,j,\ell}$ in the linear part of the exponential function and is defined as $\mu_{\ell}=\log(\frac{1}{N_{r}N_{c}}\sum_{i}^{N_{r}}\sum_{j}^{N_{c}}\frac{\bar{z}_{i,j,\ell,b}g_{i,j,\ell}}{v_{i,j,\ell}})$. As mentioned above, solving (9) requires the computation of expectations with respect to $q(\boldsymbol{B})$ which are unfortunately not available in closed form. Instead of using additional MCMC sampling to find numerical approximations (detailed in Appendix B), obtaining samples from (10) is more attractive both in terms of convergence and computational complexity than using the original full datacube which includes mixtures of background and signal photons. Indeed, (10) simply involves integrated photon counts (over the range dimension). Moreover, given the independence property of $q(\boldsymbol{B})$ among spectral bands, all the bands can be processed independently in parallel when sampling $\boldsymbol{B}$ (see Section IV-A1).

III-B Posterior distribution

Following Bayes theorem, the joint posterior distribution of the model parameters is given by

[TABLE]

where $\Psi$ denotes the set of hyperparameters $\Psi=\{\gamma_{a},\lambda_{a},\gamma_{s},\sigma^{2},\beta,\boldsymbol{K},\boldsymbol{\Theta}\}$ and $\pi(\cdot)$ is the Poisson point process reference measure. Figure 5 shows the directed acyclic graph associated with the proposed hierarchical Bayesian model.

IV Inference

In this work, we compute the same posterior statistics as in [15]: the point cloud positions and spectral signatures are estimated using the maximum-a-posteriori (MAP) estimator

[TABLE]

and the minimum mean squared error estimator is considered for $\boldsymbol{B}$

[TABLE]

As this expectation cannot be derived analytically, we propose to investigate Markov Chain Monte Carlo (MCMC) simulation methods. As in [15], we use a reversible jump MCMC algorithm that can handle the varying dimension nature of the spatial point process. This sampler generates $N_{m}$ samples of $\boldsymbol{\Phi}$ and $\boldsymbol{B}$ from the posterior distribution (11) denoted as

[TABLE]

These samples are then used to approximate the statistics of interest, i.e.,

[TABLE]

where $N_{\textrm{bi}}$ is the number of burn-in iterations.

IV-A Reversible jump MCMC moves

While other stochastic simulation algorithms for varying dimensions can be used [24], we choose an RJ-MCMC sampler, as it allows us to build proposals tailored for the MSL reconstruction problem. RJ-MCMC can be interpreted as a natural extension of the Metropolis-Hastings algorithm for problems with an a priori unknown dimensionality. Given the actual state of the chain $\boldsymbol{\theta}=\{\boldsymbol{\Phi},\boldsymbol{B}\}$ of model order $N_{\Phi}$, a random vector of auxiliary variables $\boldsymbol{u}$ is generated to create a new state $\boldsymbol{\theta}^{\prime}=\{\boldsymbol{\Phi}^{\prime},\boldsymbol{B^{\prime}}\}$ of model order $N_{\Phi^{\prime}}$, according to an appropriate deterministic function $\boldsymbol{\theta}^{\prime}=g(\boldsymbol{\theta},\boldsymbol{u})$. To ensure reversibility, an inverse mapping with auxiliary random variables $\boldsymbol{u^{\prime}}$ has to exist such that $\boldsymbol{\theta}=g^{-1}(\boldsymbol{\theta}^{\prime},\boldsymbol{u^{\prime}})$. The move $\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}^{\prime}$ is accepted or rejected with probability $\rho=\min\{1,r\left(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)\}$, where $r(\cdot,\cdot)$ is defined as

[TABLE]

and $K(\boldsymbol{\theta^{\prime}}|\boldsymbol{\theta})$ is the probability of proposing the move $\boldsymbol{\theta}\rightarrow\boldsymbol{\theta^{\prime}}$, $p(\boldsymbol{u})$ is the probability distribution of the random vector $\boldsymbol{u}$, and $\left|\frac{\partial g(\boldsymbol{\theta},\boldsymbol{u})}{\partial(\boldsymbol{\theta},\boldsymbol{u})}\right|$ is the Jacobian of the mapping $g(\cdot)$. The RJ-MCMC algorithm performs birth/death, dilation/erosion, spatial and mark shifts, and split/merge moves with probabilities $p_{\textrm{birth}}$, $1-p_{\textrm{birth}}$, $p_{\textrm{dilation}}$, $1-p_{\textrm{dilation}}$, $p_{\textrm{shift}}$, $p_{\textrm{mark}}$, $p_{\textrm{split}}$ and $1-p_{\textrm{split}}$ respectively. Due to the Markovian nature of the prior distributions, all the proposed moves are local, having a complexity proportional to the size of the neighbourhood. These moves are detailed in the following subsections. For ease of presentation, we summarize the main aspects of each move, inviting the reader to consult Appendix C for more details.

Birth and death moves

The birth move proposes a new point $\boldsymbol{\theta^{\prime}}=(\boldsymbol{c}_{N_{\Phi}+1},\boldsymbol{m}_{N_{\Phi}+1})$ uniformly at random in the 3D cube. The spectral signature of the new point is proposed by extracting a fraction $(1-u_{\ell})$ from the current value of the background level $b_{i,j,\ell}$ according to the signal-to-background-ratio (SBR) [11] $w_{\ell}$, that is for each wavelength $\ell$

[TABLE]

where $\mathcal{U}(0,1)$ denotes the uniform distribution on the interval $(0,1)$. The death move proposes the removal of a point. In contrast to the birth move, we modify the background level according to

[TABLE]

Dilation and erosion moves

Birth moves have low acceptance ratio, as the probability of randomly proposing a point within or close to the surfaces of interest is very low. However, this problem can be overcome by using the current estimation of the surface to propose in regions of high probability. The dilation move proposes a point inside the neighbourhood of an existing surface with uniform probability across all possible neighbouring positions where a point can be added. Contrary to [15], where the intensity samples are generated according to the prior distribution, the spectral signature is sampled in the same way as the birth move (16). The complementary move (named erosion), proposes to remove a point $\boldsymbol{c}_{n}$ with one or more neighbours. In this case, the background is updated in the same way as the death move.

Mark and shift moves

As in ManiPoP, the mark move updates the log-intensity of a randomly chosen point $\boldsymbol{c}_{n}$. Each wavelength is updated independently using a Gaussian proposal with variance $\delta_{m}$ as a proposal (also known as Metropolis Gaussian random walk), that is

[TABLE]

Similarly, the shift move updates the position of a uniformly chosen point using a Gaussian proposal with variance $\delta_{t}$

[TABLE]

The values of $\delta_{m}$ and $\delta_{t}$ are adjusted by cross-validation333Intensities are normalized to belong to a fixed interval across datasets. Hence, we can fix the variance of the proposal to achieve similar acceptance ratios. to yield an acceptance ratio close to $41\%$ for each move, which is the optimal value for a one dimensional Metropolis random walk, as explained in [24, Chapter 4].

Split and merge moves

Some pixels might present two points with overlapping impulse responses in depth. In such cases, a death move followed by two birth moves would happen with very low probability. Hence, as in ManiPoP, we propose a split move, which randomly picks a point $(\boldsymbol{c}_{n},\boldsymbol{m}_{n})$ and proposes two new points, $(\boldsymbol{c}^{\prime}_{k_{1}},\boldsymbol{m}^{\prime}_{k_{1}})$ and $(\boldsymbol{c}^{\prime}_{k_{2}},\boldsymbol{m}^{\prime}_{k_{2}})$. The log-intensity is proposed for each wavelength following the mapping

[TABLE]

where $\mathcal{B}(\cdot)$ denotes the beta distribution and $\eta$ is a fixed parameter. The new positions are determined according to

[TABLE]

where $\mathcal{B}$e$(\cdot)$ denotes the Bernoulli distribution, $L_{h}$ is the length in bins of the instrumental response at the wavelength with longest impulse response. The complementary move, named merge move, is performed by randomly choosing two points $\boldsymbol{c}_{k_{1}}$ and $\boldsymbol{c}_{k_{2}}$ inside the same pixel ($x_{k_{1}}=x_{k_{2}}$ and $y_{k_{1}}=y_{k_{2}}$) that satisfy the condition

[TABLE]

The merged point $(\boldsymbol{c^{\prime}}_{n},\boldsymbol{r}^{\prime}_{n})$ is obtained by the inverse mapping

[TABLE]

which preserves the total pixel intensity and weights the spatial shift of each peak according to the sum of the intensity values.

IV-A1 Sampling the background

In order to exploit the conjugacy between the Poisson likelihood and gamma priors for the background levels, we use a data augmentation scheme as in [35], which classifies each photon-detection according to the source (target(s) or background), i.e.,

[TABLE]

where $\tilde{z}_{i,j,t,\ell,n}$ are the photons at bin $t$ associated with the $n$th surface and $\tilde{z}_{i,j,t,\ell,b}$ are the ones associated with the background. In the augmented space defined by $(\tilde{z}_{i,j,t,\ell,n},\tilde{z}_{i,j,t,\ell,b})$, the background levels are sampled as follows

[TABLE]

where $\mathcal{B}(\cdot)$ denotes the Binomial distribution. The transition kernel defined by (23) produces samples of $b_{i,j,\ell}$ distributed according to the marginal posterior distribution of $\boldsymbol{B}$. In practice, we observed that only one iteration of this kernel is sufficient.

IV-B Full algorithm

We adopt a multi-resolution approach to speed up the convergence of the RJ-MCMC algorithm, in a fashion similar to [15]. The dataset is downsampled by integrating photon detections in patches of $N_{\textrm{bin}}\times N_{\textrm{bin}}$ pixels. Hence, the number of pixels is reduced by a factor of $N_{\textrm{bin}}^{2}$, meaning less points and background levels to infer with $N_{\textrm{bin}}^{2}$ times more photons per pixel. The estimated point cloud at the coarse scale is upsampled using a simple nearest neighbour algorithm and used as initialization for the next (finer) scale. In all our experiments we repeat the process for $K=3$ scales. The background hyperpriors $\boldsymbol{K}$ and $\boldsymbol{\Theta}$ are initialized with non-informative values, i.e., $k_{i,j,\ell}=0.01$ and $\theta_{i,j,\ell}=100$ for all pixels $(i,j)$ and wavelengths $\ell$. In finer scales, these hyperparameters are computed using the algorithm detailed in Section III-A3. The multi-resolution approach is finally summarized in Algorithm 1.

V Subsampling strategy

Despite not being able to design compressive measurements along the depth axis, we can still reduce the number of measurements in the two spatial (horizontal and vertical) dimensions and in the spectral dimension [16]. Given the point positions, recovering their reflectivity profile reduces to a multispectral image restoration problem using measured data corrupted by Poisson noise. While many compressive sensing strategies have been proposed for measurements under this noise assumption [36, 37, 38], MSL datasets have an additional limitation if multiple surfaces per pixel are considered: photon-detections belonging to different wavelengths cannot be integrated into a single histogram, as the reconstruction problem generally becomes highly non-convex, preventing practical reconstruction algorithms444As mentioned in the introduction, the system presented in [9] considers the integration of photons belonging to different histograms, but is limited to one surface per pixel.. Indeed, summing histograms associated with different wavelengths and including multiple peaks generates histograms with even more peaks (possibly overlapping), which makes the 3D reconstruction and the reflectivity estimation more difficult. As a consequence, we only consider subsampling of depth histograms, which incorporates all of the practical sampling limitations. Following the formulation of the observation model (2), the subsampling strategy consists of choosing the binary coefficients $\boldsymbol{G}$ for a given compression level $W/L$, with $W$ the average number of observed band per pixel. Several subsampling strategies have been proposed for different applications, such as halftoning and stippling [39], rendering, compressive spectral imaging [40, 41, 42], compressive computed tomography [43, 44], geometry processing [45], amongst others [46, 47, 48]. These approaches exploit the sampling geometry of the sensing devices to design a set of criteria. Similarly to coded aperture snapshot spectral imaging systems, the distribution of reflectivity profiles of 3D surfaces in real natural scenes suggests uniform sampling in the row, column and spectral axes. Following the design in Correa et al. [20], the coefficients are chosen according to the spatiotemporal characteristics of blue noise, which distributes the measurements in spectral and spatial dimensions as homogeneously as possible [49]. The binary cube $\boldsymbol{G}$ is obtained by minimizing the variance of (weighted) measurements per local neighbourhood, i.e.,

[TABLE]

where $\textrm{VAR}\{\cdot\}$ denotes the variance operator, $\mathcal{M}^{(i,j)}_{ss}$ denotes the set of pixels in a local neighbourhood of $(i,j)$ and $w_{i^{\prime},j^{\prime}}$ are the weights. The minimization is simplified by dividing the data in slices of $L/W$ contiguous bands and running the algorithm introduced in [20] per slice. As shown in Figure 6, the proposed strategy distributes the measurements uniformly in space, while other random strategies [16] tend to exhibit clusters, leaving some regions without measurements.

VI Experiments

To illustrate the efficacy of the proposed method, the new reconstruction algorithm is compared to other alternatives (based on the work conducted in [6]) using a synthetic dataset. Subsequently, the new subsampling scheme is compared with other random subsampling choices for a real MSL dataset. In all the experiments, the performance was measured using the following summary statistics:

• •

True detections $F_{\textrm{true}}(\tau)$: Probability of true detection as a function of the distance $\tau$, considering an estimated point as a true detection if there is another point in the ground truth point cloud in the same pixel ($x^{\textrm{true}}_{n}=x^{\textrm{est}}_{n^{\prime}}$ and $y^{\textrm{true}}_{n}=y^{\textrm{est}}_{n^{\prime}}$) such that $|t^{\textrm{true}}_{n}-t^{\textrm{est}}_{n^{\prime}}|\leq\tau$.

• •

False detections $F_{\textrm{false}}(\tau)$: Number of estimated points that cannot be assigned to a ground truth point at a distance $\tau$.

• •

Mean intensity absolute error at distance $\tau$ (IAE): Mean across all the points of the intensity absolute error $\sum_{\ell=1}^{L}|r^{\textrm{true}}_{n,\ell}-r^{\textrm{est}}_{n^{\prime},\ell}|$, normalized with respect to the total number of ground truth points. The ground truth and estimated points are coupled using the probability of detection $F_{\textrm{true}}(\tau)$. Note that if a point was falsely estimated or a ground truth point was not found, then they are considered to have resulted in an error of $\sum_{\ell=1}^{L}|r_{n,\ell}|$. The comparison is done with normalized intensity values, that is $\sum_{t=1}^{T}h_{\ell}(t)=1$ for $\ell=1,\dots,L$.

• •

Background normalized mean squared error $\textrm{NMSE}_{\boldsymbol{B}}$: Mean of the normalized squared error of the estimated background at each wavelength, i.e., $\frac{1}{L}\sum_{\ell=1}^{L}\frac{\sum_{i=1}^{N_{r}}\sum_{j=1}^{N_{r}}(b^{\textrm{true}}_{i,j,\ell}-\hat{b}_{i,j,\ell})^{2}}{\sum_{i=1}^{N_{r}}\sum_{j=1}^{N_{r}}{\left(b^{\textrm{true}}_{i,j,\ell}\right)}^{2}}$.

VI-A Synthetic data

We first assessed the performance of the proposed algorithm using a synthetic dataset created from the “Art” scene of the Middlebury dataset [50], as shown in Figure 7. The measurements were obtained by simulating the single-photon multispectral Lidar system of [17], whose bin width is 0.3 mm. The generated dataset has $N_{r}=283$ and $N_{c}=231$ pixels, $T=4500$ histogram bins and $L=4$ wavelengths (red, green, blue and yellow), where only $W=2$ wavelengths out of 4 were sampled per pixel using the coded aperture introduced in Section V. The mean number of photons per wavelength per pixel is 10, where approximately 3.4 photons are due to the background illumination. As in mono-static Lidar systems, the background levels are generated as a passive image of the scene (see Figure 9).

We compared the proposed method to a two-stage algorithm that first estimates the point positions using ManiPoP [15] (single-wavelength multiple-return state-of-the-art algorithm) by integrating the photons across wavelengths and then infers the spectral signatures with fixed point positions, similarly to the procedure suggested by Wallace et al. [6]. The resulting method, referred to as ManiPoP #1, is summarized in Algorithm 2. We also compared with ManiPoP in the strict single-wavelength setting, by choosing the most powerful wavelength and using the same total acquisition time per pixel than in the multispectral case (i.e., a per-pixel acquisition time $W=2$ longer than the one considered for each wavelength in MuSaPoP). This second alternative is referred to as ManiPoP #2.

Figure 8 illustrates $F_{\textrm{true}}(\tau)$, $F_{\textrm{false}}(\tau)$ and IAE for both methods. The proposed algorithm performs better than the other alternatives, as it finds 97.7% of the true points, whereas ManiPoP #1 and #2 only recover 95.34% and 89.6% respectively. ManiPoP #1 relies on an approximate impulse response $\tilde{h}(t)=\sum_{\ell=1}^{L}h_{\ell}(t)$, which biases the depth estimates, as the true accumulated response varies across points depending on their spectral signature555Note that this bias can be arbitrarily large depending on the variations of $h_{\ell}(t)$ across wavelengths.. The bias degrades the performance in terms of average depth absolute error (computed for true detections within a distance of 9 mm from the ground truth point). The proposed method obtains an average error of 3.9 mm, whereas the estimates by ManiPoP #1 present an average error of 5.7 mm. Despite having double acquisition time for the single-wavelength, ManiPoP #2 fails to find points corresponding to materials that have very low reflectivity in the blue wavelength (e.g., the red helicoidal structure shown in Figure 7). The proposed MuSaPoP algorithm performs slightly worse in terms of false detections, finding 3 times more false points than the competitor methods. In terms of intensity estimation, MuSaPoP obtains better results, having an asymptotic IAE of 1 photon, whereas the alternatives #1 and #2 provide IAE equal to 1.1 and 2.7. The estimated background levels are shown in Figure 9. The proposed method yields a better background NMSE (0.04) than alternatives #1 (0.14) and #2 (0.79). The improvement in background estimation over the ManiPoP alternatives can be attributed to the use of an empirical Bayes prior instead of a gamma Markov random field. The total execution time was 811 s for MuSaPoP and 294 and 348 s for alternatives #1 and #2.

VI-B Real MSL data

The proposed subsampling scheme was evaluated on a real MSL dataset [17]. The scene consists of $L=32$ wavelengths sampled at regular intervals of 10 nm from 500 nm to 810 nm, $N_{r}=N_{c}=198$ pixels and $T=4500$ histogram bins. The target is composed by a series of blocks of different types of clay and two leaves. Figure 10 shows an RGB image of the scene and the 3D reconstruction using acquisition times up to 10 ms per wavelength per pixel. We compare the blue noise codes mentioned in Section V with the random schemes introduced in [16], all yielding the same total number of measurements and acquisition time:

Random sampling without overlap: $W$ out of $L$ bands per pixel are sampled without replacement (i.e., for a given pixel, each wavelength is measured at most once). 2. 2.

Random sampling without overlap: For each wavelength, $W/L$% of the pixels are sampled without replacement. 3. 3.

Proposed sampling method: For each wavelength, $W/L$% of the pixels are chosen following the scheme presented in Section V.

The codes were evaluated for $W=1,2,4,8,16$ bands per pixel and acquisition times of 0.1, 1 and 10 ms per measurement (i.e., the histogram of one wavelength), using as ground-truth the reconstruction obtained with all the measurements and an acquisition time of 10 ms. Figure 11 shows the percentage of true detections, IAE and background NMSEs for all codes, acquisition times and numbers of sensed bands per pixel $W$. All the evaluated compressive strategies yield good results, where a small improvement can be obtained by the use of blue noise codes. In terms of total number of estimated points, the blue noise codes achieve better performance in high compression scenarios $W=1,2$ and low acquisition times (0.1 and 1 ms). For example, for an acquisition time of 1 ms, almost all points are reconstructed using blue noise codes, whereas the random codes only yield around $97\%$ of the ground-truth points. The choice of blue noise codes has a stronger impact in terms of IAE, achieving smaller IAE for all acquisition times and number of bands per pixel.

Figure 12 shows the execution time for acquisition times of 10, 1 and 0.1 ms and different numbers of sensed bands. The proposed RJ-MCMC sampler has a complexity proportional to the number of photon detections in the support of the impulse response around the 3D point being modified, whereas the background update has a complexity proportional to the total number of active histogram bins in the Lidar scene. The background extraction step required around $15\%$ of the total execution time, which could be significantly reduced if all the bands were processed in parallel instead of sequentially as it is done in the current implementation. All the experiments were performed using a Matlab R2018a implementation on a i7-3.0 GHz desktop computer (16GB RAM).

Finally, we compared the performance of MuSaPoP with the single-depth multiple-wavelength algorithm by Altmann et al. [16]. The algorithm is referred to as Depth TV and considers total variation regularizations for the background, reflectivity and depth images. Note that this method requires a (global) depth interval where all signal photons are found, which is given manually by the user. We also considered a target detection scenario (i.e., some pixels without surfaces), by removing the backplane of the scene and keeping only photons associated with background levels. In this case, we post-process the Depth TV estimates, removing points with a mean normalized intensity below 10%, which gave the best results across the evaluated datasets. In both experiments, we used the blue noise codes with $W=8$ wavelengths per pixel out of $L=32$. Figure 13 shows the 3D reconstructions obtained by MuSaPoP and Depth TV using an acquisition time of 10 ms. Figure 14 shows the performance of both algorithms in terms of true and false detections. MuSaPoP performs better in the 1 and 0.1 ms cases, whereas Depth TV obtains better depth estimates in the lowest acquisition time case (0.01 ms). However, in the 0.01 ms case without backplane, the intensity thresholding step does not remove backplane points, hence obtaining a very large number of false detections. This result illustrates the inefficiency of simple thresholding in target detection scenarios, whereas MuSaPoP includes these cases within its general formulation. Table II shows the performance of both algorithms in terms of IAE, background NMSE and execution time. The proposed method yields a better IAE than Depth TV (approximately half), as the latter tends to smooth out details within the blocks and leaves, as shown in Figure 13. Moreover, in terms of background NMSE, Depth TV fails to provide good estimates in the low-photon cases, as it only considers photon counts within the global interval without signal returns. The execution time of Depth TV was significantly higher than MuSaPoP.

VII Conclusions and future work

This paper has studied a new 3D reconstruction algorithm referred to as MuSaPoP using multispectral Lidar data, which is able to find multiple surfaces per pixel. The proposed method leads to better reconstruction quality than other alternatives, as it considers all measured wavelengths in a single observation model. While based on some ideas initially investigated in ManiPoP [15], MuSaPoP also relies on new strategies to deal with the very high dimensionality of the multispectral problem. The first novelty is the use of an empirical Bayes prior for the background levels, which speeds up significantly the RJ-MCMC algorithm. A second improvement is the adapted dilation/erosion and split/merge moves for the multispectral case, profiting from SBR estimates to increase the acceptance rate. Finally, the subsampling strategy further reduces both the algorithm’s complexity and total number of measurements, leading to faster acquisitions and reconstructions. The sparse point cloud representation of the proposed method speeds up the computations proportionally to the number of measurements, whereas other dense models [13, 14] would not achieve similar improvements.

Further work will be devoted to the design of compressive systems that are not limited to subsampling the spectral cube. Moreover, the compression can also be extended to depth information, for example using a range-gated camera [18].

Appendix A Gamma Markov random field

Most classical prior distributions (often referred to as regularization terms in the convex optimization literature) for images, such as Laplacian [26] and total variation [51], only penalize variations between neighbouring pixels, ignoring the mean intensity of the image. However, the gamma Markov random field [31], includes a penalization on the mean intensity, which promotes pixels with smaller values. This can be shown by inspecting the marginal distribution (defined as [32]),

[TABLE]

where $\alpha_{B}$ is a hyperparameter controlling the degree of smoothness and $\mathcal{M}_{B}$ denotes the set of pixels in the neighbourhood of pixel $(i,j)$. For an image of constant intensity $c$, we have $b_{i^{\prime},j^{\prime},\ell}=b_{i,j,\ell}=c$ for all pixels and spectral bands, yielding the density

[TABLE]

This dependency on the mean promotes reconstructions with lower background levels, decreasing the acceptance ratio of death and erosion moves (that propose to increase the background levels). In the case of ManiPoP, only one band is considered ($L=1$). Thus, the bias towards smaller background levels does not impact the overall reconstruction significantly. However, in the MSL case ($L\gg 1$), the reconstruction quality is reduced, hindering the use of gamma Markov random fields.

Appendix B Empirical prior for the background levels

The prior for the background levels is chosen to minimize the Kullback-Leibler divergence in (9), where the correlated model $q$ is given by (10). The minimization of (9) can be written as

[TABLE]

Considering the product of independent gamma distributions in (8), the problem reduces to

[TABLE]

for all pixels $i=1,\dots,N_{r}$ and $j=1,\dots,N_{c}$ and wavelengths $\ell=1,\dots,L$. The expected values $\mathbb{E}_{q}\{b_{i,j,\ell}\}$ and $\mathbb{E}_{q}\{\log b_{i,j,\ell}\}$ cannot be obtained in closed form for the Poisson-Gaussian model of (10). Thus, we approximate them numerically by obtaining MCMC samples of $\tilde{b}_{i,j,\ell}$. As explained in [32], off-the-shelf sampling strategies (e.g., Hamiltonian Monte Carlo [24, Chapter 9]) do not scale well with the dimension of the problem, being inefficient when applied to large multispectral Lidar datasets. Hence, we consider proposals from a Gaussian approximation of (10) (as detailed in [26]) using the perturbation optimization algorithm [52], accepting or rejecting them according to the Metropolis-Hastings rule [24, 26]. We generate $10^{3}$ samples $\{\tilde{b}^{(s)}_{i,j,\ell},s=1,\dots,10^{3}\}$ and compute the desired expected values as

[TABLE]

Finally, the values of the hyperparameters are obtained by setting $\theta_{i,j,\ell}=\mathbb{E}_{q}\{b_{i,j,\ell}\}$ and minimizing (27) with a one-dimensional Newton method.

Appendix C Complete expressions of the RJ-MCMC moves

The birth move of point $(\boldsymbol{c}_{N_{\Phi}+1},\boldsymbol{r}_{N_{\Phi}+1})$ has an acceptance ratio given by $\rho=\min\{1,r\left(\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\right)\}$ with

[TABLE]

where $C_{1}$ is defined as

[TABLE]

Similarly, the death move is accepted with probability $\rho=\min\{1,C_{1}^{-1}\}$, where the term $\frac{1}{N_{\Phi}+1}$ in the second line is replaced by $\frac{1}{N_{\Phi}}$. The dilation move of point $(\boldsymbol{c}_{N_{\Phi}+1},r_{N_{\Phi}+1})$ is accepted with probability $\rho=\min\{1,r\left(\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\right)\}$ with

[TABLE]

where $C_{2}$ is defined as

[TABLE]

A shift of the point $(\boldsymbol{c}_{n},\boldsymbol{r}_{n})$ to the new position $\boldsymbol{c}_{n}^{\prime}=[x_{n},y_{n},t^{\prime}_{n}]^{T}$ has an acceptance probability of $\rho=\min\{1,r\left(\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\right)\}$ with

[TABLE]

where

[TABLE]

A mark update randomly picks a point $(\boldsymbol{c}_{n},\boldsymbol{r}_{n})$ and proposes a new spectral signature $\boldsymbol{r}^{\prime}_{n}=\log(\boldsymbol{m}^{\prime}_{n})$. Each spectral log-intensity is accepted independently with probability $\rho=\min\{1,C_{4}\}$, where

[TABLE]

The split move from $(\boldsymbol{c}_{n}=[x_{n},y_{n},t_{n}]^{T},\boldsymbol{r}_{n})$ to $(\boldsymbol{c^{\prime}}_{k_{1}}=[x_{n},y_{n},t^{\prime}_{k_{1}}]^{T},\boldsymbol{r}^{\prime}_{k_{1}})$ and $(\boldsymbol{c^{\prime}}_{k_{2}}=[x_{n},y_{n},t^{\prime}_{k_{2}}]^{T},\boldsymbol{r}^{\prime}_{k_{2}})$ is accepted with probability $\rho=\min\{1,r\left(\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\right)\}$, where

[TABLE]

and

[TABLE]

Finally, the merge move is accepted with probability $\rho=\min\{1,C_{5}^{-1}\}$.

Acknowledgment

We would like to thank the single-photon group666http://www.single-photon.com led by Prof. G. S. Buller for providing us the real MSL dataset used in this paper.