eqfloatequation

# Deep Global Registration

Christopher Choy
Stanford University
indicates equal contribution Wei Dong∗
Carnegie Mellon University
Vladlen Koltun
Intel Labs

###### Abstract

We present Deep Global Registration, a differentiable framework for pairwise
registration of real-world 3D scans. Deep global registration is based on
three modules: a 6-dimensional convolutional network for correspondence
confidence prediction, a differentiable Weighted Procrustes algorithm for
closed-form pose estimation, and a robust gradient-based $\text{SE}(3)$
optimizer for pose refinement. Experiments demonstrate that our approach
outperforms state-of-the-art methods, both learning-based and classical, on
real-world data.

##  1 Introduction

A variety of applications, including 3D reconstruction, tracking, pose
estimation, and object detection, invoke 3D registration as part of their
operation [35, 5, 33]. To maximize the accuracy and speed of 3D registration,
researchers have developed geometric feature descriptors [22, 11, 41, 9], pose
optimization algorithms [38, 45, 24, 50], and end-to-end feature learning and
registration pipelines [41, 2].

(a) RANSAC [35]

(b) FGR [50]

(c) DCP [41]

(d) Ours

Figure 1: Pairwise registration results on the 3DMatch dataset [47]. Our
method successfully aligns a challenging 3D pair (left), while RANSAC [35],
FGR [50], and DCP [41] fail. On an easier pair (right), our method achieves
finer alignment.

In particular, recent end-to-end registration networks have proven to be
effective in relation to classical pipelines. However, these end-to-end
approaches have some drawbacks that limit their accuracy and applicability.
For example, PointNetLK [2] uses globally pooled features to encode the entire
geometry of a point cloud, which decreases spatial acuity and registration
accuracy. Deep closest point [41] makes strong assumptions on the distribution
of points and correspondences, which do not hold for partially overlapping 3D
scans.

In this work, we propose three modules for robust and accurate registration
that resolve these drawbacks: a 6-dimensional convolutional network for
correspondence confidence estimation, a differentiable Weighted Procrustes
method for scalable registration, and a robust $\text{SE}(3)$ optimizer for
fine-tuning the final alignment.

The first component is a 6-dimensional convolutional network that analyzes the
geometry of 3D correspondences and estimates their accuracy. Our approach is
inspired by a number of learning-based methods for estimating the validity of
correspondences in 2D [46, 34] and 3D [31]. These methods stack coordinates of
correspondence pairs, forming a vector
$[\mathbf{x};\mathbf{y}]\in\mathbb{R}^{2\times D}$ for each correspondence
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{D}$. Prior methods treat these $2\times
D$-dimensional vectors as a set, and apply global set processing models for
analysis. Such models largely disregard local geometric structure. Yet the
correspondences are embedded in a metric space ($\mathbb{R}^{2\times D}$) that
induces distances and neighborhood relationships. In particular, 3D
correspondences form a geometric structure in 6-dimensional space [8] and we
use a high-dimensional convolutional network to analyze the 6D structure
formed by correspondences and estimate the likelihood that a given
correspondence is correct (_i.e_., an inlier).

The second component we develop is a differentiable Weighted Procrustes
solver. The Procrustes method [15] provides a closed-form solution for rigid
registration in $\text{SE}(3)$. A differentiable version of the Procrustes
method by Wang _et al_. [41] has been used for end-to-end registration.
However, the differentiable Procrustes method passes gradients through
coordinates, which requires $O(N^{2})$ time and memory for $N$ keypoints,
limiting the number of keypoints that can be processed by the network. We use
the inlier probabilities predicted by our first module (the 6D convolutional
network) to guide the Procrustes method, thus forming a differentiable
Weighted Procrustes method. This method passes gradients through the weights
associated with correspondences rather than correspondence coordinates. The
computational complexity of the Weighted Procrustes method is linear in the
number of correspondences, allowing the registration pipeline to use dense
correspondence sets rather than sparse keypoints. This substantially increases
registration accuracy.

Our third component is a robust optimization module that fine-tunes the
alignment produced by the Weighted Procrustes solver. This optimization module
minimizes a differentiable loss via gradient descent on the continuous
$\text{SE}(3)$ representation space [52]. The optimization is fast since it
does not require neighbor search in the inner loop [48].

Experimentally, we validate the presented modules on a real-world pairwise
registration benchmark [47] and large-scale scene reconstruction datasets [18,
6, 32]. We show that our modules are robust, accurate, and fast in comparison
to both classical global registration algorithms [50, 35, 45] and recent end-
to-end approaches [31, 41, 2]. All training and experiment scripts are
available at https://github.com/chrischoy/DeepGlobalRegistration.

##  2 Related Work

We divide the related work into three categories following the stages of
standard registration pipelines that deal with real-world 3D scans: feature-
based correspondence matching, outlier filtering, and pose optimization.

Feature-based correspondence matching. The first step in many 3D registration
pipelines is feature extraction. Local and global geometric structure in 3D is
analyzed to produce high-dimensional feature descriptors, which can then be
used to establish correspondences.

Traditional hand-crafted features commonly summarize pairwise or higher-order
relationships in histograms [20, 37, 40, 36, 35]. Recent work has shifted to
learning features via deep networks [47, 22]. A number of recent methods are
based on global pooling models [12, 11, 49], while others use convolutional
networks [14, 9].

Our work is agnostic to the feature extraction mechanism. Our modules
primarily address subsequent stages of the registration pipeline and are
compatible with a wide variety of feature descriptors.

Outlier filtering. Correspondences produced by matching features are commonly
heavily contaminated by outliers. These outliers need to be filtered out for
robust alignment. A widely used family of techniques for robust model fitting
is based on RANdom SAmple Consensus (RANSAC) [38, 1, 35, 29, 19], which
iteratively samples small sets of correspondences in the hope of sampling a
subset that is free from outliers. Other algorithms are based on branch-and-
bound [45], semi-definite programming [28, 25], and maximal clique selection
[44]. These methods are accurate, but commonly require longer iterative
sampling or more expensive computation as the signal-to-noise ratio decreases.
One exception is TEASER [44], which remains effective even with high outlier
rates. Other methods use robust loss functions to reject outliers during
optimization [50, 4].

Our work uses a convolutional network to identify inliers and outliers. The
network needs only one feed-forward pass at test time and does not require
iterative optimization.

Pose optimization. Pose optimization is the final stage that minimizes an
alignment objective on filtered correspondences. Iterative Closest Points
(ICP) [3] and Fast Global Registration (FGR) [50] use second-order
optimization to optimize poses. Makadia  _et al_. [26] propose an iterative
procedure to minimize correlation scores. Maken  _et al_. [27] propose to
accelerate this process by stochastic gradient descent.

Recent end-to-end frameworks combine feature learning and pose optimization.
Aoki _et al_. [2] combine PointNet global features with an iterative pose
optimization method [24]. Wang _et al_. [41, 42] train graph neural network
features by backpropagating through pose optimization.

We further advance this line of work. In particular, our Weighted Procrustes
method reduces the complexity of optimization from quadratic to linear and
enables the use of dense correspondences for highly accurate registration of
real-world scans.

##  3 Deep Global Registration

3D reconstruction systems typically take a sequence of partial 3D scans as
input and recover a complete 3D model of the scene. These partial scans are
scene fragments, as shown in Fig. 1. In order to reconstruct the scene,
reconstruction systems often begin by aligning pairs of fragments [6]. This
stage is known as pairwise registration. The accuracy and robustness of
pairwise registration are critical and often determine the accuracy of the
final reconstruction.

Our pairwise registration pipeline begins by extracting pointwise features.
These are matched to form a set of putative correspondences. We then use a
high-dimensional convolutional network (ConvNet) to estimate the veracity of
each correspondence. Lastly, we use a Weighted Procrustes method to align 3D
scans given correspondences with associated likelihood weights, and refine the
result by optimizing a robust objective.

The following notation will be used throughout the paper. We consider two
point clouds, $X=[\mathbf{x}_{1},...,\mathbf{x}_{N_{x}}]\in\mathbb{R}^{3\times
N_{x}}$ and $Y=[\mathbf{y}_{1},...,\mathbf{y}_{N_{y}}]\in\mathbb{R}^{3\times
N_{y}}$, with $N_{x}$ and $N_{y}$ points respectively, where
$\mathbf{x}_{i},\mathbf{y}_{j}\in\mathbb{R}^{3}$. A correspondence between
$\mathbf{x}_{i}$ and $\mathbf{y}_{j}$ is denoted as
$\mathbf{x}_{i}\leftrightarrow\mathbf{y}_{j}$ or $(i,j)$.

###  3.1 Feature Extraction

To prepare for registration, we extract pointwise features that summarize
geometric context in the form of vectors in metric feature space. Our pipeline
is compatible with many feature descriptors. We use Fully Convolutional
Geometric Features (FCGF) [9], which have recently been shown to be both
discriminative and fast. FCGF are also compact, with dimensionality as low as
16 to 32, which supports rapid neighbor search in feature space.

###  3.2 Correspondence Confidence Prediction

Given the features
$\mathcal{F}_{x}=\\{\mathbf{f}_{\mathbf{x}_{1}},...,\mathbf{f}_{\mathbf{x}_{N_{x}}}\\}$
and
$\mathcal{F}_{y}=\\{\mathbf{f}_{\mathbf{y}_{1}},...,\mathbf{f}_{\mathbf{y}_{N_{x}}}\\}$
of two 3D scans, we use the nearest neighbor in the feature space to generate
a set of putative correspondences or matches
$\mathcal{M}=\\{(i,\mathop{\mathrm{argmin}}_{j}\|\mathbf{f}_{\mathbf{x}_{i}}-\mathbf{f}_{\mathbf{y}_{j}}\|)|i\in[1,...,N_{x}]\\}$.
This procedure is deterministic and can be hand-crafted to filter out noisy
correspondences with ratio or reciprocity tests [50]. However, we propose to
learn this heuristic filtering process through a convolutional network that
learns to analyze the underlying geometric structure of the correspondence
set.

