# Automatic Dialect Density Estimation for African American English

###### Abstract

In this paper, we explore automatic prediction of dialect density of the
African American English (AAE) dialect, where dialect density is defined as
the percentage of words in an utterance that contain characteristics of the
non-standard dialect. We investigate several acoustic and language modeling
features, including the commonly used X-vector representation and ComParE
feature set, in addition to information extracted from ASR transcripts of the
audio files and prosodic information. To address issues of limited labeled
data, we use a weakly supervised model to project prosodic and X-vector
features into low-dimensional task-relevant representations. An XGBoost model
is then used to predict the speaker’s dialect density from these features and
show which are most significant during inference. We evaluate the utility of
these features both alone and in combination for the given task. This work,
which does not rely on hand-labeled transcripts, is performed on audio
segments from the CORAAL database. We show a significant correlation between
our predicted and ground truth dialect density measures for AAE speech in this
database and propose this work as a tool for explaining and mitigating bias in
speech technology.

Index Terms: Dialect Density, Dialect Identification, African American
English, Fairness in ASR, Low-Resource Tasks

##  1 Introduction

This paper111This work is supported in part by the NSF. presents a novel
framework for automatically assessing a speaker’s dialect density. Dialect
density is the degree to which speech includes aspects of a dialect separate
from the mainstream dialect of the language [1]. A commonly used dialect
density measure is the number of non-standard dialect tokens that appear in an
utterance divided by the number of words in the utterance. Previous literature
[2] shows that high-dialect speakers often face more socioeconomic
discrimination and bias in early education than their low-dialect
counterparts. The work in [3] demonstrates that the performance of ASR systems
trained to recognize General American English (GAE) degrades when performing
recognition on high dialect density African American English (AAE). As
dialectal changes are composed of both phonological and morphosyntactic
aspects [4], dialect-heavy speech can present challenges to both the acoustic
model and the language model (LM) of an ASR system. There is currently a lack
of training data needed to train ASR systems to recognize AAE speech with the
same fidelity they display for GAE speech, and no investigation on low-
resource methods has been done in this area before. As a result, bias against
speakers of AAE and other dialects can propagate through spoken language
systems, creating inequitable outcomes in the technology.

To alleviate this issue, one could create an ASR system that first predicts
the dialect density of its user given a small amount of speech and
subsequently adjusts its model or model hyperparameters to best recognize
speech of the target speaker’s dialect. This strategy has been applied with
wide success to accented speech recognition [5, 6, 7]. However, several
differences exist between accent and dialect identification that make it
difficult to directly apply previous work to dialect density prediction:
First, accents are generally considered to be a binary characteristic. That
is, each speaker either does or does not have a given accent. Dialect density
is considered to be a continuum in which some speakers of a dialect possess
more dialectal attributes than others. As a result, previous work in accent-
robust ASR typically only considers accent identification as a binary
classification problem and does not lend itself well to identifying the degree
of difference from mainstream speakers. Second, most accent identification
systems do not explicitly consider grammar and within-language diction as
relevant to their task. However, dialects can be composed of a variety of
changes in pronunciation, grammar, and diction which make them difficult to
characterize by looking at any one feature of speech individually. In
addition, speakers may also consciously change the amount of dialect they use
depending on social context (translanguaging or code switching) [8], making it
difficult to assign a consistent dialect density measure to a long portion of
a speaker’s speech as one would with a type of accent.

For these reasons, new methods are needed in order for machines to accurately
infer dialect density. Although a primary use of this technology is in
dialect-robust ASR, dialect density estimation could also be beneficial in
data mining speech of a particular dialect, ensuring fair training of speech-
language technologies, and evaluating human bias against speakers of dialects.
In this work, we propose a new method of automatic dialect density estimation
for African American English from short utterances in a low-resource task. We
first analyze the correlation of several commonly-used and novel speech
features in predicting dialect density (Sec. 2). We then utilize a combination
of the most effective of those features in order to create a model that best
estimates the dialect density of a target speaker (Sec. 3). We follow our
results with a discussion (Sec. 4) and conclusions (Sec. 5) on the usefulness
of the features in predicting dialect density.

##  2 Methods

###  2.1 Data

We use speech samples from the Corpus of Regional African American Language
(CORAAL) [9]. This database contains spoken interactions between an
interviewer and an interviewee who speaks a regional variant of AAE. The set
of speakers range in age from approximately 20 to 80 years old and contains
roughly equal numbers of male and female identifying participants. The
original sampling frequency is 44.1kHz. We down-sampled the audio files to
16kHz in all experiments. We used the following numbers of speakers with
regional AAE dialects from five US cities: 22 speakers from Washington DC
(DCB), 10 speakers from Princeville, NC (PRV), 11 speakers from Rochester, NY
(ROC), 10 speakers from Lower East Side Manhattan, NY (LES), and 12 speakers
from Valdosta, GA (VLD). For each speaker, several utterances with good audio
quality ranging in length from 5sec to 1min were selected, and their dialect
densities were scored by hand as ground truths. The subset of speakers from
DCB, PRV, and ROC were chosen because their dialect densities were annotated
by the authors of [3]. All speakers from VLD and LES were used, and their
dialect densities were scored by the authors of this paper. This results in a
total of approximately 3 hours of dialect density-scored utterances from 65
speakers. This data set was split into speaker independent sets for which 70%
was used in training, 15% was used in validation, and 15% was used in testing.

###  2.2 Feature Extraction

From each utterance, we extracted the following six feature sets:

Wav2Vec2.0 Transcripts: For the first three feature sets, we generated ASR
transcripts using a pretrained Wav2Vec2.0 model [10] trained on the 960hr
LibriSpeech database [11]. While these transcripts contained errors and
misrepresented the out-of vocabulary (OOV) words, we implicitly attempted to
utilize consistent errors and accurate portions of the transcripts to identify
useful phonetic and grammatical information.

1\. ASR Output Character Combination Frequency: The frequency of each sequence
of two characters (bigram) in the transcript was counted and used as a
feature. The Wav2Vec2.0 model can output 31 different characters leading this
feature to be a 961 x 1 vector which can be thought of as the flattened 31 x
31 matrix in which the element in row i and column j is the number of times
character i was followed directly by character j in the generated transcript
for the given utterance. We hypothesize that this feature will capture
consonant clusters that commonly occur in a particular dialect.

2\. ASR Output Character Duration: From the output logits of the Wav2Vec2.0
model, the average duration of each output character was computed. We
hypothesize that this feature will be useful in determining which sounds are
more or less frequently spoken or stressed by speakers of a particular
dialect.

3\. ASR Output Language Modeling: In addition to the previously mentioned
features, we were interested in how neural language modeling techniques could
be applied to automatically generated transcripts of speech in order to
predict dialect density. We noticed that, of the most commonly noted features
of AAE [1, 3], language differences relating to the tense, collocation, and
negation of verbs (eg. absence of copula, negative concord, generalization of
“is” and “was” to use with plural and second person subjects, etc.) were
especially prevalent. This led us to pay particular attention to verbs. First,
the verbs in each utterance were found using a pre-trained FLAIR part-of-
speech tagging model [12, 13]. We then used the Fisher corpus [14], consisting
largely of GAE conversations, to train word- and character-based LSTM language
models, which provide probability distributions over the next word or
character in an utterance given the history. To measure mismatch of verbs in
the GAE training data and AAE testing data, we then extracted the verb OOV
rate (using the word-based model vocabulary) and the average verb surprisal
[15] (using the character-based model) for each utterance, where the surprisal
of the $i$-th word ($S(w_{i})$) is calculated from the letter LM as:

$\displaystyle S(w_{i})$ $\displaystyle=$
$\displaystyle-\log(p(w_{i}|w_{i-1},\ldots,w_{1}))$ $\displaystyle=$
$\displaystyle-\sum_{j=1}^{l_{i}}\log(p(c_{i,j}|c_{i,j-1},\ldots,c_{1,1}))$

where $l_{i}$ is the length of the word $w_{i}$ with characters $c_{i},_{j}$.
We also calculate the overall utterance perplexity from the character-based
model (char_ppl), the average surprisal for all words, and the ratio of
average verb surprisal to average overall surprisal. Since the LM is trained
on GAE, word choices more characteristic of AAE will have high surprisal.

