# Phase-contrast THz-CT for non-destructive testing

Peter Fosodeder  1,* Simon Hubmer  2 Alexander Ploier  2,3 Ronny Ramlau  2,4
Sandrine van Frank1 and Christian Rankl1  1Research Center for Non Destructive
Testing GmbH (RECENDT), Altenbergerstraße 69, A-4040 Linz, Austria
2Johann Radon Institute Linz, Altenbergerstraße 69, A-4040 Linz, Austria
3Doctoral Program Computational Mathematics, Johannes Kepler University Linz,
Altenbergerstraße 69, A-4040 Linz, Austria
4Institute of Industrial Mathematics, Johannes Kepler University Linz,
Altenbergerstraße 69, A-4040 Linz<EMAIL_ADDRESS>
###### Abstract

A new approach for image reconstruction in THz computed tomography (THz-CT) is
presented. Based on a geometrical optics model containing the THz signal
amplitude and phase, a novel algorithm for extracting an average phase from
the measured THz signals is derived. Applying the algorithm results in a
phase-contrast sinogram, which is further used for image reconstruction. For
experimental validation, a fast THz time-domain spectrometer (THz-TDS) in
transmission geometry is employed, enabling CT measurements within several
minutes. Quantitative evaluation of reconstructed 3D printed plastic profiles
reveals the potential of our approach in non-destructive testing of plastic
profiles.

††journal: oe††articletype: Research Article

##  1 Introduction

Over the past years, THz technology has evolved into a widely recognized
technology for industrial non-destructive testing (NDT), its most prominent
application being the thickness control of thin coating layers [1, 2, 3].
Equally important and well-established is the application of THz spectroscopy
[4, 5, 6] in chemical process control [7] or even in security applications
[8]. In this work, we concentrate on one of the less established NDT methods,
namely THz computed tomography (THz-CT). The following paragraphs contain a
short motivation for the use of THz-CT, a description of the state-of-the-art,
current challenges and a description of our approach for extending the
applicability of THz-CT.

Computed tomography (CT) is an imaging technique which uses multiple
transmission measurements from different angles to produce tomographic (cross-
sectional) images. Most prominently, X-radiation [9] or ultrasound [10, 11] is
used to probe the sample. We are particularly interested in using THz
radiation for CT, since it provides industrially relevant advantages over the
complementary NDT methods mentioned above. THz radiation is completely
harmless to humans and does not require extensive shielding, which is
necessary e.g. in X-Ray based systems. In addition, THz measurements are
operated contact-less, without the need for coupling media, which are required
in ultrasound testing. In the past, limitations concerning the measurement
speed represented a major obstacle for the use of THz systems in industrial
environments. However, the current development of fast and compact measurement
systems [12] provides a suitable solution to this challenge. We propose a
potential application of THz-CT for in-line quality control of plastic
profiles. Applied within such a major industrial sector, the system
demonstrated in this work could provide immediate feedback regarding the
quality of produced plastic profiles and therefore effectively reduce the
amount of produced plastic waste in case of production failure.

Due to demands for high output power or the early limitations of measurement
speed of pulsed systems, existing THz-CT approaches are often based on
continuous wave (cw) source [13, 14] and detector technology. However, one
major drawback regarding the use of cw systems is their low operating
frequency below several hundred GHz, resulting in a low diffraction limited
resolution. Furthermore, conventional cw systems are typically restricted to
imaging based on the transmitted intensity only. On the other hand, approaches
based on pulsed technology [15, 16, 17, 18] provide access to the THz pulse
phase and allow imaging with frequencies up to several THz, thus providing
more comprehensive information about the sample.

Regarding the reconstruction of cross-sectional images, most of the existing
approaches, including the ones reported in the literature (cw and pulsed)
cited above, utilize methods developed for X-Ray CT. These methods are
typically based on the Radon imaging model [19] and thus do not account for
the physical effects of electromagnetic wave propagation like scattering,
diffraction and refraction. Approaches to reconstruct images based on more
generally valid imaging models such as Snell’s law of refraction [20] or
Helmholtz’s equation [21] have been reported to improve image reconstruction.
This however comes at the cost of having to introduce preliminary knowledge
about the sample geometry and/or a significantly increased computational
effort.

In this paper, we present an approach for THz-CT based on a fast THz time-
domain spectrometer. Based on geometrical optics, we formulate the imaging
model and present two methods for calculating a sinogram from the measured THz
signals. Analogously to X-Ray CT, the first method utilizes amplitude
information only, while the second method exploits the fact that in THz-TDS
the signal is sampled with sub-picosecond resolution, which allows us to
extract a pulse phase sinogram. Cross-sectional images are reconstructed from
such sinograms using different image reconstruction approaches, namely
filtered backprojection, Tikhonov regularization and Landweber iteration [19,
22, 23]. The performance of our approach was tested by measuring 3D-printed
plastic profiles and quantitatively comparing the outcome of the above
procedure with the known geometry.

The goal of this article is to demonstrate the robustness, accuracy and
therefore the potential applicability of our THz-CT approach to in-line NDT of
extruded plastic profiles.

##  2 Materials and methods

###  2.1 The THz-TDS setup

Our setup is based on the TeraFlash smart [12] THz-TDS system (TOPTICA,
Munich). This system uses two photoconductive antennas (PCA) and two
synchronized femtosecond lasers, operating at a center wavelength of 1560 nm
for generating and detecting pulsed THz radiation.

As illustrated in Fig. 1, the THz pulses are guided and focused by 4 off-axis
parabolic mirrors (OPM), and the sample is mounted in the focal plane. Sample
positioning for the CT measurement is performed by one translation and one
rotation stage. For imaging in transmission it is crucial to find a reasonable
trade-off between depth of field and focal spot-size, therefore the effective
focal length of the focusing OPM was dimensioned with 101.6 mm. A Gaussian
focal spot with a FWHM of approximately 1 mm (measured by the knife-edge
technique) was achieved while still preserving enough depth of field for
sample diameters of several centimeters.

Fig. 1 also illustrates two THz signals, one reference measurement in air and
one measurement on a 3D printed plate of Polypropylene (PP) with a thickness
of 2 mm. Typically, the bandwidth of a single shot THz reference signal
exceeds 3.5 THz at a measurement rate of 1600 signals per second. The total
length of one THz signal is 152.3 ps, sampled with a resolution of $\approx$
49 fs.

Figure 1: Schematic drawing of the measurement setup. The THz radiation is
generated by a PCA (Tx). Two Off-Axis Parabolic Mirrors (OPM) are used to form
a focused THz beam. After interacting with an object, the THz beam is guided
to the detector PCA (Rx) by two OPMs again. Sample movement is performed by a
motorized translation and rotation stage. Exemplarily, a measured reference
signal through air and a measured signal through a 2 mm plate of PP is shown.

###  2.2 Sample preparation

The samples for this study were 3D printed from natural polypropylene (PP) at
a temperature of 240 °C. Due to minor visible structural inhomogenities from
the 3D printers layer resolution of 150 µm, a THz spectroscopy experiment on a
printed plate of PP was performed. The results in Fig. 2 were calculated,
using a frequency-domain model-based parameter estimation [5, 6]. Despite the
visible surface structure of the printed samples, no significant absorption
features compared to other THz spectroscopy experiments performed on PP were
found [4].

Figure 2: Material parameters of a 3D printed plate of PP (thickness = 2 mm).
The absorption coefficient $\alpha$ and the refractive index $n$ were
calculated from a THz-TDS measurement in transmission. No significant
absorption features caused by potential structural inhomogenities were found.

###  2.3 Computed tomography (CT) measurements

In conventional X-Ray CT measurements, detector arrays are used in order to
minimize acquisition time. Since our THz-TDS setup employs two PCAs,
effectively acting as a static point-like emitter/detector pair, we perform
continuous line scans by moving the sample through the THz focus, at typical
velocities of 150 mm/s. After each linescan, the sample is rotated by 1° until
the full 360° scan is completed. Given the typical linescan velocity of 150
mm/s and a measurement rate of 1600 THz signals per second, a spatial
discretization smaller than 0.1 mm and an angular discretization of 1° could
be achieved in a total acquisition time of less than 5 minutes per scan.

###  2.4 Data processing algorithms

In order to improve the signal-to-noise ratio (SNR) of the acquired THz
signals, a digital FIR low-pass filter [24] with a passband frequency of 2.5
THz and a stopband frequency of 3.5 THz (attenuation: -80 dB) was applied to
the raw data.

Consecutively acquired THz signals were filtered by a moving average filter
with a window size of five signals. This effectively corresponds to a spatial
low-pass filtering of the dataset and therefore further reduces the SNR.

Due to the continuous motion of the sample, the acquired measurement positions
are available on non-equidistantly distributed positions only. After
calculating the sinogram values on the corresponding positions (see section
3.1), the values were interpolated on an equidistant grid with a spatial
discretization of 0.1 mm and an angular discretization of 1°.

After image reconstruction, the retrieved images were thresholded in order to
remove values smaller than zero. These values typically appear due to
numerical reasons and are nonphysical.

##  3 The THz imaging model

The basis for image reconstruction is an accurate and ideally compact forward
model, that will be derived within this section. In its most general form, the
propagation of electromagnetic radiation through an object is described by the
macroscopic Maxwell equations [25]. In this framework, the common phenomena of
wave propagation, such as refraction and diffraction, are considered.
Recently, efforts have been taken in order to perform THz CT image
reconstruction within this framework [21].

Since our work is specifically dedicated to 3D imaging of plastic profiles,
some assumptions will be introduced in the following, leading to a more
compact imaging model for our specific problem. The basis of our imaging model
is the well-known regime of geometrical optics, where a light beam is treated
as a mathematical ray, undergoing reflection, refraction and absorption by the
sample. In this regime, the material-light interaction is commonly assumed to
be linear and isotropic. Furthermore, the appearance of multiple reflections
within a sample is neglected, since they result only in minor contributions to
the measured signal.

Since the samples in this work represent a special case of sparsely populated
objects, we introduce a significant simplification to the imaging model by
neglecting diffraction and refraction in order to preserve a reasonably
complicated framework for computationally efficient image reconstruction.

Using the assumptions introduced above, the remaining radiation-matter
interaction can be fully described by the complex and frequency dependent
refractive index of the material

$\widetilde{n}\left(\omega\right)=n+i\kappa\,,$ (1)

where $n$ denotes the real refractive index and $\kappa$ denotes the
absorption index. With the notation in Eq. (1), the measured THz signal, after
travelling a distance $d$ through the object, can be written as

