LoGG3D-Net: Locally Guided Global Descriptor Learning
for 3D Place Recognition

The task of point cloud based retrieval for place recognition is generally formulated as follows. Given a point cloud $\mathcal{P}\in\mathbb{R}^{N\times 4}$ representing varying number of $N$ points with associated $x,y,z$ and intensity, a mapping function $\Phi:\mathcal{P}\rightarrow g\in\mathbb{R}^{d^{\prime}}$ is developed that represents the point cloud with a fixed-size global descriptor $g$. End-to-end learning of global descriptors is formulated as follows:

Problem: Given a training set $\{\mathcal{S}_{i}\}=\{(\mathcal{P}_{i},\mathtt{x}_{i})\}$ consisting of pairs of point cloud $\mathcal{P}_{i}$ with associated geo-location $\mathtt{x}_{i}$, find the parameters $\theta$ of the mapping function $\Phi_{\theta}:\mathcal{P}\rightarrow g\in\mathbb{R}^{d^{\prime}}$ , such that for any subset of samples $\{\mathcal{S}_{a},\mathcal{S}_{i},\mathcal{S}_{j}\}$,

$$\mathcal{D}(\mathtt{x_{a}},\mathtt{x_{i}})\leq\mathcal{D}(\mathtt{x_{a}},\mathtt{x_{j}})\implies\lVert g_{a}-g_{i}\rVert\leq\lVert g_{a}-g_{j}\rVert,$$

(1)

where, $\mathcal{D}$ represents geometric distance and $\lVert.\rVert$ represents a distance in the feature space (typically ${L}_{2}$).

Learning the parameters $\theta$ that address the above problem is generally performed using a metric learning setting by applying a loss function $\mathcal{L}_{g}$ that acts on the global descriptors extracted from a tuple $\mathcal{T}$ of training samples. We use the subscript ‘$g$’ in $\mathcal{L}_{g}$ to highlight that the loss is typically only applied on the global descriptors, i.e., at the scene-level.

We note that generally $\Phi$ can be decomposed into two functions such that,

$$\Phi\equiv\varphi\circ\phi,$$

(2)

where, $\phi:\mathcal{P}\rightarrow\{f\}$ extracts local features $f\in\mathbb{R}^{d}$, and $\varphi:\{f\}\rightarrow g$ aggregates the local features into a single global descriptor $g$. Under this setting, optimizing $\theta$ with a training signal applied only to the global descriptor seems limited, i.e., there are desirable properties of the intermediate representations $\{f\}$ (local features) that cannot be fully enforced by $\mathcal{L}_{g}$.

Hypothesis: Given two point clouds from nearby locations (i.e., with considerable overlap of points), enforcing consistency of local features of corresponding points will result in a more repeatable global descriptors after aggregation.

We define ‘corresponding points’ as two points from different point clouds that are nearby when represented with respect to a global coordinate frame. ‘Consistency of local features’ implies that local features of corresponding points should be nearby in the embedding space. Towards testing the above hypothesis, we introduce an additional training signal to enforce the features $\{f\}$ of corresponding points to be consistent in the embedding space.

For the local feature extractor $\phi$, we use a Sparse U-Net style backbone [22] which performs sparse point-voxel convolution. The backbone consists of two branches: a voxel-based branch and a point-based branch. The voxel-based branch learns local neighborhood information at varying receptive field sizes using convolutional layers, and the point-based branch learns high-resolution information that may be lost in the voxelized branch using MLP layers. The information across the two branches are fused at intermediate steps to capture complementary information from the input 3D point cloud.

Given a point cloud $\mathcal{P}=\{p_{i}\}$ and the set of point features $\{f(p_{i})\}$ where $f(p_{i})\in\mathbb{R}^{d}$, the second-order pooling $F^{O_{2}}$ of the features is defined as,

$$F^{O_{2}}=\{F^{O_{2}}_{xy}\},\quad F^{O_{2}}_{xy}=\max_{p_{i}\in\mathcal{P}}\ f_{xy}^{o_{2}}(p_{i}),$$

(4)

where $F^{O_{2}}$ is a matrix with elements $F^{O_{2}}_{xy}(1\leq x,y\leq d)$ and $f^{o_{2}}(p_{i})=f(p_{i})f(p_{i})^{T}\in\mathbb{R}^{d\times d}$ is the outer product of the point feature with itself. This accounts to taking the element-wise maximum of the second-order features of all points in the point cloud.

In order to make the scene descriptor matrix $F^{O_{2}}$ more discriminative, we use Eigen-value Power Normalization (ePN) [20, 19, 21]. Given the singular value decomposition $F^{O_{2}}=U\lambda V^{T}$ the ePN result $F^{O_{2}}_{\alpha}$ is obtained by raising each of the singular values by a power of $\alpha$ as follows,

$$F^{O_{2}}_{\alpha}=U\hat{\lambda}V^{T},\quad\hat{\lambda}=\textit{diag}(\lambda^{\alpha}_{1,1},..,\lambda^{\alpha}_{d,d}),$$

