3D Surface Reconstruction from Voxel-based Lidar Data

TL;DR

This paper presents a novel 3D surface reconstruction algorithm from voxel-based Lidar data using an adaptive TSDF approach, improving mesh accuracy and density for autonomous vehicle navigation.

Contribution

Introduces an adaptive neighborhood kernel based on Gaussian confidence for TSDF voxel representation, enhancing surface reconstruction quality.

Findings

Effective on synthetic and real datasets

Outperforms existing surface reconstruction methods

Balances mesh density and accuracy

Abstract

To achieve fully autonomous navigation, vehicles need to compute an accurate model of their direct surrounding. In this paper, a 3D surface reconstruction algorithm from heterogeneous density 3D data is presented. The proposed method is based on a TSDF voxel-based representation, where an adaptive neighborhood kernel sourced on a Gaussian confidence evaluation is introduced. This enables to keep a good trade-off between the density of the reconstructed mesh and its accuracy. Experimental evaluations carried on both synthetic (CARLA) and real (KITTI) 3D data show a good performance compared to a state of the art method used for surface reconstruction.

Tables

Table 1. TABLE I: Performance on synthetic data

		Average Error (m)			Hausdorff Distance (m)
Method		P$\to$GT	GT$\to$P	Sym	P$\to$GT	GT$\to$P	Sym
IMLS [18]		$0.37$	$0.09$	$0.23$	$4.54$	$8.32$	6.43
Ours AN+GC	${𝐤}_{\text{max}}=\mathrm{𝟓}$	$0.14$	$0.13$	0.14	$1.39$	$20.49$	$10.94$
	$k_{\text{max}}=3$	$0.09$	$0.26$	$0.17$	$0.69$	$30.84$	$15.77$
Ours AN	$k_{\text{max}}=5$	$0.30$	$0.12$	$0.21$	$1.69$	$20.36$	$11.02$
	$k_{\text{max}}=3$	$0.14$	$0.25$	$0.19$	$0.87$	$30.83$	$15.85$
Ours CN+GC	$k=5$	$0.15$	$0.16$	$0.15$	$1.38$	$20.50$	$10.94$
	$k=3$	$0.09$	$0.27$	$0.18$	$0.69$	$30.85$	$15.77$
	$k=1$	$0.03$	$3.44$	$1.73$	$0.24$	$65.28$	$32.76$
Ours CN	$k=5$	$0.30$	$0.14$	$0.22$	$1.69$	$20.36$	$11.03$
	$k=3$	$0.14$	$0.26$	$0.20$	$0.87$	$30.83$	$15.85$
	$k=1$	$0.03$	$3.44$	$1.73$	$0.24$	$65.28$	$32.76$

Equations

$$|\mathcal{H}|\ge N_{min},$$