$E\left(d,t\right)=\int\limits_{-\infty}^{+\infty}\hat{E}^{\mathrm{ref}}\left(\omega\right)\exp\left[i\left(\omega
t-\frac{\omega}{c}\left(\widetilde{n}\left(\omega\right)-n_{0}\right)d\right)\right]\mathrm{d}\omega\,.$
(2)

Here $c$ refers to the speed of light in a vacuum, $n_{0}$ denotes the real
refractive index of air and $\hat{E}^{\mathrm{ref}}$ is the Fourier transform
of the reference signal $E^{\mathrm{ref}}$ measured in air.

The complex refractive index $\widetilde{n}$ is commonly approximated by its
Taylor series

$\widetilde{n}\left(\omega\right)=\widetilde{n}\left(\omega_{0}\right)+\mathrm{\frac{d}{d\omega}}\widetilde{n}\left(\omega\right)|_{\omega_{0}}\left(\omega-\omega_{0}\right)+\mathcal{O}\left(\omega^{2}\right)\,.$
(3)

A comparison with Fig. 2 shows that $n$ is only marginally varying over the
frequency bandwidth of our measurement system and is therefore sufficiently
approximated by a constant 0-th order term
$\overline{n}=n\left(\omega_{0}\right)$.

We will further rewrite the absorption index $\kappa\left(\omega\right)$ by
using the absorption coefficient
$\alpha\left(\omega\right)=2\kappa\left(\omega\right)\frac{\omega}{c}$. Fig. 2
shows that the frequency dependence of the absorption coefficient is low in
contrast to other dielectric materials (e.g. ABS, PET [4]). Low material
parameter dispersion typically implies that the THz pulse shape is not
significantly altered after propagating through the material. In Fig. 1 a
comparison of a reference pulse and a pulse after travelling through 2 mm of
PP is shown, supporting the statement that the influence of dispersion is
negligibly low. We therefore introduce the constant absorption coefficient
$\overline{\alpha}=\alpha\left(\omega_{0}\right)$ as a further simplification
to our model.

Inserting the average quantities $\overline{n}$ and $\overline{\alpha}$ into
Eq. (2) leads to

$E\left(d,t\right)=\int\limits_{-\infty}^{+\infty}\hat{E}^{\mathrm{ref}}\left(\omega\right)\exp\left[-i\frac{\omega}{c}\left(\overline{n}-n_{0}\right)d\right]\exp\left(-\frac{\overline{\alpha}}{2}d\right)\exp\left(i\omega
t\right)\mathrm{d}\omega\,.$ (4)

The approximations introduced with $\overline{n}$ and $\overline{\alpha}$
allow us to transform our model back to the time-domain analytically. A more
general approach would be to proceed with the frequency domain model in Eq.
(4) in order to account for the frequency dependent material parameters.

We leave this as subject to future work and continue by performing the inverse
Fourier transform in equation 4. This leads to the time-domain representation
of the model

$E\left(d,t\right)=E^{\mathrm{ref}}\left(t-\frac{\overline{n}-n_{0}}{c}d\right)\exp\left(-\frac{\overline{\alpha}}{2}d\right)\,.$
(5)

Eq. (5) describes the interaction of a mathematical THz ray with a material of
thickness $d$ and material parameters $\overline{n}$ and $\overline{\alpha}$.
One can see that the interaction between sample and THz pulse is only
expressed by a time-delay of the THz pulse with respect to the reference
pulse, and an absorption term according to Lambert-Beers law.

In order to describe the CT measurement procedure (see section 2.3), we
introduce the discrete sample positions $s_{i}$ and projection angles
$\theta_{j}$ for each point measurement, and the density function
$f\left(x\right)$ representing the spatially dependent absorption coefficient.
Analogously to our previous work [26], we now introduce the Radon transform
$R$ [19] which is defined as the line integral along the line $L_{i,j}$

$\left(Rf\right)\left(s_{i},\theta_{j}\right)=\int_{L_{i,j}}f\left(x\right)\mathrm{d}x\,,$
(6)

where the line $L_{i,j}$ is characterized by the sample position $s_{i}$ and
the projection angle $\theta_{j}$.

In our specific use-case of extruded plastic profiles, the samples consist of
only one material and air. The density function $f$ in Eq. (6) is therefore
piecewise constant with values either 0 or $\overline{\alpha}$. For a given
$\left(s_{i},\theta_{j}\right)$, the Radon transform of $f$ is therefore
proportional to the sample thickness $d_{i,j}$ along the line $L_{i,j}$

$\left(Rf\right)\left(s_{i},\theta_{j}\right)=\bar{\alpha}d_{i,j}\,.$ (7)

We can use Eq. (7) to rewrite Eq. (5) for the time domain THz-CT model

$E_{i,j}\left(t\right)=E^{\mathrm{ref}}\left(t-\frac{\bar{n}-n_{0}}{c\overline{\alpha}}\left(Rf\right)\left(s_{i},\theta_{j}\right)\right)\exp\left(-\frac{1}{2}\left(Rf\right)\left(s_{i},\theta_{j}\right)\right)\,,$
(8)

where $E_{i,j}$ denotes the final measurement signal in time-domain.

In the following, we deduce two different approaches for mapping each
measurement signal $E_{i,j}(t)$ to a scalar sinogram value $g_{i,j}$ that acts
as the basis for the image reconstruction approaches shown in section 4.

###  3.1 Pulse energy sinogram mapping

In the conventional approach of X-Ray CT, information about the sample is
retrieved from the transmitted intensity of X-Ray radiation $I_{i,j}$ with
respect to a reference intensity $I^{\mathrm{ref}}$. The sinogram values for
X-Ray CT are given according to
$g_{i,j}^{X}=-\ln\left[I_{i,j}/I^{\mathrm{ref}}\right]$.

Analogously, the energy of a THz pulse is defined as
$\int_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}\,\mathrm{d}t$. By
calculating the energy from Eq. (8) we find

$g_{i,j}^{E}:=-\ln\left[\frac{\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}\,\mathrm{d}t}{\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t\right)^{2}\,\mathrm{d}t}\right]=\left(Rf\right)\left(s_{i},\theta_{j}\right)\,,$
(9)

with $g_{i,j}^{E}$ denoting the pulse energy sinogram values.

###  3.2 Pulse phase sinogram mapping

The approach presented above extracts information about the sample based on
the signal amplitude. Since the measurement setup is designed to measure the
time evolution of the electric field of the THz pulse, it is crucial to find a
way to access the phase information as well. According to our model in Eq.
(8), the phase of the THz pulse is given by the time-shift term
$\frac{\bar{n}^{\prime}-n_{0}}{c\overline{\alpha}}\left(Rf\right)\left(s_{i},\theta_{j}\right)$
.

Standard signal processing offers a straightforward way to calculate the time-
shift, for instance by means of peak-finder algorithms and thresholding.
However, the resulting sinogram is non-continuous, which can lead to artifacts
in the reconstructed image. Furthermore, this approach does not offer a robust
way of dealing with THz signals containing multiple peaks. In practice, these
signals typically appear due to the finite cross-section of the THz beam
(FWHM: 1 mm) travelling through differently thick parts of the sample
simultaneously. We have recently discussed this effect in [26] by means of a
more complex, non-linear imaging model.

In this work, we therefore proceed with our more compact model of the
mathematical THz ray and extract the time-shift term in Eq. (8)

$\triangle
t:=\frac{\bar{n}-n_{0}}{c\overline{\alpha}}\left(Rf\right)\left(s_{i},\theta_{j}\right)\,.$
(10)

Calculating the quantity
$\int_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}t\,\mathrm{d}t$ from Eq. (8)
by substituting $t^{\prime}=t-\triangle t$ results in

$\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}t\,\mathrm{d}t=\left[\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t^{\prime}\right)^{2}t^{\prime}\,\mathrm{d}t^{\prime}+\triangle
t\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t^{\prime}\right)^{2}\,\mathrm{d}t^{\prime}\right]\exp\left[\left(Rf\right)\left(s_{i},\theta_{j}\right)\right]\,.$
(11)

The remaining integrals in Eq. (11) are independent of $f$ and calculated from
the reference measurement. From Eq. (9) we find that dividing Eq. (11) by
$\int_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}\,\mathrm{d}t$ and
rearranging terms results in

$\triangle
t=\frac{\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}t\,\mathrm{d}t}{\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}\,\mathrm{d}t}-\frac{\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t\right)^{2}t\,\mathrm{d}t}{\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t\right)^{2}\,\mathrm{d}t}\,.$
(12)

Reinserting the definition Eq. (10) back into Eq. (12) and rearranging terms
leads to the pulse phase sinogram values

$g_{i,j}^{P}:=\frac{c\overline{\alpha}}{\bar{n}-n_{0}}\left[\frac{\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}t\,\mathrm{d}t}{\int\limits_{-\infty}^{+\infty}E_{i,j}\left(t\right)^{2}\,\mathrm{d}t}-\frac{\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t\right)^{2}t\,\mathrm{d}t}{\int\limits_{-\infty}^{+\infty}E^{\mathrm{ref}}\left(t\right)^{2}\,\mathrm{d}t}\right]=\left(Rf\right)\left(s_{i},\theta_{j}\right)\,.$
(13)

The latter term in Eq. (12) is calculated from the reference measurement and
is therefore constant, while the material parameters cause a linear scaling of
the sinogram values.

The advantage of this procedure over standard peak finding algorithms is its
robustness and the absence of tuning parameters (e.g. thresholds).
Furthermore, integrating over the THz signals has a time-domain averaging
effect, leading to a continuous sinogram. In addition, it automatically takes
into account the case of multiple peaks appearing in one THz signal, as
discussed above.

##  4 Image reconstruction approaches

In this section, we shortly introduce the different image reconstruction
approaches used throughout this work. Reconstructing based on both the pulse
energy and the pulse phase mappings in Eq. (9) and (13) amount to solving a
linear system

$(Rf)(s_{i},\theta_{j})=g_{i,j}\,,$ (14)

for the unknown density $f$, with right hand sides $g_{i,j}=g_{i,j}^{E}$ and
$g_{i,j}=g_{i,j}^{P}$, respectively. Since Eq. (14) is nothing else than a
discretized version of the continuous Radon transform equation [19, 27]

$(Rf)(s,\theta)=g(s,\theta)\,,$ (15)

standard reconstruction methods for this equation can be adapted to solve our
THz-CT problem. Perhaps the most well-known of these methods is filtered back-
projection, which is based on the classic Radon inversion formula [19, 27].
The filtering step in this method acts as a regularization, which is necessary
since the problem is in general ill-posed and unstable with respect to
unavoidable measurement noise in the data. Another popular approach for
solving Eq. (15) is Tikhonov regularisation [23], which determines a stable
approximation of $f$ by minimizing