We first provide a 1-dimensional analogy to explain the geometry of
correspondences. Let $A$ be a set of 1-dimensional points $A=\\{0,1,2,3,4\\}$
and $B$ be another such set $B=\\{10,11,12,13,14\\}$. Here $B$ is a
translation of $A$: $B=\\{a_{i}+10|a_{i}\in A\\}$. If an algorithm returns a
set of possible correspondences
$\\{(0,10),(1,11),(2,12),(3,13),(4,14),(0,14),(4,10)\\}$, then the set of
correct correspondences (inliers) will form a line (first 5 pairs), whereas
incorrect correspondences (outliers) will form random noise outside the line
(last 2 pairs). If we extend this to 3D scans and pointclouds, we can also
represent a 3D correspondence $\mathbf{x}_{i}\leftrightarrow\mathbf{y}_{j}$ as
a point in 6-dimensional space
$[\mathbf{x}_{i}^{T},\mathbf{y}_{j}^{T}]^{T}\in\mathbb{R}^{6}$. The inlier
correspondences will be distributed on a lower-dimensional surface in this 6D
space, determined by the geometry of the 3D input. We denote
$\mathcal{P}=\\{(i,j)|\;\|T^{*}(\mathbf{x}_{i})-\mathbf{y}_{j}\|<\tau,(i,j)\in\mathcal{M}\\}$
as a set of inliers or a set of correspondences $(i,j)$ that align accurately
up to the threshold $\tau$ under the ground truth transformation $T^{*}$.
Meanwhile, the outliers $\mathcal{N}=\mathcal{P}^{C}\cap\mathcal{M}$ will be
scattered outside the surface $\mathcal{P}$. To identify the inliers, we use a
convolutional network. Such networks have been proven effective in related
dense prediction tasks, such as 3D point cloud segmentation [16, 7]. The
convolutional network in our setting is in 6-dimensional space [8]. The
network predicts a likelihood for each correspondence, which is a point in 6D
space $[\mathbf{x}_{i}^{T},\mathbf{y}_{j}^{T}]^{T}$. The prediction is
interpreted as the likelihood that the correspondence is true: an inlier.

Note that the convolution operator is translation invariant, thus our 6D
ConvNet will generate the same output regardless of the absolute position of
inputs in 3D. We use a similar network architecture to Choy et al. [9] to
create a 6D convolutional network with skip connections within the spatial
resolution across the network. The architecture of the 6D ConvNet is shown in
Fig. 2. During training, we use the binary cross-entropy loss between the
likelihood prediction that a correspondence $(i,j)$ is an inlier,
$p_{(i,j)}\in[0,1]$, and the ground-truth correspondences $\mathcal{P}$ to
optimize the network parameters:

$\small
L_{\text{bce}}(\mathcal{M},T^{*})=\frac{1}{|\mathcal{M}|}\bigg{(}\sum_{(i,j)\in\mathcal{P}}\log
p_{(i,j)}+\sum_{(i,j)\in\mathcal{N}}\log p_{(i,j)}^{C}\bigg{)},$ (1)

where $p^{C}=1-p$ and $|\mathcal{M}|$ is the cardinality of the set of
putative correspondences.

Figure 2: 6-dimensional convolutional network architecture for inlier
likelihood prediction (Sec. 3.2). The network has a U-net structure with
residual blocks between strided convolutions. Best viewed on the screen.

###  3.3 Weighted Procrustes for $\text{SE}(3)$

The inlier likelihood estimated by the 6D ConvNet provides a weight for each
correspondence. The original Procrustes method [15] minimizes the mean squared
error between corresponding points
$\frac{1}{N}\sum_{(i,j)\in\mathcal{M}}\|\mathbf{x}_{i}-\mathbf{y}_{j}\|^{2}$
and thus gives equal weight to all correspondences. In contrast, we minimize a
weighted mean squared error
$\sum_{(i,j)\in\mathcal{M}}w_{(i,j)}\|\mathbf{x}_{i}-\mathbf{y}_{j}\|^{2}$.
This change allows us to pass gradients through the weights, rather than
through the position [41], and enables the optimization to scale to dense
correspondence sets.

Formally, Weighted Procrustes analysis minimizes:

$\displaystyle e^{2}$ $\displaystyle=e^{2}(R,\mathbf{t};\mathbf{w},X,Y)$ (2)
$\displaystyle=\sum_{(i,j)\in\mathcal{M}}\tilde{w}_{(i,j)}(\mathbf{y}_{j}-(R\mathbf{x}_{i}+\mathbf{t}))^{2}$
(3)
$\displaystyle=\text{Tr}\left((Y-RX-\textbf{t}\textbf{1}^{T})W(Y-RX-\textbf{t}\textbf{1}^{T})^{T}\right),$
(4)

where $\mathbf{1}=(1,...,1)^{T}$,
$X=[\mathbf{x}_{1},...,\mathbf{x}_{|\mathcal{M}|}]$, and
$Y=[\mathbf{y}_{J_{1}},...,\mathbf{y}_{J_{|\mathcal{M}|}}]$. $J$ is a list of
indices that defines the correspondences
$\mathbf{x}_{i}\leftrightarrow\mathbf{y}_{J_{i}}$.
$\mathbf{w}=[w_{1},\cdots,w_{|\mathcal{M}|}]$ is the weight vector and
$\tilde{\mathbf{w}}=[\tilde{w}_{1},\cdots,\tilde{w}_{|\mathcal{M}|}]\triangleq\frac{\phi(\mathbf{w})}{||\phi(\mathbf{w})||_{1}}$
denotes the normalized weight after a nonlinear transformation $\phi$ that
applies heuristic prefiltering. $W=\text{diag}(\tilde{\mathbf{w}})$ forms the
diagonal weight matrix.

###### Theorem 1

: The $R$ and t that minimize the squared error
$e^{2}(R,\textbf{t})=\sum_{(i,j)}w_{(i,j)}(\textbf{y}_{j}-R\textbf{x}_{i}-\textbf{t})^{2}$
are $\hat{\textbf{t}}=(Y-RX)W\mathbf{1}$ and $\hat{R}=USV^{T}$ where $U\Sigma
V^{T}=\text{SVD}(\Sigma_{xy})$, $\Sigma_{xy}=YKWKX^{T}$,
$K=I-\sqrt{\tilde{\mathbf{w}}}\sqrt{\tilde{\mathbf{w}}}^{T}$, and
$S=\text{diag}\left(1,\cdots,1,\text{det}(U)\text{det}(V)\right)$.

Sketch of proof. First, we differentiate $e^{2}$ w.r.t. $\mathbf{t}$ and
equate the partial derivative to 0. This gives us
$\hat{\mathbf{t}}=(Y-RX)W\mathbf{1}$. Next, we substitute
$X=KX+X\sqrt{\tilde{\mathbf{w}}}\sqrt{\tilde{\mathbf{w}}}^{T}$ on Eq. 4 and do
the same for $Y$. Then, we substitute $\mathbf{t}=\hat{\mathbf{t}}$ into the
squares. This yields $e^{2}=\text{Tr}\left((Y-RX)KWK^{T}(Y-RX)^{T}\right)$ and
expanding the term results in two squares plus
$-2\text{Tr}(YKWK^{T}X^{T}R^{T})$. We maximize the last negative term whose
maximum is the sum of all singular values, which leads to $\hat{R}$. The full
derivation is in the supplement. $\square$

We can easily extend the above theorem to incorporate a scaling factor
$c\in\mathbb{R}^{+}$, or anisotropic scaling for tasks such as scan-to-CAD
registration, but in this paper we assume that partial scans of the same scene
have the same scale.

The Weighted Procrustes method generates rotation $\hat{R}$ and translation
$\hat{\mathbf{t}}$ as outputs that depend on the weight vector $\mathbf{w}$.
In our current implementation, $\hat{R}$ and $\hat{\mathbf{t}}$ are directly
sent to the robust registration module in Section 4 as an initial pose.
However, we briefly demonstrate that they can also be embedded in an end-to-
end registration pipeline, since Weighted Procrustes is differentiable. From a
top-level loss function $L$ of $\hat{R}$ and $\hat{\mathbf{t}}$, we can pass
the gradient through the closed-form solver, and update parameters in
downstream modules:

$\displaystyle\frac{\partial}{\partial\mathbf{w}}L(\hat{R},\hat{\mathbf{t}})=$
$\displaystyle\frac{\partial
L(\hat{R},\hat{\mathbf{t}})}{\partial\hat{R}}\frac{\partial\hat{R}}{\partial\hat{\mathbf{w}}}+\frac{\partial
L(\hat{R},\hat{\mathbf{t}})}{\partial\hat{\mathbf{t}}}\frac{\partial\hat{\mathbf{t}}(\hat{R},\hat{\mathbf{w}})}{\partial\hat{\mathbf{w}}},$
(5)

where $L(\hat{R},\hat{\mathbf{t}})$ can be defined as the combination of
differentiable rotation error (RE) and translation error (TE) between
predictions $\hat{R},\hat{\mathbf{t}}$ and ground-truth
$R^{*},\mathbf{t}^{*}$:

$\displaystyle L_{\text{rot}}(\hat{R})=$
$\displaystyle\arccos{\frac{\text{Tr}(\hat{R}^{T}R^{*})-1}{2}},$ (6)
$\displaystyle L_{\text{trans}}(\hat{\mathbf{t}})=$
$\displaystyle||\hat{\mathbf{t}}-\mathbf{t}^{*}||_{2}^{2},$ (7)

or the Forbenius norm of relative transformation matrices defined in [2, 41].
The final loss is the weighted sum of $L_{\text{rot}}$, $L_{\text{trans}}$,
and $L_{\text{bce}}$.

##  4 Robust Registration

In this section, we propose a fine-tuning module that minimizes a robust loss
function of choice to improve the registration accuracy. We use a gradient-
based method to refine poses, where a continuous representation [52] for
rotations is adopted to remove discontinuities and construct a smooth
optimization space. This module initializes the pose from the prediction of
the Weighted Procrustes method. During iterative optimization, unlike Maken
_et al_. [27], who find the nearest neighbor per point at each gradient step,
we rely on the correspondence likelihoods from the 6D ConvNet, which is
estimated only once per initialization.

In addition, our framework naturally offers a failure detection mechanism. In
practice, Weighted Procrustes may generate numerically unstable solutions when
the number of valid correspondences is insufficient due to small overlaps or
noisy correspondences between input scans. By computing the ratio of the sum
of the filtered weights to the total number of correspondences, _i.e_.
$\sum_{i}\phi(w_{i})/|\mathcal{M}|$, we can easily approximate the fraction of
valid correspondences and predict whether an alignment may be unstable. When
this fraction is low, we resort to a more time-consuming but accurate
registration algorithm such as RANSAC [38, 1, 35] or a branch-and-bound method
[45] to find a numerically stable solution. In other words, we can detect when
our system might fail before it returns a result and fall back to a more
accurate but time-consuming algorithm, unlike previous end-to-end methods that
use globally pooled latent features [2] or a singly stochastic matrix [41] –
such latent representations are more difficult to interpret.

###  4.1 SE(3) Representation and Initialization

We use the 6D representation of 3D rotation proposed by Zhou  _et al_. [52],
rather than Euler angles or quaternions. The new representation uses 6
parameters $\mathbf{a}_{1},\mathbf{a}_{2}\in\mathbb{R}^{3}$ and can be
transformed into a $3\times 3$ orthogonal matrix by