$$\overrightarrow{n}=\{\begin{array}{ll}\overrightarrow{e_3}& \text{if e3⋅(sp−Hμ)>0},\\ -\overrightarrow{e_3}& \text{otherwise}.\end{array}$$

$$\mathcal{N}^{k}_{\text{PDF}}(w^{k}\;|\;\mathcal{H}^{k}_{\mu},\Sigma)\geq\tau\,;\,\,\Sigma={\scriptsize\left[{\begin{array}[]{cc}\lambda_{1}&0\\ 0&\lambda_{2}\end{array}}\right]}\,.$$

$$\text{TSDF}(v)=\overrightarrow{n^{k^{\prime }}}\cdot (v-{\mathcal{H}}_{\mu }^{k^{\prime }}).$$

Full text

I Introduction

A robust and accurate model of the environment is crucial for autonomous vehicles. In fact, imprecise representations of the vehicle’s surrounding may lead to unexpected situations that could endanger the passengers. Although many different geometrical representations have been proposed by both the robotics [1] and the graphics [2] communities, creating an accurate 3D model of the surroundings from mounted sensors remains a challenge.

In recent years, the constant evolution of Lidar sensors enables to obtain rich and accurate information, even in large and complex scenes. However, the size of the input data, along with the noise, misaligned scans and density variation, make the task a very challenging one. To overcome this, coarse 3D representations are often used and restrictive hypotheses are applied (i.e. planar assumptions for road detection). While these representations can be sufficient to perform path planning or obstacle avoidance, complex maneuvers might require a more accurate description of the scene’s geometry. This can be specially convenient for terrain mapping or off road navigation.

In this paper, we present an algorithm capable to perform a fine and accurate 3D surface reconstruction of the environment from depth sensors. From the statistics of the input point cloud sampled into a voxel grid, local approximations of the surface are performed using an adaptive neighborhood capable to cope with the heterogeneous density of the input scan. Then, a truncated signed distance field (TSDF) is estimated to obtain a continuous mesh that maintains a high level of detail and density in areas close to the vehicle. Our output mesh can be of special interest for both the robotics and the graphics community to perform different tasks, such as terrain traversability assessment or physical modeling.

II Related work

The 3D geometry of a scene described from range sensors can be represented in many different ways. The choice of the representation is often related to its final purpose: visualization, mapping, localization, terrain traversability, among many others.

Some methods use directly the 3D set of points received from the input sensor, which might be useful for visualization or obstacle detection tasks. However, the level of details of the representation depends on the amount of data used, that rapidly becomes prohibitive for outdoor scenes. Conversely, other works propose a regularly sampled grid as introduced in [3], where occupancy information is stored into each cell. This enables to handle big amounts of data more efficiently and reduce memory needs, specially by using recursive structures such as octrees [4, 5]. These approaches have become widely used for terrain traversability assessment, mapping and visualization. However, their discrete nature does not enable a continuous representation, which might be desirable for other tasks such as physical modeling.

Alternatively, the graphics community has explored different methods to create a triangular mesh from the 3D points of the scanned surface. These methods have a wide range of applications on completely different fields as described in [2]. For simplicity, we distinguish between explicit and implicit methods.

Explicit methods: These techniques are either parametric or triangulated with respect to how they represent the reconstructed surface. For some applications, the continuity of the reconstruction is not required. Some researchers propose to fit local primitive shapes to the input point cloud and obtain a segmented mesh with simple primitives for the surface description [6, 7]. These approaches are not suitable for reconstructing complex scenes.

Other methods create local descriptions of the surface by a set of unoriented discs (-surfels-) calculated from the points distribution inside defined neighborhoods. Surfels have been used for traversability assessment [8], and more recently to perform simultaneous localization and mapping (SLAM) in outdoor urban environments [9].

Implicit methods: For other applications such as physical modeling or detailed terrain traversability, a more accurate and continuous representation might be preferred. To obtain this, some approaches commonly define a Truncated Signed Distance Field (TSDF) to represent the surface implicitly by a gradient field [2]. This implicit representation needs to be post-processed by meshification algorithms (e.g. Marching Cubes [10]) to obtain the final mesh.

In [11], range information across different viewpoints are integrated to average a TSDF from where the scalar field is obtained. By using this technique, 3D modeling has been performed in both small indoor [12] and large scale outdoor scenes [13, 14]. These methods typically require a large number of viewpoints to output a dense reconstruction and are susceptible to outliers.

Other implicit methods perform local approximations of the surface from where the TSDF is obtained. Poisson reconstruction [15, 16] is a well-known technique for creating watertight surfaces from oriented point samples. Point-based methods such as Hoppe’s [17] or Implicit Moving Least Squares (IMLS) [18, 19, 20] locally fit the data to a lower degree polynomial by using a projection operator. For such methods, variable density inputs might represent a challenge since the data is treated equally along the complete space.

In contrast to most methods presented above, that are unable to accommodate to the heterogeneous density of the input data while keeping the accuracy of the reconstruction, we introduce a method that is able to handle this variable density while keeping a good trade-off between the accuracy and the density of the outputted surface.

III Methodology

Let us consider an input point cloud $\mathcal{P}$ obtained from any range sensor at a known viewpoint and sampled by a voxel grid, our aim is to reconstruct the underlying 3D mesh. From the statistical distributions of the points in each voxel (Sec. III-A), we first approximate local planar surfaces (Sec. III-B), and then compute the implicit surface representation (Sec. III-C) that encodes the distances to the local planar surfaces. To accommodate to the input data heterogeneous density, as well as to gain robustness to noise, we use an adaptive neighborhood kernel, resulting in a denser and smoother reconstruction. The overall overview of our methodology is shown in Fig. 2.

III-A Voxel representation

We benefit of the work of [21] to update efficiently a regular voxel-wise representation of the point cloud $\mathcal{P}~{}\rightarrow~{}\{V^{1},V^{2},...,V^{n}\}$, with voxel size $\alpha$. In addition to the number of points $|V|$, each voxel $V$ stores the 3D statistical distribution of the points lying inside, that is: the mean $V_{\mu}$ and the covariance $V_{\sigma}$. This enables a rich compact representation of the points inside each voxel, while being significantly lighter than storing all the points. The statistical distribution is computed incrementally upon insertion of new points.

III-B Explicit local surfaces

It has been shown in [8, 9], that complex environments are efficiently approximated as a set of primitive local surfaces. This is especially suitable for mobile robotics, as robots usually evolve in well structured environments. Following this observation, we compute local planar surfaces using the 3D statistical distribution described in the previous section. While a naive implementation would fit planar surfaces to each voxel independently, this would inherently lead to a noisy reconstruction since some voxels may have very few points, if not none, due to the heterogeneous density of Lidar data.

In our pipeline we propose using an adaptive neighborhood definition where local surfaces are estimated from multi-scale voxels statistics. Our neighborhood definition presents two main advantages: a) it increases the statistical robustness which improves large planar surface estimation (e.g. ground, walls), b) it counterbalances the lack of local data due to low density or occlusion.