$||Rf-g||^{2}+\beta||f||^{2}\,,$ (16)

where $\beta$ denotes a regularization parameter. The minimizer can be
computed explicitly, in our discrete case as the solution of a linear matrix-
vector equation. Alternatively, a common iterative approach for solving (15)
is given by Landweber iteration [23], which is defined by

$f_{k+1}=f_{k}+\gamma R^{*}(g-Rf_{k})\,,$

where $\gamma$ denotes a step-size parameter. Here, the choice of stopping
index acts as regularization.

For our THz-CT problem we focus on the above three reconstruction methods,
i.e., on filtered back-projection, Tikhonov regularization, and Landweber
iteration. For further details on these methods, as well as for other possible
reconstruction methods we refer to [23, 27].

##  5 Results

In order to compare the outcomes of the two different sinogram mapping methods
and reconstruction approaches presented above, we comparatively walk through
and analyze the results obtained for the sample in Fig. 3. Further on, more
reconstructed sample cross-sections, retrieved with the most promising
combination of sinogram mapping and reconstruction approach, are shown below
and in the additionally provided supplementary material (see Supplement 1).

Figure 3: 3D printed sample for THz CT measurements: 3D representation of a
circular sample with the letters "THZ" inscribed (left). Cross-section of the
sample (right). Cross-hatched areas represent solid walls of PP.

###  5.1 Qualitative verification

The purpose of the following measurements is to compare different methods of
sinogram mapping. After choosing the most efficient method, different image
reconstruction approaches are qualitatively evaluated. For this, the sample in
Fig. 3 was scanned over a range of 0 - 359° and the raw measurement data was
preprocessed according to section 2.3.

In section 3.1 and 3.2, two different approaches for extracting information
from THz signals are discussed. For both approaches, the calculated sinograms
are shown in Fig. 4 (a, c). The reconstructed cross-sections in Fig. 4 (b, d)
were obtained by applying the classical filtered backprojection algorithm to
the sinograms. Those image reconstructions were implemented using the built-in
MATLAB [28] function iradon() with default settings. The calculations were
performed on a standard desktop PC (Intel i7-4790 CPU with 16 GB RAM) with a
typical computation time of one second per reconstruction.

Figure 4: Comparison of sinogram mapping methods. (a) Pulse energy sinogram,
(b) FBP of pulse energy sinogram, (c) Pulse phase sinogram, (d) FBP of pulse
phase sinogram. The reconstructed image (d) obtained by pulse phase sinogram
mapping (c) shows a significantly better representation of the sample geometry
than conventional pulse energy based FBP (b).

Using the pulse phase sinogram given in Fig. 4 (c), the different image
reconstruction approaches discussed in section 4 were tested in order to give
a qualitative comparison between the filtered back-projection, Tikhonov
regularization and Landweber iterations. The latter two were implemented in
MATLAB using the AIR TOOLS II toolbox [29] for assembling the Radon transform
matrix. Fig. 5 shows the results of Tikhonov regularization (left) and
Landweber iteration (right). These were obtained with the choice $\beta=3700$
for the Tikhonov regularisation parameter and $\gamma=3.8\cdot 10^{-6}$ for
the step-size in Landweber iteration, which was stopped after $70$ iterations.
These parameters were optimized manually for minimal noise and maximum image
sharpness. The total computation time was in the range of one minute on the
same desktop PC as mentioned above. As no effort was made to increase
computational efficiency, there is still potential to speed up the
reconstruction algorithm.

Figure 5: Comparison of image reconstruction approaches: (a) Tikhonov
regularization, (b) Landweber iteration. Both reconstruction approaches yield
similar results, with Landweber iteration (b) resulting in a more homogeneous
representation of the letters THZ.

###  5.2 Quantitative evaluation of wall thicknesses

For quantitative validation of our THz-CT measurements, we designed and
measured two further samples with more ordinary geometries. Due to the sample
design, scattering artifacts in the reconstruction are to be expected only due
to 3D-printed inhomogenities (e.g. layer stacks, joints).

While the geometry of the first object in Fig. 6 (top) is parameterized by a
spiral with linearly increasing radius, the second object (bottom) is highly
symmetrical and consists of planar elements only.

Both images were reconstructed using Landweber iteration with the same
parameters as in section 5.1. As a reference measurement for both objects, the
nominal wall thickness (2 mm) was verified within minor deviations of $\pm
0.1\,\mathrm{mm}$ by a mechanical caliper. The contours of the reconstructed
objects were determined by calculating the iso-line (blue line) at the half
maximum value of the reconstructed image.

The wall thickness of the spiral was measured at discrete positions indicated
by green lines. A variing thickness between 2.0 - 2.6 mm was found.

Due to the linear geometry of the second object, we applied the linear Hough
transform (LHT) to the contour of the object in order to identify a valid
boundary along each sample wall. As a constraint for the LHT, only vertical
and horizontal lines with a minimum length of 20 mm were accepted as valid
object boundaries. The determined thicknesses for the three horizontal and
vertical wall segments are marked in Fig. 6 (b).

Figure 6: Quantitative evaluation of reconstructed images: Sample cross-
sections (left), Quantitative evaluation of a curved sample (a), Quantitative
evaluation of a planar sample (b). The blue hatched area in the sample cross-
section represents the sample walls. Blue lines in the reconstructed image
represent the calculated sample contour and green lines mark the identified
sample boundaries using the Hough transform. Sample a) shows a gradually
increasing deviation of the wall thickness towards its center of rotation.
Overall good agreement within the accuracy of the 3d printer was achieved for
sample b).

##  6 Discussion

We begin this section with a qualitative comparison of the sinogram mapping
methods in Fig. 4. In both cases (b, d), the cylinder around the letters THZ
is well identifiable. Since the cylinder is 3D printed from multiple circles
next to and on top of each other, an inhomogenious joint/stack of layers is
produced. This joint can be identified in both reconstructed images at the
position
$\left(y|z\right)\approx\left(55\,\mathrm{mm\,}|\,45\,\mathrm{mm}\right)$.

Generally, the pulse energy sinogram (a) and reconstruction (b) yields more
pronounced edges and inhomogenities than the pulse phase sinogram. We presume
that this is caused by the strong influence of scattering effects on small
sample structures and edges, effectively reducing the detected signal
amplitude. On the other hand, the pulse phase mapping approach eliminates the
influence of the absolute signal amplitude by normalizing each measured signal
with respect to its energy (see Eq. 13). This property effectively suppresses
the presence of large peaks in the reconstructed image (d), and leads to a
better representation of the sample geometry in general.

Furthermore, it is worth noting that the thickness of the letters in Fig. 4
(d) is significantly larger than the cylinders wall thickness, despite both
being printed with a thickness of $2.0\pm 0.1\,\mathrm{mm}$ (measured with a
mechanical caliper). We assume that this is a consequence of refraction in the
cylinder, and point to further experimental evidence given below.

In general, for sparse sample structures with sharp edges we prefer the method
of pulse phase mapping for geometry identification. For more specific
applications, as for example detection of small inhomogenities in the
dimension of the THz wavelength, pulse energy mapping might be superior.

In addition to the filtered backprojection image reconstruction method, we
have qualitatively shown that Tikhonov regularization and Landweber iterations
yield similar results (see Fig. 5). Especially within the letters THZ we
observe a much more homogeneous structure and less artifacts for Tikhonov
regularization and Landweber iteration. After manual parameter optimization,
both reconstruction approaches seem to outperform the classical approach of
filtered backprojection.

We have further examined two more samples for quantitative testing of the
accuracy of the proposed THz-CT approach. In Fig. 6 (a) we have observed a
locally increasing wall thickness of up to 2.6 mm. The tendency indicates that
the amount of surrounding sample structure layers artificially increases the
wall thickness in the reconstructed image. Throughout this work, this effect
was consistently and explicitly observed for all samples with curved
surrounding structures (see Supplement 1). Thus, we strongly assume that this
artifact is predominantly caused by the wave propagation effect of refraction.
According to geometrical optics, refraction on parallel, planar surfaces only
applies a spatial shift to a ray, while refraction on a curved surface
entirely changes the direction of the ray, thus introducing a larger
inaccuracy of our forward model. So far, we have not accounted for this effect
in our model. Nevertheless, the spiral structure is qualitatively well
identifiable and we therefore plan to use the current approach as the basis
for further improvements in the image reconstruction of curved objects. In the
industrially more relevant case of a planar profile (see Fig. 6, b), excellent
agreement with the nominal wall thickness of $2.0\pm 0.1\,\mathrm{mm}$ was
achieved. The high accuracy of this result shows that our THz-CT approach is
capable of detecting structural features down to industrially relevant
dimensions.

Previously reported results obtained by THz-CT often involve samples with a
very simple geometry (e.g. filled cylinders [30]). Measurements performed on
slightly more complex samples [18] often show unnatural structural artifacts
(e.g. non-uniform wall thicknesses, unnatural edges), besides the effect of
broadened walls, which is also reported in this work. Other reported systems
[17] are very slow, resulting in long data acquisition times or poor image
quality and thus are not useable in an industrial environment. To our
knowledge, the THz-CT system presented here outperforms the systems reported
in the literature in terms of data acquisition time, image quality and
quantitative accuracy versus sample complexity, and robustness for different
sample geometries.

##  7 Conclusion

In this work, we have presented two computationally efficient and robust
approaches for performing THz-CT measurements on plastic profiles by means of
a fast THz-TDS system. While the first approach, in the style of classical
X-Ray CT, utilizes the amplitude information of the measured signals for image
reconstruction, the second approach makes use of the availability of phase
information in the THz-TDS measurement. We have experimentally shown that THz-
CT based on the signal amplitude is prone to wave scattering artifacts.
Structures with sizes in the dimension of the wavelength or sharp edges lead
to extended peaks in the reconstructed image, potentially overshadowing other
structures. In contrast, the influence of scattering is mostly eliminated in
the second approach, leading to overall better representations of the sample
geometry.

We have further compared the performance of multiple image reconstruction
methods, namely filtered backprojection, Tikhonov regularization, and
Landweber iteration. It was empirically found that Tikhonov regularization and
Landweber iteration yield more homogeneous images than filtered
backprojection.