4\. CompareE16 Features: The widely used ComParE16 features [16] were
extracted from the audio segments using the OpenSmile toolkit [17]. This set
includes pitch, energy, spectral, cepstral coefficients (MFCCs) and voicing
related frame-level features which are referred to as low-level descriptors
(LLDs). It also includes the zero crossing rate, jitter, shimmer, the
harmonic-to-noise ratio (HNR), spectral harmonicity and psychoacoustic
spectral sharpness. In total, this feature set contains 6373 features
resulting from the computation of various statistics, polynomial regression
coefficients, and transformations calculated over the low-level descriptor
contours.

Weak Supervision: To create the following two feature sets, we employed a
weakly-supervised learning technique. We noticed that the five cities used
from the CORAAL database have widely varying average dialect densities, with
the averages from PRV and VLD being much higher than those from ROC and LES,
and with the DCB average in between. Therefore, we believed that an
utterance’s city of origin could serve as a weak label in a preliminary step
before dialect density estimation. We gathered the set of utterances from the
entirety of the 200hr CORAAL database from the five cities of interest that
matched the following criteria: 1) Contained at least 10 words to have enough
speech to estimate dialect density, 2) Contained no interruptions from the
interviewer 3) Were not contained in the set of dialect density-scored
utterances. We then used shallow neural networks to map larger input feature
vectors into 5-dimensional vectors for which the ith element represents the
probability that the utterance was spoken by a speaker from the ith city in
the database. This step is intended to project larger sources of information
into smaller features vectors which contain only relevant dialect information.
The idea is that training a model to classify diverse utterances by region
would prompt it to learn region-specific information such as dialectal traits
without the need to label the dialect density of all of the utterances in the
training set. The output 5-dimensional vector is then used as the
representative feature.

5\. X-Vector: The popular X-vector was incorporated to capture speaker-
specific information [18]. These 512-dimensional neural network-generated
embeddings contain speaker-specific information that may relate to dialect. As
described earlier, the 512-dimensional vectors were projected into
5-dimensional feature vectors using the fully connected network shown in
Figure 1. This network achieved a validation accuracy of 72.6%.

6\. Prosodic Embedding: Inspired by [19], four pitch and energy features were
extracted across time from the utterances: F0 (extracted with Praat [20]), the
total energy in the frame, the energy in the spectrum below 1kHz, and the
energy in the spectrum above 1kHz. These features were then normalized and
used as the input to a CNN (as shown in Figure 1) that was trained to predict
the region of origin of the speaker. This forces the CNN to classify region
specific information from only the prosodic information contained in the
speaker’s changes in pitch and energy. This CNN achieved a validation accuracy
of 70.7%. The output probability vector was then used as the final prosodic
embedding.

Figure 1: The architecture for the fully connected (FC) network used to
project the X-vectors (left) and the CNN used to project the prosodic
information (right). The inputs to the CNN are the pitch (F0) and three energy
contours of the utterance. The output of both networks is a vector whose
elements represent the probability of the speaker belonging to each of the
cities used from the CORAAL database. Figure 2: Overview of the features used
in the proposed work.

###  2.3 Model Training

First, one distinct XGBoost model [21] was trained for each of the six feature
sets. This boosted decision tree model has the advantage of allowing us to
easily measure the impact of the input features on the output value for
explainability. Each of the six models was trained to predict dialect density
scores from one of the given input feature sets. Then the correlation between
the predicted dialect density labels and actual dialect density labels was
calculated. We chose correlation as the performance metric because human-
performed dialect density assessments are subject to possibly high inter-rater
variability within the ranges of their scores [22], and so evaluation methods
that rely heavily on the absolute value of the dialect density may be subject
to measurement noise. However, raters do tend to assign higher or lower scores
to the same speakers, and so we expect correlation between predicted and
ground truth scores to be meaningful. As some features may only correlate with
phonological aspects or only correlate with morphosyntactic aspects of dialect
density, we train each model to predict each of the three types of dialect
density scores:

$DDMphon=\frac{N_{ph}}{N}\qquad DDMgram=\frac{N_{ms}}{N}$
$DDM=\frac{N_{ph}+N_{ms}}{N}$

where $N_{ph}$ is the number of phonological AAE tokens in the utterance,
$N_{ms}$ is the number of morphosyntactic AAE tokens in the utterance, and N
is the total number of words in the utterance. The average dialect density for
each city used is shown in Table 1.

| DDMphon | DDMgram | DDM
---|---|---|---
DCB | 0.083 | 0.004 | 0.088
ROC | 0.041 | 0.006 | 0.047
PRV | 0.166 | 0.028 | 0.194
LES | 0.018 | 0.025 | 0.042
VLD | 0.122 | 0.029 | 0.141
Table 1:  Average dialect density by city for each of the dialect density
measures shown.

Finally, we used the set of all features as the input to the XGBoost model, as
shown in Figure 2. As the ComParE feature set was large, only the most
impactful 10 ComParE features were used in the combined feature set.

##  3 Results

Table 2 gives the Pearson Correlation of the predicted dialect density measure
with the ground truth labels for the test set for an XGBoost model trained on
the listed feature sets. We also include the SHAP value plots [23] which give
the relative importance of each feature to the model during prediction.
Figures 2 and 3 give the SHAP value plots for the models trained on all
features for predicting DDMphon and DDMgram, respectively. As the DDMphon term
dominiates the total dialect density measure, the SHAP value plot for DDM is
nearly identical to that of DDMphon. In these plots, the Wav2Vec2.0 Char Comb
features are listed as $char1\\_char2$ (e.g. N_space is the frequency of the
“N” character being followed by a space character), and the Wav2Vec2.0 Char
Dur features are listed as the character whose duration was used as the input
feature from the letters A-Z, period, apostrophe, space, or silence (sil)
characters.

In order to demonstrate the reliability of our results, we also perform random
hold out on the highest performing features. Here, we randomly select speaker-
independent train and test split (80% train, 20% test) from the data 200 times
and report the average scores over all runs in Table 3.

Correlation | DDMphon | DDMgram | DDM
---|---|---|---
| Wav2Vec2.0
---
Char Dur.
0.382 | -0.013 | 0.359
| Wav2Vec2.0
---
Char Comb
0.303 | 0.124 | 0.503
Wav2Vec2.0 LM | 0.520 | 0.108 | 0.637
X-vector | 0.404 | 0.392 | 0.369
ComParE | 0.102 | 0.189 | 0.443
Prosody | 0.029 | 0.376 | 0.008
All features | 0.552 | 0.430 | 0.718
Table 2:  Pearson Correlation between actual and predicted dialect density measures for each of the three metrics: Only the phonological component of the dialect density (DDMphon), only the morphosyntactic component of the dialect density measure (DDMgram), and the entire dialect density measure (DDM). The results for the model trained on six feature sets individually as well as the model trained on the combination of all of the features are shown. Correlation | DDMphon | DDMgram | DDM
---|---|---|---
| Wav2Vec2.0
---
Char Comb
0.339 | 0.126 | 0.495
Wav2Vec2.0 LM | 0.502 | 0.173 | 0.629
All features | 0.569 | 0.385 | 0.678
Table 3:  Average Pearson Correlation between actual and predicted dialect
density measures for each of the three DDMs over 200 iterations of Random Hold
Out validation. Figure 3: SHAP value plot for the XGBoost model trained to
predict DDMphon from the set of all features. The features are listed from top
to bottom in order of significance. Figure 4: SHAP value plot for the XGBoost
model trained to predict DDMgram from the set of all features. The features
are listed from top to bottom in order of significance.

##  4 Discussion