III-B1 Neighborhood definition

We define the multi-scale neighborhood at the vertices location rather than voxels, since implicit surfaces will later be estimated at each vertex of the voxel representation. Let’s consider $v$ a vertex from the voxel-grid representation, its 8 adjacent voxels make up the first neighborhood level denoted $\mathcal{H}^{1}(v)$. Subsequently, the union of $\mathcal{H}^{1}(v)$ and the voxels adjacent to $\mathcal{H}^{1}(v)$ make up the neighborhood $\mathcal{H}^{2}(v)$. More formally, the neighborhood $\mathcal{H}^{k}(v)$ is composed of the $(2k)^{3}$ nearest voxels surrounding $v$. A two-dimensional illustration of $\mathcal{H}^{k}(v)$ at levels $k=1$ and $k=2$ is presented in Fig. 5.

Following [8], for a given neighborhood $\mathcal{H}$ (indices are dropped for clarity) we obtain the cardinal $|\mathcal{H}|$, statistical mean $\mathcal{H}_{\mu}$, and covariance $\mathcal{H}_{\sigma}$ from the merging of statistical data of all voxels in $\mathcal{H}$.

III-B2 Planar estimation

Having obtained the local statistics, we use the covariance $\mathcal{H}_{\sigma}$ to estimate the local planar surface of $\mathcal{H}$ through a Principal Component Analysis (PCA) if the following equation is satisfied:

[TABLE]

where $N_{min}$ is a hyper parameter. Suppose ($\overrightarrow{e_{1}},\overrightarrow{e_{2}},\overrightarrow{e_{3}}$) the eigen vectors, and ($\lambda_{1},\lambda_{2},\lambda_{3}$) the eigen values with $\lambda_{1}~{}\geq~{}\lambda_{2}\geq~{}\lambda_{3}$, we define $\overrightarrow{e_{3}}$ (the least dominant axis) as the unoriented normal of the planar estimation of the surface at neighborhood $\mathcal{H}$. Since normals need to be consistently oriented, the normal $\overrightarrow{n}$ of the plane is oriented towards the sensor pose $s_{p}$ as follows:

[TABLE]

where $\cdot$ stands for the dot product. We denote $\Pi$ the plane formed by the pair of normal and center $(\overrightarrow{n},\mathcal{H}_{\mu})$. An example of the local planar surfaces estimation is visible in Fig. 3.

III-C Implicit global surface

To reconstruct the global continuous surface, we compute the Truncated Signed Distance Field (TSDF) for each vertex $v\subseteq V$ close to the scanned surface.

To cope with the varying density of points in the point cloud, we first compute an optimal neighborhood level $k^{\prime}$ at each vertex, and then estimate the TSDF value given this optimal neighborhood.

III-C1 Adaptive neighborhood

A naive implementation of our neighborhood definition would use a constant $k$ throughout the scene. However, large $k$ values will over smooth high density regions, while small values of $k$ will lead to noisy estimation in low density regions. Instead, we compute the optimal neighborhood level $k^{\prime}$ for each vertex $v$, given the probability of $v$ to belong to a multivariate Gaussian distribution projected on $\Pi^{k^{\prime}}$ (the neighborhood planar estimation).

For each vertex, we calculate the projection $w^{k}$ of $v$ onto the plane $\Pi^{k}(v)$ and evaluate the likelihood of this projection to belong to the Gaussian $\mathcal{N}^{k}$, where $\mathcal{N}^{k}$ is the 2D planar-Gaussian of the statistical distribution of $\mathcal{H}^{k}$ projected onto $\Pi^{k}$, as illustrated in Fig. 4. The optimal neighborhood level $k^{\prime}$ is defined as the smallest level for which the projection of $v$ onto $\Pi^{k}$ has a probability to belong to the Gaussian $\mathcal{N}^{k}$ greater than $\tau$ (a hyper parameter). Formally, it is the smallest integer for which the probability density function $\mathcal{N}_{\text{PDF}}$ satisfies