As a last point, the quantitative performance of our THz-CT system was
evaluated using the example of a curved and a planar sample (see Fig. 6).
Local deviations were found for the special case of the curved sample. The
planar sample could be reconstructed correctly within an accuracy of $\pm
0.1\,\mathrm{mm}$. Due to the successful measurement on the planar sample and
further collected experimental data (see Supplement 1), we conclude that the
artificially increased wall thickness observed for curved samples is a
consequence of refraction. However, the quantitative accuracy of the
reconstruction on the planar profile proves that our THz-CT approach is a
valid tool for NDT on plastic profiles.

Funding Österreichische Forschungsförderungsgesellschaft (871974); European
Commission (777222); Austrian Science Fund (F6805-N36).

Acknowledgments Part of this work was performed under the scope of the COMET
program within the research project "Photonic Sensing for Smarter Processes".
This program is promoted by BMK, BMDW, the federal state of Upper Austria and
the federal state of Styria, represented by SFG. This project was supported by
the strategic economic- and research program "Innovative Upper Austria 2020"
of the province of Upper Austria. "This work was supported by the Austrian
Science Fund (FWF) in the project F6805-N36 within the Special Research
Programme SFB F68: “Tomography Across the Scales”

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

Supplemental document See Supplement 1 for supporting content.

## References

  * [1] S. Krimi, J. Klier, J. Jonuscheit, G. v. Freymann, R. Urbansky, and R. Beigang, “Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology,” Applied Physics Letters 109, 021105 (2016).
  * [2] W. Tu, S. Zhong, Y. Shen, and A. Incecik, “Nondestructive testing of marine protective coatings using terahertz waves with stationary wavelet transform,” Ocean Engineering 111, 582–592 (2016).
  * [3] K. Su, Y. Shen, and J. A. Zeitler, “Terahertz sensor for non-contact thickness and quality measurement of automobile paints of varying complexity,” IEEE Transactions on Terahertz Science and Technology 4, 432–439 (2014).
  * [4] Y. Jin, G. Kim, and S. Jeon, “Terahertz dielectric properties of polymers,” Journal of the Korean Physical Society 49, 513–517 (2006).
  * [5] I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with thz time-domain spectroscopy,” Opt. Express 15, 4335–4350 (2007).
  * [6] M. Scheller, C. Jansen, and M. Koch, “Analyzing sub-100-µm samples with transmission terahertz time domain spectroscopy,” Optics Communications 282, 1304–1306 (2009).
  * [7] C. J. Strachan, T. Rades, D. A. Newnham, K. C. Gordon, M. Pepper, and P. F. Taday, “Using terahertz pulsed spectroscopy to study crystallinity of pharmaceutical materials,” Chemical Physics Letters 390, 20–24 (2004).
  * [8] Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Applied Physics Letters 86, 241116 (2005).
  * [9] S. C. Garcea, Y. Wang, and P.J.Withers, “X-ray computed tomography of polymer composites,” Composites Science and Technology 156, 305–319 (2018).
  * [10] H. Choi and J. S. Popovics, “NDE application of ultrasonic tomography to a full-scale concrete structure,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 62, 1076–1085 (2015).
  * [11] J. Wiskin, B. Malik, D. Borup, N. Pirshafiey, and J. Klock, “Full wave 3D inverse scattering transmission ultrasound tomography in the presence of high contrast,” Sci. Rep. 10, 20166 (2020).
  * [12] R. J. B. Dietz, N. Vieweg, T. Puppe, A. Zach, B. Globisch, T. Göbel, P. Leisching, and M. Schell, “All fiber-coupled THz-TDS system with kHz measurement rate based on electronically controlled optical sampling,” Optics Letters 39, 6482–5 (2014).
  * [13] B. Recur, J. P. Guillet, L. Bassel, C. Fragnol, I. Manek-Hönninger, J. Delagnes, W. Benharbone, P. Desbarats, J. Domenger, and P. Mounaix, “Terahertz radiation for tomographic inspection,” Optical Engineering 51, 1–8 (2012).
  * [14] G. Trichopoulos and K. Sertel, “Broadband terahertz computed tomography using a 5k-pixel real-time THz camera,” Journal of Infrared, Millimeter, and Terahertz Waves 36 (2015).
  * [15] J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (THz-TDS),” Journal of Applied Physics 124, 231101 (2018).
  * [16] B. Ferguson, S. Wang, D. Gray, D. Abbot, and X.-C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002).
  * [17] E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Optics Communications 283, 2050–2055 (2010).
  * [18] S. R. Tripathi, Y. S., K. Murate, K. Imayama, and K. Kawase, “Terahertz wave three-dimensional computed tomography based on injection-seeded terahertz wave parametric emitter and detector,” Opt. Express 24, 6433–6440 (2016).
  * [19] F. Natterer, _The Mathematics of Computerized Tomography_ (Society for Industrial and Applied Mathematics, 2001).
  * [20] J. Tepe, T. Schuster, and B. Littau, “A modified algebraic reconstruction technique taking refraction into account with an application in terahertz tomography,” Inverse Problems in Science and Engineering (2016).
  * [21] A. Wald and T. Schuster, “Tomographic terahertz imaging using sequential subspace optimization,” in _New Trends in Parameter Identification for Mathematical Models,_ B. Hofmann, A. Leitão, and J. Zubelli, eds. (Birkhäuser, 2018).
  * [22] M. L. Louarn, C. Vérinaud, V. Korkiakoski, N. Hubin, and E. Marchetti, “Adaptive optics simulations for the European Extremely Large Telescope,” in _in Advances in Adaptive Optics II no. 6272 in Proc. SPIE,_ (2006), p. 627234.
  * [23] H. W. Engl, M. Hanke, and A. Neubauer, _Regularization of inverse problems_ (Kluwer Academic Publishers, 1996).
  * [24] B. A. Shenoi, _Finite Impulse Response Filters_ (John Wiley & Sons, Ltd, 2005), chap. 5, pp. 249–302.
  * [25] D. J. Griffiths, _Introduction to electrodynamics_ (Pearson, 2014), 4th ed.
  * [26] S. Hubmer, A. Ploier, R. Ramlau, P. Fosodeder, and S. van Frank, “A mathematical approach towards THz tomography for non-destructive imaging,” submitted, available from: https://arxiv.org/abs/2010.14938 (2021).
  * [27] A. K. Louis, _Inverse und schlecht gestellte Probleme_ , Teubner Studienbücher Mathematik (Vieweg+Teubner Verlag, 1989).
  * [28] MATLAB, _9.1.0.441655 (R2016b)_ (The MathWorks Inc., Natick, Massachusetts, 2016).
  * [29] P. C. Hansen and J. Jorgensen, “Air tools ii: algebraic iterative reconstruction methods, improved implementation,” Numerical Algorithms 79 (2017).
  * [30] S. Mukherjee and J. Federici, “Study of structural defects inside natural cork by pulsed terahertz tomography,” in _2011 International Conference on Infrared, Millimeter, and Terahertz Waves,_ (2011).


# Efficient Map Sparsification Based on 2D and 3D Discretized Grids

Xiaoyu Zhang Yun-Hui Liu
T Stone Robotics Institute, Chinese University of Hong Kong
Hong Kong Centre for Logistics Robotics
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Localization in a pre-built map is a basic technique for robot autonomous
navigation. Existing mapping and localization methods commonly work well in
small-scale environments. As a map grows larger, however, more memory is
required and localization becomes inefficient. To solve these problems, map
sparsification becomes a practical necessity to acquire a subset of the
original map for localization. Previous map sparsification methods add a
quadratic term in mixed-integer programming to enforce a uniform distribution
of selected landmarks, which requires high memory capacity and heavy
computation. In this paper, we formulate map sparsification in an efficient
linear form and select uniformly distributed landmarks based on 2D discretized
grids. Furthermore, to reduce the influence of different spatial distributions
between the mapping and query sequences, which is not considered in previous
methods, we also introduce a space constraint term based on 3D discretized
grids. The exhaustive experiments in different datasets demonstrate the
superiority of the proposed methods in both efficiency and localization
performance. The relevant codes will be released at
https://github.com/fishmarch/SLAM_Map_Compression.

##  1 Introduction

To realize autonomous navigation for robots, localization in a pre-built map
is a basic technique. Lots of algorithms have been proposed for mapping and
localization using different sensors, including camera [20] and lidar [27].
These algorithms commonly work well in some small-scale environments now. When
applied in large-scale environments or in long-term, however, new challenges
appear and need to be settled for practical applications. Some of these
problems are time and memory consuming.

When using cameras, visual simultaneous localization and mapping (SLAM) is a
commonly used method to build maps. In visual SLAM, redundant features are
extracted from images and then constructed as landmarks in a map, such that
camera poses can be tracked robustly and accurately. These redundant landmarks
promise good localization results in small-scale environments. As more and
more images are received when working in large-scale environments or in long-
term, however, memory consumption is increasing unboundedly. Localizing in
such large maps will also be more time-consuming. These problems are
especially severe for some low-cost robots.

Figure 1: Localization results in a compact map. The original map is
constructed from sequence 00 of the KITTI dataset [14]. The original map
consists of $141$K landmarks, indicated by gray points; while the compact map
consists of only $5$K landmarks, indicated by blue points. $96.15\%$ of the
query images are localized successfully in the compact map, indicted by green
poses.

Actually, not all landmarks are necessary for robots to localize in pre-built
maps. In theory, even only 4 matched landmarks can determine a camera pose
using EPnP [16] (more landmarks are commonly used for robust and accurate
estimation). This reveals that maps can be compressed and still retain
comparable performance for localization. Map compression can be classified
into two types: descriptor compression [18, 22, 5] and landmark sparsification
[17, 23]. The research of this paper falls in the latter one, which is to find
a subset of an original map while maintaining comparable localization
performance. A subset map is called compact map in this paper. For example, as
indicated in Fig. 1, only $3.91\%$ of the original landmarks are selected as
the compact map, in which more than $96\%$ of the query images can still be
localized successfully.

To select an optimal subset for localization, map sparsification is related to
a $K$-cover problem [17], which means the number of landmarks in a compact map
is minimized while keeping the number of associated landmarks in each image
larger than a threshold (i.e., $K$). To solve the $K$-cover problem, it can be
formulated as mixed-integer linear programming, through which an optimal
subset is obtained. The original formulation only considers the number of
landmarks for localization, while their distribution also affects localization
performance. Therefore, some works design and add quadratic terms, formulating
mixed-integer quadratic programming to enforce a more uniform distribution of
selected landmarks [23, 11]. However, these quadratic terms slow down the
optimization speed heavily. The required high memory capacity and heavy
computation become severe limitations of these map sparsification methods. For
example, in our experiments, the mixed-integer quadratic programming methods
cannot be used for the maps containing more than 55K landmarks, because all
computer memory has been consumed.