Looking at the individual features, we note that the Wav2Vec2.0 Character
Combinations and Wav2Vec2.0 LM features were especially effective in
estimating dialect density. Many of the character combinations appear to
relate to word initial and word final sounds (eg. N_space (frequency of an “N”
followed by a space character in the ASR transcripts), F_space (frequency of
an “F” followed by a space), and space_U (frequency of a space followed by a
“U”)). This is in line with observations that AAE includes dropping of word
final nasals and glides and simplification of word initial and word final
consonant clusters. Character perplexity (char_ppl) from the language modeling
features was the most impactful feature in estimating all three DDM scores.
This feature is particularly useful in providing an objective distance metric
between the GAE of the Fisher Corpus and the ASR transcripts of the target
dialect speech which, unlike WER, does not require ground truth transcripts or
suffer as heavily in the presence of OOV words. The features derived through
weakly supervised embedding (projected X-vector and Prosody embedding) have
the most significant correlation with DDMgram. This may indicate that learning
grammar from audio files or imperfect transcripts requires larger amounts of
data which our method of weak supervision allows us to utilize. In general,
the ComParE features using Auditory Rasta filtering proved to be most useful.
The RASTA-style filtered auditory spectrum is inspired by psychoacoustics and
has been shown to capture context-dependent information useful in speech
recognition [24].

As Figure 3 shows, the combination of the five most impactful features in
predicting DDMphon was: character perplexity (char_ppl), mean rising slope of
the Rasta-filtered auditory spectrum, the frequency of an “N” character
followed by a space character in the ASR transcripts (N_space), the standard
deviation of distances between peaks in the Rasta-filtered auditory spectrum,
and the PRV component of the projected X-vector. As Figure 4 shows, the five
most impactful features in estimating DDMgram are character perplexity,
duration of sounds predicted to be silence or unintelligible by Wav2Vec2.0
(sil), the PRV component of the prosody embedding, the ROC component of the
projected X-vector, and the frequency of an “F” followed by an “A” in the ASR
transcripts (F_A). The frequency of F_A as a feature may due to a raising of
the vowel following “F” in several words such as “fell” or “fire” as is seen
in some dialects of the US South.

##  5 Conclusions

In this work, we are able to predict dialect density with a high correlation
($~{}68\%$) between the predicted and groundtruth labels as well as explain a
possible combination of factors that make up dialectal differences in ASR. The
results may improve when other factors are considered such as finer aspects of
prosody or within-speaker variability. The system may also improve with more
accurate ASR transcripts. Overall, we have created a novel framework from
which to later explore these factors. Future work includes further exploring
the role of prosody in dialect, applying this framework to other dialects, and
applying this work in dialect-robust ASR systems.

## References

  * [1] J. A. Washington, L. Branum-Martin, C. Sun, and R. Lee-James, “The impact of dialect density on the growth of language and reading in African American children.” Language, speech, and hearing services in schools, 49(2), pp. 232-247. 2018.
  * [2] C. Puranik, L. Branum-Martin, and J. A. Washington, J.A, “The relation between dialect density and the codevelopment of writing and reading in African American children.” Child Development, 91(4), e866–e822. PMID: 31612998
  * [3] A. Koenecke, A Nam, E. Lake, J. Nudell, M. Quartey, Z. Mengesha, C. Toups, J. R. Rickford, D. Jurafsky, and S. Goel, “Racial disparities in automated speech recognition.” Proceedings of the National Academy of Sciences Apr 2020, 117 (14) 7684-7689; DOI: 10.1073/pnas.1915768117
  * [4] S. Lanehart and Ayesha M. Malik. “Language use in African American communities: An introduction.” The Oxford handbook of African American language, ed. by Sonja Lanehart, pp. 1–19. Oxford: Oxford University Press. 2015.
  * [5] S. Li, B. Ouyang, D. Liao, S. Xia, L. Li and Q. Hong, “End-To-End Multi-Accent Speech Recognition with Unsupervised Accent Modelling,” ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 6418-6422, doi: 10.1109/ICASSP39728.2021.9414833.
  * [6] T. Tan, Y. Lu, R. Ma, S. Zhu, J. Guo and Y. Qian, “AISpeech-SJTU ASR System for the Accented English Speech Recognition Challenge,” ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 6413-6417, doi: 10.1109/ICASSP39728.2021.9414471.
  * [7] R. Huang, J. H. L. Hansen and P. Angkititrakul, “Dialect/Accent Classification Using Unrestricted Audio,“ in IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 2, pp. 453-464, Feb. 2007, doi: 10.1109/TASL.2006.881695.
  * [8] O. García, and L. Wei, (2014). “Translanguaging and education. In Translanguaging: Language, bilingualism and education,” pp. 63-77. Palgrave Pivot, London, 2014.
  * [9] Tyler Kendall and Charlie Farrington, “The Corpus of Regional African American Language. Version 2020.05.” Eugene, OR: The Online Resources for African American Language Project. 2020 http://oraal.uoregon.edu/coraal.
  * [10] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, Michael Auli, “Wav2Vec2.0: A Framework for Self-Supervised Learning of Speech Representations.” NeurIPS 2020
  * [11] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur. Librispeech: an asr corpus based on public domain audio books.” In Proc. of ICASSP, pages 5206–5210. IEEE, 2015.
  * [12] A. Akbik, D. Blythe and R. Vollgraf, “Contextual String Embeddings for Sequence Labeling,” in Proceedings of the 27th International Conference on Computational Linguistics, 2018.
  * [13] A. Akbik, T. Bergmann, D. Blythe, K. Rasul, S. Schweter and R. Vollgraf, “FLAIR: An easy-to-use framework for state-of-the-art NLP,” 2019 Annual Conference of the North American Chapter of the Association for Computational Linguistics (Demonstrations), 2019.
  * [14] C. Cieri, D. Graff, O. Kimball, D. Miller, and K. Walker. “Fisher English Training Speech Part 1 Transcripts, LDC2004T19.” Web Download. Philadelphia: Linguistic Data Consortium, 2004.
  * [15] John Hale, “A probabilistic earley parser as a psycholinguistic model”. In Proceedings of the second meeting of the North American Chapter of the Association for Computational Linguistics on Language technologies, pages 1–8. Association for Computational Linguistics, 2001.
  * [16] B. Schuller, S. Steidl, A. Batliner, J. Hirschberg, J. K. Burgoon, A. Baird, A. Elkins, Y. Zhang, E. Coutinho, and K. Evanini, “The INTERSPEECH 2016 Computational Paralinguistics Challenge: Deception, Sincerity & Native Language,” in Proc. Interspeech, 2016, pp. 2001–2005.
  * [17] F. Eyben, M. Wollmer, and B. Schuller, “Opensmile: the Munich Versatile and Fast Open-Source Audio Feature Extractor,” in Proc. ACM Multimedia, Firenze, Italy, 2010, pp. 1459–1462.
  * [18] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey and S. Khudanpur, “X-Vectors: Robust DNN Embeddings for Speaker Recognition,” 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 5329-5333, doi: 10.1109/ICASSP.2018.8461375.
  * [19] Trang Tran, Shubham Toshniwal, Mohit Bansal, Kevin Gimpel, Karen Livescu, and Mari Ostendorf. 2018. Parsing Speech: A Neural Approach to Integrating Lexical and Acoustic-Prosodic Information. In Proc. Conf. North American Chapter Assoc. for Computational Linguistics (NAACL)
  * [20] P. Boersma, “Praat, a system for doing phonetics by computer,” Glot International 5:9/10, pp. 341-345.
  * [21] Tianqi Chen and Carlos Guestrin, “XGBoost: A Scalable Tree Boosting System.” Krishnapuram, Balaji; Shah, Mohak; Smola, Alexander J.; Aggarwal, Charu C.; Shen, Dou; Rastogi, Rajeev (eds.). In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016. ACM. pp. 785–794. arXiv:1603.02754. doi:10.1145/2939672.2939785.
  * [22] Cynthia Puranik, Lee Branum-Martin, and Julie Washington, “The Relation Between Dialect Density and the Codevelopment of Writing and Reading in African American Children.” Child Development. 2019. 91. 10.1111/cdev.13318.
  * [23] S. M. Lundberg and S. I. Lee, “A unified approach to interpreting model predictions,” in Proc. Adv. Neural Inf. Process. Syst., 2017, pp. 4768–4777.
  * [24] H. Hermansky and N. Morgan, “RASTA processing of speech,” in IEEE Transactions on Speech and Audio Processing, vol. 2, no. 4, pp. 578-589, Oct. 1994, doi: 10.1109/89.326616.


# SLAMesh: Real-time LiDAR Simultaneous Localization and Meshing