(5)

where, $\alpha=0.5$. The matrix $F^{O_{2}}_{\alpha}$ is flattened and normalized to obtain the final global descriptor vector $g\in\mathbb{R}^{d^{2}}$. To incorporate Eq. (5) into our end-to-end pipeline, we utilize differentiable SVD as introduced in [26], and its PyTorch implementation.

After the aggregation of point features into a global descriptor using second-order pooling, the scene-level loss $\mathcal{L}_{g}$ is applied to a tuple $\mathcal{T}$ of training samples. We use the quadruplet loss [27] where a tuple is denoted as $\mathcal{T}_{i}=(\mathcal{S}_{a},\{\mathcal{S}_{p}\},\{\mathcal{S}_{n}\},\mathcal{S}_{n^{*}})$, where $\mathcal{S}_{a}$ is the anchor sample, {$\mathcal{S}_{p}$} a set of positives (such that $\mathcal{D}(\mathtt{x_{a}},\mathtt{x_{p}})<\tau_{p}$), {$\mathcal{S}_{n}$} a set of negatives (such that $\mathcal{D}(\mathtt{x_{a}},\mathtt{x_{n}})>\tau_{n}$), and $\mathcal{S}_{n^{*}}$ is sampled such that it is not a positive to the query nor to all previous negatives (such that $\mathcal{D}(\mathtt{x_{n^{*}}},\mathtt{x})>\tau_{p}\,\forall\,\mathtt{x}=\{\mathtt{x_{a}},\mathtt{x_{n}}\}$).

In each tuple we first find the hardest positive sample,

$$\mathcal{P}_{hp}=\underset{\mathcal{P}_{p^{i}}\in\left\{\mathcal{P}_{p}\right\}}{\rm max}\,\,||g(\mathcal{P}_{a})-g(\mathcal{P}_{p^{i}})||_{2},$$

(6)

The quadruplet loss is then defined as:

$$\begin{split}\mathcal{L}_{g}&\!=\sum_{i}^{\mathcal{N}}\Bigg{\{}\left[||g(\mathcal{P}_{a})\!-\!g(\mathcal{P}_{hp})||_{2}^{2}\!-\!||g(\mathcal{P}_{a})\!-\!g(\mathcal{P}_{n^{i}})||_{2}^{2}\!+\!\alpha\right]_{+}\\ &\!+\!\left[||g(\mathcal{P}_{a})\!-\!g(\mathcal{P}_{hp})||_{2}^{2}\!-\!||g(\mathcal{P}_{n^{*}})\!-\!g(\mathcal{P}_{n^{i}})||_{2}^{2}\!+\!\beta\right]_{+}\Bigg{\}},\end{split}$$

(7)

where, $\alpha$ and $\beta$ are constant margins and $\mathcal{N}$ is the number of sampled negatives.

Our network is jointly optimized by a weighted sum of the global scene-level loss and the local consistency loss described as:

$$\mathcal{L}=\mathcal{L}_{g}+\omega\cdot\mathcal{L}_{lc}$$

(8)

where, $\omega$ is a scalar hyperparameter term.

The proposed network is implemented using the PyTorch framework and trained on 12 Nvidia Tesla P100-16GB GPUs using $\mathtt{DistributedDataParallel}$. The TorchSparse library [22] is used for sparse convolutions. During training, the ground plane is first removed using RANSAC plane fitting followed by down-sampling using a voxel grid filter of $10cm$. Finally, input point clouds are limited to a maximum of $35K$ points. To reduce overfitting, we apply the following data augmentations for training. Random point jitter is introduced using Gaussian noise sampled from $\mathcal{N}(\mu=0,\sigma=0.01)$ clipped at $0.03m$. Each point cloud is also randomly rotated about the $z$-axis by an angle between $\pm 180^{\circ}$. Note that ground plane removal is not used during evaluation to speed up inference time as it does not affect evaluation performance of our proposed method¹¹1The additional features from the ground plane are still well-formed (property in Fig. 3), and after aggregation, the global descriptors still remain discriminative due to the robustness of ePN Eq. (5) to feature burstiness [20]..

For a fair comparison with PointNetVLAD [7], we set the same global descriptor dimension ($d^{2}=256$), hence the dimension of the local features is set to $d=16$. In the local consistency loss $\mathcal{L}_{lc}$, the margins $m_{p}$ and $m_{n}$ are set to 0.1 and 2.0 respectively, and $\lambda{{}_{n}}$ is set to 0.5. The quadruplet loss margins are set to $\alpha=0.5$ and $\beta=0.3$. The distances for sampling positive and negative point cloud pairs are set to $\tau_{p}=3m$ and $\tau_{n}=20m$. For $\mathcal{L}_{g}$ we use 2 positives, 9 negatives and 1 other negative. We train our model using the Adam optimizer with an initial learning rate of 0.001 and a multi-step scheduler to drop the learning rate by a factor of 10 after 10 epochs and train until convergence for a maximum of 24 hours.