To select uniformly distributed landmarks and at the same time maintain the
computation efficiency, we keep map sparsification formulation in a linear
form in this paper. We firstly discretize images into 2D fix-sized grids. Then
for all observed landmarks, we can find which cells they fall in. Therefore,
more occupied cells reflect a more uniform distribution of landmarks. This can
be formulated in a linear form easily. In this way, uniformly distributed
landmarks are selected efficiently, and thus the localization performance in
compact maps will be better.

Another severe limitation is that all of past works assume the spatial
distribution of the query sequence is close to that of the mapping sequence.
Then landmarks are selected only based on their association with the images of
the mapping sequence. However, this assumption cannot be guaranteed in real
robotic applications. The perspective difference between query and mapping
sequences may cause the localization in compact maps to fail unexpectedly.

To ensure more query images from whole 3D space can be localized successfully
in compact maps, we propose to select landmarks based on not only their
association with mapping images but also their visibility in 3D space. The
visibility region of a landmark is defined based on its viewing angle and
distance. The 3D space is also discretized into a 3D discretized grid. For
each 3D cell, all visible landmarks are collected, and a constraint on the
number of visible landmarks is added into map sparsification formulation. In
this way, landmarks are selected to maintain localization performance for
query images from the whole space.

In summary, the contributions of this paper are as follows: 1) We propose an
efficient map sparsification method formulating uniform landmark distribution
in a linear form to keep the computation efficiency. 2) We propose to perform
map sparsification involving the visibility of landmarks to achieve better
localization results for the query images from the whole space. 3) We conduct
exhaustive experiments in different datasets and compare with other state-of-
the-art methods, showing the effectiveness and superiority of our methods for
different kinds of query sequences.

##  2 Related Work

Map compression or sparsification has received many research attention in SLAM
and structure from motion (SfM) because of the burden of maintaining and
utilizing a large map. Different methods are also proposed to solve this
problem. In [18, 22, 5], the descriptors are quantized to compact forms to
reduce memory consumption. Mera-Trujillo et al. [19] propose a constrained
quadratic program method and design an efficient solver to select landmarks.
Contreras et al. [8] develop a map compression method based on traveled
trajectories, and they also propose an online system performing SLAM and map
compression simultaneously [9]. Some SLAM systems [26, 28, 7] are also
designed to select and build a compact map online, but much more landmarks are
still remained compared with the offline methods. Recently, some learning-
based map sparsification methods have also appeared and achieved good results
[25, 4], but they require high-cost GPUs to train and struggle to be applied
in different unseen scenes.

Map sparsification is more commonly regarded as a $K$-cover problem. Since
$K$-cover is an NP-hard problem, Li et al. [17] propose to use a greedy
algorithm to select the landmark observed by the largest number of not-yet-
fully covered frames. Cao et al. [3] also use a greedy algorithm to select
landmarks, considering both the coverage and distinctiveness. In [6], the
original $K$-cover method is improved by predicting the parameter $K$.
Similarly, the landmarks are scored and sampled based on the observation
count, covariance, and descriptor stability in [12]. Moreover, the authors
come up with more scoring factors to rank landmarks for selecting in [13].
These methods can provide compact maps for localization, but depend on a
careful design of selecting methods and cannot get the optimal solution.

Another $K$-cover based map sparsification solution is to be formulated as a
mixed-integer programming problem, which is introduced in detail in [23]. The
linear formulation does not consider the distribution of the selected
landmarks but can be solved efficiently, and therefore it is adopted in some
later works [13, 24]. The mixed-integer quadratic programming formulation
encourages the uniformly distributed landmarks but requires a high capacity of
memory and heavy computation [23] and becomes impractical for large maps. To
improve the computation efficiency, Dymczyk et al. [11] have to divide the
whole map into several small parts and perform map sparsification separately.
Camposeco et al. [2] keep the linear form of map sparsification but divide
each image into $q$ uniformly-sized cells and let each cell be covered by
$K/q$ landmarks. This method seems reasonable, but we find in the experiments
that the number of corresponding selected landmarks in some images is much
less than $K$ actually, and thus the localization results get worse.
Similarly, the proposed method also makes use of discretized images, but the
image is still considered as a whole and thus avoids the above problems.

All of the above methods are based on the existing association between
landmarks and images, thus landmarks are selected only to ensure localization
performance of these images. For an arbitrary query image with a different
viewing perspective, localization may fail unexpectedly, especially when the
compression ratio is high. To solve this problem, the proposed method
considers landmark visibility in the 3D space and adds constraints to preserve
enough visible landmarks for the whole space in compact maps.

##  3 Map Sparsification Formulation

The proposed method is designed mainly for the maps built from feature-based
visual SLAM systems (e.g., ORB-SLAM [20]). Receiving sequences of images, the
feature points observed by keyframes are constructed as landmarks. Thus, an
original map $\mathcal{M}$ consists of landmarks
$\left\\{\boldsymbol{p}_{i}\right\\}_{i=1,2,\dots,N}$, keyframes
$\left\\{\boldsymbol{F}_{j}\right\\}_{j=1,2,\dots,M}$, and their association.
The map sparsification aims to decrease the number of landmarks. We firstly
introduce commonly used formulations of map sparsification in past works [23].

Using mixed-integer programming, map sparsification can be formulated as:

$\begin{array}[]{cl}\underset{\mathbf{x},\boldsymbol{\xi}}{\operatorname{minimize}}&\mathbf{q}^{\top}\mathbf{x}+\lambda\mathbf{1}^{\top}\boldsymbol{\xi}\\\
&\mathbf{Ax}\geq K\mathbf{1}-\boldsymbol{\xi}\\\ \operatorname{subject\
to}&\mathbf{x}\in\left\\{0,1\right\\}^{N}\\\
&\boldsymbol{\xi}\in\left\\{0,\mathbb{Z}_{+}\right\\}^{M}\end{array}$ (1)

where $\mathbf{x}$ is a binary vector whose $i^{\operatorname{th}}$ element
indicates whether the landmark $\boldsymbol{p}_{i}$ is selected or not.
$\mathbf{q}$ is a weight vector to select landmarks. $\mathbf{A}$ is a binary
matrix describing the association between landmarks and keyframes, i.e., the
element $\mathbf{A}_{ij}$ indicates whether the landmark $\boldsymbol{p}_{j}$
is visible in the keyframe $\boldsymbol{F}_{i}$. $K$ is the desired minimum
number of visible landmarks in each keyframe and can be used to adjust
compression ratio. To handle the situation when the total number of associated
landmarks of a keyframe is less than $K$, a slack variable $\boldsymbol{\xi}$
is used to relax the hard constraints, and $\lambda$ determines the hardness
of the constraints.

Therefore, Eq. 1 describes the mixed-integer linear programming method for map
sparsification. The distribution of selected landmarks is not considered in
this method. To involve the distribution of landmarks, quadratic terms are
commonly used:

$\begin{array}[]{cl}\underset{\mathbf{x},\boldsymbol{\xi}}{\operatorname{minimize}}&\frac{1}{2}\mathbf{x}^{\top}\mathbf{Q}\mathbf{x}+\mathbf{q}^{\top}\mathbf{x}+\lambda\mathbf{1}^{\top}\boldsymbol{\xi}\\\
&\mathbf{Ax}\geq K\mathbf{1}-\boldsymbol{\xi}\\\ \operatorname{subject\
to}&\mathbf{x}\in\left\\{0,1\right\\}^{N}\\\
&\boldsymbol{\xi}\in\left\\{0,\mathbb{Z}_{+}\right\\}^{M}\end{array}$ (2)

where $\mathbf{Q}$ is a symmetric matrix. The element $\mathbf{Q}_{ij}$
indicates the relation between $\boldsymbol{p}_{i}$ and $\boldsymbol{p}_{j}$.
For example, it can be the number of common observations for each pair of
landmarks. Therefore, the landmarks that are observed together with others
frequently would be discarded. As a result, selected landmarks are distributed
more uniformly. Eq. 2 describes the commonly used map sparsification
formulation in past works, which considers both the number and distribution of
landmarks for localization.

##  4 Map Sparsification Based on 2D Grids

Uniformly distributed landmarks commonly provide better localization results.
Because of the quadratic terms, however, the computation complexity increases
largely. Besides, Eq. 2 encourages selection of rarely observed landmarks,
which may have larger errors and be detrimental for localization [11].

To solve these problems, we design a linear term to encourage a uniform
distribution of selected landmarks. As illustrated in Fig. 2, the image is
discretized into fixed-size grids (e.g., $C\times{R}$ cells). Each cell has
two stats: _empty_ or _occupied_. Based on the camera model, projected
positions of landmarks can be computed easily. Then a cell is _occupied_ if at
least one landmark is projected into it; otherwise, it is _empty_. Therefore,
to acquire uniformly distributed landmarks, we encourage selecting the
landmarks that make more cells to be _occupied_. This concept can be easily
formulated into Eq. 1, which becomes:

$\begin{array}[]{cl}\underset{\mathbf{x},\boldsymbol{\xi},\boldsymbol{\phi}}{\operatorname{minimize}}&\mathbf{q}^{\top}\mathbf{x}+\lambda_{1}\mathbf{1}^{\top}\boldsymbol{\xi}+\lambda_{2}\mathbf{1}^{\top}\boldsymbol{\phi}\\\
&\mathbf{Ax}\geq K\mathbf{1}-\boldsymbol{\xi}\\\
&\mathbf{Bx}\geq\mathbf{1}-\boldsymbol{\phi}\\\ \operatorname{subject\
to}&\mathbf{x}\in\left\\{0,1\right\\}^{N}\\\
&\boldsymbol{\xi}\in\left\\{0,\mathbb{Z}_{+}\right\\}^{M}\\\
&\boldsymbol{\phi}\in\left\\{0,1\right\\}^{M\times C\times R}\\\ \end{array}$
(3)

where $\mathbf{B}$ is a binary matrix describing the association between the
landmarks and all cells, i.e., the element $\mathbf{B}_{ij}$ indicates whether
the projected position of the landmark $\boldsymbol{p}_{j}$ falls within the
cell $\boldsymbol{c}_{i}$. We encourage every cell to keep at least one
corresponding landmark and use the slack variable $\boldsymbol{\phi}$ to relax
this constraint. An example is illustrated in Fig. 2, the dotted points can be
discarded without changing the number of _occupied_ cells.

Figure 2: Illustration of discretized images and associated landmarks. The
_occupied_ cells are drawn in deeper color. The red points represent the
projected position of associated landmarks. Dotted points can be discarded
without changing occupancy.