$f\left(\begin{bmatrix}|&|\\\ \mathbf{a}_{1}&\mathbf{a}_{2}\\\
|&|\end{bmatrix}\right)=\begin{bmatrix}|&|&|\\\
\mathbf{b}_{1}&\mathbf{b}_{2}&\mathbf{b}_{3}\\\ |&|&|\end{bmatrix},$ (8)

where $\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3}\in\mathbb{R}^{3}$ are
$\mathbf{b}_{1}=N(\mathbf{a}_{1})$,
$\mathbf{b}_{2}=N(\mathbf{a}_{2}-(\mathbf{b}_{1}\cdot\mathbf{a}_{2})\mathbf{b}_{1})$,
and $\mathbf{b}_{3}=\mathbf{b}_{1}\times\mathbf{b}_{2}$, and $N(\cdot)$
denotes L2 normalization. Thus, the final representation that we use is
$\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{t}$ which are equivalent to
$R,\mathbf{t}$ using Eq. 8.

To initialize $\mathbf{a}_{1},\mathbf{a}_{2}$, we simply use the first two
columns of the rotation matrix $R$, i.e., $\mathbf{b}_{1}$, $\mathbf{b}_{2}$.
For convenience, we define $f^{-1}$ as $f^{-1}(f(R))=R$ though this inverse
function is not unique as there are infinitely many choices of
$\mathbf{a}_{1},\mathbf{a}_{2}$ that map to the same $R$.

###  4.2 Energy Minimization

We use a robust loss function to fine-tune the registration between predicted
inlier correspondences. The general form of the energy function is

$E(R,\mathbf{t})=\sum_{i=1}^{n}\phi(w_{(i,J_{i})})L(\mathbf{y}_{J_{i}},R\mathbf{x}_{i}+\mathbf{t}),$
(9)

where $\tilde{w}_{i}$ and $J_{i}$ are defined as in Eq. 3 and $\phi(\cdot)$ is
a prefiltering function. In the experiments, we use $\phi(w)=I[w>\tau]w$,
which clips weights below $\tau$ elementwise as neural network outputs bounded
logit scores. $L(\mathbf{x},\mathbf{y})$ is a pointwise loss function between
$\mathbf{x}$ and $\mathbf{y}$; we use the Huber loss in our implementation.
The energy function is parameterized by $R$ and $\mathbf{t}$ which in turn are
represented as $\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{t}$. We can apply first-
order optimization algorithms such as SGD, Adam, etc. to minimize the energy
function, but higher-order optimizers are also applicable since the number of
parameters is small. The complete algorithm is described in Alg. 1.

Input: $X\in\mathbb{R}^{n\times 3},Y\in\mathbb{R}^{m\times 3}$

Output: $R\in\text{SO}(3),\mathbf{t}\in\mathbb{R}^{3\times 1}$

$\mathcal{F}_{x}\leftarrow\text{Feature}(X)$

// § 3.1

1 $\mathcal{F}_{y}\leftarrow\text{Feature}(Y)$

$J_{x\rightarrow
y}\leftarrow\text{NearestNeighbor}(\mathcal{F}_{x},\mathcal{F}_{y})$

// § 3.2

2 $\mathcal{M}\leftarrow\\{(i,J_{x\rightarrow y,i})\mid i\in[1,...,n]\\}$

3 $\mathbf{w}\leftarrow\text{InlierProbability}(\mathcal{M})$

4

5if  _$\mathbb{E}_{i}\phi(w_{i}) <\tau_{s}$_ then

return SafeGuardRegistration$(X,Y)$

// §4

6

7 else

$\hat{R},\hat{\mathbf{t}}\leftarrow\mathop{\mathrm{argmin}}_{R,\mathbf{t}}e^{2}(R,\mathbf{t};\mathbf{w},X,Y)$

// § 3.3

$\mathbf{a}\leftarrow f^{-1}({\hat{R}}),\mathbf{t}\leftarrow\hat{\mathbf{t}}$

// § 4.1

8 while  _not converging_ do

9
$\ell\leftarrow\sum_{(i,j)\in\mathcal{M}}\phi(w_{(i,j)})L(Y_{j},f(\mathbf{a})X_{i}+\mathbf{t})$

10
$\mathbf{a}\leftarrow\text{Update}(\mathbf{a},\frac{\partial}{\partial\mathbf{a}}\ell(\mathbf{a},\mathbf{t}))$

11
$\mathbf{t}\leftarrow\text{Update}(\mathbf{t},\frac{\partial}{\partial\mathbf{t}}\ell(\mathbf{a},\mathbf{t}))$

12 return $f(\mathbf{a}),\mathbf{t}$

Algorithm 1 Deep Global Registration

##  5 Experiments

We analyze the proposed model in two registration scenarios: pairwise
registration where we estimate an $\text{SE}(3)$ transformation between two 3D
scans or fragments, and multi-way registration which generates a final
reconstruction and camera poses for all fragments that are globally
consistent. Here, pairwise registration serves as a critical module in multi-
way registration.

For pairwise registration, we use the 3DMatch benchmark [47] which consists of
3D point cloud pairs from various real-world scenes with ground truth
transformations estimated from RGB-D reconstruction pipelines [17, 10]. We
follow the train/test split and the standard procedure to generate pairs with
at least 30% overlap for training and testing [12, 11, 9]. For multi-way
registration, we use the simulated Augmented ICL-NUIM dataset [6, 18] for
quantitative trajectory results, and Indoor LiDAR RGB-D dataset [32] and
Stanford RGB-D dataset [6] for qualitative registration visualizations. Note
in this experiment we use networks trained on the 3DMatch training set and do
not fine-tune on the other datasets. This illustrates the generalization
abilities of our models. Lastly, we use KITTI LIDAR scans [13] for outdoor
pairwise registration. As the official registration splits do not have labels
for pairwise registration, we follow Choy  _et al_. [9] to create pairwise
registration train/val/test splits.

For all indoor experiments, we use 5cm voxel downsampling [35, 51], which
randomly subsamples a single point within each 5cm voxel to generate point
clouds with uniform density. For safeguard registration, we use RANSAC and the
safeguard threshold $\tau_{s}=0.05$, which translates to 5% of the
correspondences should be valid. We train learning-based state-of-the-art
models and our network on the training split of the 3DMatch benchmark. During
training, we augment data by applying random rotations varying from $-$180 to
180 degrees around a random axis. Ground-truth pointwise correspondences are
found using nearest neighbor search in 3D space. We train the 6-dimensional
ConvNet on a single Titan XP with batch size 4. SGD is used with an initial
learning rate $10^{-1}$ and an exponential learning rate decay factor 0.99.

###  5.1 Pairwise Registration

Figure 3: Global registration results of our method on all 8 different test
scenes in 3DMatch [47]. Best viewed in color. Figure 4: Overall pairwise
registration recall (y-axis) on the 3DMatch benchmark with varying rotation
(left image) and translation (right image) error thresholds (x-axis). Our
approach outperforms baseline methods for all thresholds while being
$6.5\times$ faster than the most accurate baseline.

In this section, we report the registration results on the test set of the
3DMatch benchmark [47], which contains 8 different scenes as depicted in Fig.
3. We measure translation error (TE) defined in Eq. 6, rotation error (RE)
defined in Eq. 7, and recall. Recall is the ratio of successful pairwise
registrations and we define a registration to be successful if its rotation
error and translation error are smaller than predefined thresholds. Average TE
and RE are computed only on these successfully registered pairs since failed
registrations return poses that can be drastically different from the ground
truth, making the error metrics unreliable.

Figure 5: Analysis of 3DMatch registration results per scene. Row 1: recall
rate (higher is better). Row 2-3: TE and RE measured on successfully
registered pairs (lower is better). Our method is consistently better on all
scenes, which were not seen during training. Note: a missing bar corresponds
to zero successful alignments in a scene. Table 1: Row 1-6: registration
results of our method and classical global registration methods on point
clouds voxelized with 5cm voxel size. Our method outperforms RANSAC and FGR
while being as fast as FGR. Row 7-10: results of ICP variants. Row 11 - 12:
results of learning-based methods. The learning-based methods generally fail
on real-world scans. Time includes feature extraction.

| Recall | TE (cm) | RE (deg) | Time (s)
---|---|---|---|---
Ours w/o safeguard | 85.2% | 7.73 | 2.58 | 0.70
Ours |  91.3% | 7.34 | 2.43 | 1.21
FGR [50] | 42.7% | 10.6 | 4.08 | 0.31
RANSAC-2M [35] | 66.1% | 8.85 | 3.00 | 1.39
RANSAC-4M | 70.7% | 9.16 | 2.95 | 2.32
RANSAC-8M | 74.9% | 8.96 | 2.92 | 4.55
Go-ICP [45] | 22.9% | 14.7 | 5.38 | 771.0
Super4PCS [29] | 21.6% | 14.1 | 5.25 | 4.55
ICP (P2Point) [51] | 6.04% | 18.1 | 8.25 | 0.25
ICP (P2Plane) [51] | 6.59% | 15.2 | 6.61 | 0.27
DCP [41] | 3.22% | 21.4 | 8.42 | 0.07
PointNetLK [2] | 1.61% | 21.3 | 8.04 | 0.12

Table 2: ATE (cm) error on the Augmented ICL-NUIM dataset with simulated depth
noise. For InfiniTAM, the loop closure module is disabled since it fails in
all scenes. For BAD-SLAM, the loop closure module only succeeds in ‘Living
room 2’.

|  ElasticFusion [43] |  InfiniTAM [21] |  BAD-SLAM [39] |  Multi-way + FGR [50] |  Multi-way + RANSAC [51] | Multi-way + Ours
---|---|---|---|---|---|---
Living room 1 | 66.61 | 46.07 | fail | 78.97 | 110.9 | 21.06
Living room 2 | 24.33 | 73.64 | 40.41 | 24.91 | 19.33 | 21.88
Office 1 | 13.04 | 113.8 | 18.53 | 14.96 | 14.42 | 15.76
Office 2 | 35.02 | 105.2 | 26.34 | 21.05 | 17.31 | 11.56
Avg. Rank | 3 | 5 | 5 | 3.5 | 2.5 | 2

(a) Real-world: Apartment

(b) Real-world: Boardroom

(c) Synthetic: Office

(d) Real-world: Copyroom

(e) Real-world: Loft

(f) Synthetic: Livingroom

Figure 6: Fragment registrations on [6, 32]. From left to right: FGR [50],
RANSAC [35], Ours. Row 1-3: our method succeeds on scenes with small overlaps
or ambiguous geometry structures while other methods fail. Row 4-6: by
combining Weighted Procrustes and gradient-based refinement, our method
outputs more accurate registrations in one pass, leading to better aligned
details.