[TABLE]

In practice, we compute the optimal $k^{\prime}$ iteratively starting at level 1 and stops when the above equation is satisfied as it is shown in Fig. 5. To avoid exponential computation time, $k$ is bounded between $[1,k_{max}]$.

III-C2 TSDF computation

We now compute the TSDF value for each vertex of the voxel grid, from the optimal local neighborhood level, while accounting for the normal estimation at the corresponding level. In other words, we compute $\text{TSDF}(v)$ such that:

[TABLE]

Our adaptive neighborhood level selection can efficiently handle varying local density and fill gaps between missing data, such as those between two adjacent Lidar layers. Subsequently, the condition on the probability density function given in Eq. 3 avoids the extension of the surface at its borders which is visible in qualitative results.

After TSDF computation, we use the popular marching cubes [10] to extract the mesh from the zero-crossing level of the gradient field. Fig. 3 shows the reconstructed mesh with constant neighborhood ($k=1$).

IV Experimental results

We evaluate our proposal on synthetic data from CARLA [22] and real data from KITTI [23], and compare our performance against the IMLS baseline [18]. For both synthetic and real data we used 100 frames equi-sampled from either an urban-like sequence for CARLA, or the residential sequence 0018 for the public KITTI dataset.

IV-A Methodology

Unless stated otherwise, our hyper-parameters remain unchanged in all experiments, with: $\alpha=0.2m$, $\tau=0.2$, $N_{min}=10$, and $k_{\text{max}}=5$. For fair comparison, the spherical neighborhood radius of IMLS is set to $\alpha\times k_{\text{max}}=1m$, and its k-nearest neighbor search to $N_{min}$. The noise parameter of IMLS is set to $h=\nicefrac{{1}}{{3}}\,(\alpha\times k_{\text{max}})=0.33$.

Noteworthy, it is impractical to have a real mesh ground-truth and KITTI can only be used for qualitative evaluation. Conversely, we use synthetic data from CARLA to evaluate the benefit of each of our contributions. The setup in CARLA replicates a top-mounted Velodyne HDL64E, similarly to the real KITTI setup. We also simulate a collocated noise-free Lidar with abnormally high resolution (316 layers), which serves as ground-truth for the reconstruction. Consequently, we frame the evaluation as a set-to-set distance problem, and measure the quality as the distance of the predicted mesh vertices (P) to the ground-truth (GT) set.

We use two metrics derived from the literature: the Average Error (AE), and the Haussdorf Distance (HD). The former measures the average distance error from one point to its nearest point in the other set, that is: $\displaystyle{}AE_{P\rightarrow{}GT}=\sum_{a\in{}P}\frac{1}{|P|}\min_{b\in{}GT}|a-b|$. The Haussdorf distance [24] is classically used for point set distances and gives a sense of the largest minimum error, that is: $\displaystyle{}\text{HD}_{P\rightarrow{}GT}=\max_{a\in{}P}\min_{b\in{}GT}|a-b|$. As each of the two metrics are directed, we also report the symmetrical metric, as the average of both directed metrics. For Haussdorf, $\text{HD}_{\text{sym}}=0.5(\text{HD}_{P\rightarrow{}GT}+\text{HD}_{GT\rightarrow{}P})$. We chose not to use the Chamfer distance used in machine learning [25] because it isn’t a metric-scale and is thus harder to interpret intuitively.

IV-B Ablation study

We assess the average performance on 100 frames from CARLA and subsequently report quantitative performance in Table I and qualitative results in Fig. 6. In the former, our pipeline is compared to IMLS [18].

To evaluate the importance of our contributions, we compare the benefits of our Adaptive Neighborhood (AN, Sec. III-C1) with or without the Gaussian Confidence (GC, Eq. 3) and also with a naive Constant Neighborhood (CN) approach. More precisely, the neighborhood $\mathcal{H}^{k^{\prime}}$ considered for the plane estimation (cf. Sec III-B2) is different:

AN+GC: $\mathcal{H}^{k^{\prime}}$ is the minimum neighborhood among $\left\{\mathcal{H}^{1},\cdots,\mathcal{H}^{k_{max}}\right\}$ that satisfies Eq. 1 and 3.

AN: $\mathcal{H}^{k^{\prime}}$ is the minimum neighborhood among $\left\{\mathcal{H}^{1},\cdots,\mathcal{H}^{k_{max}}\right\}$ that satisfies Eq. 1.