Eq. 3 formulates map sparsification as a mixed-integer linear programming
problem and thus can be solved efficiently. The selected landmarks are
directly encouraged to spread uniformly in all corresponding keyframes. Thus
better localization performance will be achieved.

##  5 Map Sparsification In 3D Space

From the previous section, it is noticed that map sparsification is formulated
based on the association between landmarks and keyframes. This association can
be easily constructed from feature matching during mapping. All past works are
based on such association and assume the spatial distribution of the query
sequence will be close to that of the mapping sequence. Thus selected
landmarks can be matched again in query images. Obviously, this assumption
cannot be promised in real robotic applications. Some selected landmarks would
not be matched in query images because of different viewing perspectives, and
thus, the localization may fail unexpectedly.

In this section, we will introduce our proposed map sparsification method that
involves landmark visibility in 3D space and thus improves localization
performance for query images from the whole space.

###  5.1 Landmark Visibility

Feature points (e.g., ORB [21]) are commonly extracted and constructed as
landmarks in visual SLAM. One of the basic requirements of used features is
robust matching from different viewing perspectives. In other words, it is
assumed that a valid feature point is visible and can be matched within a
certain region. Specifically, we assume a feature point can be matched when
the viewing angle and distance are within certain thresholds. Such specific
regions can be approximated from mapping process. Similarly defined visibility
is also used for efficient feature matching [20]. Deng et al. [10] also
predict such feature matching to ensure localization performance in an active
SLAM framework.

For a landmark $\boldsymbol{p}$ observed by multiple keyframes, an average
viewing direction ${\bar{\boldsymbol{v}}^{p}}$ can be computed from all
viewing directions, pointing from $\boldsymbol{p}$ to the optical centers of
these keyframes. Therefore, the first condition for $\boldsymbol{p}$ to be
visible is that the viewing angle is within the threshold,

$\langle\boldsymbol{\bar{v}}^{p},\boldsymbol{v}^{i}\rangle<\theta_{th}^{p}$
(4)

where $\boldsymbol{v}^{i}$ is the viewing direction for a query position. In
our experiments, $\theta_{th}^{p}$ is set according to the largest viewing
angle during mapping.

Besides, an image pyramid is commonly employed for matching features from
different distance. Therefore, the second condition for $\boldsymbol{p}$ to be
visible is that the viewing distance $d_{i}^{p}$ from a query position is
within the following range,

$d_{\min}^{p}<d_{i}^{p}<d_{\max}^{p}$ (5)

where $d_{\min}^{p}$ and $d_{\max}^{p}$ are also set according to the smallest
and largest viewing distance during mapping.

Figure 3: 3D cells in which the landmark $\boldsymbol{p}$ is visible.

Therefore, a landmark is assumed to be visible when meeting the conditions Eq.
4 and Eq. 5. The visible region forms a truncated spherical cone, as indicated
in Fig. 3.

###  5.2 Map Sparsification Based on 3D Grid

For map sparsification methods described in Sec. 3, enough visible landmarks
are preserved for every keyframe based on the constraint $\mathbf{Ax}\geq
K\mathbf{1}-\boldsymbol{\xi}$. Therefore, the compact map can be used to
localize the query images that are close to the mapping sequence as expected.
However, for query images taken from different viewing perspectives, which is
common in practical robotic applications, the number of matched visible
landmarks may decrease, and thus the localization may fail unexpectedly.

To solve the above problem, we propose a map sparsification method considering
landmark visibility in 3D space. We firstly discretize the whole space into a
3D grid. The corresponding 3D cells from which it is visible are collected for
each landmark, as indicated in Fig. 3. Therefore, each cell stores all visible
landmarks from that position. To decrease computation, we define a threshold
$K_{2}$, and the 3D cell is _valid_ only if the number of corresponding
visible landmarks is larger than $K_{2}$.

In the proposed map sparsification method, landmarks are selected to keep as
many _valid_ 3D cells as possible. Therefore, an extra constraint is
introduced into the formulation for each valid 3D cell:

$\begin{array}[]{cl}\underset{\mathbf{x},\boldsymbol{\xi},\boldsymbol{\phi},\boldsymbol{\varphi}}{\operatorname{minimize}}&\mathbf{q}^{\top}\mathbf{x}+\lambda_{1}\mathbf{1}^{\top}\boldsymbol{\xi}+\lambda_{2}\mathbf{1}^{\top}\boldsymbol{\phi}+\lambda_{3}\mathbf{1}^{\top}\boldsymbol{\varphi}\\\
&\mathbf{Ax}\geq K_{1}\mathbf{1}-\boldsymbol{\xi}\\\
&\mathbf{Bx}\geq\mathbf{1}-\boldsymbol{\phi}\\\ &\mathbf{Cx}\geq
K_{2}\mathbf{1}-\boldsymbol{\varphi}\\\ \operatorname{subject\
to}&\mathbf{x}\in\left\\{0,1\right\\}^{N}\\\
&\boldsymbol{\xi}\in\left\\{0,\mathbb{Z}_{+}\right\\}^{M}\\\
&\boldsymbol{\phi}\in\left\\{0,1\right\\}^{M\times C\times R}\\\
&\boldsymbol{\varphi}\in\left\\{0,\mathbb{Z}_{+}\right\\}^{S}\end{array}$ (6)

where $\mathbf{C}$ is a binary matrix describing the association between
landmarks and valid 3D cells, i.e., the element $\mathbf{C}_{ij}$ indicates
whether the landmark $\boldsymbol{p}_{j}$ is visible for the 3D cell
$\boldsymbol{g}_{i}$. As before, the slack variable $\boldsymbol{\varphi}$ is
used to relax the constraints.

Therefore, the proposed map sparsification method aims to preserve enough
visible landmarks not only for mapping keyframes but also for the whole 3D
space. As a result, more query images from the whole space can be localized
successfully.

##  6 Experiments

The proposed map sparsification based on 2D grids (i.e., Eq. 3) is firstly
compared with previous works to show its superiority in performance and
efficiency for localizing query images close to mapping sequences. Then we
demonstrate the effectiveness of map sparsification involving the visibility
of landmarks in 3D space (i.e., Eq. 6) and its superiority when query images
are not close to mapping sequences. More experimental results can be found in
the supplementary material, including visualization results and the comparison
between original and compact maps in terms of memory consumption and
localization efficiency.

###  6.1 Experimental Setup

All experiments are tested on a laptop with an i7-7700HQ 2.80 Hz CPU and 16GB
RAM. The mixed-integer programming problem is solved using
Gurobi111https://www.gurobi.com.

The proposed methods are tested on three commonly used SLAM datasets: ICL-NUIM
[15], EuRoC [1], and KITTI [14]. These datasets cover indoor and outdoor,
small-scale and large-scale environments; the number of landmarks in original
maps changes from thousands to more than one hundred thousand. The original
map is constructed by running a state-of-the-art visual SLAM system (i.e.,
ORB-SLAM2 [20]). Query images are localized in the map based on the re-
localization module in ORB-SLAM2, in which localization is successful when
enough inliers are matched. Different map sparsification methods are compared
in terms of localization rate and run-time. Localization rate is the ratio of
the number of successfully localized images to the number of all query images
in one sequence.

The weight vector $\mathbf{q}$ in all methods is set based on the number of
matched times for each landmark, just like in [23]. The proposed method based
only on 2D grids (i.e., Eq. 3) is denoted as Ours-2D, while the method based
on 2D and 3D grids (i.e., Eq. 6) is denoted as Ours-3D. The proposed methods
are compared with several state-of-the-art $K$-cover-based map sparsification
methods. The map sparsification using mixed-integer linear programming (i.e.,
Eq. 1) is denoted as LP. For the methods using mixed-integer quadratic
programming (i.e., Eq. 2), the weight matrix $\mathbf{Q}$ is set according to
the number of co-observation times[23] or average projected distance[11], and
these two methods are denoted as QP1 and QP2 respectively. The map
sparsification based on divided images [2] is denoted as DI.

###  6.2 Map Sparsification Based on 2D Grids