We compare our methods with various classical methods [50, 35, 45] and state-
of-the-art learning based methods [41, 42, 2, 31]. All the experiments are
evaluated on an Intel i7-7700 CPU and a GTX 1080Ti graphics card except for
Go-ICP [45] tested on an Intel i7-5820K CPU. In Table 1, we measure recall
with the TE threshold 30cm which is typical for indoor scene relocalization
[30], and RE threshold 15 degrees which is practical for partially overlapping
scans from our experiments. In Fig. 4, we plot the sensitivity of recall on
both thresholds by changing one threshold and setting the other to infinity.
Fig. 5 includes detailed statistics on separate test scenes. Our system
outperforms all the baselines on recall by a large margin and achieves the
lowest translation and rotation error consistently on most scenes.

Classical methods. To compare with classical methods, we evaluate point-to-
point ICP, Point-to-plane ICP, RANSAC [35], and FGR [50], all implemented in
Open3D [51]. In addition, we test the open-source Python bindings of Go-ICP
[45] and Super4PCS [29]. For RANSAC and FGR, we extract FPFH from voxel-
downsampled point clouds. The results are shown in Table 1.

ICP variants mostly fail as the dataset contains challenging 3D scan sequences
with small overlap and large camera viewpoint change. Super4PCS, a sampling-
based algorithm, performs similarly to Go-ICP, an ICP variant with branch-and-
bound search.

Feature-based methods, FGR and RANSAC, perform better. When aligning 5cm-
voxel-downsampled point clouds, RANSAC achieves recall as high as 70%, while
FGR reaches 40%. Table 1 also shows that increasing the number of RANSAC
iterations by a factor of 2 only improves performance marginally. Note that
our method is about twice as fast as RANSAC with 2M iterations while achieving
higher recall and registration accuracy.

Learning-based methods. We use 3DRegNet [31], Deep Closest Point (DCP) [41],
PRNet [42], and PointNetLK [2] as our baselines. We train all the baselines on
3DMatch with the same setup and data augmentation as ours for all experiments.

For 3DRegNet, we follow the setup outlined in [31], except that we do not
manually filter outliers with ground truth, and train and test with the
standard realistic setup. We find that the registration loss of 3DRegNet does
not converge during training and the rotation and translation errors are
consistently above 30 degrees and 1m during test.

We train Deep Closest Point (DCP) with 1024 randomly sampled points for each
point cloud for 150 epochs [41]. We initialize the network with the pretrained
weights provided by the authors. Although the training loss converges, DCP
fails to achieve reasonable performance for point clouds with partial overlap.
DCP uses a singly stochastic matrix to find correspondences, but this
formulation assumes that all points in point cloud $X$ have at least one
corresponding point in the convex hull of point cloud $Y$. This assumption
fails when some points in $X$ have no corresponding points in $Y$, as is the
case for partially overlapping fragments. We also tried to train PRNet [42] on
our setup, but failed to get reasonable results due to random crashes and
high-variance training losses.

Lastly, we fine-tune PointNetLK [2] on 3DMatch for 400 epochs, starting from
the pretrained weights provided by the authors. PointNetLK uses a single
feature that is globally pooled for each point cloud and regresses the
relative pose between objects, and we suspect that a globally pooled feature
fails to capture complex scenes such as 3DMatch.

In conclusion, while working well on object-centric synthetic datasets,
current end-to-end registration approaches fail on real-world data. Unlike
synthetic data, real 3D point cloud pairs contain multiple objects, partial
scans, self-occlusion, substantial noise, and may have only a small degree of
overlap between scans.

###  5.2 Multi-way Registration

Multi-way registration for RGB-D scans proceeds via multiple stages. First,
the pipeline estimates the camera pose via off-the-shelf odometry and
integrates multiple 3D scans to reduce noise and generate accurate 3D
fragments of a scene. Next, a pairwise registration algorithm roughly aligns
all fragments, followed by multi-way registration [6] which optimizes fragment
poses with robust pose graph optimization [23].

We use a popular open-source implementation of this registration pipeline [51]
and replace the pairwise registration stage in the pipeline with our proposed
modules. Note that we use the networks trained on the 3DMatch training set and
test on the multi-way registration datasets [18, 32, 6]; this demonstrates
cross-dataset generalization.

We test the modified pipeline on the Augmented ICL-NUIM dataset [6, 18] for
quantitative trajectory results, and Indoor LiDAR RGB-D dataset [32] and
Stanford RGB-D dataset [6] for qualitative registration visualizations. We
measure the absolute trajectory error (ATE) on the Augmented ICL-NUIM dataset
with simulated depth noise. As shown in Table 2, compared to state-of-the-art
online SLAM [43, 21, 39] and offline reconstruction methods [50], our approach
yields consistently low error across scenes.

For qualitative results, we compare pairwise fragment registration on these
scenes against FGR and RANSAC in Fig. 6. Full scene reconstruction results are
shown in the supplement.

###  5.3 Outdoor LIDAR Registration

We use outdoor LIDAR scans from the KITTI dataset [13] for registration,
following [9]. The registration split of Choy  _et al_. [9] uses GPS-IMU to
create pairs that are at least 10m apart and generated ground-truth
transformation using GPS followed by ICP to fix errors in GPU readings. We use
FCGF features [9] trained on the training set of the registration split to
find the correspondences and trained the 6D ConvNet for inlier confidence
prediction similar to how we trained the system for indoor registration. We
use voxel size 30cm for downsampling point clouds for all experiments.
Registration results are reported in Tab. 3 and visualized in Fig. 7.

(a) KITTI registration test pair 1

(b) KITTI registration test pair 2

(c) KITTI registration test pair 3

Figure 7: KITTI registration results. Pairs of scans are at least 10m apart.
Table 3: Registration on the KITTI test split [13, 9]. We use thresholds of
0.6m and 5 degrees. ‘Ours + ICP’ refers to our method followed by ICP for
fine-grained pose adjustment. The runtime includes feature extraction.

| Recall | TE (cm) | RE (deg) | Time (s)
---|---|---|---|---
FGR [50] | 0.2% | 40.7 | 1.02 | 1.42
RANSAC [35] | 34.2% | 25.9 | 1.39 | 1.37
FCGF [9] | 98.2% | 10.2 | 0.33 | 6.38
Ours | 96.9% | 21.7 | 0.34 | 2.29
Ours + ICP | 98.0% | 3.46 | 0.14 | 2.51

##  6 Conclusion

We presented Deep Global Registration, a learning-based framework that
robustly and accurately aligns real-world 3D scans. To achieve this, we used a
6D convolutional network for inlier detection, a differentiable Weighted
Procrustes algorithm for scalable registration, and a gradient-based optimizer
for pose refinement. Experiments show that our approach outperforms both
classical and learning-based registration methods, and can serve as a ready-
to-use plugin to replace alternative registration methods in off-the-shelf
scene reconstruction pipelines.