Jianyuan Ruan1, Bo Li2, Yibo Wang1, and Yuxiang Sun1,∗ This work was supported
by the HK PolyU Start-up Fund under Grant P0034801, which is awarded to Dr.
Yuxiang Sun.1The Hong Kong Polytechnic University, Hung Hom, Hong Kong
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>sun.yuxiang@outlook.com).2Zhejiang University of Science
and Technology, Hangzhou, China (e-mail: 120052@zust.edu.cn).∗Corresponding
author: Yuxiang Sun.

###### Abstract

Most current LiDAR simultaneous localization and mapping (SLAM) systems build
maps in point clouds, which are sparse when zoomed in, even though they seem
dense to human eyes. Dense maps are essential for robotic applications, such
as map-based navigation. Due to the low memory cost, mesh has become an
attractive dense model for mapping in recent years. However, existing methods
usually produce mesh maps by using an offline post-processing step to generate
mesh maps. This two-step pipeline does not allow these methods to use the
built mesh maps online and to enable localization and meshing to benefit each
other. To solve this problem, we propose the first CPU-only real-time LiDAR
SLAM system that can simultaneously build a mesh map and perform localization
against the mesh map. A novel and direct meshing strategy with Gaussian
process reconstruction realizes the fast building, registration, and updating
of mesh maps. We perform experiments on several public datasets. The results
show that our SLAM system can run at around $40$Hz. The localization and
meshing accuracy also outperforms the state-of-the-art methods, including the
TSDF map and Poisson reconstruction. Our code and video demos are available
at: https://github.com/lab-sun/SLAMesh.

##  I Introduction

Simultaneous localization and mapping (SLAM) aims to estimate sensor poses and
reconstruct traversed environments at the same time. 3-D LiDAR SLAM has
attracted significant attention in the robotics research community in recent
years. Usually, LiDAR SLAM produces a point-cloud map, which is a basic
representation for 3-D environment models [1]. With such map type,
registration and fusion of new range data can be straightforwardly operated
[2]. Unfortunately, a point-cloud map requires extra data structures, such as
KD-tree, to speed up the query of the unordered points. The memory cost also
grows significantly with the increased scale or density. Moreover, although
point-cloud map seems dense to human eyes, it would still be sparse when
zoomed in. Due to the sparsity, it usually requires post-processing operations
before the use for navigation, such as voxelization, to produce dense voxel
maps.

Figure 1: Left: Mesh maps produced by SLAMesh, Puma [3] and Voxblox [4] on the
Mai City dataset [3]. A zoomed view in each map highlights the connections of
the mesh edges. The mesh map of Puma is water-tight but with relatively
complex multi-layers. The mesh maps produced by Voxblox and our SLAMesh are of
high orderliness, and the vertices are evenly distributed. Right: Mesh maps
built by our SLAMesh on the KITTI sequences 00, 07, and 09 demonstrate the
accuracy and scalability of SLAMesh in large-scale environments.

In this work, we consider a dense map as a kind of representation that models
3-D objects in terms of surfaces or volumes. Under this definition, there
exist several types of dense maps. For instance, occupancy grid maps [5] and
truncated signed distance function (TSDF) maps [6] divide the 3-D space into
voxels, making them applicable for robot navigation [4]. However, the
discretization introduces the scalability problem and limits the localization
accuracy to the grid size. The normal distribution transform (NDT) [7][8] and
surfel maps [9][10][11] can also be seen as dense maps. They use parameterized
representations to describe local structures. Parameterization ensures memory
efficiency but limits the ability to describe fine details unless a high
resolution is used. Besides, gaps exist around the NDT and surfel map
elements.

Mesh is a 3-D solid surface comprised of faces, edges, and vertices. In the
field of 3-D modeling, the triangular mesh has become the dominating
representation because it is simple and can approximate most complex 3-D
structures. Unlike grid maps, mesh maps alleviate the discretization problem
and can model smooth surfaces for robotics applications. Mesh can also depict
manifold structures compared to surfel and NDT maps. In addition, mesh map is
memory efficient, scalable, and retains the topological information [12][13].

However, real-time 3-D LiDAR mesh mapping is still hard to realize now,
especially in large-scale environments. The main difficulty is that mesh
building and updating are usually time-consuming. A widely-adopted solution to
this problem is first focusing on point-cloud mapping and then performing
meshing separately using the point-cloud maps. This two-step solution could
not enable localization and meshing to benefit each other [14].

To address this issue, this work designs a simultaneous localization and
meshing system named SLAMesh. The system can run in real time with only a CPU.
In SLAMesh, we make the vertices uniformly distributed, enabling mesh
building, registration, and updating more efficiently. Fig. 1 shows such a
mesh map constructed by the Gaussian process (GP) [15] reconstruction. Our
previous works [16][17] have shown the benefits of utilizing GP in SLAM,
especially when the range data is sparse. However, they produce point-cloud
maps without exporting the surfaces model. In addition, they cannot achieve
real-time performance when using sensors with denser point clouds, like
64-beam LiDAR. This work exploits GP for mesh map building and is an upgraded
version upon GP-SLAM+ [17]. The main contributions of this work are summarized
as follows:

  * •

We propose a novel meshing strategy based on GP reconstruction and vertices
connection, which allows fast building, query, and updating of mesh maps.

  * •

We design a point-to-mesh registration method. Together with constraints
combination and multi-thread implementation, we ensure the efficiency and
accuracy of localization and mapping with the mesh map.

  * •

We develop a dense, real-time, and scalable open-source LiDAR SLAM system upon
mesh map and demonstrate the advantages through experiments.

##  II Related Work

Indoor dense mapping is more investigated than the outdoor case. A widely-
adopted idea to build mesh is fusing depth data into a TSDF map and then
extracting the implicit mesh on demand through the Marching cubes algorithm
[18] like Voxblox [4] and Kinectfusion [19]. On the contrary, some methods
produce mesh maps explicitly. Ryde et al. [20] voxelized point clouds and
regressed a plane inside each voxel to form polygons. Piazza et al. [21] used
the Delaunay triangulation [22] on the sparse feature points to build and
update the mesh. Schreiberhuber et al. [23] connected the depth point in the
camera image and then stitched the local mesh patches together. Our work
explicitly produces mesh maps based on GP. Different from [23], our SLAMesh
tackles unordered point clouds and builds mesh maps in the world frame rather
than merging mesh patches in different image frames.

Recently, reliable 3-D LiDARs have enabled outdoor mapping. The outdoor
environments are more challenging than indoor ones because the scenarios are
several magnitudes larger, vehicle movements are fast, and the point cloud is
sparse in remote areas. Suma [9] and Suma++ [10] provide accurate odometry,
but the surfel map is relatively cluttered. LiTAMIN2 [8] shows the efficiency
advantage of NDT registration with a map consisting of ellipsoids. Based on
the TSDF map, Kühner et al. [24] modified the ray-tracking model for LiDAR and
included an offline refinement process. Roldão et al. [25] combined explicit
and implicit meshing approaches for better accuracy and mesh density. Those
methods [24][25] mainly focus on mapping function and can not provide real-
time odometry. Puma [3] and our SLAMesh both uses mesh maps. However, Puma
utilizes ray-casting for registration and needs to rebuild the mesh map
repetitively via the Poisson reconstruction [26] rather than updating it. Lin
et al. [27] also investigated online meshing in LiDAR SLAM. They use a
separate SLAM system and a direct connecting meshing strategy upon point-cloud
maps. Our method reconstructs unordered points to integrate the mesh map into
the whole system.

The ray-tracking model in the above TSDF-based methods neglects the fact that
local surfaces are continuous. Conversely, GP builds a continuous model
accompanied by uncertainty. This property attracts people to use GP for
recovering occupancy field [28] and TSDF field [29], which are continuous
naturally. Different from [28][29], we use GP to predict the Euclidean
coordinate of surfaces. Moreover, our system can also apply to large-scale
environments.

##  III The Proposed Approach

Figure 2: The general diagram of our SLAMesh. The system has three parts:
meshing, registration, and mesh management. (a) The raw point cloud is
transformed into the world frame $\\{W\\}$, down-sampled, and assigned into
voxel cells; (b) GP models the local surface (light purple) of each cell to
get uniformly distributed vertices. The color of vertices indicates their
uncertainty (warmer means higher), which increases with its distance to the
raw points; (c) We build mesh (light blue) by connecting adjacent vertices
directly; (d) Point-to-mesh registration aligns the reconstructed current scan
(red points) to the mesh map. Fast matching is feasible based on the locations
of vertices. The normal is smoothed from the surrounding mesh. GPed scan means
scan after the GP reconstruction. Finally, we only need to adjust the 1-D
predictions of vertices and their uncertainty to maintain the mesh map.