Table 1: Comparison of different map sparsification methods in terms of localization rate. Num denotes the number of landmarks in the mam, Rate denotes the localization rate in the corresponding map. The highest localization rates in compact maps are labelled in bold. ‘-’ means the method cannot work because of memory limit.  Dataset | Original Map | Ours-2D | LP [23] | QP1 [23] | QP2 [11] | DI [2] [t]
---|---|---|---|---|---|---
Num |  Rate |  Num |  Rate |  Num |  Rate |  Num |  Rate |  Num |  Rate |  Num |  Rate [b]
ICL-NUIM [15] |  liv0 | 4061 |  98.34% |  519 |  79.97% |  519 |  79.51% |  519 |  79.64% |  519 |  78.65% |  734 |  77.39% [t]
liv1 | 4112 |  99.89% |  445 |  87.47% |  445 |  86.23% |  445 |  87.37% |  445 |  85.51% |  516 |  55.59%
liv2 | 6435 |  100.0% |  312 |  91.36% |  310 |  89.77% |  310 |  90.57% |  310 |  91.02% |  833 |  89.89%
liv3 | 6173 |  93.38% |  532 |  76.21% |  531 |  73.71% |  531 |  75.32% |  531 |  73.55% |  987 |  61.85%
off0 | 5875 |  99.93% |  398 |  95.62% |  398 |  95.42% |  398 |  95.76% |  398 |  95.76% |  791 |  84.88%
off1 | 5992 |  93.57% |  485 |  80.00% |  485 |  77.93% |  485 |  79.07% |  484 |  76.89% |  736 |  37.20%
off2 | 6838 |  100.0% |  299 |  92.27% |  298 |  91.02% |  298 |  90.34% |  298 |  91.36% |  759 |  90.45%
off3 | 5367 |  99.52% |  396 |  83.47% |  395 |  83.06% |  395 |  78.31% |  396 |  83.15% |  430 |  51.05% [b]
EuRoC [1] |  MH1 | 12506 |  99.04% |  2978 |  97.33% |  2976 |  97.03% |  2976 |  96.76% |  2977 |  96.54% |  4210 |  94.56% [t]
MH2 | 11306 |  99.73% |  2546 |  98.17% |  2544 |  96.43% |  2547 |  98.33% |  2547 |  98.30% |  3665 |  93.30%
MH3 | 15234 |  97.34% |  2105 |  95.32% |  2104 |  95.32% |  2107 |  95.44% |  2105 |  95.17% |  3720 |  91.30%
MH4 | 17490 |  99.69% |  744 |  97.37% |  743 |  97.32% |  742 |  97.27% |  742 |  97.01% |  2493 |  96.05%
MH5 | 16956 |  99.77% |  1116 |  98.87% |  1116 |  98.69% |  1116 |  98.60% |  1116 |  98.60% |  2933 |  96.08%
V101 | 6864 |  99.03% |  461 |  86.04% |  460 |  84.92% |  461 |  85.52% |  461 |  86.35% |  1011 |  79.14%
V102 | 9878 |  98.12% |  507 |  85.93% |  507 |  84.86% |  506 |  85.87% |  506 |  85.49% |  1255 |  84.17%
V103 | 11518 |  93.22% |  1140 |  83.76% |  1140 |  83.91% |  1140 |  83.09% |  1140 |  83.05% |  1661 |  77.89%
V201 | 7346 |  94.69% |  565 |  79.64% |  564 |  78.75% |  564 |  78.62% |  564 |  78.97% |  990 |  78.75%
V202 | 11791 |  99.48% |  1254 |  93.25% |  1254 |  92.86% |  1254 |  93.07% |  1254 |  93.42% |  2285 |  90.39%
V203 | 14557 |  87.51% |  1796 |  81.48% |  1795 |  81.27% |  1796 |  81.11% |  1795 |  81.48% |  2384 |  75.45% [b]
KITTI [14] |  00 | 141673 |  100.0% |  5536 |  96.15% |  5524 |  95.68% |  - |  - |  - |  - |  18542 |  91.15% [t]
01 | 89443 |  91.01% |  4218 |  76.02% |  4216 |  76.02% |  - |  - |  - |  - |  10204 |  71.84%
02 | 182499 |  99.96% |  7660 |  97.47% |  7661 |  97.15% |  - |  - |  - |  - |  25133 |  86.76%
03 | 24284 |  100.0% |  837 |  96.00% |  836 |  95.63% |  837 |  96.13% |  836 |  95.76% |  2952 |  88.76%
04 | 16824 |  100.0% |  583 |  99.63% |  585 |  99.26% |  572 |  99.63% |  572 |  99.63% |  2019 |  97.41%
05 | 74158 |  100.0% |  2737 |  97.46% |  2738 |  97.25% |  - |  - |  - |  - |  9655 |  90.62%
06 | 42522 |  100.0% |  1879 |  98.27% |  1878 |  98.09% |  1878 |  98.18% |  1877 |  98.18% |  5635 |  95.82%
07 | 27036 |  100.0% |  1242 |  92.28% |  1240 |  91.28% |  1240 |  91.73% |  1240 |  92.55% |  3517 |  92.01%
08 | 119323 |  100.0% |  5903 |  94.87% |  5899 |  94.82% |  - |  - |  - |  - |  16872 |  83.03%
09 | 55984 |  99.94% |  2945 |  95.79% |  2944 |  95.60% |  2942 |  95.85% |  2943 |  95.79% |  8224 |  88.37%
10 | 30933 |  99.92% |  1729 |  97.09% |  1728 |  96.67% |  1728 |  96.83% |  1729 |  96.42% |  4743 |  86.34% [b]

Since previous methods assume the spatial distribution of the query sequence
is close to that of the mapping sequence, one sequence is firstly used to
build an original map and then the images from the same mapping sequence are
re-localized in compact maps. In this subsection, the proposed map
sparsification method is only based on 2D discretized grids (i.e., Eq. 3).
Tab. 1 compares the localization rates in different maps obtained by different
methods. An example of the selected sparse landmarks is illustrated by blue
points in Fig. 1.

The experimental results demonstrate that the number of landmarks can be
reduced largely while the localization is still successful for most query
images. For example, only $4.37\%$ landmarks are selected from the map of
sequence off2 and $3.46\%$ for sequence 04, but more than $90\%$ of the query
images can still be localized successfully. The compression ratio is mainly
controlled by the threshold $K$ in Eq. 1, Eq. 2, and Eq. 3, and we use $K=50$
in our experiments. With a larger $K$, more landmarks will be selected in
compact maps, and localization rate will be higher.

For more sequences, our proposed method achieves the highest localization rate
(labelled in bold), especially compared with the method only considering the
number of landmarks for map sparsification (i.e., LP). Our method also gets
slightly better results than other methods that consider landmark distribution
(i.e., QP1, QP2, and DI). In QP1, the weight matrix $\mathbf{Q}$ is set
according to the number of co-observation times between landmarks and as a
result, rarely observed landmarks are preferred, which may be detrimental for
localization [11]. To overcome this limitation, QP2 sets the weight matrix
$\mathbf{Q}$ based on the average distance between projected positions. But
instead of uniform distribution, QP2 actually prefers to select landmarks that
are projected far from the image center. Unlike these quadratic programming-
based methods, our method formulates uniform distribution in a more clear and
simple way, achieving better results.

DI gets unexpected bad performance in these testings. DI divides images into
$q$ cells and forces each cell is covered by at least $K/q$ landmarks.
Therefore, DI also achieves map sparsification in a linear programming form
and considers both the number and distribution of landmarks. But we find in
these testings that the slack variable $\lambda$ in Eq. 1 needs to be small
for DI to get a high compression ratio. As a result, the number of selected
landmarks is much less than $K$ for some images, where localization would
fail. Although based on similar discretized grids, our method still considers
an image as a whole and thus avoids this problem.

Table 2: Comparison of different map sparsification methods in terms of run-time and consumed memory. $N$ denotes the number of landmarks, $N_{k}$ denotes the number of keyframes. Seq. | Original Map | Run-time (s) [t]
---|---|---
$N$ | $N_{k}$ |  Ours |  LP |  QP1 |  QP2 |  DI [b]
liv0 | 4.1K | 107 |  0.2 |  0.1 |  13.6 |  12.5 |  0.2 [t]
off0 | 5.9K | 133 |  0.3 |  0.1 |  27.9 |  18.9 |  0.2
V103 | 11.5K | 216 |  1.1 |  0.3 |  59.8 |  36.3 |  0.6
MH5 | 17.0K | 369 |  4.5 |  2.2 |  107 |  61.5 |  1.2
07 | 27.0K | 258 |  2.6 |  1.1 |  127 |  95.5 |  0.7
09 | 56.0K | 603 |  7.7 |  9.9 |  754 |  491 |  1.8
08 | 119K | 1269 |  33.2 |  19.2 |  - |  - |  4.9
02 | 182K | 1794 |  69.6 |  63.2 |  - |  - |  7.9
| $N$ | $N_{k}$ | Consumed memory (MB) [t]
liv0 | 4.1K | 107 |  35 |  14 |  555 |  476 |  33 [t]
off0 | 5.9K | 133 |  61 |  27 |  842 |  708 |  21
V103 | 11.5K | 216 |  128 |  56 |  3048 |  1367 |  83
MH5 | 17.0K | 369 |  255 |  107 |  3532 |  1656 |  77
07 | 27.0K | 258 |  141 |  62 |  6315 |  5766 |  100
09 | 56.0K | 603 |  676 |  230 |  13005 |  10952 |  211
08 | 119K | 1269 |  1162 |  637 |  - |  - |  512
02 | 182K | 1794 |  3494 |  1502 |  - |  - |  1187 [b]

Another superiority of our method is computational efficiency. The comparison
of run-time and memory consumption on different sequences is shown in Tab. 2.
The run-time for map sparsification depends on the number of landmarks,
keyframes and associations among them. As the number of landmarks or keyframes
increases, the run-time will commonly increase quickly. It is clear that the
methods using quadratic programming are much slower, especially for the maps
containing a large number of landmarks (e.g., sequence 09). In addition to
solving more complex optimization problems, QP1 and QP2 also need much more
time to set the weight matrix $\mathbf{Q}$, since the relation between every
pair of landmarks needs to be found and this is not recorded in original maps.
In our method, both the matrix $\mathbf{A}$ and $\mathbf{B}$ in Eq. 3 are set
according to the relation between landmarks and keyframes, which can be
retrieved from original maps directly and efficiently. In addition, since
formulated in a linear form, our method can also be solved efficiently and is
only slightly slower than LP.

Similarly for memory consumption, ours consumes much less memory than
quadratic programming. QP1 and QP2 cannot work for some sequences (e.g., 00)
on our device with 16 GB memory. The testing for sequence 09 consumes nearly
all of the memory when using QP1 and QP2, in which the original map contains
56K landmarks.

Therefore, if the mapping sequence has covered most of the query space for
localization, the proposed map sparsification based on 2D discretized grids
can be used to decrease the number of landmarks to save memory and improve
localization efficiency.

###  6.3 Map Sparsification Based on 2D and 3D Grids

Figure 4: Comparison of map sparsification results. The original map is
constructed from sequence V201. There are 3401 blue points indicating valid 3D
cells of the compact map acquired by Ours-2D; while there are 6468 green
points indicating valid 3D cells of the compact map acquired by Ours-3D.

We have demonstrated that compact maps with much fewer landmarks achieve good
localization performance when the spatial distribution of the query sequence
is close to that of the mapping sequence. In this subsection, we continue to
test the localization performance for query images from different
distributions.

In ICL-NUIM [15] and EuRoC [1] datasets, different sequences are recorded in
the same environments. Therefore, one sequence is used to build original maps
and then query images of _other_ sequences are also localized in compact maps.
The comparison of experimental results is presented in Table 3. For example,
sequence off0 is used to build the map, and then sequence off0, off2 and off3
are used as query sequences to be localized in the original and compact maps.
Not all sequences of these two datasets are used in this experiment since the
successful localization is very low even in original maps for some sequences.