## References

  * [1] Dror Aiger, Niloy J Mitra, and Daniel Cohen-Or.  4-points congruent sets for robust pairwise surface registration.  ACM Transactions on Graphics, 2008.
  * [2] Yasuhiro Aoki, Hunter Goforth, Rangaprasad Arun Srivatsan, and Simon Lucey.  PointNetLK: Robust & efficient point cloud registration using PointNet.  In CVPR, 2019.
  * [3] Paul J Bes, Neil D McKay, et al.  A method for registration of 3-d shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992\.
  * [4] Sofien Bouaziz, Andrea Tagliasacchi, and Mark Pauly.  Sparse iterative closest point.  In Eurographics, 2013.
  * [5] Qin Cai, David Gallup, Cha Zhang, and Zhengyou Zhang.  3D deformable face tracking with a commodity depth camera.  In ECCV, 2010.
  * [6] Sungjoon Choi, Qian-Yi Zhou, and Vladlen Koltun.  Robust reconstruction of indoor scenes.  In CVPR, 2015.
  * [7] Christopher Choy, JunYoung Gwak, and Silvio Savarese.  4D spatio-temporal ConvNets: Minkowski convolutional neural networks.  In CVPR, 2019.
  * [8] Christopher Choy, Junha Lee, René Ranftl, Jaesik Park, and Vladlen Koltun.  High-dimensional convolutional networks for geometric pattern recognition.  In CVPR, 2020.
  * [9] Christopher Choy, Jaesik Park, and Vladlen Koltun.  Fully convolutional geometric features.  In ICCV, 2019.
  * [10] Angela Dai, Matthias Nießner, Michael Zollöfer, Shahram Izadi, and Christian Theobalt.  BundleFusion: Real-time globally consistent 3D reconstruction using on-the-fly surface re-integration.  ACM Transactions on Graphics, 2017.
  * [11] Haowen Deng, Tolga Birdal, and Slobodan Ilic.  PPF-FoldNet: Unsupervised learning of rotation invariant 3D local descriptors.  In ECCV, 2018.
  * [12] Haowen Deng, Tolga Birdal, and Slobodan Ilic.  PPFNet: Global context aware local features for robust 3D point matching.  In CVPR, 2018.
  * [13] Andreas Geiger, Philip Lenz, and Raquel Urtasun.  Are we ready for autonomous driving? the kitti vision benchmark suite.  In CVPR, 2012.
  * [14] Zan Gojcic, Caifa Zhou, Jan Dirk Wegner, and Andreas Wieser.  The perfect match: 3D point cloud matching with smoothed densities.  In CVPR, 2019.
  * [15] John C Gower.  Generalized procrustes analysis.  Psychometrika, 1975.
  * [16] Benjamin Graham, Martin Engelcke, and Laurens van der Maaten.  3D semantic segmentation with submanifold sparse convolutional networks.  In CVPR, 2018.
  * [17] Maciej Halber and Thomas Funkhouser.  Fine-to-coarse global registration of RGB-D scans.  In CVPR, 2017.
  * [18] Ankur Handa, Thomas Whelan, John McDonald, and Andrew J. Davison.  A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM.  In ICRA, 2014.
  * [19] Dirk Holz, Alexandru E Ichim, Federico Tombari, Radu B Rusu, and Sven Behnke.  Registration with the point cloud library: A modular framework for aligning in 3-D.  RA-L, 2015.
  * [20] Andrew E. Johnson and Martial Hebert.  Using spin images for efficient object recognition in cluttered 3D scenes.  Pattern Analysis and Machine Intelligence, 1999.
  * [21] Olaf Kähler, V. A. Prisacariu, and David W. Murray.  Real-time large-scale dense 3D reconstruction with loop closure.  In ECCV, 2016.
  * [22] Marc Khoury, Qian-Yi Zhou, and Vladlen Koltun.  Learning compact geometric features.  In ICCV, 2017.
  * [23] Rainer Kümmerle, Giorgio Grisetti, Hauke Strasdat, Kurt Konolige, and Wolfram Burgard.  g2o: A general framework for graph optimization.  In ICRA, 2011.
  * [24] Bruce D Lucas and Takeo Kanade.  An iterative image registration technique with an application to stereo vision.  In DARPA Image Understanding Workshop, 1981.
  * [25] João Maciel and Joao Paulo Costeira.  A global solution to sparse correspondence problems.  Pattern Analysis and Machine Intelligence, 2003.
  * [26] Ameesh Makadia, Alexander Patterson, and Kostas Daniilidis.  Fully automatic registration of 3D point clouds.  In CVPR, 2006.
  * [27] Fahira Afzal Maken, Fabio Ramos, and Lionel Ott.  Speeding up iterative closest point using stochastic gradient descent.  In ICRA, 2019.
  * [28] Haggai Maron, Nadav Dym, Itay Kezurer, Shahar Kovalsky, and Yaron Lipman.  Point registration via efficient convex relaxation.  ACM Transactions on Graphics, 2016.
  * [29] Nicolas Mellado, Dror Aiger, and Niloy J Mitra.  Super 4PCS fast global pointcloud registration via smart indexing.  Computer Graphics Forum, 2014.
  * [30] Montiel J. M. M. Mur-Artal, Raúl and Juan D. Tardós.  ORB-SLAM: A versatile and accurate monocular SLAM system.  IEEE Transactions on Robotics, 2015.
  * [31] G. Dias Pais, Pedro Miraldo, Srikumar Ramalingam, Jacinto C. Nascimento, Venu Madhav Govindu, and Rama Chellappa.  3DRegNet: A deep neural network for 3D point registration.  arXiv, 2019.
  * [32] Jaesik Park, Qian-Yi Zhou, and Vladlen Koltun.  Colored point cloud registration revisited.  In ICCV, 2017.
  * [33] Chen Qian, Xiao Sun, Yichen Wei, Xiaoou Tang, and Jian Sun.  Realtime and robust hand tracking from depth.  In CVPR, 2014.
  * [34] René Ranftl and Vladlen Koltun.  Deep fundamental matrix estimation.  In ECCV, 2018.
  * [35] R. B. Rusu, N. Blodow, and M. Beetz.  Fast point feature histograms (FPFH) for 3D registration.  In ICRA, 2009.
  * [36] R. B. Rusu, N. Blodow, Z. C. Marton, and M. Beetz.  Aligning point cloud views using persistent feature histograms.  In IROS, 2008.
  * [37] S. Salti, F. Tombari, and L. di Stefano.  SHOT: Unique signatures of histograms for surface and texture description.  CVIU, 2014.
  * [38] Ruwen Schnabel, Roland Wahl, and Reinhard Klein.  Efficient RANSAC for point-cloud shape detection.  Computer Graphics Forum, 2007.
  * [39] Thomas Schops, Torsten Sattler, and Marc Pollefeys.  BAD SLAM: Bundle adjusted direct RGB-D SLAM.  In CVPR, 2019.
  * [40] F. Tombari, S. Salti, and L. Di Stefano.  Unique shape context for 3D data description.  In ACM Workshop on 3D Object Retrieval, 2010.
  * [41] Yue Wang and Justin M. Solomon.  Deep closest point: Learning representations for point cloud registration.  In ICCV, 2019.
  * [42] Yue Wang and Justin M. Solomon.  PRNet: Self-supervised learning for partial-to-partial registration.  In NeurIPS, 2019.
  * [43] Thomas Whelan, Stefan Leutenegger, Renato F. Salas-Moreno, Ben Glocker, and Andrew J. Davison.  ElasticFusion: Dense SLAM without a pose graph.  In Robotics: Science and Systems, 2015.
  * [44] Heng Yang and Luca Carlone.  A polynomial-time solution for robust registration with extreme outlier rates.  In Robotics: Science and Systems, 2019.
  * [45] Jiaolong Yang, Hongdong Li, Dylan Campbell, and Yunde Jia.  Go-ICP: A globally optimal solution to 3D ICP point-set registration.  Pattern Analysis and Machine Intelligence, 2015.
  * [46] Kwang Moo Yi, Eduard Trulls, Yuki Ono, Vincent Lepetit, Mathieu Salzmann, and Pascal Fua.  Learning to find good correspondences.  In CVPR, 2018.
  * [47] Andy Zeng, Shuran Song, Matthias Nießner, Matthew Fisher, Jianxiong Xiao, and Thomas Funkhouser.  3DMatch: Learning local geometric descriptors from RGB-D reconstructions.  In CVPR, 2017.
  * [48] Zhengyou Zhang.  Iterative point matching for registration of free-form curves and surfaces.  IJCV, 1994.
  * [49] Chen Zhao, Zhiguo Cao, Chi Li, Xin Li, and Jiaqi Yang.  NM-Net: Mining reliable neighbors for robust feature correspondences.  In CVPR, 2019.
  * [50] Qian-Yi Zhou, Jaesik Park, and Vladlen Koltun.  Fast global registration.  In ECCV, 2016.
  * [51] Qian-Yi Zhou, Jaesik Park, and Vladlen Koltun.  Open3D: A modern library for 3D data processing.  arXiv, 2018.
  * [52] Yi Zhou, Connelly Barnes, Jingwan Lu, Jimei Yang, and Hao Li.  On the continuity of rotation representations in neural networks.  In CVPR, 2019.


# Classifying basins of attraction using the basin entropy

Alvar Daza  Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento
de Física, Universidad Rey Juan Carlos, Móstoles, Madrid, Tulipán s/n, 28933,
Spain<EMAIL_ADDRESS>Alexandre Wagemakers  Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Móstoles, Madrid, Tulipán s/n, 28933, Spain  Miguel A.F. Sanjuán  Nonlinear
Dynamics, Chaos and Complex Systems Group, Departamento de Física, Universidad
Rey Juan Carlos, Móstoles, Madrid, Tulipán s/n, 28933, Spain

###### Abstract

A basin of attraction represents the set of initial conditions leading to a
specific asymptotic state of a given dynamical system. Here, we provide a
classification of the most common basins found in nonlinear dynamics with the
help of the basin entropy. We have also found interesting connections between
the basin entropy and other measures used to characterize the unpredictability
associated to the basins of attraction, such as the uncertainty exponent, the
lacunarity or other different parameters related to the Wada property.

##  I Introduction

Multistability is a common phenomenon in nature [1]. From a deterministic
point of view, it implies the existence of different asymptotic states in
phase space for different choices of initial conditions. The set of all the
initial conditions leading to one particular fate, such as an attractor or an
exit of the system, is called a basin of attraction or an escape basin [2].
Since the early days of chaos theory, basins have been the object of extensive
study. For example, one of the many facets of dissipative chaos is the fractal
nature of the basin boundaries [3, 4]. Along the last decades, a plethora of
different kinds of basins with special characteristics has been reported:
locally connected and disconnected [3], Wada [5], riddled [6, 7], intermingled
[8], etc. At the same time, a variety of different tools has arisen to
describe and characterize them: the uncertainty exponent [9], the lacunarity
[10], the basin stability [11], the Wada parameter [12], the Wada index [13],
etc.

Among them, the basin entropy [14] integrates some previous measures in a
natural way. It is extremely useful to characterize the unpredictability
associated to the basins in different scientific contexts [15, 16, 17, 18,
19]. The present paper analyzes the different kinds of basins (Wada, riddled,
intermingled, etc.) under the viewpoint of the basin entropy. From this
inquiry, a classification emerges where iconic types of basins maximize one
aspect of this measure. Needless to say, the basin taxonomy proposed here is
just one of the many that can be chosen attending to different criteria.
Nevertheless, it is striking to observe how well the paradigmatic basins of
attraction fit into this classification. Probably, the main contribution of
the present work is to provide a framework to understand better the
unpredictability associated to the different types of basins. Such framework
is based on the basin entropy, and goes beyond of the mere quantification of
the unpredictability: it lays the ground to a deeper understanding of concepts
such as fractality and smoothness, Wada boundaries, riddled basins, symmetry
of the basins and more.

The paper is organized as follows. In Sec. II, we review the basics of the
basin entropy, its definition and interpretation. Sections III-V are devoted
to analyze the role of each ingredient of the basin entropy and its relation
to some of the most common types of basins. We also dissect the relation of
the basin entropy with other relevant measurements for each case. Finally, in
Sec. VI, we summarize our results providing a rational classification of the
basins according to their unpredictability, as measured by the basin entropy.

##  II Basics of the basin entropy

Since its appearance [20], the basin entropy has been prolifically used in
different contexts and with different aims. For instance, it is an efficient
alternative to the box-counting algorithm for estimating the Hausdorff
dimension [21]. It also has been successfully applied to the detection of KAM
islands in conservative systems [15]. Furthermore, it has been applied to
investigate features of spatial patterns emerging in networks with competing
species, using a model inspired in the rock-paper-scissors game [17].

Among other purposes, the basin entropy has been a response to solve the
vagueness of some affirmations concerning the unpredictability associated to
either basins of attraction or exit basins. For example, Wada basins were
typically assumed to be more unpredictable than fractal basins without the
Wada property [22]. Certainly, such kind of statements could have been
received a more rigorous consideration. Precisely, these considerations lead
to a central question: how can we quantify the intuitive notion of the
unpredictability associated to the basins of attraction? The basin entropy
provides some answers.

Here, we will sketch the definition of the basin entropy and some of the main
results explained in [20]. The basins of attraction can be considered as a
coloring of the phase space where each initial condition is matched to a color
corresponding to an attractor. For the definition of the basin entropy, we use
a covering of size $\varepsilon$ of the basins. The proportion of colors in
each $\varepsilon$-ball defines the probability associated to each color
$p_{i}$. Then, the Gibbs entropy of each ball can be computed as
$S_{i}=-\sum_{i=1}^{N_{A}}p_{i}\log{p_{i}}$, where $N_{A}$ is the number of
different colors in the basin. The basin entropy is defined as the average of
the Gibbs entropy for all balls $S_{b}=\langle S_{i}\rangle$.

Although computing the basin entropy requires no extra assumptions but just
applying the instructions of the previous paragraph, we can delve deeper into
its meaning with the following analysis. Using the definition of the box-
counting dimension, the number of $\varepsilon$-balls covering a phase space
of dimension $D$ grows as $N=n\varepsilon^{D}$. This number, for the covering
of a basin boundary of dimension $D_{k}$, grows as
$N_{k}=n_{k}\varepsilon^{D_{k}}$. Also, we can assume that the basins are
equally distributed around their boundaries with identical probabilities
$p_{k}$. The entropy of an $\varepsilon$-ball in such a boundary is
$S_{i}=-\log~{}m_{k}$, where $m_{k}\geq 2$ refers to the number of basins
separated by the boundary. The sum of all $S_{i}$ for the boundary $k$ is
equal to $N_{k}\log~{}m_{k}$. Bringing everything back together, we obtain the
basin entropy