CN+GC: $\mathcal{H}^{k^{\prime}}$ with $k^{\prime}=k$, is considered only if Eq. 1 and 3 are satisfied.

CN: $\mathcal{H}^{k^{\prime}}$ with $k^{\prime}=k$, is considered only if Eq. 1 is satisfied.

From Table I, the accuracy on the reconstruction directly affects the P$\rightarrow$GT distances while the GT$\rightarrow$P distances are mostly influenced by the reconstruction density.

As expected, the benefit of our adaptive neighborhood strategy is noticeable when comparing AN+GC and CN+GC. Not only AN+GC exhibits higher accuracy (lower P$\rightarrow$GT) but it also increases the reconstruction density (lower GT$\rightarrow$P). When using our Gaussian Confidence (GC), there is a slight density loss (lower GT$\rightarrow$P) but a significant accuracy increase. Qualitatively from Fig. 6, our adaptive neighborhood helps to maintain the details of the surface by not over-smoothing the data, while the Gaussian Confidence strategy avoids to extend the surface, keeping its structure and generating fewer artifacts.

Overall, we observe that our complete pipeline (AN+GC with $k_{\text{max}=5}$) keeps a good trade-off between accuracy and density, which is shown by the best result obtained $AE_{\text{sym}}$ ($0.14m$ vs. $0.23m$ for IMLS). With larger neighborhood (bigger $k_{\text{max}}$), the density of the reconstruction increases but at the expense of a lower accuracy which is intuitive as there is a need to extrapolate more the data. Since IMLS performs a denser reconstruction, lower GT$\rightarrow$P distances are obtained, which explains the best results on the symmetric Hausdorff Distance ($10.94m$ vs. $6.43m$ for IMLS). While our method is less dense, our predicted mesh is significantly more precise ($1.39m$ vs. $4.54m$ for IMLS).

IV-C Synthetic data

To further evaluate the performance, we report the average error as a function of the sensor distance in Fig.8a, and the cumulative delta error that indicates the percentage of vertices of the output meshes having an error lower than a given value in Fig. 8b.

For both metrics, the accuracy of our reconstruction is significantly higher than the one obtained by IMLS, as the average error of our reconstruction is often 50% lower (Fig. 8a). The percentage of vertices below a given error exhibits the same behavior (Fig. 8b), with a significant advantage for our method. Furthermore, almost $80\%$ of the vertices of our mesh have an error lower than $0.2m$, while $40\%$ of vertices lie below the same threshold with IMLS.

Finally, a qualitative comparison on CARLA is shown in Fig. 7. Again, one can see that IMLS outputs a more dense reconstruction of the scene but also that the method tends to extend all surfaces, generating artifacts and inaccuracies in the reconstruction. This can be observed in the circled areas where lampposts, traffic signs and walls are abnormally enlarged by IMLS. On the other hand, our method is able to keep these details on the reconstruction and maintains a high density in areas close to the vehicle, which confirms the quantitative results obtained in Fig. 8.

IV-D Real data

Qualitative results on KITTI residential sequence 0018 are shown in Fig. 9. As for the results presented on the synthetic data, IMLS outputs a denser reconstruction but extends all surfaces out of the borders and generates artifacts. This can be observed at the circled areas on the rightmost images. Conversely, our method outputs a more accurate reconstruction, keeping the structure of the scanned surface and completing the surface at no data points areas close to the sensor. Even though our method outputs a less dense reconstruction, the density is still adequate for further robotics tasks such as trajectory planning or obstacle avoidance. Further quantitative results might be obtained by randomly sub-sampling the input data and calculating the reconstruction error over the non sub-sampled mesh.

V Conclusion

In this paper we presented a surface reconstruction method based on a light voxel representation that benefits of statistical information. Our pipeline uses an adaptive neighborhood strategy coupled with a Gaussian confidence estimation, to best estimate local surfaces and the implicit surface estimation (TSDFs).

The performance of our method was demonstrated on both synthetic (CARLA) and real data (KITTI), significantly outperforming the classical IMLS though sparser reconstruction. While the approach is computationally too expensive for online applications, it can be useful for terrain mapping or off road terrain analysis.

For mobile robots, an accurate representation of the environment is a prerequisite. Our method contributes to improve the reconstruction accuracy, which can be useful for mapping and localization applications.

While our results show that local planar surfaces are sufficient for accurate reconstruction, we intend to extend this research to a more complex surface definition, for example with polynomials.