###  III-A Approach Overview

The motivation of this work is to build, utilize, and maintain a mesh map in a
LiDAR SLAM system. Fig. 2 illustrates a general overview. The system mainly
consists of three components: meshing, registration, and mesh management.
Firstly, each new LiDAR scan is transformed into the world frame $\\{W\\}$ as
$\mathcal{S}^{raw}$ using the initial guess from a constant-velocity model.
The following operations are also performed in $\\{W\\}$. Then, points are
assigned into voxel cells. GP reconstructs the local surface inside each cell
and obtains vertices $\boldsymbol{v}_{i}$, which are connected to form a mesh.
In the registration component, a point-to-mesh registration is designed to
align the reconstructed current scan with the built mesh map $\mathcal{M}$.
Finally, the mesh map is updated iteratively.

###  III-B Meshing Strategy

As aforementioned, building and updating mesh is time-consuming. To tackle
this problem, we adopt a reconstructing and connection strategy to facilitate
the following pipeline so that the whole system can run in real-time. As shown
in Fig. 2, GP recovers the local surface from noisy and sparse point clouds
inside voxels. Then, the vertices are interpolation results of the surface.
Two coordinates of the 3-D vertices are evenly located (named as locations),
and the other one (named as predictions) has a continuous value domain. The
locations serve as the indices to enable fast query in constant time. The
value domain of predictions is continuous to avoid discretization-induced
accuracy loss.

Here the GP process is described (a more detailed description of GP can be
found in [16][17]). The bold lowercase letters represent vectors, and the
uppercase letters represent matrices. The input to GP is a subset of the raw
point cloud
$\mathcal{S}_{k}^{raw}=\left\\{\boldsymbol{p}_{i}=(x_{i},y_{i},z_{i},\sigma_{in}^{2}),i=1,....,n_{i}\right\\}$
containing $n_{i}$ points, where $k$ represents the $k$-th cell
$\mathcal{C}_{k}$ in the current scan, $\sigma_{in}^{2}$ is the isotropic
variance of input noise. The output is a layer containing $n_{j}$ vertices
$\mathcal{L}_{k}^{z}=\left\\{\boldsymbol{v}_{j}=(x_{j},y_{j},z_{j},\sigma_{j}^{2}),j=1,....,n_{j}\right\\}$
with uncertainty $\sigma_{j}^{2}$. The superscript $z$ means that GP predicts
the coordinate $z$ as
$\boldsymbol{\widetilde{f}}=\left\\{z_{j},j=1,....,n_{j}\right\\}$, and it is
omitted when the coordinate can be any one in $x$, $y$, or $z$. We denote the
locations of input and output point set as $\boldsymbol{i}$ and
$\boldsymbol{j}$. Given the input observation $\boldsymbol{f}$, the
predictions $\boldsymbol{\widetilde{f}}$ also follow the Gaussian
distribution. The expectation of $\boldsymbol{\widetilde{f}}$ is:

$\boldsymbol{\widetilde{f}}=\boldsymbol{k}_{\boldsymbol{ij}}^{T}\left(\sigma_{in}^{2}\boldsymbol{I}+\boldsymbol{K}_{\boldsymbol{ii}}\right)^{-1}\boldsymbol{f}.$
(1)

The uncertainty of predictions is its variance:

$\boldsymbol{\sigma_{j}}^{2}=\boldsymbol{k}_{\boldsymbol{jj}}-\boldsymbol{k}_{\boldsymbol{ij}}^{\boldsymbol{T}}\left(\sigma_{in}^{2}\boldsymbol{I}+\boldsymbol{K}_{\boldsymbol{ii}}\right)^{-1}\boldsymbol{k}_{\boldsymbol{ij}},$
(2)

where $\boldsymbol{k}_{\boldsymbol{jj}}$, $\boldsymbol{k}_{\boldsymbol{ij}}$,
and $\boldsymbol{K}_{\boldsymbol{i}\boldsymbol{i}}$ represent different
combinations of the kernel function of locations. For instance,

$k\left(i,j\right)=\exp\left(-\kappa\left|i-j\right|\right),$ (3)

where $k\left(i,j\right)$ is a scale value, and $\kappa$ is a constant
($\kappa=1$ in our algorithm). This exponential kernel function can represent
a local smooth surface in a 2-D manifold. Thus, one cell can contains more
mesh layers with different GP functions (up to $3$ if all coordinates are
interpolated) when modeling a complicated structure [17].

A vertex with uncertainty $\sigma_{j}^{2}$ under a certain threshold
$\sigma_{match}^{2}$ is accurate enough or valid. Triangle mesh faces are
built by connecting adjacent or diagonal vertices in the 2-D space of
locations. A mesh face is valid if all vertices are valid. Similar to the
Delaunay triangulation, which establishes edges in a 2-D space, our method can
prevent sliver mesh faces.

###  III-C Point-to-Mesh Registration

As aforementioned, we perform localization and meshing simultaneously so that
meshing can benefit localization. To this end, one intuitive idea is to treat
vertices as points or extract points from faces and then employ traditional
point-cloud registration for pose estimation. However, this idea neglects the
normal information of the mesh faces. Puma [3] indicates that a point-to-mesh
error could improve accuracy. Unlike the ray-casting data association in Puma,
our SLAMesh establishes correspondences based on locations.

For a vertex $\boldsymbol{v}_{p}$ in the reconstructed current scan
$\mathcal{S}^{gp}$, we first query the sub-mesh layer located in the same or
adjacent cells (see Fig. 3), and then find the vertex $\boldsymbol{v}_{q}$
sharing the same locations (see Fig. 2(d)). Data association is established
between $\boldsymbol{v}_{p}$ and the valid faces that contain
$\boldsymbol{v}_{q}$. Finding the $n_{q}$ neighboring vertices along the edges
in a mesh is quick. The normals of those faces that contain
$\boldsymbol{v}_{q}$ are averaged to get a smoothed normal
$\overline{\boldsymbol{n}}$ in case the surface is rugged:

$\overline{\boldsymbol{n}_{q}}=\frac{\sum_{q=1}^{n_{q}-2}(\boldsymbol{v}_{q}-\boldsymbol{v}_{q+1})\times(\boldsymbol{v}_{q}-\boldsymbol{v}_{q+2})}{\sum_{q=1}^{n_{q}-2}\left\|(\boldsymbol{v}_{q}-\boldsymbol{v}_{q+1})\times(\boldsymbol{v}_{q}-\boldsymbol{v}_{q+2})\right\|},$
(4)

where $\left\|\cdot\right\|$ is the 2-norm. We denote the number of
correspondences in one cell as $n_{q}$ and the number of overlapped cells as
$K$. The point $\boldsymbol{v}_{p}^{\prime}$ and the point-to-mesh residual
after applying the transformation $\boldsymbol{T}\in SE3$ including rotation
$\boldsymbol{R}$ and translation $\boldsymbol{t}$ on $\boldsymbol{v}_{p}$ are:

$\left[\begin{array}[]{c}\boldsymbol{v}_{p}^{\prime}\\\ 1\\\
\end{array}\right]=\boldsymbol{T}\left[\begin{array}[]{c}\boldsymbol{v}_{p}\\\
1\\\
\end{array}\right]=\left[\begin{array}[]{cc}\boldsymbol{R}&\boldsymbol{t}\\\
\boldsymbol{0}&1\\\
\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{v}_{p}\\\ 1\\\
\end{array}\right],$ (5)
$e_{p}=\overline{\boldsymbol{n}_{q}^{T}}(\boldsymbol{v}_{p}^{\prime}-\boldsymbol{v}_{q}),$
(6)

where the superscript $T$ means matrix transpose. The optimal relative
transformation $\boldsymbol{T}$ is computed in the optimization problem where
point-to-mesh residuals in all $K$ overlapped layers $\mathcal{L}_{k}$ are
joined:

${\boldsymbol{T}}=\underset{\boldsymbol{T}}{{argmin}}\sum_{k=1}^{K}\sum_{\boldsymbol{v}_{p}\in\mathcal{L}_{k},p=1}^{n_{p}}e_{p}.$
(7)

We use the Levenberg–Marquardt algorithm in Ceres solver111http://ceres-
solver.org/ to solve this non-linear least square problem. The analytical
Jacobin matrix is derived to accelerate the solving process:

$\boldsymbol{J}=\frac{\partial e_{p}}{\partial\boldsymbol{T}}=\frac{\partial
e_{p}}{\partial\boldsymbol{v}_{p}^{\prime}}\frac{\partial\boldsymbol{v}_{p}^{\prime}}{\partial\boldsymbol{T}},$
(8) $\frac{\partial
e_{p}}{\partial\boldsymbol{v}_{p}^{\prime}}=\overline{\boldsymbol{n}_{p}^{T}},\hskip
15.00002pt\frac{\partial\boldsymbol{v}_{p}^{\prime}}{\partial\boldsymbol{t}}=\boldsymbol{I},\hskip
15.00002pt\frac{\partial\boldsymbol{v}_{p}^{\prime}}{\partial\boldsymbol{R}}=-\left(\boldsymbol{R}\boldsymbol{v}_{p}+\boldsymbol{t}\right)_{\times},$
(9)

where the symbol $\left(\cdot\right)_{\times}$ represents the corresponding
skew-symmetric matrix of a vector. The mesh may contain many faces so the
optimization problem can be huge. We combine the residuals inside each layer
to be one by averaging before optimization to speed up the optimization
process.

Figure 3: The data association can be established across cells alongside the
predicted coordinates (green cells). The query length $b$ can be decreased as
registration tends to converge. From (a) to (c), $b$ is decreased from $2$ to
$0$.

Storing data in voxel cells induces discontinuities of data association on the
border. If we only perform registration on overlapped cells, when the vehicle
moves fast, the convergence basins of registration would be too narrow. Thus,
we allow cross-cell overlapping by querying cells alongside the axis of
predictions for each layer as shown in Fig. 3. The length of the query is
decreased when the registration tends to converge.

###  III-D Mesh Management and Multithreading

Mesh maps are more challenging to maintain than point-cloud maps and grid-
based maps. Points and grids can be independent of their neighborhoods. On the
contrary, mesh retains the topological structure with vertices connected to
each other. The map update process should maintain existing connections and
avoid sliver triangles. Typical solutions include re-meshing with Delaunay
triangulation [21], maintaining an implicit filed similar to Voxblox [4], or
repeating the meshing process like Puma [3]. SLAM requires a fast update of
mesh. In our SLAMesh, the locations of vertices are fixed, so only the 1-D
predictions need to be updated. Given the previous $t$ data below the
threshold $\sigma_{update}^{2}$, the new prediction can be solved by the
iterative least square:

${\widetilde{f}_{t_{k}}}={\sum_{t=1}^{t_{k}}({\widetilde{f}_{t}}\sigma_{t}^{2})}/{\sum_{t=1}^{t_{k}}(\sigma_{t}^{2}),\text{if
$\sigma_{t}^{2}<\sigma_{update}^{2}$}}.$ (10)

The cells are stored in a hash map, which ensures a linear complexity for
inserting, deleting, and querying. Also, this structure is flexible because it
can grow incrementally. An advantage of our cell-based map shared with VGICP
[30] is the compatibility for parallel processing. Since each cell is
independent in these steps, we can accelerate the GP reconstruction and mesh
publishing processes with multithreading.

##  IV Experimental Results and Discussions

###  IV-A Experimental Setup

We compare SLAMesh with several state-of-the-art SLAM systems or meshing
tools. Algorithms are run on a PC with a 3.6 GHz 8-core Intel i7-11700KF CPU.
In SLAMesh, the number of predictions in each cell $n_{j}=6^{2}$. The cell
size is $1.5$m so that the resolution of locations is identical to the voxel
size in Voxblox when evaluating mesh. In other tests, the cell size is $1.6$m.
The variance of the input noise is $\sigma_{in}^{2}=0.02$. The thresholds
$\sigma_{match}^{2}=0.5$ and $\sigma_{update}^{2}=1.0$. These parameters are
tuned manually. $8$ threads are allocated in multithreading. We use the Mesh
Tools [12] for visualization. Other methods are run in their default
configurations.

Figure 4: Mesh maps compared with the ground-truth model on the Mai City
dataset [3]. The left and right columns illustrate Precision and Recall
values, respectively. Darker colors indicate worse performance and vice versa.
Voxblox recovers good results but fails to build thin structures like trees.
Puma presents less accuracy on the edges of structures or areas where the pose
drifts. Our SLAMesh presents both good Recall and Precision.

###  IV-B Meshing Evaluation

The first experiment evaluates the meshing quality in terms of accuracy and
completeness. We use the Mai City dataset [3] built on the CARLA simulation
environment with ground truth of maps. In this dataset, the sensor is a
simulated 64-beam Velodyne HDL-64E LiDAR. The vehicle drives for $99$m and
produces $100$ frames of LiDAR scans. The ground truth of the environment is a
dense point cloud scanned by a high-resolution sensor. We evaluate the
combination of Voxblox with A-LOAM222https://github.com/HKUST-Aerial-
Robotics/A-LOAM, Puma, and SLAMesh. Because this work aims to build a SLAM
system rather than an offline mapping tool, and the ground truth pose is
usually unavailable in real environments, we use the poses estimated in each
pipeline or an extra odometry system (i.e., using A-LOAM for Voxblox).

Fig. 1 presents the mesh map produced by each method. All maps reconstruct the
main structures well. In the zoomed views, we can observe the characteristic
of each strategy. The vertices in SLAMesh and Voxblox are evenly organized due
to the reconstruction or voxelization. The mesh of Puma is water-tight with
Poisson reconstruction but is relatively complex with multi-layer phenomenons.
We believe that this is because SLAMesh and Voxblox both iteratively fuse
observation into surfaces, while Puma accumulates several scans and then
conducts a meshing upon them.

The quantitative evaluation uses the standard point-cloud-based metrics:
Precision, Recall, and F1-score [31]. We densely and uniformly sample the mesh
maps to form a point cloud, and then compare it with the ground-truth point
cloud. The quantitative results are displayed in Table. I, where the distance
threshold $d$ is set to $0.3$m. The performance of our SLAMesh is best,
followed by Voxblox+A-LOAM and Puma. Fig. 4 visualizes the Precision and
Recall results. The values decrease from lighter to darker. Voxblox erases
tiny or thin objects like trees (the black areas in the Recall sub-figures)
with its ray-tracing method. Puma accumulates errors along the trajectory. The
water-tight assumption makes the top edge of walls slightly warped. Our
SLAMesh can recover both large and small structures with high accuracy.

###  IV-C Odometry Evaluation