$\displaystyle
S_{b}=\sum\limits_{k=1}^{k_{max}}\dfrac{n_{k}}{n}\varepsilon^{\alpha_{k}}\log(m_{k}),$
(1)

where $\alpha_{k}=D-D_{k}$ is the uncertainty exponent of the $k$ boundary. As
explained in [20], this expression reveals the main ingredients that
contribute to the basin entropy and therefore to the uncertainty associated to
the long-term prediction of the associated system. In the following, we show
the role of each of these ingredients in our attempt to classify the most
common types of basins of attraction found in nonlinear dynamics.

##  III Uncertainty exponent: smoothness of the boundaries

The first obvious classification of basin boundaries is between smooth and
fractal. This distinction can be analyzed in terms of the uncertainty on a
$\varepsilon$ neighborhood of an initial condition. If all the initial
conditions in the neighborhood lead to the same attractor, this initial
condition is called certain. If at least two initial conditions lead to
distinct attractors, this is an uncertain initial condition. For an average
over the phase space, the ratio $f$ of uncertain to the certain initial
conditions is a function of $\varepsilon$. The usual interpretation of
$f(\varepsilon)$ is a measure of the ignorance on the final state of initial
conditions in the boundary. Smooth basins are well behaved in the sense that,
if the initial uncertainty $\varepsilon$ is decreased, $f$ decreases
proportionally. This is probably the expected behavior, but this is not what
happens with fractal boundaries. In this case, the ratio $f$ decreases as
$f\sim\varepsilon^{D-d}=\varepsilon^{\alpha}$, being $D$ the topological
dimension of the phase space and $d$ the box-counting dimension of the basin
boundary. As a consequence, fractal boundaries have $0\leq\alpha<1$, so one
might decrease the uncertainty of the initial conditions $\varepsilon$ by a
half (e.g. doubling the resolution of the grid of points), and still get a
ratio of uncertain initial conditions $f$ very similar to the original one.
Indeed, in the extreme cases when $\alpha=0$, $f$ remains constant for all the
values of $\varepsilon$. This is exactly what happens for the case of riddled
basins [6].

Therefore, the uncertainty exponent $\alpha$ measures the rate of growth of
the space occupied by the boundary at different magnification scales. This
parameter $\alpha$ does not only allows us to discern between smooth and
fractal boundaries, but also conveys information about how the boundary fills
the phase space. The closer $\alpha$ is to zero, the closer the boundary is to
fill completely the phase space. However, as we shall see later, the
uncertainty exponent fails to grasp other relevant aspects of the basin
boundaries.

The uncertainty exponent of each boundary appears explicitly in Eq. 1, so the
basin entropy also captures the scaling of the boundaries. However, it can
also provide an efficient test of fractality when the basin is computed with
at a single scale. The solution comes as a test [23] able to distinguish
between smooth boundaries and fractal basins. This test, based on the basin
entropy, cannot give information about the rate of growth of the boundary, but
can reveal the existence of structures visible at a given resolution. The
principle of detection is to compute the basin entropy only at the boundary as
accurately as possible with a ball of radius $\varepsilon$. If the value lies
within a theoretical interval, the boundary is smooth. Otherwise, it means
that it has structures below the scale $\varepsilon$.

##  IV Lacunarity: connectedness of the basins

The factor $n_{k}/n$ of Eq. 1 can be identified with the size of the
boundaries, which is related to the concept of lacunarity [24], although it is
not quite the same. To illustrate the meaning of the lacunarity and its role
in the basin entropy we can think of a simple example. Assume we have two
different systems, both with two attractors $N_{A}=2$ and both with smooth
basin boundaries, so that $D_{k}=1$. The basin entropy can still be different
if, for example, the basin of an attractor is just a black disk or if it
consists of $N$ disks with the same total area. As long as we use disks or
other smooth shapes, we would have the same uncertainty exponent $\alpha=1$.
Furthermore, if we carefully keep the total area of all the disks unaltered,
the basin stability [11] would remain constant as well. Having just one disk
with area $\pi R^{2}$ means having a boundary of length $l=2\pi R$, while
having $N$ disks with total area $\pi R^{2}$ implies a perimeter $l_{i}=2\pi
R/\sqrt{N}$ for each disk. The total length of the boundary in this latter
case can grow indefinitely as we keep increasing the number of circles
$l_{T}=\sum_{i=1}^{N}l_{i}=2\pi R\sqrt{N}$. This difference in the uncertainty
has its origin in their different lacunarity.

Figure 1: Basins of attraction of the map $z_{n+1}=z_{n}^{2}+c$ for parameters
$c=-0.4849+0.4472i$ (left) and $-0.4548+0.7688i$ (right). The two figures have
fractal boundaries with the same uncertainty exponent $\alpha=0.8$.

For fractal boundaries, the lacunarity can be very different from one basin to
another. To investigate its role more in detail we will focus on basins with
the same uncertainty exponent. For example, the Fatou sets in Fig. 1(a) and
Fig. 1(b) have this property, but one is a connected set and the other one is
disconnected. In Fig. 2(c), we use a colormap to show the uncertainty exponent
of Fatou sets with different values of the complex parameter $c$, which
produces a picture that resembles the Mandelbrot set, as expected. We can
observe that the same colors are inside and outside the Mandelbrot set,
meaning that the same values of the uncertainty exponent $\alpha$ can be found
for connected and disconnected sets, alike. Thus, the uncertainty exponent
provides no information about the connectivity of the basins, as stated for
example in [25].

Figure 2:  (a) Uncertainty dimension calculation for both Fatou sets in Fig.
1. Both of them have the same uncertainty exponent (parallel curves) but their
lacunarity is different. (b) Heat map of the uncertainty exponent for
different values of the complex parameter $c$ for the Fatou sets. Although the
uncertainty exponent is typically between 0 and 1, here it can take larger
values because we are dealing with Cantor-dust structures, as explained in the
text. (c) Same as (b), but here the color represents the basin entropy. While
for (b) many points inside and outside the Mandelbrot set have the same color
(same uncertainty exponent for connected and disconnected sets, alike), the
basin entropy is clearly different for points that belong to the Mandelbrot
set. In panels (b) and (c), the values corresponding to the Fatou sets of Fig.
1-(a) and (b) are represented as a circle $(c=-0.4849+0.4472i)$ and a cross
($c=-0.4548+0.7688i$).

The Fatou sets depicted in Fig. 1 have the same uncertainty exponent but
different lacunarity. The sparse structure of the disconnected set indicates a
higher lacunarity than the connected set. Figure 2(a) represents $f$ versus
$\varepsilon$ for the two basins in a log-log scale. Since the slope of the
linear fit estimates the uncertainty exponent, the parallel plots reveals the
same fractal dimension for the two examples. The lacunarity is sometimes
interpreted as the intercept of the linear fit on the vertical axis. This
explains why the plot corresponding to the disconnected Fatou with higher
lacunarity (in red) is above the plot of the connected set (in blue). The
basin entropy is larger for the connected basin than for the disconnected
basin for all scales.

Although the lacunarity of the disconnected set is higher, the probability of
landing in one of its islands is very small. Indeed, from a strictly
mathematical perspective, the basin entropy of all the disconnected Fatou sets
should be equal to zero, because they are zero Lebesgue measure sets (two-
dimensional Cantor dusts) and the probability of randomly landing exactly on a
point belonging to them is exactly zero. However, if we represent these
disconnected Fatou sets with finite resolution, such as in Fig. 1(b), their
basin entropy is different from zero and it conveys some information about
their structure. For this same reason, the values of the basin entropy for the
connected Fatou sets are much larger than for the disconnected ones, as can be
clearly seen in Fig. 2(c).

In the previous example, both basins had the same uncertainty exponent
although there was a topological distinction between connected and
disconnected boundaries. In Fig. 3, the two basins from different dynamical
systems have an uncertainty exponent $\alpha=0$ 111In this article, $\alpha=0$
is considered a sufficient condition to decide if a basin is riddled, see [6,
7, 26, 4] for further discussions on this topic.. For both boundaries, the
disconnected set fills all the phase space. Nonetheless, there are striking
differences in the texture of both basins.

Figure 3:  In (a), we reproduce the basins published in [26] that have the
riddled property. In (b) we display the basins of the forced damped pendulum
$\ddot{x}+0.2\dot{x}+\sin(x)=1.3636\sin 0.5t$, which is also riddled. While
the two basins have the same uncertainty exponent, their fine structure
clearly differs as can be observed from the figures.

The basin entropy is an average of the Gibbs entropy for boxes all over the
phase space. The restriction of this average to boxes on the boundary is
called the boundary basin entropy $S_{bb}$. For riddled basins, the boundary
fills the phase space and both averages are the same: $S_{b}=S_{bb}$. Still,
the distribution of basins inside each box may vary depending on the location
of the box. Here we have to interpret the lacunarity in a wider sense that
also includes the distribution of the basins inside each box. Figure 3(a) was
first described in Ref. [26] where the authors have proved the riddled nature
of the basins. The computation of the basins of the forced damped pendulum in
Fig. 3(b) reveals riddled basins with an aspect very similar to white noise.
The basin entropy obtained for Fig. 3(a) is $S_{b}\simeq 0.4$ and
$S_{b}\simeq\log 2$ for Fig. 3(b). This difference in the basin entropy
between the two figures reflects the fine structure of the basins that cannot
be grasped with the fractal dimension. The uncertainty is lower in one case
due to the asymmetry in the distribution of the basins of each attractor.
Also, the uncertainty is maximum in the case of the forced pendulum since the
value of the basin entropy is maximum (for two attractors).

##  V Number of attractors and Wada property

The previous sections account for situations with two attractors, but in
multistable systems a large number of attractors with their corresponding
basins can coexist [27]. Furthermore, fractal geometry allows for three or
more basins to share a common boundary. Such situation is often referred to as
Wada basins, and many works have been devoted to identify, prove and
characterize the Wada property [12, 28, 29, 30]. Its interest resides
precisely in the particular uncertainty of these basins: any deviation of an
initial condition lying in a Wada boundary can evolve into any of the three or
more different fates of the system. This unique feature provokes an increase
of the basin entropy since clearly the uncertainty grows when the number of
different attractors increases.