We evaluate the proposed method on two public LiDAR datasets (KITTI, MulRan), both of which were collected from a moving vehicle in multiple dynamic urban environments. Note that these datasets are collected using different LiDAR sensors (Velodyne HDL-64E in KITTI, Ouster OS1-64 in MulRan) and in different countries (Germany, Korea).

We compute the cosine similarity between the global descriptors of each query with a database of global descriptors of previously seen point clouds in each sequence. Previous entries adjacent to the query by less than $t_{r}$ time difference are excluded from the search to avoid matching to the same instance. For $t_{r}$ we use $30s$ and $90s$ for KITTI and MulRan, respectively. Methods are compared using the Precision-Recall curve and its scalar metric the maximum $F1$ score ($F1_{max}$). The 3m, 20m thresholds are used to classify true positives and false positives respectively, as done in [14].

We evaluate the effect of inclusion of the local consistency loss through an ablation study on selected sequences of the MulRan dataset. We train on sequences DCC1, Riverside1 and evaluate performance on DCC2 and Riverside2. Table I summarizes the $F1_{max}$ for each test sequence by varying the weight of the local consistency loss $\omega$ in Eq. (8). It is evident that the inclusion on the local consistency loss leads to an improvement in place recognition performance. $\omega=0.1$ leads to an improvement of $26\%$ mean $F1_{max}$ with respect to the baseline while $\omega=1.0$ leads to an improvement of $61\%$ leading to the best performance. All the following experiments are carried out with $\omega=1.0$.

A qualitative depiction of the effect of $\mathcal{L}_{lc}$ is shown in Fig. 3 between two point clouds extracted from nearby locations. The left half shows the alignment of the point cloud after ground plane removal to estimate point correspondences during training. The right half shows the two point clouds separately with each point colored based on the t-SNE embedding of the local features extracted using our pre-trained model (note that point correspondence information is not used during inference). The visualization clearly highlights that distinct regions in a single point cloud are in distinct regions in the feature space. Additionally, across point clouds, corresponding regions have similar point features implying that the local features extracted from our method are repeatable.

We compare LoGG3D-Net with the state-of-the-art handcrafted method ScanContext²²2https://github.com/irapkaist/scancontext [5], the recently proposed hybrid method Locus³³3https://github.com/csiro-robotics/locus [14], and the popular end-to-end method PointNetVLAD [7]. For ScanContext we use the python version of the code provided by the authors. We use the 20x60 descriptor size and use ring-key search to find the top-10 candidates for descriptor matching. For PointNetVLAD we use a PyTorch-based re-implementation of the original tensorflow implementation⁴⁴4https://github.com/mikacuy/pointnetvlad.

The results are summarized in Table II. On the KITTI dataset, Locus remains the highest performing method with a mean $F1_{max}$ score $1\%$ higher than ours. On the MulRan dataset we obtain the best mean $F1_{max}$ score which is $5\%$ higher than the next best performing method, i.e., ScanContext. It should be noted that the local feature extractor of Locus was trained on the KITTI sequences 05 and 06 thus leading to its extremely high performance on KITTI and relatively low performance on MulRan.

Precision-Recall plots for the sequences KITTI 02 and DCC 03 are depicted in Fig. 4. KITTI 02 contains several repetitive environments and a single revisit to an intersection in the opposite direction. DCC 03 contains long traversals of revisits from the reverse direction. All these scenarios prove challenging for most methods while the adverse effect on the performance of LoGG3D-Net is not as pronounced.

The feature dimensions of Locus and ScanContext give these methods an unfair advantage in terms of representation power. The ScanContext descriptor is 20$\times$60 (a total of 1200 floating point numbers) and the Locus descriptor is 4096 dimensions. Both LoGG3D-Net and PointNetVLAD have descriptor sizes of 256 which are compact and scale well to large databases essential for real-time robotic applications.

The computation time for pre-processing, description and querying is demonstrated in  Table III. All experiments in this section are run on a system with an 8 core Intel i7-9700 processor with 32GB RAM and a single Nvidia RTX 2080Ti GPU. Since the time for retrieval increases with the size of the database, we present the average retrieval time for all methods on MulRan DCC1 which consists of 5541 point clouds. ScanContext uses ring-key retrieval to find the top-10 candidates for full descriptor distance calculation.

The results show that Locus has very high pre-processing time due to ground-plane removal and a high description time (due to around 30-80 sequential forward passes through the segment feature extraction network which is not parallelized). ScanContext has very high retrieval time even after limiting the number of ring-key candidates to just 10. The need for cosine similarity computation on each column shifted variant of the descriptor increases the descriptors query time. PointNetVLAD is the most efficient method for descriptor extraction due to its light network architecture. However, we note that PointNetVLAD uses a considerable amount of pre-processing to first remove the ground plane and then iteratively downsample the point cloud to exactly 4096 points which makes its total time longer than ours.

Our method has the lowest pre-processing and retrieval times allowing it to run at real-time ($\sim{10Hz}$), enabling its integration into SLAM systems as a loop closure detection module.