TABLE I: Quantitative comparative results ($\%$) of meshing. The bold font indicates best results. Metric | Voxblox[4]+A-LOAM | Puma [3] | SLAMesh (Ours)
---|---|---|---
Precision | 69.50 | 43.00 | 74.96
Recall | 85.91 | 62.90 | 86.09
F1-score | 76.83 | 51.08 | 80.14
TABLE II: Quantitative results of odometry accuracy on the KITTI dataset. The first and the second row are the Relative Translation Error (%) and Rotation Error (deg/100m), respectively. The red bold and the blue fonts indicate the best and the second best results, respectively. Comb. means constraints combination, and P2Mesh means point-to-mesh registration. Map type | Method | 00 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | Mean
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Point Cloud | LOAM [1] | 0.78 | 1.43 | 0.92 | 0.86 | 0.71 | 0.57 | 0.65 | 0.63 | 1.12 | 0.77 | 0.79 | 0.839
- | - | - | - | - | - | - | - | - | - | - | -
A-LOAM | 0.97 | 2.75 | 4.91 | 1.22 | 1.35 | 0.63 | 0.61 | 0.48 | 1.17 | 1.11 | 1.58 | 1.525
0.74 | 0.96 | 2.70 | 1.06 | 0.67 | 0.51 | 0.59 | 0.66 | 0.76 | 0.75 | 0.82 | 0.928
Surfel | Suma [9] | 0.70 | 1.70 | 1.10 | 0.70 | 0.40 | 0.50 | 0.40 | 0.40 | 1.00 | 0.50 | 0.70 | 0.736
0.30 | 0.30 | 0.40 | 0.50 | 0.30 | 0.20 | 0.20 | 0.30 | 0.40 | 0.30 | 0.30 | 0.318
Suma++ [10] | 0.64 | 1.60 | 1.00 | 0.67 | 0.37 | 0.40 | 0.46 | 0.34 | 1.10 | 0.47 | 0.66 | 0.701
0.22 | 0.46 | 0.37 | 0.46 | 0.26 | 0.20 | 0.21 | 0.19 | 0.35 | 0.23 | 0.28 | 0.294
NDT | Litamin2 [8] | 0.70 | 2.10 | 0.98 | 0.96 | 1.05 | 0.45 | 0.59 | 0.44 | 0.95 | 0.69 | 0.80 | 0.883
0.28 | 0.46 | 0.32 | 0.48 | 0.52 | 0.25 | 0.34 | 0.32 | 0.29 | 0.40 | 0.47 | 0.375
Mesh | Puma [3] | 1.46 | 3.38 | 1.86 | 1.60 | 1.63 | 1.20 | 0.88 | 0.72 | 1.44 | 1.51 | 1.38 | 1.551
0.68 | 1.00 | 0.72 | 1.10 | 0.92 | 0.61 | 0.42 | 0.55 | 0.61 | 0.66 | 0.84 | 0.737
SLAMesh (Ours) | Full | 0.77 | 1.25 | 0.77 | 0.64 | 0.50 | 0.52 | 0.53 | 0.36 | 0.87 | 0.57 | 0.65 | 0.676
0.35 | 0.30 | 0.30 | 0.43 | 0.13 | 0.30 | 0.22 | 0.23 | 0.27 | 0.25 | 0.42 | 0.291
w/o Comb. | 0.96 | 1.70 | 1.04 | 0.81 | 0.67 | 0.51 | 0.64 | 0.46 | 0.86 | 0.74 | 1.02 | 0.857
0.41 | 0.41 | 0.39 | 0.55 | 0.30 | 0.29 | 0.30 | 0.28 | 0.30 | 0.32 | 0.44 | 0.363
w/o P2Mesh | 0.70 | 1.33 | 0.86 | 0.59 | 0.69 | 0.49 | 0.37 | 0.43 | 1.14 | 0.77 | 0.93 | 0.756
0.43 | 0.29 | 0.38 | 0.42 | 0.36 | 0.34 | 0.13 | 0.40 | 0.42 | 0.29 | 0.57 | 0.366
Figure 5: Trajectories estimated by our SLAMesh on the KITTI $00\sim 10$
sequences (solid blue line) compared with the ground truth (dash gray line).
Most trajectories are well consistent with the ground truth even without loop
closures.

We also quantitatively evaluate our pose estimation performance in the widely-
used KITTI [32] odometry benchmark. The KITTI dataset provides point clouds
from a Velodyne HDL-64E LiDAR and ground-truth poses. The results on the
sequences $00\sim 10$ (including urban, country, and highway environments) are
shown in Fig. 5. Our estimated trajectories present high consistency compared
with ground truth.

Several state-of-the-art methods with different types of maps are tested. Suma
[9] is a surfel-based method, and Suma++ [10] enhances the odometry accuracy
with dynamic object removal. LiTAMIN2 [8] maintains a voxelized NDT map.
A-LOAM is an implementation of LOAM [1]. Puma [3] is also a mesh-based method.
Tab. II shows the quantitative evaluation results with the standard relative
translation error and rotation error [32]. The results of these methods are
directly imported from their published papers. Our SLAMesh achieves superior
performance with $0.676\%$ average translation error and $0.291$deg/100m
rotation error, which outperforms those point cloud, NDT, surfel, and mesh-
based methods. Notice that we can achieve this without any loop closure. Also,
as reported in [9][10], Suma needs GPU, and the surfel map is relatively
cluttered.

Fig. 1 also displays the online mesh maps built by SLAMesh in the sequences
00, 07, and 09. The maps exhibit fine alignments when the vehicle travels back
to the start. The maximum travel length of these sequences is $5.07$km, which
demonstrates the scalability of SLAMesh. An interesting observation is that
when we allow the current scan to register with revisited areas, the map and
the whole trajectory show better consistency, but when evaluating the relative
error of odometry in KITTI, disabling this implicit loop closure enables
higher accuracy. We believe this phenomenon is common in methods using cell-
organized maps like [30].

Figure 6: Memory consumption of maps with the growing number of frames on the
KITTI sequence 07. The raw point cloud can occupy several GB of disk space.
The three mesh-based methods have lower memory costs. However, the mesh maps
produced by Puma are more complex and heavier than those produced by Voxblox
and SLAMesh. Figure 7: Time cost per frame on the KITTI sequence 07. SLAMesh
can run at 40Hz. The Voxblox+A-LOAM pipeline [4] achieves approximate real-
time performance, while Puma cannot.

###  IV-D Memory and Computational Efficiency

The map size from each method in the KITTI sequence 07 is reported in Fig. 6.
The raw point-cloud map proliferates to several GBs, which would be too heavy
for online tasks. Three mesh-based methods demonstrate their superiority in
memory efficiency. Voxblox and SLAMesh consume fewer memory resources than
Puma. We conjecture the reason is that Puma does not iteratively fuse point
clouds.

Fig. 7 displays the time cost per frame on the KITTI sequence 07. The time
cost of SLAMesh does not increase with the scale of the scenario due to the
cell-organized map. Note that the LiDAR runs at $10$Hz in this dataset.
However, Puma costs $4.7$s per frame on average. For the first 30 frames, Puma
uses iterative closest point (ICP) rather than mesh. The A-LOAM+Voxblox
pipeline costs $129$ms per frame ($7.8$Hz), which is approximately real-time.
Our SLAMesh can run at $40$Hz and surpasses other methods.

The time costs of each module are further analyzed in Fig. 8. They are
generally categorized into three groups. The first part, pre-process, includes
down-sampling and GP reconstruction in SLAMesh, feature extraction in A-LOAM,
and normal computation in Puma. The second part is registration, including
data association and optimization. The last part is map maintenance, including
mesh map update and publishing. In Puma, the Poisson reconstruction is the
main burden. In SLAMesh and A-LOAM, registration both iterates twice. The main
computational cost of SLAMesh is the pre-process (account for $63\%$), while
for A-LOAM, the registration and map update cost $87\%$ time. This observation
shows our different strategies for using reconstruction to structure data in
advance. Voxblox is an efficient meshing tool due to its fast-integrate
strategy. The max range of the sensor is cut into $50$m rather than $100$m,
which also benefits its efficiency. However, the pipeline needs to maintain
two maps, a point cloud map with KD-tree and a TSDF map, which costs more
computation and memory resources.

###  IV-E Ablation Study

We perform an ablation study to investigate the contribution of each technique
in our system. The parameter in SLAMesh is the same for all the variants. Tab.
II shows the results. If SLAMesh disables the point-to-mesh error metric, the
translation and rotation errors increase. This demonstrates that the pose
estimation could be more accurate given more geometric information, like the
normal of surfaces. The constraints combination also reduces the errors. The
reason is that individual point-to-mesh error metric may not be robust enough
in unstructured areas, and the constraints combination could suppress the
influence of outliers.

Figure 8: The time cost breakdown for each module in different methods. The
main burdens in SLAMesh, Puma, and A-LOAM are GP reconstruction, mesh
building, and registration, respectively. The Voxblox+A-LOAM pipeline needs to
maintain two maps. The constraints combination (Comb.) and multithreading
(Thread.) notably accelerate the GP-reconstruction and optimization modules in
our SLAMesh. Note that Puma uses the bottom x-axis, and the other methods use
the top x-axis. The figure is best viewed in color.