The merging method [28], one of the algorithms used to detect Wada basins,
provides some insights into the relation of the number of attractors and the
basin entropy. In Fig. 4(a), the points on the top and its corresponding
fitting line (both in red), show the basin entropy for the Newton fractal with
ten roots, i.e., ten different basins. Following the procedure described in
[28], we have merged the basins one by one, so initially we had ten basins and
at the end we only had two basins. A possible result of this merging process
is shown in the basins of Fig. 4(b). Interestingly, Wada boundaries remain
unaltered by this merging procedure, so the only difference between the merged
versions of the basins is the number of attractors $N_{A}$ (the number of
different colors of the picture). This is clearly seen in Fig. 4(a), where in
the log-log plot of the basin entropy $S_{b}$ versus the box size
$\varepsilon$, the lines for the different merged basins have the same slope
and they only differ by their intercept. Also, going from the bottom line of
Fig. 4(a) (corresponding to two attractors) to the top line (corresponding to
ten attractors), we can see how the values of the basin entropy tend to a
limiting value as the number of attractors increases. Given that the basins of
the Newton fractal have a single boundary $k=1$, we have that Eq. 1 simply
becomes $S_{b}=(n_{b}/n)\varepsilon^{\alpha}\log m$. This result is in
agreement with Eq. 5 of [13], so that the unpredictability of the system
increases more when it changes from $2$ to $3$ attractors than when it changes
from $9$ to $10$.

Figure 4:  (a) Basin entropy calculation for the Newton fractal with $10$
attractors. Following the merging method, the number of attractors $N_{A}$
varies from $2$ (bottom) to $10$ (top). The slope is always the same because
the boundaries remain the same after the merging operation, but the slope
increases with the number of attractors converging geometrically. (b) Newton
fractal for $10$ attractors, where we have merged the basins creating a binary
picture ($N_{A}=2$). (c) $S_{b}$ (in blue) and $S_{bb}$ (in green) calculation
for different values of the energy in the Hénon-Heiles system. As the critical
energy $E=1/6$ is approached (vertical line), both values converge to
$S_{b}=S_{bb}=\log 3$. We have also plotted the Wada index multiplied by $\log
N_{A}/N_{A}$ as black circles, showing that for Wada boundaries Eq. 3 is
fulfilled.

Another important aspect of these calculations is related to symmetry. In this
case, the symmetry of the basins of the Newton fractal implies the invariance
of the basin entropy for any merging operation. However, this is not true in
general. Even though the boundaries remain unaltered by the merging process,
the basins do change. Since the basin entropy measures the unpredictability
associated to the basins, if they are different, the value of the basin
entropy will be different too.

Thus, the number of basins separated by the boundaries plays an important role
in the basin entropy, but so does another important parameter: the Wada index
[13]. This number accounts for the distribution of the colors inside the
boundary boxes, so that it is maximal for equiprobable situations. Briefly,
the Wada index $w$ for a given box of size $s$ in pixels can be computed as:

$w=\left\\{\begin{array}[]{l}0,~{}m<3\\\
-\frac{m}{\log(m)}\sum_{k=1}^{m}p_{k}\log(p_{k}),~{}m\geq
3,\end{array}\right.$ (2)

where m is the number of different attractors inside the box. The Wada index
is clearly connected to the boundary basin entropy of Eq. 1, which is the
basin entropy computed only on the boxes of the boundary. Indeed, for Wada
basins we obtain

$S_{bb}=-\dfrac{logN_{A}}{N_{A}}W,$ (3)

where $W$ is the average of the Wada index $w$ over the boxes covering the
boundary.

Figure 5: Interpretation of the basin entropy as a function of two of the
important ingredients described in Eq. 1. The basin entropy is decomposed
artificially into three orthogonal components. Although these components are
not independent, this picture is useful to visualize the classification
between basins and their associated unpredictability. In (a) we have a heatmap
representing the value of the basin entropy on the plane lacunarity-$\alpha$.
The coloring is completely arbitrary, but it illustrates the correlation
between the basin entropy and these two quantities. The disconnected Fatou set
is not included on this plot since its basin entropy varies depending on the
resolution chosen for the representation of the points of the set. This may
lead to unintuitive situations where basins with higher lacunarity have a
lower basin entropy. In (b) we have extended the plane to a three dimensional
cube to include the number of attractors in the analysis. The examples
described here, such as the Wada basins obtained from the Hénon-Heiles system,
have been placed on their corresponding place.

This simple relation, which only holds for Wada basins, has been numerically
verified for the Hénon-Heiles Hamiltonian. This paradigmatic 2D open
Hamiltonian system is defined by
$H=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})+\frac{1}{2}(x^{2}+y^{2})+x^{2}y-\frac{1}{3}y^{3}$.
For an energy above its critical value, it has three exits and the
corresponding escape basins show the Wada property [22]. Varying the total
energy of the system, we can find hyperbolic or non-hyperbolic regimes and the
uncertainty exponent $\alpha$ changes accordingly, but in all cases the escape
basins present the Wada property [15]. In Fig. 4(c), we depict the computation
of the basin entropy $S_{b}$ (blue line) and the boundary basin entropy
$S_{bb}$ (green line) for different values of the energy $E$ in the Hénon-
Heiles system. First, the Eq. 3 holds numerically, since the dots calculated
as $-\dfrac{logN_{A}}{N_{A}}W$ perfectly match the boundary basin entropy
$S_{bb}$. Second, as the energy approaches its critical value $E_{c}=1/6$
(vertical dashed line), both the basin entropy $S_{b}$ and the boundary basin
entropy $S_{bb}$ tend to their maximum possible value, in this system $\log 3$
(horizontal dashed line). This result is connected to the total fractalization
of the basins taking place at the critical value of the energy [31]. At the
critical value, the boundary occupies all the phase space, so that all the
boxes are in the boundary and therefore $S_{b}=S_{bb}$. Also, given the
$2\pi/3$ rotational symmetry of the system, the three basins are equiprobable.
It maximizes the possible uncertainty for $N_{A}=3$, defined by
$S_{b}=S_{bb}=\log 3$. This example underscores the importance of the Wada
property and the Wada index in the unpredictability of a dynamical system.

##  VI Discussion

Along the previous lines we have explained the contributions of different
ingredients to the basin entropy and, consequently, to the unpredictability
associated to basins of attraction in multistable systems. The interesting
point is that it opens a way to classify basins of attraction. We can identify
some of the most notable types of basins as cases where these ingredients take
extreme values. For example, basins with smooth boundaries have $\alpha=1$,
which is the maximum possible value for the uncertainty exponent and therefore
they are the most predictable in this respect.

On the other side of the spectrum we find riddled basins, which have the
minimum possible value of $\alpha=0$. They are the most unpredictable basins
in the sense of how the ratio of uncertain initial conditions grows with
resolution. Also, riddled basins have large values for the lacunarity, given
their Cantor-like disconnected structure.

Concerning the number of attractors, the Wada property maximizes the number of
basins separated by one boundary. Nevertheless, the distribution of the
different probabilities plays a crucial role too. This can be quantified by
the Wada index, and its value often depends on the symmetry of the system.
Wada boundaries are fractal boundaries so they always have $\alpha<1$, but
there is a case where boundaries separate all possible basins (like Wada
basins) and also has $\alpha=0$ (like riddled basins). This kind of basins is
called intermingled [8, 32, 33] and they maximize the unpredictability in two
different ways at the same time. Although the Wada index may vary,
intermingled basins are the most unpredictable basins according to their basin
entropy.

The basin entropy can be used to detect these extreme cases using the criteria
listed in Table 1 for each case. These necessary conditions on $S_{b}$,
$S_{bb}$ and $W$ can help to identify these special basins. However, the basin
entropy remains a classification tool since it sets a hierarchy among the
basins from the point of view of the unpredictability in the initial
conditions space.

Basin type | Ingredients values | Basin entropy criterion
---|---|---
|  $\alpha$ |  lacunarity |  $N_{A}$ |
Smooth | 1 | - | - |  $S_{bb}=0.439$ [23]
Fractal | $\neq 1$ | - | - |  $S_{bb}\neq 0.439$ [23]
Wada | $<1$ | - | $\geq 3$ | $S_{bb}=-\dfrac{logN_{A}}{N_{A}}W$
Riddled | 0 | max | 2 | $S_{b}=S_{bb}$
Intermingled | 0 | max | $\geq 3$ | $S_{b}=S_{bb}=-\dfrac{logN_{A}}{N_{A}}W$
Table 1:  Summary of the classification of the different basins according to
the basin entropy. The main ingredients present in the basin entropy are
listed with their typical values. The column on the right expresses a
necessary condition on the basin entropy and boundary basin entropy to detect
the type of basin involved. This table summarizes the extreme cases, but the
basin entropy can also compare two basins from the perspective of the
unpredictability in the phase space.

The classification can be visualized graphically in Fig. 5(a) and 5(b) where
the basin entropy has been decomposed artificially on three independent axis.
The examples presented in this work have been represented and located on this
graph. Although the lacunarity and the uncertainty exponent are dependent of
each other, this figure is illustrative of the main contributions of these
ingredients to the basin entropy.

## Conflict of interest

The authors declare that they have no conflict of interest.

## CRediT authorship contribution statement

Alvar Daza: Investigation, Visualization, Software, Formal analysis, Writing -
original draft. Alexandre Wagemakers: Investigation, Visualization, Software,
Formal analysis, Writing - review & editing. Miguel A.F. Sanjuán: Supervision,
Conceptualization, Investigation, Formal analysis, Writing - review & editing,
Funding acquisition.

##  VII Acknowledgments

The authors acknowledge support from the Spanish State Research Agency (AEI)
and the European Regional Development Fund (ERDF, EU) under project
PID2019-105554GB-I00.