Table 3: Comparison of different map sparsification methods in terms of localization rate. Num denotes the number of landmarks in the map, Rate denotes the localization rate in the corresponding map. The sequence for mapping is indicated in bold and the images from several different sequences are localized in the map. The highest localization rates in compact maps are labelled in bold. | off0 | off2 |  liv2
---|---|---|---
Map | Num | Rate | Num | Rate | Num | Rate [t]
off0 | off2 | off3 | off0 | off2 | liv1 | liv2 [b]
Original | 5875 | 99.93% | 66.70% | 52.26% | 6838 | 65.65% | 100.00% | 6435 | 55.80% | 100.00% [t]
Ours-2D | 398 | 95.62% | 27.95% | 7.74% | 299 | 4.31% | 92.27% | 312 | 2.48% | 91.36%
Ours-3D | 399 | 95.49% | 28.52% | 7.98% | 301 | 6.17% | 91.70% | 314 | 4.04% | 90.68%
| MH1 | MH4 |  MH5
Map | Num | Rate | Num | Rate | Num | Rate [t]
MH1 | MH2 | MH3 | MH4 | MH5 | MH4 | MH5 [b]
Original | 12506 | 99.04% | 87.10% | 52.53% | 17490 | 99.69% | 84.20% | 16956 | 82.34% | 99.77% [t]
Ours-2D | 2978 | 97.33% | 73.67% | 45.76% | 744 | 97.37% | 63.77% | 1116 | 74.75% | 98.87%
Ours-3D | 2982 | 96.78% | 73.97% | 46.18% | 745 | 97.06% | 64.27% | 1121 | 75.71% | 98.74%
| MH2 | V101 |  V102
Map | Num | Rate | Num | Rate | Num | Rate [t]
MH1 | MH2 | MH3 | V101 | V102 | V101 | V102 [b]
Original | 11306 | 91.78% | 99.73% | 51.81% | 6864 | 99.03% | 61.56% | 9878 | 72.25% | 98.12% [t]
Ours-2D | 2546 | 74.80% | 98.17% | 35.50% | 461 | 86.04% | 32.60% | 507 | 34.61% | 85.93%
Ours-3D | 2550 | 77.41% | 97.77% | 38.69% | 462 | 85.79% | 34.61% | 506 | 37.67% | 85.99%
| MH3 | V201 |  V202
Map | Num | Rate | Num | Rate | Num | Rate [t]
MH1 | MH2 | MH3 | V201 | V202 | V201 | V202
Original | 15234 | 59.93% | 52.80% | 97.34% | 7346 | 94.69% | 52.77% | 11791 | 75.22% | 99.48%
Ours-2D | 2105 | 46.85% | 37.23% | 95.32% | 565 | 79.64% | 19.48% | 1254 | 49.82% | 93.25%
Ours-3D | 2108 | 46.94% | 37.60% | 95.25% | 569 | 80.00% | 22.77% | 1256 | 50.04% | 93.46%

As shown in Tab. 3, compact maps commonly achieve good localization
performance for query images from mapping sequences, while localization rates
may decrease largely for other query sequences. Because the map sparsification
is mainly selecting landmarks to maintain the observations in keyframes of the
mapping sequences, and thus discarding other redundant landmarks which may be
useful for localizing query images from the whole space. Our proposed map
sparsification method in Eq. 6 aims to add such extra constraints for the
whole space. In most testings, ours-3D achieves higher localization rates for
the query sequences that are different from mapping sequences. For the query
images from the mapping sequence, ours-3D also achieves comparable
localization performance with other map sparsification methods. Other map
sparsification methods achieve similar results as Ours-2D, which are included
in the supplementary material.

A comparison of compact maps acquired by ours-2D and ours-3D is illustrated in
Fig. 4. The blue points indicate the valid 3D cells of ours-2D, while the
green points indicate the valid 3D cells of ours-3D. The numbers of selected
landmarks are similar using these two methods. But it is clear that much more
valid 3D cells are kept using ours-3D and they are spreading over a larger
space. As a result, more query images from the whole 3D space can be localized
successfully.

Table 4: Comparison of average run-time (s) in different scenes. Scene | Resolution (m) | Ours-2D | Ours-3D [t]
---|---|---|---
liv | 0.15 | 0.25 | 1.00 [t]
off | 0.15 | 0.20 | 1.19
V1 | 0.20 | 0.99 | 8.06
V2 | 0.20 | 1.24 | 10.36
MH | 0.40 | 3.62 | 11.34 [b]

Since Ours-3D is also formulated in a linear form, it can also be solved
efficiently. Tab. 4 shows the average run-time in different scenes. The run-
time of Ours-3D is related to the resolution of the 3D grids. In the
experiment, we choose different resolutions for different scenes and thousands
of valid cells are used for map sparsification. The run-time of Ours-3D does
not increase much compared with Ours-2D, and Ours-3D is still much faster than
QP1 and QP2.

Therefore, if the mapping sequence has not covered the whole query space and
query images may be taken with different viewing perspectives, the map
sparsification involving the constraints from the 3D discretized grid is a
better choice to select landmarks.

##  7 Conclusion

In this paper, we propose two novel terms for efficient map sparsification.
The first term is to enforce a uniform distribution of selected landmarks
based on 2D discretized grids. The second one adds a space constraint based on
landmark visibility to weaken the assumption that the spatial distribution of
the query sequence is close to that of the mapping sequence. Both two terms
are formulated in efficient linear forms, and thus avoid heavy computation.
The proposed terms can be chosen and set according to different conditions of
query sequences. The exhaustive experiments have demonstrated the
effectiveness and superiority of the proposed method, especially the
computation efficiency compared with previous QP methods.

Acknowledgments: This work is supported by the CUHK T Sone Robotics Institute,
and the InnoHK of the Government of Hong Kong via the Hong Kong Centre for
Logistics Robotics.

## References

  * [1] Michael Burri, Janosch Nikolic, Pascal Gohl, Thomas Schneider, Joern Rehder, Sammy Omari, Markus W Achtelik, and Roland Siegwart.  The EuRoC micro aerial vehicle datasets.  The International Journal of Robotics Research, 35(10):1157–1163, 2016.
  * [2] Federico Camposeco, Andrea Cohen, Marc Pollefeys, and Torsten Sattler.  Hybrid scene compression for visual localization.  In CVPR, pages 7645–7654, 2019.
  * [3] Song Cao and Noah Snavely.  Minimal scene descriptions from structure from motion models.  In CVPR, pages 461–468, 2014.
  * [4] M.-F. Chang, Y Zhao, R Shah, J J Engel, M Kaess, and S Lucey.  Long-term visual map sparsification with heterogeneous GNN.  In CVPR, pages 2406–2415, 2022.
  * [5] Wentao Cheng, Weisi Lin, Kan Chen, and Xinfeng Zhang.  Cascaded parallel filtering for memory-efficient image-based localization.  In ICCV, pages 1032–1041, 2019.
  * [6] Wentao Cheng, Weisi Lin, Xinfeng Zhang, Michael Goesele, and Ming-Ting Sun.  A data-driven point cloud simplification framework for city-scale image-based localization.  TIP, 26(1):262–275, 2017.
  * [7] Siddharth Choudhary, Vadim Indelman, Henrik I Christensen, and Frank Dellaert.  Information-based reduced landmark SLAM.  In ICRA, pages 4620–4627, 2015.
  * [8] Luis Contreras and Walterio Mayol-Cuevas.  Trajectory-driven point cloud compression techniques for visual SLAM.  In IROS, pages 133–140, 2015.
  * [9] Luis Contreras and Walterio Mayol-Cuevas.  O-POCO: Online point cloud compression mapping for visual odometry and SLAM.  In ICRA, pages 4509–4514, 2017.
  * [10] Xinke Deng, Zixu Zhang, Avishai Sintov, Jing Huang, and Timothy Bretl.  Feature-constrained active visual SLAM for mobile robot navigation.  In ICRA, pages 7233–7238, 2018.
  * [11] Marcin Dymczyk, Simon Lynen, Michael Bosse, and Roland Siegwart.  Keep it brief: Scalable creation of compressed localization maps.  In IROS, pages 2536–2542, 2015.
  * [12] Marcin Dymczyk, Simon Lynen, Titus Cieslewski, Michael Bosse, Roland Siegwart, and Paul Furgale.  The gist of maps-summarizing experience for lifelong localization.  In ICRA, pages 2767–2773, 2015.
  * [13] Marcin Dymczyk, Thomas Schneider, Igor Gilitschenski, Roland Siegwart, and Elena Stumm.  Erasing bad memories: Agent-side summarization for long-term mapping.  In IROS, pages 4572–4579, 2016.
  * [14] Andreas Geiger, Philip Lenz, and Raquel Urtasun.  Are we ready for autonomous driving? The KITTI vision benchmark suite.  In CVPR, pages 3354–3361, 2012.
  * [15] Ankur Handa, Thomas Whelan, John McDonald, and Andrew J Davison.  A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM.  In ICRA, pages 1524–1531, 2014.
  * [16] Vincent Lepetit, Francesc Moreno-Noguer, and Pascal Fua.  Epnp: An accurate o (n) solution to the pnp problem.  IJCV, 81(2):155–166, 2009.
  * [17] Yunpeng Li, Noah Snavely, and Daniel P Huttenlocher.  Location recognition using prioritized feature matching.  In ECCV, pages 791–804. Springer, 2010.
  * [18] Simon Lynen, Torsten Sattler, Michael Bosse, Joel A Hesch, Marc Pollefeys, and Roland Siegwart.  Get out of my lab: Large-scale, real-time visual-inertial localization.  In RSS, 2015.
  * [19] Marcela Mera-Trujillo, Benjamin Smith, and Victor Fragoso.  Efficient scene compression for visual-based localization.  In 3DV, pages 1–10, 2020.
  * [20] Raul Mur-Artal and Juan D Tardos.  ORB-SLAM2: An open-source SLAM system for monocular, stereo, and RGB-D cameras.  IEEE Transactions on Robotics, 33(5):1255–1262, 2017.
  * [21] Ethan Rublee, Vincent Rabaud, Kurt Konolige, and Gary Bradski.  ORB: An efficient alternative to SIFT or SURF.  In ICCV, pages 2564–2571, 2011.
  * [22] Torsten Sattler, Michal Havlena, Filip Radenovic, Konrad Schindler, and Marc Pollefeys.  Hyperpoints and fine vocabularies for large-scale location recognition.  In ICCV, pages 2102–2110, 2015.
  * [23] Hyun Soo Park, Yu Wang, Eriko Nurvitadhi, James C Hoe, Yaser Sheikh, and Mei Chen.  3D point cloud reduction using mixed-integer quadratic programming.  In CVPR Workshops, 2013.
  * [24] Dominik Van Opdenbosch, Tamay Aykut, Nicolas Alt, and Eckehard Steinbach.  Efficient map compression for collaborative visual SLAM.  In IEEE Winter Conference on Applications of Computer Vision (WACV), pages 992–1000, 2018.
  * [25] Luwei Yang, Rakesh Shrestha, Wenbo Li, Shuaicheng Liu, Guofeng Zhang, Zhaopeng Cui, and Ping Tan.  SceneSqueezer: Learning to compress scene for camera relocalization.  In CVPR, pages 8259–8268, 2022.
  * [26] Guangcong Zhang and Patricio A Vela.  Good features to track for visual SLAM.  In CVPR, jun 2015.
  * [27] Ji Zhang and Sanjiv Singh.  LOAM: Lidar odometry and mapping in real-time.  In RSS, 2014.
  * [28] Yipu Zhao and Patricio A Vela.  Good feature matching: Toward accurate, robust VO/VSLAM with low latency.  IEEE Transactions on Robotics, 36(3):657–675, 2020.