Constraints combination and multithreading also play essential roles in
efficiency improvement. In Fig. 8, the multithreading reduces the processing
time of SLAMesh from $57.1$ms to $23.8$ms per scan. In detail, multithreading
makes the GP reconstruction three times faster and the mesh publishing $2.5$
times faster. We note that Puma and Voxblox also leverage the multi-cores of
the CPU in practice, which guarantees a fair comparison. The constraints
combination saves about $95\%$ time cost in optimization as the average number
of constraints is reduced from $9,000$ to $540$.

We conduct extensive tests with a set of parameters and find that our method
has a wide tolerance for parameters. As the cell size grows to $3$m, the frame
rate of SLAMesh can increase to about $80$Hz with 0.88% relative translation
error. The reason is that we use an un-parameterized model inside each voxel
rather than a simple element.

##  V Conclusions

We proposed a novel real-time localization and meshing system based on the
reconstruction and connection strategy. Uniformly populated mesh vertices are
the main property of the map, which lays a good foundation for mesh building,
data association, and mesh update. Experiments showed that SLAMesh presents
high odometry accuracy, efficiency, and mesh quality compared to state-of-the-
art methods.

## References

  * [1] Ji Zhang and Sanjiv Singh.  Low-drift and real-time lidar odometry and mapping.  Autonomous Robots, 41(2):401–416, 2017.
  * [2] Paul J Besl and Neil D McKay.  Method for registration of 3-D shapes.  In Sensor fusion IV: control paradigms and data structures, volume 1611, pages 586–606. International Society for Optics and Photonics, 1992\.
  * [3] Ignacio Vizzo, Xieyuanli Chen, Nived Chebrolu, Jens Behley, and Cyrill Stachniss.  Poisson surface reconstruction for lidar odometry and mapping.  In Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA), 2021.
  * [4] Helen Oleynikova, Zachary Taylor, Marius Fehr, Roland Siegwart, and Juan Nieto.  Voxblox: Incremental 3d euclidean signed distance fields for on-board mav planning.  In 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 1366–1373. IEEE, 2017.
  * [5] Armin Hornung, Kai M Wurm, Maren Bennewitz, Cyrill Stachniss, and Wolfram Burgard.  Octomap: An efficient probabilistic 3d mapping framework based on octrees.  Autonomous robots, 34(3):189–206, 2013.
  * [6] Brian Curless and Marc Levoy.  A volumetric method for building complex models from range images.  In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pages 303–312, 1996.
  * [7] Martin Magnusson, Achim Lilienthal, and Tom Duckett.  Scan registration for autonomous mining vehicles using 3D-NDT.  Journal of Field Robotics, 24(10):803–827, 2007.
  * [8] Masashi Yokozuka, Kenji Koide, Shuji Oishi, and Atsuhiko Banno.  Litamin2: Ultra light lidar-based slam using geometric approximation applied with kl-divergence.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 11619–11625. IEEE, 2021.
  * [9] Jens Behley and Cyrill Stachniss.  Efficient surfel-based slam using 3d laser range data in urban environments.  In Robotics: Science and Systems, volume 2018, 2018.
  * [10] Xieyuanli Chen, Andres Milioto, Emanuele Palazzolo, Philippe Giguere, Jens Behley, and Cyrill Stachniss.  Suma++: Efficient lidar-based semantic slam.  In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 4530–4537. IEEE, 2019.
  * [11] Chanoh Park, Peyman Moghadam, Jason L. Williams, Soohwan Kim, Sridha Sridharan, and Clinton Fookes.  Elasticity meets continuous-time: Map-centric dense 3d lidar slam.  IEEE Transactions on Robotics, 38(2):978–997, 2022.
  * [12] Sebastian Pütz, Thomas Wiemann, and Joachim Hertzberg.  The mesh tools package–introducing annotated 3d triangle maps in ros.  Robotics and Autonomous Systems, 138:103688, 2021.
  * [13] Xieyuanli Chen, Ignacio Vizzo, Thomas Läbe, Jens Behley, and Cyrill Stachniss.  Range image-based lidar localization for autonomous vehicles.  arXiv preprint arXiv:2105.12121, 2021.
  * [14] Hugh Durrant-Whyte and Tim Bailey.  Simultaneous localization and mapping: part i.  IEEE robotics & automation magazine, 13(2):99–110, 2006.
  * [15] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer School on Machine Learning, pages 63–71. Springer, 2003\.
  * [16] Bo Li, Yingqiang Wang, Yu Zhang, Wenjie Zhao, Jianyuan Ruan, and Ping Li.  GP-SLAM: laser-based SLAM approach based on regionalized Gaussian process map reconstruction.  Autonomous Robots, 44(6):947–967, 2020.
  * [17] Jianyuan Ruan, Bo Li, Yinqiang Wang, and Zhou Fang.  Gp-slam+: real-time 3d lidar slam based on improved regionalized gaussian process map reconstruction.  In 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 5171–5178. IEEE, 2020.
  * [18] William E Lorensen and Harvey E Cline.  Marching cubes: A high resolution 3d surface construction algorithm.  ACM siggraph computer graphics, 21(4):163–169, 1987.
  * [19] Richard A Newcombe, Shahram Izadi, Otmar Hilliges, David Molyneaux, David Kim, Andrew J Davison, Pushmeet Kohi, Jamie Shotton, Steve Hodges, and Andrew Fitzgibbon.  Kinectfusion: Real-time dense surface mapping and tracking.  In 2011 10th IEEE international symposium on mixed and augmented reality, pages 127–136. IEEE, 2011.
  * [20] Julian Ryde, Vikas Dhiman, and Robert Platt.  Voxel planes: Rapid visualization and meshification of point cloud ensembles.  In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 3731–3737. IEEE, 2013.
  * [21] Enrico Piazza, Andrea Romanoni, and Matteo Matteucci.  Real-time cpu-based large-scale three-dimensional mesh reconstruction.  IEEE Robotics and Automation Letters, 3(3):1584–1591, 2018.
  * [22] Der-Tsai Lee and Bruce J Schachter.  Two algorithms for constructing a delaunay triangulation.  International Journal of Computer & Information Sciences, 9(3):219–242, 1980.
  * [23] Simon Schreiberhuber, Johann Prankl, Timothy Patten, and Markus Vincze.  Scalablefusion: High-resolution mesh-based real-time 3d reconstruction.  In 2019 International Conference on Robotics and Automation (ICRA), pages 140–146. IEEE, 2019.
  * [24] Tilman Kühner and Julius Kümmerle.  Large-scale volumetric scene reconstruction using lidar.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 6261–6267. IEEE, 2020.
  * [25] Luis Roldão, Raoul de Charette, and Anne Verroust-Blondet.  3d surface reconstruction from voxel-based lidar data.  In 2019 IEEE Intelligent Transportation Systems Conference (ITSC), pages 2681–2686. IEEE, 2019.
  * [26] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe.  Poisson surface reconstruction.  In Proceedings of the fourth Eurographics symposium on Geometry processing, volume 7, page 0, 2006.
  * [27] Jiarong Lin, Chongjiang Yuan, Yixi Cai, Haotian Li, Yuying Zou, Xiaoping Hong, and Fu Zhang.  Immesh: An immediate lidar localization and meshing framework.  arXiv preprint arXiv:2301.05206, 2023.
  * [28] Soohwan Kim and Jonghyuk Kim.  GPmap: A Unified Framework for Robotic Mapping Based on Sparse Gaussian Processes, pages 319–332.  Springer International Publishing, Cham, 2015.
  * [29] Bhoram Lee, Clark Zhang, Zonghao Huang, and Daniel D Lee.  Online continuous mapping using gaussian process implicit surfaces.  In 2019 International Conference on Robotics and Automation (ICRA), pages 6884–6890. IEEE, 2019.
  * [30] Kenji Koide, Masashi Yokozuka, Shuji Oishi, and Atsuhiko Banno.  Voxelized gicp for fast and accurate 3d point cloud registration.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 11054–11059. IEEE, 2021.
  * [31] Arno Knapitsch, Jaesik Park, Qian-Yi Zhou, and Vladlen Koltun.  Tanks and temples: Benchmarking large-scale scene reconstruction.  ACM Transactions on Graphics (ToG), 36(4):1–13, 2017.
  * [32] Andreas Geiger, Philip Lenz, Christoph Stiller, and Raquel Urtasun.  Vision meets robotics: The kitti dataset.  The International Journal of Robotics Research, 32(11):1231–1237, 2013.