## References

  * Pisarchik and Feudel [2014] Pisarchik AN, Feudel U.  Control of multistability.  Phys Rep 2014;540:167–218.
  * Nusse and Yorke [1996] Nusse HE, Yorke JA.  Basins of attraction.  Science 1996;271:1376–1380.
  * McDonald et al. [1985] McDonald SW, Grebogi C, Ott E, Yorke JA.  Fractal basin boundaries.  Physica D 1985;17:125–153.
  * Aguirre et al. [2009] Aguirre J, Viana RL, Sanjuán MAF.  Fractal structures in nonlinear dynamics.  Rev Mod Phys 2009;81:333–386.
  * Kennedy and Yorke [1991] Kennedy J, Yorke JA.  Basins of Wada.  Physica D 1991;51:213–225.
  * Alexander et al. [1992] Alexander J, Yorke JA, You Z, Kan I.  Riddled basins.  Int J Bifurcat Chaos 1992;02:795–813.
  * Ott et al. [1993] Ott E, Sommerer JC, Alexander JC, Kan I, Yorke JA.  Scaling behavior of chaotic systems with riddled basins.  Phys Rev Lett 1993;71:4134–4137.
  * Lai and Grebogi [1995] Lai YC, Grebogi C.  Intermingled basins and two-state on-off intermittency.  Phys Rev E 1995;52:R3313.
  * Grebogi et al. [1983] Grebogi C, McDonald SW, Ott E, Yorke JA.  Final state sensitivity: An obstruction to predictability.  Phys Lett A 1983;99:415–418.
  * Mandelbrot [1995] Mandelbrot BB.  Measures of fractal lacunarity: Minkowski content and alternatives.  In: Fractal geometry and stochastics. Springer; 1995, p. 15–42.
  * Menck et al. [2013] Menck PJ, Heitzig J, Marwan N, Kurths J.  How basin stability complements the linear-stability paradigm.  Nat Phys 2013;9:89–92.
  * Daza et al. [2015] Daza A, Wagemakers A, Sanjuán MAF, Yorke JA.  Testing for Basins of Wada.  Sci Rep 2015;5:16579.
  * Saunoriene et al. [2021] Saunoriene L, Ragulskis M, Cao J, Sanjuán MAF.  Wada index based on the weighted and truncated shannon entropy.  Nonlinear Dynamics 2021;104:739–751.
  * Daza et al. [2016] Daza A, Wagemakers A, Georgeot B, Guéry-Odelin D, Sanjuán MAF.  Basin entropy: a new tool to analyze uncertainty in dynamical systems.  Sci Rep 2016;6:1–10.
  * Nieto et al. [2020] Nieto AR, Zotos EE, Seoane JM, Sanjuán MAF.  Measuring the transition between nonhyperbolic and hyperbolic regimes in open hamiltonian systems.  Nonlinear Dynamics 2020;99:3029–3039.
  * Halekotte and Feudel [2020] Halekotte L, Feudel U.  Minimal fatal shocks in multistable complex networks.  Sci Rep 2020;10:1–13.
  * Mugnaine et al. [2019] Mugnaine M, Andrade F, Szezech Jr J, Bazeia D.  Basin entropy behavior in a cyclic model of the rock-paper-scissors type.  EPL 2019;125:58003.
  * Daza et al. [2017] Daza A, Georgeot B, Guéry-Odelin D, Wagemakers A, Sanjuán MAF.  Chaotic dynamics and fractal structures in experiments with cold atoms.  Phys Rev A 2017;95:013629.
  * Zotos and Suraj [2018] Zotos EE, Suraj MS.  Basins of attraction of equilibrium points in the planar circular restricted five-body problem.  Astrophys Space Sci 2018;363:20.
  * Daza et al. [2018a] Daza A, Wagemakers A, Georgeot B, Guéry-Odelin D, Sanjuán MAF.  Basin entropy, a measure of final state unpredictability and its application to the chaotic scattering of cold atoms.  In: Chaotic, Fractional, and Complex Dynamics: New Insights and Perspectives. Springer; 2018a, p. 9–34.
  * Gusso and de Mello [2021] Gusso A, de Mello LE.  Fractal dimension of basin boundaries calculated using the basin entropy.  Chaos Solitons Fractals 2021;153:111532.
  * Aguirre et al. [2001] Aguirre J, Vallejo JC, Sanjuán MAF.  Wada basins and chaotic invariant sets in the Hénon-Heiles system.  Phys Rev E 2001;64:066208.
  * Puy et al. [2021] Puy A, Daza A, Wagemakers A, Sanjuán MAF.  A test for fractal boundaries based on the basin entropy.  Commun Nonlinear Sci Numer Simul 2021;95:105588.
  * Mandelbrot [1982] Mandelbrot B.  The fractal geometry of nature.  WH Freeman New York; 1982.
  * Falconer [2004] Falconer K.  Fractal geometry: mathematical foundations and applications.  John Wiley & Sons; 2004.
  * Ott et al. [1994] Ott E, Alexander J, Kan I, Sommerer JC, Yorke JA.  The transition to chaotic attractors with riddled basins.  Physica D 1994;76:384–410.
  * Feudel et al. [1996] Feudel U, Grebogi C, Hunt BR, Yorke JA.  Map with more than 100 coexisting low-period periodic attractors.  Physical Review E 1996;54(1):71.
  * Daza et al. [2018b] Daza A, Wagemakers A, Sanjuán MAF.  Ascertaining when a basin is wada: the merging method.  Sci Rep 2018b;8:1–8.
  * Wagemakers et al. [2021] Wagemakers A, Daza A, Sanjuán MAF.  How to detect wada basins.  Discrete Continuous Dyn Syst Ser B 2021;26:717–739.
  * Wagemakers et al. [2020] Wagemakers A, Daza A, Sanjuán MAF.  The saddle-straddle method to test for wada basins.  Commun Nonlinear Sci Num Sim 2020;84:105167.
  * Aguirre and Sanjuán [2003] Aguirre J, Sanjuán MAF.  Limit of small exits in open hamiltonian systems.  Phys Rev E 2003;67:056201.
  * Sommerer and Ott [1996] Sommerer JC, Ott E.  Intermingled basins of attraction: uncomputability in a simple physical system.  Phys Lett A 1996;214:243–251.
  * Alexander et al. [1996] Alexander JC, Hunt BR, Kan I, Yorke JA.  Intermingled basins for the triangle map.  Ergod Theory and Dyn Syst 1996;16:651–662.


date22022022

# A micro-CT of a human skull

Thomas Kirchner
Institut für Physik
Martin-Luther-Universität Halle-Wittenberg, Halle (Saale), Germany
<EMAIL_ADDRESS>
(; date)

###### Abstract

This is a data descriptor for an X-ray microtomography of a human skull. The
data is available at doi:10.5281/zenodo.6108435, open data and adheres to the
FAIR [1] principles.

_Keywords_ micro-CT $\cdot$ X-ray microtomography $\cdot$ _ex vivo_ $\cdot$
human $\cdot$ Homo sapiens $\cdot$ skull

##  1 Data set

The data set includes the X-ray projection images (the measured raw data), a
3D reconstruction in approximate Hounsfield units (HU) and a semi-manual bone
segmentation.

Figure 1: Raw data example – one of the 2501 2D cone beam X-ray projection
images.

### Raw data

‘halle_skull_2D_projections.zip‘ contains 2501 projection TIFF images,
recorded by a cone beam micro-CT as described in section 3.

### Reconstructed 3D volume

‘halle_skull.nrrd‘ is a nearly raw raster data (NRRD) [2] file containing the
reconstructed 3D micro-CT. The file was produced by

  1. 1.

reconstructing a volume with the proprietary software CT Pro (Metris, Tring,
UK) using the settings defined in the file
‘etc/halle_skull_3D_reconstruction.xtekVolume‘, then

  2. 2.

converting the CT Pro output to NRRD using ImageJ [3], and then

  3. 3.

scaling the reconstruction to approximate HU.

HU in a voxel are defined by
$\mathrm{HU}=1000\times(\mu_{\mathrm{voxel}}-\mu_{\mathrm{water}})/(\mu_{\mathrm{water}}-\mu_{\mathrm{air}})$.
The absorption values $\mu_{\mathrm{air}}$ for air and $\mu_{\mathrm{water}}$
for water were manually segmented from three embedded water samples and
surrounding air. Note that these HU are only approximated – as this is an ad
hoc calibrated cone beam CT, the resulting HU values should not be interpreted
as quantitative [4].

The finished 3D volume has a isotropic spacing of $0.125\text{\,}\mathrm{mm}$
with 1150 $\times$ 1700 $\times$ 1400 16-bit signed integer voxels.

### Bone segmentation

‘halle_skull_segmented.nrrd‘ is a semi-manual segmentation of the skull bone
based on ‘etc/threshold.py‘ which combines a fixed threshold and a Sauvola
binarization [5]. Incorrectly segmented voxels, such as reconstruction
artifacts, were cleaned up from the segmentation by hand using MITK [6].

##  2 Subject

The subject is an _ex vivo_ skull of a body donor, a 70 year old male. He
signed a declaration of consent at lifetime for the use of his body for
scientific purposes. The skull was provided by the Institute of Anatomy and
Cell Biology, Martin-Luther-University Halle-Wittenberg, Halle (Saale),
Germany. Other than the calcifications around a small trepanation in the right
parietal bone, this is a typical human skull. There is one coronal cut, which
was necessary during preparation of the skull. The skull was imaged without
jaw bone.

##  3 Measurement

X-ray tomography was performed with an industrial cone beam micro-CT (XT H
225, Nikon Metrology) with $150\text{\,}\mathrm{kV}$,
$400\text{\,}\mathrm{\SIUnitSymbolMicro A}$ and using a
$1\text{\,}\mathrm{mm}$ copper filter. The detector had 1750 $\times$ 2000
pixels with $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ isotropic spacing.
This and other metadata of the acquisition can be found in the files
‘etc/halle_skull.xtekct‘ and ‘etc/halle_skull.ctprofile.xml‘. 2501 angles were
measured, each by averaging 8 exposures with $708\text{\,}\mathrm{ms}$
exposure times. This covered a $360.349\text{\,}\mathrm{\SIUnitSymbolDegree}$
rotation – exact angles are listed in the file
‘halle_skull_2D_projections/_ctdata.txt‘.

The skull was placed on a low-density polyethylene packing foam together with
three water filled $1.5\text{\,}\mathrm{mL}$ polypropylene micro-centrifuge
tubes. The cut-off top of the skull was held in position by
$48\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick polypropylene tape.

## Acknowledgments

Thanks to Prof. Heike Kielstein of the Institute of Anatomy and Cell Biology
of the Martin-Luther-University Halle-Wittenberg, Halle (Saale), Germany, for
providing the sample. Thanks to Dr. John Maximilian Köhne of the Department of
Soil Physics, Helmholtz-Zentrum für Umweltforschung UFZ, Halle (Saale),
Germany, for his essential contributions in the micro-CT data acquisition.
Thanks to Prof. Hans-Jörg Vogel of the Department of Soil Physics, Helmholtz-
Zentrum für Umweltforschung UFZ, Halle (Saale), Germany, for providing access
to the micro-CT of his department. Thanks to Prof. Jan Laufer of the
Fachgruppe Medizinische Physik, Institut für Physik, Martin-Luther-Universität
Halle-Wittenberg, Halle (Saale), Germany, for hosting this research.

This work was funded by the German Research Foundation (DFG, Deutsche
Forschungsgemeinschaft) under grant number 471755457.

## References

  * [1] M. D. Wilkinson et al.  _The FAIR Guiding Principles for scientific data management and stewardship_. Scientific data, 2016, doi:10.1038/sdata.2016.18
  * [2] Nearly raw raster data (NRRD) filetype,  http://teem.sourceforge.net/nrrd/
  * [3] ImageJ image processing software,  https://imagej.net/software/imagej/
  * [4] P. Mah, T. E. Reeves, and W. D. McDavid.  _Deriving Hounsfield units using grey levels in cone beam computed tomography_. Dentomaxillofacial Radiology, 2010. doi:10.1259/dmfr/19603304.
  * [5] J. Sauvola and M. Pietikainen.  _Adaptive document image binarization_. Pattern Recognition, 2000. doi:10.1016/S0031-3203(99)00055-2
  * [6] The Medical Imaging Interaction Toolkit (MITK),  https://mitk.org