# An Indoor Localization Dataset and Data Collection Framework with High
Precision Position Annotation

F. Serhan Daniş  A. Teoman Naskali  A. Taylan Cemgil  Cem Ersoy  Dept. of
Computer Engineering, Galatasaray University, Istanbul, Turkey  Dept. of
Computer Engineering, Boğaziçi University, Istanbul, Turkey

###### Abstract

We introduce a novel technique and an associated high resolution dataset that
aims to precisely evaluate wireless signal based indoor positioning
algorithms. The technique implements an augmented reality (AR) based
positioning system that is used to annotate the wireless signal parameter data
samples with high precision position data. We track the position of a
practical and low cost navigable setup of cameras and a Bluetooth Low Energy
(BLE) beacon in an area decorated with AR markers. We maximize the performance
of the AR-based localization by using a redundant number of markers. Video
streams captured by the cameras are subjected to a series of marker
recognition, subset selection and filtering operations to yield highly precise
pose estimations. Our results show that we can reduce the positional error of
the AR localization system to a rate under 0.05 meters. The position data are
then used to annotate the BLE data that are captured simultaneously by the
sensors stationed in the environment, hence, constructing a wireless signal
data set with the ground truth, which allows a wireless signal based
localization system to be evaluated accurately.

###### keywords:

Data collection , indoor positioning , augmented reality , data annotation ,
localization evaluation

††journal: Pervasive and Mobile Computing

##  1 Introduction

The majority of the world population has started to dwell in cities and
urbanized areas, where people spend most of their time indoors. As the
buildings and structures become ever more complex, the need to use navigation
tools gains importance. However, outdoor navigation technologies like GPS do
not work indoors. Therefore, there is a need to use other means of
localization indoors.

Wireless technologies are preferred for indoor localization as they are
globally available and affordable. Whereas the localization algorithms have
been improved with the ubiquitous Wi-Fi technology for indoors [1], due to its
low energy consumption, price and size, Bluetooth Low Energy (BLE) is
particularly suited to this domain [2]. BLE beacons or sensors are utilized to
track objects or people, and most hand held mobile devices can easily
incorporate this technology. However, the nature of the wireless signals
render the localization task rather difficult as the signals are easily
reflected or absorbed by obstacles. To calibrate and optimize wireless based
localization systems, ground truth data must be collected. More importantly,
to measure the real performance of these systems, we require trajectories
labeled with ground truth positions alongside the wireless signal data, rather
than separately measured reference points. This ground truth collection
process can be rather tedious and usually labor intensive [3].

In the literature, there exist several valuable datasets formed by wireless
signal parameters and corresponding ground truth information, however, such
data collections do not reflect the true nature of the trajectories. Instead,
the sets are formed by signal parameters collected with sparsely located
reference points [4], or by modifying the paths artificially making the
carrier remaining still on predefined reference points [5]. Moreover, the
ground truth notion may not even belong to the positions. In the work of
Girolami et al. [6], the ground truth of social interactions are reported with
the BLE received signal strength indicator (RSSI) data. Likewise, in the work
of Sansano-Sansano et al. [7], the ground truth is in terms of the user gate
speed.

High resolution ground truth for indoor positioning is rare because of the
difficulties in collecting accurate position data on every point the object of
interest passes. This difficulty prevents the researchers from setting a
generalized benchmark for indoor positioning systems [8]. The positioning
community lacks a clear standard for evaluating positioning systems. Multiple
attempts have been made in order to fill this gap to validate the accuracy of
an indoor positioning system. De la Osa et al. propose a method to define a
ground truth for real life purposes [9]. The technique bases on creating
predefined paths with checkpoints which are used as indicators when the tester
passes over them. The ground truth trajectory is then estimated by
interpolated points between checkpoints. Another one is the skeleton path
technique in which the user follows a predetermined path with a constant speed
and the positions are sampled. Clearly the techniques based on walking are
prone to human errors, incorrectly assuming constant paces and true timings
[10]. Moreover, the positions are estimated using inaccurate dead-reckoning
data based on step counting [11]. According to Adler et al. [12], this
predefined path technique, or namely office walk, is the most preferred ground
truth collection method. This inaccurate method is also the cheapest and the
most precise ground truth data collection method so far used to validate
indoor positioning systems.

Our ground truth data collection technique provides a set of timestamped
wireless signal parameters. Whereas the technique can be easily modified to
collect and annotate other wireless signal based parameters, we work with the
prototypical BLE RSSI data. The RSSI data are annotated both in temporal and
spatial domains accurately with the corresponding ground truth positions using
a vision-based localization system. We employ the markers designed for AR
applications, to annotate the RSSI data with ground truth positions.

AR is another technology that performs exceptionally well in indoor locations
when based on special markers. Whereas AR applications are initially developed
for adding artificial figures into real world views, a camera with proper
configurations can detect the precise positions and orientations of its
special markers relative to the camera. This information is inversely used to
calculate the location of the camera if the marker is stationary, so they are
preferred in indoor navigation when we aim to track the viewer [13]. Although
they function well, camera based devices are usually more expensive than BLE
devices and they require the active movement of the camera to determine the
location, which is not ideal if objects or people need to be passively
tracked. Accordingly, La Delfa et al. compare different marker technologies
for environmental suitability [14]. Xavier et al. utilize markers for SLAM
[15], whereas Ujkani et al. perform sensor calibration with AR markers in a
robotic cell achieving accuracies of 5-10 cm [16].

In the wireless positioning community, AR positioning techniques are used
hybridly to enhance the performance of the former systems as wireless signals
are prone to yield high error rates [17, 18]. Amongst these studies, Byrne et
al. propose a methodology for collecting RSSI information and accelerometer
data using markers placed on the ground and cameras mounted on people to
collect the ground truth [19].

We propose a framework that makes use of an AR-based localization system to
collect precisely annotated RSSI data in both time and space. The system
estimates positions of a specially designed camera-beacon bundle that captures
video and RSSI data simultaneously, and is shown to implement a very accurate
positioning with an error rate under 0.05 m and it can function alongside the
BLE based localization without causing extra interference. With its low error
rate, this method provides ground truth positions to evaluate or verify a pure
wireless signal based localization system. After the postprocessing, we
possess RSSI data points, each of which is annotated with its precise 3D
position. Moreover, the proposed technique has a potential to greatly speed up
the calibration procedure and decrease or fully automatize the workload of
fingerprinting and radio frequency map generation.

As an end product, we also provide a publicly available dataset collected
using the proposed method. The dataset includes position annotated RSSI data
of several trajectories alongside the preprocessed probabilistic radio
frequency map estimations in fingerprint or grid formats 111The dataset is
publicly available with DOI: 10.34740/KAGGLE/DS/1662453, and the framework is
public in Github as “Position-Annotated-BLE-RSSI-Dataset”..

In Section 2, the methods for data selection, filtering and synchronization
that make up the position annotated wireless data collection framework are
given. We describe our experimental setup and dataset specifications in
Section 3. The results supporting our proposed framework are reported in
Section 4 and we conclude with future insights in Section 5.

##  2 Data Annotation Methodology

Data annotation methodology comprises of four main consecutive sections as
seen in Figure 1. With the camera calibration matrices, premeasured poses
(positions and orientations) of the AR markers, timestamped RSSI readings and
accompanying video frames in hand, we want to annotate each RSSI reading on a
trajectory with a precise position estimated using the videos and marker
information.

Figure 1: Position annotation on RSSI data: we first (i) estimate the camera
poses in the video frames using the marker poses and the camera calibration
matrices, and then (ii) filter the likely poses by eliminating highly
irrelevant estimations. The likely poses are (iii) smoothed using a properly
designed Kalman filter to find a final pose for each frame, which is finally
(iv) synchronized with the simultaneous RSSI data in time to obtain a position
annotated RSSI dataset for each experiment in the area. The arrow weights
denote the data load coarsely.

Before going through the data annotation procedure, we first present the
mathematical notations, definitions and necessary conversions in 3D space in
Section 2.1. In Section 2.2 we employ the ArUco marker detector to estimate
marker poses relative to the camera(s) [20, 21]. With the defined mathematical
foundations, estimated marker poses relative to the cameras are combined with
the marker poses relative to the world to obtain the pose of the camera-beacon
bundle in the world frame. In Section 2.3, we apply pruning techniques to
eliminate the misleading estimations in a batch of pose estimations, which we
observe to be very noisy due to erroneous marker detections. We apply a Kalman
filter specially tailored for this problem to obtain the most relevant pose
given a batch of poses in Section 2.4. Finally, in Section 2.5 these smoothed
precise camera poses are synchronized in time with the RSSI readings.

###  2.1 Notations, definitions and conversions

Throughout this article, vectors are represented by bold lowercase letters
(e.g. $\bm{\tau}$ or $\bm{\theta}$), while bold capital letters represent
matrices (e.g. $\mathbf{R}$). The variables, $m$, $c$ and $w$ denote the
entities, a marker, the camera-beacon bundle and the world, respectively. We
represent the position (translation) and the rotation (orientation) vectors
with ${}^{i}\bm{\tau}_{jk}$ and ${}^{i}\bm{\theta}_{j}$ respectively. The
superscript on the left, $i$, is for the frame the vector is defined in, and
the subscript on the right, $j$, is for the entity the vector belongs to. For
the position vectors, a double subscript $jk$ shows the start and the end
points. If the start position is the frame origin itself, we omit the first
entity, ${}^{j}\bm{\tau}_{jk}\equiv{}^{j}\bm{\tau}_{k}$. For example,
${}^{c}\bm{\theta}_{m}$ defines the rotation vector of the marker $m$, in the
camera frame, $c$, and ${}^{w}\bm{\tau}_{cm}$, specifies a vector from the
camera $c$ to the marker $m$ in the world frame $w$.

We use the rotation matrices to convert vectors from a predefined coordinate
frame to another. A position vector defined in the frame $i$ is multiplied
with the rotation matrix, ${}^{k}_{i}\mathbf{R}$, from the left to be
converted into the coordinate frame $k$. Note that with a proper
multiplication, the coordinate frames in the diagonals are eliminated, ($i$ in
(1)).

$\displaystyle{}^{k}\bm{\tau}_{j}={}^{k}_{i}\mathbf{R}~{}{}^{i}\bm{\tau}_{j}$
(1)

Similarly, if we have a rotation matrix from the frame $i$ to $j$ and another
one from $j$ to $k$, we can obtain transitively the rotation matrix from $i$
to $k$ by multiplying them as in (2). Note the elimination of the frame $j$.

$\displaystyle{}^{k}_{i}\mathbf{R}={}^{k}_{j}\mathbf{R}~{}{}^{j}_{i}\mathbf{R}$
(2)

Since rotation matrices have a rank of three, a rotation vector is a
convenient and compact representation of a rotation matrix. Conversion from a
rotation vector to a rotation matrix is straightforward with the use of
Rodrigues’ formula (3), as follows:

$\displaystyle\begin{split}\mathbf{R}&=Rodrigues(\bm{\theta})\\\
&=\mathbf{I}\cos(\theta)+(1-\cos(\theta))\bm{r}\bm{r}^{T}+\sin(\theta)\left[\begin{matrix}0&-r_{z}&r_{y}\\\
r_{z}&0&-r_{z}\\\ -r_{y}&r_{x}&0\end{matrix}\right]\end{split}$ (3)

where $\theta=\|\bm{\theta}\|$, is the norm of the rotation vector,
$\bm{r}=\bm{\theta}/\theta$, the normalized rotation vector, and $r_{i}$ are
the individual orientation components for the axes in 3D. If we have a
rotation vector of the entity $i$ defined in the coordinate frame $j$, we can
obtain the associated rotation matrix from $i$ to $j$, using the Rodrigues’
rotation formula [22]:

$\displaystyle{}^{j}_{i}\mathbf{R}$
$\displaystyle=Rodrigues({}^{j}\bm{\theta}_{i})$ (4)

The inverse of the Rodrigues’ formula is also available, but with a more
complicated setting:

$\displaystyle\begin{split}\bm{\theta}&=Rodrigues^{-1}(\mathbf{R})\\\
&=\left\\{\begin{matrix}[l]\bm{\theta}=0,&\text{if }s=0\text{ and }d=1\\\
\bm{\theta}=S_{1/2}\left(\frac{\mathbf{v}}{\|\mathbf{v}\|}\pi\right),&\text{if
}s=0\text{ and }d=-1\\\ \bm{\theta}=\frac{\bm{\rho}}{s}\theta,&\text{if
}\sin\theta\neq\ 0\\\ \end{matrix}\right.\end{split}$ (5)

where $\mathbf{A}=\frac{\mathbf{R}-\mathbf{R}^{T}}{2}$,
$\bm{\rho}=[a_{32}~{}a_{13}~{}a_{21}]^{T}$, $s=\|\bm{\rho}\|$,
$d=(\text{tr}(\mathbf{R}-1))/2$, $\theta=\arctan_{2}(s,d)$, $\mathbf{v}$ is a
nonzero column of $\mathbf{R}+\mathbf{I}$ and
$S_{1/2}(\bm{\theta})=-\bm{\theta},\text{if }\|\bm{\theta}\|=\pi\text{ or
}S_{1/2}(\bm{\theta})=\bm{\theta},\text{ otherwise}$. Both of these formulas
can be easily employed through the OpenCV library. If we have a rotation
matrix from the coordinate frame $i$ to the frame $j$, we can obtain the
associated rotation vector of the entity $i$ relative to the coordinate frame
$j$, using the inverse of the Rodrigues’ rotation formula:

$\displaystyle{}^{j}\bm{\theta}_{i}$
$\displaystyle=Rodrigues^{-1}({}^{j}_{i}\mathbf{R})$ (6)

We denote the positions of stationary markers as ${}^{w}\bm{\tau}_{m}$, and
their orientations as ${}^{w}\bm{\theta}_{m}$. When markers are detected and
recognized in a video frame, the pose estimator supplied by the ArUco suite
estimates the position vectors pointing to the center of the marker,
${}^{c}\bm{\tau}_{m}$, and rotation vector, ${}^{c}\bm{\theta}_{m}$, relative
to the camera coordinate frame.

###  2.2 Camera pose estimation

We are interested in extracting the poses of the camera-beacon bundle, which
we will denote (${}^{w}\bm{\tau}_{c}$, ${}^{w}\bm{\theta}_{c}$). The ArUco
suite provides the marker detector routines. With properly calibrated cameras,
we can estimate the poses of the markers relative to the camera coordinate
system, denoted as (${}^{c}\bm{\tau}_{m}$, ${}^{c}\bm{\theta}_{m}$). We relate
the latter information to the former using the premeasured marker poses
relative to the world frame, (${}^{w}\bm{\tau}_{m}$, ${}^{w}\bm{\theta}_{m}$)
(see Section 3.2 for details on markers). For the positions part, we do:

$\displaystyle\begin{split}{}^{w}\bm{\tau}_{c}&={}^{w}\bm{\tau}_{m}-{}^{w}\bm{\tau}_{cm}\\\
&={}^{w}\bm{\tau}_{m}-{}^{w}_{m}\mathbf{R}~{}{}^{m}_{c}\mathbf{R}~{}{}^{c}\bm{\tau}_{cm}\\\
&={}^{w}\bm{\tau}_{m}-{}^{w}_{m}\mathbf{R}~{}\left({}^{c}_{m}\mathbf{R}\right)^{T}~{}{}^{c}\bm{\tau}_{m}\end{split}$
(7)

where ${}^{w}_{m}\mathbf{R}=Rodrigues({}^{w}\bm{\theta}_{m})$ and
${}^{c}_{m}\mathbf{R}=Rodrigues({}^{c}\bm{\theta}_{m})$. The computation
consists of a mere redefinition of the position vector from the camera to the
marker in the world coordinate frame, ${}^{w}\bm{\tau}_{cm}$, and its
subtraction from the marker position in the same coordinate frame,
${}^{w}\bm{\tau}_{m}$ as seen Figure 2.

Figure 2: The position vectors and the coordinate frames: $w$, $c$ and $m$;
for the world, the camera and the marker frames respectively. The axes, $x$,
$y$ and $z$ are coded with red, green and blue arrows respectively. We are to
compute the cyan vector, given the orange vectors.

We can also easily obtain the rotation vector using the rotation matrix from
the world frame to the camera frame, ${}^{w}_{c}\mathbf{R}$. Whereas this
matrix can be used to convert a vector from camera coordinate system to the
world coordinate system, it also defines how the camera is oriented within the
world frame. The rotation matrix was already computed as an intermediate step
of the position computation (7). Employing the inverse of the Rodrigues’
formula, we convert the rotation matrix into the rotation vector of the camera
in the world frame (8) as follows:

$\displaystyle\begin{split}{}^{w}\bm{\theta}_{c}&=Rodrigues^{-1}({}^{w}_{c}\mathbf{R})\\\
&=Rodrigues^{-1}\left({}^{w}_{m}\mathbf{R}~{}\left({}^{c}_{m}\mathbf{R}\right)^{T}\right)\end{split}$
(8)

The conversions and computations above are repeated for each detected marker
in an image frame. Since we construct our detection setup aiming to capture
several markers at a time, we expect to estimate multiple camera pose vectors
at each frame, (${}^{w}\bm{\tau}_{c}$, ${}^{w}\bm{\theta}_{c}$), namely the
raw poses.

###  2.3 Pose elimination

The variances of the camera pose estimations are observed to be large due to
noisy marker pose detections caused by minor errors in manual measurements,
faulty camera matrices and coefficients, motion blur or varying lighting
conditions. The combination of these factors lead to inaccurate rotation
matrices, which in turn accumulate to varying positions for the camera. Before
directly applying a filtering algorithm on the raw poses, we employ two
elimination strategies (see Figure 3).

Figure 3: Position elimination strategies: close marker selection (left) and
outlier elimination (right). Black points are selected for the next stage,
while the red points are eliminated.

At each video frame, we expect to recognize multiple markers. We use two
cameras to expand the field of view, so if one camera misses the markers, the
other one captures possibly more than one marker. Moreover, to make sure that
there exist multiple marker readings, we process $F$ video frames in one shot.
The traversed distance of a batch, $d_{F}$, for the camera-beacon bundle
moving at a constant speed of $v$ can be computed with $d_{F}=vF/(2fps)$,
where $2fps$ is the frame capture speed of a two camera system. For example,
with $F=15$, $fps=60$, and navigating at $v=0.35~{}m/s$, each batch of frames
corresponds to about $d_{F}=0.044$ meters.

We observe that camera poses from distant markers can be more misleading. This
defect is due to the number of pixels of the detected marker image. As the
distance increases, the lack of pixels lead to larger measurement errors.
Therefore, as long as multiple poses are estimated in a time window of $F$
frames, the poses based on distant markers may be discarded. We set a
closeness parameter, $C$, that takes into account the closest $C$ marker
readings. This is a quick technique that sorts the detected markers in a frame
with respect to their distances to the camera and keeps only the least $C$
distant poses.

Another elimination is achieved by an outlier detection test. In a set of pose
estimations, we eliminate the outliers with respect to the distances to a
consensus point [23]. The parameter $U$, is defined to determine the number of
remaining poses after consecutive outlier detection iterations. This technique
is merely a slow backward subset selection method that is repeated in a greedy
manner.

Employing the elimination methods with proper values for the parameters, $C$
and/or $U$, we prune the raw poses into a smaller set of poses, namely the
likely poses.

###  2.4 Kalman filter for smoothing

Even though some marker detections are eliminated, we still receive erroneous
estimations. Error fluctuations can be controlled by employing appropriate
filtering mechanisms. The current section is linked to the previous section
through pose vector estimates of the camera-beacon bundle, which will be our
observations in the filtering design:
$\mathbf{y}_{t}=({}^{w}\bm{\tau}_{c},{}^{w}\bm{\theta}_{c})_{t}$. We model
this new problem as a state space model, in which the observed poses,
$\mathbf{y}_{t}$, are assumed to be generated noisily from the latent poses,
$\mathbf{x}_{t}$. The latent pose vector of the tracked object,
$\mathbf{x}_{t}$ at time $t$, is a function of the pose vector at a previous
time, $\mathbf{x}_{t-1}$. This function is also called the motion model that
stands for the dynamics of the object.

The transition densities are assumed not to follow explicit motion models.
Instead, we assume that $\mathbf{x}_{t}$ resides close to $\mathbf{x}_{t-1}$.
Observation densities are designed similarly. $\mathbf{y}_{t}$ should not be
far away from $\mathbf{x}_{t}$. With a prior distribution at the initial time,
$p(\mathbf{x}_{0})$, we are to find the filtering distribution
$p(\mathbf{x}_{t}|\mathbf{y}_{1:t})$. These assumptions can be realized by
forming transition and measurement densities with Gaussians (9), whose noise
covariance matrices will be denoted as $\mathbf{R}_{t}$ and $\mathbf{Q}_{t}$
respectively.

$\displaystyle\begin{split}\mathbf{x}_{0}&\sim p(\mathbf{x}_{0})\\\
\mathbf{x}_{t}&\sim\mathcal{N}(\mathbf{x}_{t};\mathbf{x}_{t-1},\mathbf{R}_{t})\\\
\mathbf{y}_{t}&\sim\mathcal{N}(\mathbf{y}_{t};\mathbf{x}_{t},\mathbf{Q}_{t})\end{split}$
(9)

The dimensions of the latent space and the measurement space are assumed to be
independent, so we use diagonal matrices. The values of the covariance
matrices are chosen to compensate for unexpected high errors in measurements
[24, 25]. The exact inference for the filtering distribution becomes
tractable, and can be solved with the standard Kalman filter:

Predict: $\displaystyle\hat{\mathbf{x}}_{t}$
$\displaystyle=\mathbf{A}\mathbf{x}_{t-1}$ $\displaystyle\hat{}\mathbf{P}_{t}$
$\displaystyle=\mathbf{A}\mathbf{P}_{t-1}\mathbf{A}^{T}+\mathbf{R}_{t}$
Update: $\displaystyle\mathbf{K}_{t}$
$\displaystyle=\hat{}\mathbf{P}_{t}\mathbf{H}^{T}(\mathbf{H}\hat{}\mathbf{P}_{t}\mathbf{H}^{T}+\mathbf{Q}_{t})^{-1}$
$\displaystyle\mathbf{x}_{t}$
$\displaystyle=\hat{}\mathbf{x}_{t}+\mathbf{K}_{t}(\mathbf{y}_{t}-\mathbf{H}\hat{}\mathbf{x}_{t})$
$\displaystyle\mathbf{P}_{t}$
$\displaystyle=(\mathbf{I}_{d}-\mathbf{K}_{t}\mathbf{H})\hat{}\mathbf{P}_{t}$

where transition and observation matrices, $\mathbf{A}$ and $\mathbf{H}$, are
chosen as identity matrices ($\mathbf{I}_{d}$), $\mathbf{P}_{t}$ are the
estimated covariance matrices, and $\mathbf{K}_{t}$ are the optimal Kalman
gains. We obtain $\mathbf{x}_{t}$ as the pose estimations, namely the final
poses.

###  2.5 RSSI data synchronization

After the filtering operation, timestamped pose data are obtained, but there
still exist problems in synchronizing the timestamped RSSI data to the
timestamped position data.

(a) Detecting covering visually. (b) Detecting covering in RSSI readings.

Figure 4: Synchronization of video frames with RSSI data in time. We apply
artificial disturbances on the videos and the RSSI data by our hands. We
observe three stages: the bundle is uncovered (white); the hands are hold on
the bundle (gray) and the transition stage where the hands are moving to cover
or release the bundle (light gray).

The first problem is that the cameras cannot be synchronized to a precise
universal clock. We overcome this problem by adding artificial disturbances
that affect both the camera views and the RSSI data so that a synchronization
may be possible. After various unsuccessful Faraday cage-like attempts, we
found that a pair of hands would be practical to give disturbances to
electromagnetic signals and visual data. We cover and release the camera-
beacon bundle with two hands at the beginning of each experiment. This hand
movement can be clearly seen in the video frames as shown in Figure 4(a).
Because body parts absorb electromagnetic signals, we observe that this method
creates a significant reduction in the RSSI data. The artificial disturbance
can easily be detected in RSSI readings as seen in Figure 4(b). In the
figures, we divide this hand movement into four color-coded stages. We first
cover the camera-beacon bundle quickly (light gray). When the bundle is
covered with hands, we are in the gray region where the cameras’ views are
largely blocked and we observe that strong RSSI readings coarsely between -85
to -65 decibels are lost because of the absorption due to the hands. Then we
release (light gray) the bundle as quickly as possible to reach the white
region where the camera views are no longer blocked and the RSSI readings
return to their normal course. We select the timestamp of the second light
gray region that points to the change point between the white region and the
gray region. We modify the data streams to start from this selected timestamp,
and we are left with synchronized video frames and RSSI data.

Whereas the data streams are set to begin at the same time, RSSI timestamps
and video timestamps do not match exactly. Secondly and obviously, we have to
unify the two timestamps. RSSI data is far less frequent than the pose data (2
Hz vs 10-30 Hz). To coincide the two data types in time, we perform a linear
interpolation of the video frame timestamps at the RSSI timestamps, and we
label each RSSI timestamp with the linearly interpolated pose data as its
ground truth. After the time synchronization stage, we have a dataset of RSSI
data points annotated with the position it is captured on. Two of the
trajectories from our dataset are displayed in Figure 5.

Figure 5: RSSI data overlaid on the position data: time flows from red bars to
black bars.

##  3 Data Collection Setup

The aim of the data collection setup is to exploit as many sensors as possible
and build a BLE sensor network that covers the whole area at a minimal cost,
and finally build a dataset with the highest information within. We use a set
of BLE beacons and adapters, low profile computers, motion cameras and visual
markers installed in an indoor area of size $20.66\times 17.64$ m2. Intended
originally as an office space for fellow graduate students, the area also
includes columns, walls and furniture that cause multipath interference on the
wireless signals (see Figure 6).

Figure 6: Sample images from the area with markers, extension cables, BLE
adapters, computers and other obstacles.

(a) Illustration of the sensor distribution: yellow rectangles are the
computers, colored circles the sensors and the gray lines are the USB cables
that link the sensors to the computers.

(b) The navigable platform (left); the beacon-camera bundle of a POI beacon
and two GoPro Hero 4 motion cameras (top middle); a Raspberry Pi 3 Model B+
computer (down middle); Logilink Bluetooth adapter (top right); and a KLPRO
KLLMZ60 laser distance meter (down right).

Figure 7: Hardware setup

###  3.1 Hardware setup

We consider that an efficient configuration of the sensor positions in terms
of data quality can be achieved by putting them on the walls surrounding the
test area. However, there appear difficulties in both powering the sensors and
reading data from them with more than 10 meters of the extension cables.
Because of the cable length constraints, in our design, each two sensors are
controlled by one computer. To cover the whole area as equidistantly as
possible, we install 8 Logilink Bluetooth 4.0 adapters on the walls, and use
the 4 onboard Bluetooth embedded chips of the computers, making 12 Bluetooth
sensors (Sensor-$ij$) in total as seen in Figure 7(a). The computers are
placed on the cabinets to ensure that their embedded sensors have a standard
elevation of one meter (see Figure 7(b)). The USB cables are attached to the
ceiling to minimize possible signal interferences. The adapters are attached
on the walls elevated about two meters from the floor. In this configuration
the extension cables are tested to supply enough power to the adapters. The
computers responsible for the data capture system are four Raspberry Pi 3
Model B+ single board computers (RPi-$i$), which are based on Broadcom
BCM2837B0, Cortex-A53 (ARMv8) 64-bit quad-core processors. They run Raspbian
Stretch operating system, on which the Robot Operating System (ROS) Melodic
middleware is installed [26]. The computers scan the BLE packets emitted by
the transmitters and publish them into the local area network (see Figure
7(b)).

We use one sample of POI Beacons as the BLE signal transmitter, which is
constructed upon Nordic Semiconductor multi-protocol SoC and 32-bit Arm Cortex
M0 processor with 256 kB flash and 16 kB on chip. POI Beacons use the
Bluetooth 4.0 BLE protocol and iBeacon protocol, and have a signal range of 70
meters. The beacons send two BLE identity packets per second and have a
running time up to 4 years with two standard AA batteries (see Figure 7(b)).

A 3D printed custom case, which we call “the camera-beacon bundle”, holds the
BLE beacon and two GoPro Hero 4 motion cameras together (see Figure 7(b)). The
cameras are capable of capturing images with a resolution of $1920\times 1080$
at 60 fps. The special design of the custom case lets the cameras stay in such
a formation that they are placed back to back, and their lenses are aligned
facing opposite directions. One camera stays upside down with respect to the
other one. The design also minimizes the translations between camera
coordinate frames and the beacon, provides a large view angle in order to
detect as many markers as possible at an instant, and maximizes the
functionality of the cameras by allowing reaching both of the buttons. The
cameras and the BLE beacon are deliberately put in the same place as we are to
collect BLE data originating from the beacon and simultaneous videos that will
be used to implicitly estimate precise positions of the mobile beacon. By this
setup, we build a data collection infrastructure and then by processing the
videos captured by the cameras, we estimate the positions of the camera-beacon
bundle, so that the RSSI readings by the BLE beacon can be labeled with those
positions. The labels will serve as ground truth for the future pure BLE based
localization purposes.

The bundle is attached onto a tripod which in turn is mounted on a wheeled
cabinet that allows smooth navigation in the test area (see Figure 7(b)). The
beacon height is set to be 2 meters to obtain a clear view of the test area in
terms of both visual data and radio frequency data. The tripod functionality
also allows us to switch to different alignments of the camera-beacon bundle,
from which we prefer horizontal alignment to maximize the view angle, hence,
the number of detectable markers, and to minimize the interference on the
beacon due to the cameras.

The distances are rigorously measured with a KLPRO KLLMZ60 laser distance
meter with a working range of 60 m and an accuracy of $\pm 0.002$ meters (see
Figure 7(b)).

(a) ArUco marker samples (top) and the calibration chart (bottom). (b) The
markers in the test area as black dots and arrows.

Figure 8: ArUco markers and their distribution in the test area.

###  3.2 Markers and camera calibration

We choose 30 square markers from the dictionary of $6\times 6$ matrix markers
of the ArUco marker suite (Figure 8(a)), and print them on hard paper with
edges of 0.3 meters. The markers are attached onto the walls or columns where
they can be viewed from almost any angle in the area (Figure 8(b)). We record
the label and the pose (of the center) of each marker with respect to a
predetermined world coordinate frame (${}^{w}\bm{\tau}_{m}$,
${}^{w}\bm{\theta}_{m}$). The measurement of the position vector,
${}^{w}\bm{\tau}_{m}$, is straightforward, that is, we manually measure the
distances to the walls of reference and to the floor with the help of the
laser distance meter, but the computation of a rotation vector,
${}^{w}\bm{\theta}_{m}$, may be challenging. To ease this process, we install
the markers on the walls or columns in such a way that rotations between the
marker coordinates and the world coordinates can be easily computed by
multiples of $\pi/2$. Specifically the $y$-axis of the markers are always
aligned with the $z$-axis of the world coordinate system. To compute a
rotation vector of a marker in world coordinates, we first compute the
required rotation around the $x$-axis, and then around its $z$-axis. For
example, for the marker $m$ in Figure 2, we first rotate it around its
$x$-axis by $\pi/2$, and then around its $z$-axis by $\pi/2$ again to coincide
the marker coordinate frame with the world coordinate frame. The rotation
vector of the marker $m$ in the world coordinate system is then recorded as
${}^{w}\bm{\theta}_{m}=[\pi/2,0,\pi/2]$.

The ArUco suite provides calibration procedures as well as generating the
marker images and detecting them in still images. Using several views of the
Charuco calibration chart [27] (Figure 8(a)), the camera matrix and the
corresponding distortion coefficients are tuned. A highly accurate camera
matrix is mandatory for precise pose detection of the markers.

The current positioning of two cameras allows us to have an almost spherical
world view with the fish eye mode of the cameras, but we prefer the planar
view mode of cameras to the fish eye mode for three reasons: (i) the
calibration procedure is challenging for the fish eye view, (ii) the fish eye
distortions have to be corrected before running the AR Marker detection
algorithms as they have to be fed with planar markers, and (iii) we do not
need to see ceilings and with the planar view mode, the cameras already have a
large enough angle of view.

###  3.3 Handling setup errors

The errors in marker based localization are due to multiple factors that can
be divided into two sets: human errors and environmental errors. Human errors
include mostly errors in measuring the marker poses, reference points and
unstable navigation. We can also put the camera calibration based errors into
this set, which usually lack the necessary views for the calibration process.
Environmental errors are due mostly to ambient light fluctuations.

The markers are installed with the following factors taken into consideration
in order to minimize the errors due to human factors: (i) From any position
with any orientation in the area, at least one marker should be detectable.
Markers should be above the objects that can cause visual occlusion. (ii)
Markers are attached in the places illuminated uniformly and sufficiently.
(iii) Markers are placed with the help of a jig to ensure that they are at a
fixed height and orientation to facilitate the height and rotation
measurements.

We minimize the errors due to camera calibration by using several views of the
Charuco calibration chart [27] (Figure 8(a)), by which the camera matrix and
the corresponding distortion coefficients are tuned properly. To evaluate the
calibration and detection performance before running the main set of
experiments, we collect a preliminary set of still images. For this, the
navigable platform is put in the area in front of each marker at different
angles and distances, but the office has a raised floor that is prone to
flexing while people are stationary or mobile in the proximity of the image
capture device. The tilting or oscillation of the image capture device can
have a small effect on the captured position information. For calibration and
verification purposes, the authors waited until the capture platform was
stabilised and image capture was triggered wirelessly to obtain the best
results. Sample images can be seen in Figure 6. By adjusting the camera matrix
and the parameters to their best values, we achieve an empirical minimum mean
squared distance error of 0.05 meters for all markers installed in the area.
As the recent studies on BLE based localization techniques report localization
errors no less than 1 meters for practical setups [28, 29], with our
controlled error rate at 0.05 meters, marker based localization is highly
better than any other BLE or WiFi based localization method. Thus this
technique can be used as the ground truth to evaluate the RF RSSI based
localization performance.

The installation of the AR markers is a tedious and time consuming work, but
this procedure can be eased by considering the following points: (i) Instead
of sticking the markers on the walls, fixing them on portable stands at
predetermined elevations will reduce the installation time significantly.
Polycarbonate stands should be preferred to minimize wireless signal
attenuation and interference due to these artificially introduced materials.
The portable marker stands can be removed from the area quickly after the data
collection stage. (ii) Since markers can be recognized from any angle,
choosing common orientations like the multiples of $\frac{\pi}{2}$ will help
greatly. With standard elevations and easy orientation labeling, the
experimenter is left with measuring the marker positions in 2D coordinates.
(iii) We use two motion cameras not to miss any marker at any time interval,
but this can be regarded as overkill. A single camera will also be sufficient
to detect the markers, but the focus of the camera should be static, as
variable focal length will disrupt the marker pose estimations. (iv) Instead
of mobile beacons, one can use a mobile sensor (i.e. computer or smartphone)
to capture the BLE packets emitted by the beacons that furnish the area.
Several other sensors can be introduced into the setup, leading to more data
that can improve the position estimation accuracy.

###  3.4 Data collection software setup

For sharing data among multiple computers connected through Ethernet, we use
the ROS middleware. ROS middleware is a collection of open-source frameworks
that are intended to be used in robot software development. Besides the pure
robotics services like low-level device control, hardware abstraction, package
management and common robotics functionality; the middleware provides a
message-passing service between processes running at different devices. The
messages of sensor data, control, state, planning, actuator information are
shared among processes, which are called nodes. The nodes can publish or
subscribe to the topics where the messages are streamed. The middleware is
developed to be platform and language independent.

Figure 9: Data collection and communication architecture: A ROS master
computer (i) runs a master node that orchestrates and takes the role of a
system access point for the other nodes. Dedicated RPi devices (ii) run a
publisher node for each sensor to publish the sensor data into the system. Any
other device (iii) can be used for live data visualization, collection and
higher level processing.

In this work, the middleware is used for message-passing, visualizing and
recording data. However, with its high and easy portability degree, the
infrastructure allows plugging any node in real-time. The proposed
architecture is open to improvement in the domain of wireless sensor networks.
Our ROS architecture, as seen in Figure 9, has three main node types: we have
a master node (i) for orchestration and an entry point for other nodes. We
prefer to dedicate another computer, RPi-m, that runs the master node. This
node is responsible for triggering the publishers in the network. The network
can easily be expanded with other nodes. A node that wants to enter into the
ROS network has to present itself to this master node. Each computer, RPi-$i$,
runs three publisher nodes, Publisher-$ij$, (ii), which are responsible for
scanning the corresponding Bluetooth devices, Sensor-$ij$, generate RSSI
parameters from incoming Bluetooth packets, and publish these parameters as
ros messages. The published messages include a timestamp, MAC address of the
source beacon, MAC address of the sensor that captures the Bluetooth packet
and the RSSI value in decibels. Finally, we have a visualizer (iii) that
subscribes to the RSSI topic messages and displays them in live windows with
respect to time. The visualizer is also responsible for recording the data in
log files for each sensor. The master, publisher and visualizer nodes are
running on Python-3 scripts. The visualizers also have a graphical user
interface based on PyQt5 wrappers (see Figure 10).

Figure 10: RSSI data flow on ROS (left) and the live data visualizer (right).
Figure 11: AR localization toolkit for estimating, visualizing, synchronizing
and logging data.

To ease the synchronization of video and RSSI data, we develop and utilize a
custom toolkit written with PyQt5. A snapshot of the toolkit can be seen in
Figure 11. The toolkit is used for performing several tasks required for
calibration, pose estimation and synchronization processes: such as displaying
video frames from two cameras and navigate through the frames freely;
detecting markers given the proper calibration parameters and displaying them
in video frames; estimating raw positions given the markers; applying
different elimination strategies; toggling Kalman filter and tuning its noise
parameters on the fly; displaying raw, eliminated or smoothed poses on the map
with their related markers; taking snapshots from the map (most of the figures
including maps in this article are generated by this toolkit); marking the
needed timestamps for synchronization; displaying logged or predetermined
trajectories; and logging the obtained poses into data files.

###  3.5 Data collection experiments

We collect two types of data simultaneously: timestamped RSSI data captured by
BLE sensors, and position estimations through the videos captured by the
cameras (see Figure 12(a)).

(a) Detected markers in forward and backward camera views (left) and
corresponding raw pose estimations (right). (b) Predetermined paths of the
experiments.

Figure 12: Position data collection and experiments

We prepare five experimental predetermined paths for testing the performance
of our methodology: three distinct straight paths, a rectangular path and a
zigzagging one. The camera-beacon bundle is driven along these paths trying to
respect a linear and smooth motion, simulating human movement. The
predetermined paths are also given in Figure 12(b). The data collection
processes are varied by changing the pace on the straight paths and adding
rotation on the rectangular and zigzagging paths, which makes a dataset of 15
trajectories.

The collected data are processed by varying parameters for the elimination
strategies, the closest selection, $C$, and the outlier elimination, $U$; and
with the measurement error matrix of Kalman filter,
$\mathbf{Q}=q\mathbf{I}_{d}$. We predefine the transition noise matrix as
$\mathbf{R}_{t}=0.01\mathbf{I}_{d}$, which favors the pose estimations being
closer to the previous ones. Some of the estimated paths are given in Figure
13.

Figure 13: Sample paths of collected data: Different straight paths, bendy
paths with and without point rotations, and paths at different paces.

There is a reason that we collect data from more complicated paths along with
the pure straight paths. The rectangular paths investigate an expected
phenomenon whether the beginning and the end of the trajectories point the
same positions on the map, closing the loop. We observe that these same
positions (and even orientations) are successfully captured. Moreover, we
witness that even a point rotation may complicate the path on its own, because
of the relative position of the camera-beacon bundle with respect to the
rotation axis. However, we see that the selected parameters of the Kalman
filter can compensate for this jumping behavior successfully. In our
perspective, a complicated path refers to different things: this may be
related to paths with turns and covering more space in the map, but we think
that the pace makes things more complicated than choosing bendy paths, so we
investigate this phenomenon and add the related results at the end of Section
4.

###  3.6 Position annotated dataset

The dataset provides the trajectory data files that include the position
annotated RSSI data collected while navigating the special platform in the
test area. There are 15 different trajectory files with varying lengths, forms
and paces. Each sample in a file comes with a timestamp, a transmitter
address, a sensor address, an RSSI reading, a ground truth position in 3D and
its associated rotation matrix. We also supply the map of the test area, the
files for the configuration of the sensors and various parameters that are
required to convert the measurements in meters into pixels, and vice versa.
For the completeness of an indoor localization project, we publish auxiliary
data of raw fingerprints, their processed histograms, two sets of estimated
probabilistic radio frequency maps in the grid form and sets of occupancy maps
[30].

##  4 Results

We hypothesize that raw poses of the camera-beacon bundle display two types of
faults, by visualizing raw poses given in Figure 12(a): (i) they may be highly
erroneous and positioned very far from the true position, because of highly
incorrect and inverted axes of occluded or very distant markers (e.g. two
points in the lower part of the map), or (ii) they may fluctuate around a true
point because of pixel-wise problems with smaller errors (e.g. the ones around
the dominant cluster). The former types of errors can be controlled by
discarding the estimations due to the farthest markers, and then an outlier
targeted pruning can be employed on the remaining poses for the latter cases.

Figure 14: Deviations from the predetermined linear paths of raw poses and
isolated elimination strategies: $C$ for close marker selection and $U$ for
outlier elimination.

To evaluate the estimated trajectories, we report how the estimations deviate
from the predetermined linear paths, knowing that data are collected by the
experimenter that follows straight segments as given in Figure 12(b). We first
compare the isolated versions of the two elimination strategies with the
errors of raw poses (without elimination). According to Figure 14 and in
accordance with our hypothesis, it is obvious that a pose elimination strategy
is mandatory, and between these isolated elimination strategies, the outlier
elimination strategy with $U=15$ performs better than others, giving an error
median of 0.045. Moreover, the error distributions for the selection strategy
display undesired skewed distributions with high numbers of outliers. However,
even though a pure close marker selection strategy performs significantly
worse than a pure outlier elimination strategy, we still prefer to begin with
a selection strategy because of its speed and ability to discard probable
systematic outliers due to misestimated rotation axes.

Figure 15: Deviations from the predetermined linear paths of the elimination
strategy combinations with different parameters: $C$ for close marker
selection and $U$ for outlier elimination.

In Figure 15, we compare different combinations of the pose elimination
strategies with the isolated selection strategy with $U=15$. Whereas there is
not a significantly best pose elimination combination, a combination of a
selection strategy with $C=20$ and a succeeding outlier elimination with
$U=10$ distinguishes among the experiments performed, with a median and mean
deviations of 0.043 m and 0.068 m, respectively. Usage of a close marker
selection is also preferable because a direct application of the greedy
outlier elimination strategy is computationally intensive, so it would be also
intuitive to precede it with the quicker close marker selection strategy,
since the far markers tend to yield misleading estimations. Aligned with this
preference, we select a combination of the pose elimination strategies with
$C=20$, $U=10$, which corresponds to quickly reducing the number of poses from
about 100-150 down to 20 with the close marker selection strategy, and then
applying the outlier elimination method on a set of 20 poses to obtain a final
set of 10 poses.

Figure 16: Deviations from the predetermined linear path after the Kalman
filter with different noise covariance multipliers, preceded by a combination
of selection and outlier detection strategies.

On the selected pose estimations after the elimination techniques, we employ a
Kalman filter to attain smoother trajectories. Figure 16 summarize the
parameters for the noise covariance multiplier, $q$, and the deviations after
employing the Kalman filter on the preprocessed pose estimations. We decide
that the deviation medians can be reduced down to 0.034 meters and the means
under 0.05 meters with a covariance matrix of $\mathbf{Q}=2.0\mathbf{I}_{d}$.

Figure 17: Deviations from the predetermined linear path after the Kalman
filter with different covariance multipliers and lag penalty application.

Smoothing with the Kalman filter misleads to inaccurate trajectories, because
we discard time up to this point. We measure and compare the deviations of the
estimations from the predefined ideal paths. Smoother paths yield lower
deviations, however, when the Kalman filter is applied on the pose
estimations, dragging is also possible if the noise covariance multiplier,
$q$, of the Kalman filter is chosen to be too high. If this covariance is
excessively high the estimated trajectory is surely smooth, yielding very low
deviances, but fails to estimate the real positions with respect to time. This
dragging can be evaluated by measuring the distance between a real known end
point of a segment and the estimated points at the same timestamp. To avoid
high dragging situations, we penalize the experiment by accumulating these
drag distances as errors on the deviations, and we obtain Figure 17. The
parameters with $q>0.2$ are highly penalized, meaning that estimated paths do
not lead to accurate localization, and with $q<0.1$, less smoothing is
applied. We conclude that choosing $q\in[0.1,0.2]$ for these experiments
results with the best performances.

(a) (b) (c) (d)

Figure 18: Experiments with different strategies with black dots denoting the
estimated positions and green segments the predefined paths: (a) raw position
estimations, (b) estimations after only the Kalman filter, (c) only
elimination strategies, (d) and the Kalman filter preceded by elimination
strategies.

One may propose to apply a Kalman filter directly on the raw position
estimations (Figure 18(a)) without any prior elimination strategy (see Figure
18(b)), however, as the raw poses are highly erroneous because of mainly
misdetection of marker orientations, the Kalman filter is unable to cope with
raw poses unaidedly. An elimination strategy is mandatory. The positions after
the elimination strategies still draw fluctuating trajectories as seen in
Figure 18(c), but the selected position estimations are distributed around the
true paths. These fluctuations may be eliminated with online filtering
mechanisms. Figure 18(d) is the result of first applying the aforementioned
elimination strategies on the positions estimated by the ArUco marker
detection algorithm, and then smoothing the intermediate results using the
Kalman filter with proper noise covariance parameters. We see that the green
predetermined path can also be visually matched with the filtered positions.

We consider a final situation where Kalman parameters are also affected by the
speed of the navigation. We also test the same experiment at three different
paces. The results in Figure 19 show that choosing $q=0.2$ and $q=0.1$ for
slow navigations (at about 0.3 m/s) and for quicker navigations (at about 0.4
m/s) respectively will yield better performances.

Figure 19: Deviations from the predetermined linear path after employing the
Kalman filter with different noise covariance multipliers and lag penalty
application for the experiments at varying paces.

##  5 Conclusion

This work proposes techniques for collecting high precision position data
using AR markers and annotating BLE RSSI with precise positions in both time
and space domain. We decorate a 364 square meter office area with AR markers
and BLE sensing equipment to collect simultaneous video and RSSI data using a
low cost and practical custom setup. The videos are used to estimate the noisy
positions, which are then eliminated and filtered into smoother trajectories,
highly correlated with the ideal paths. With the aforementioned processing and
synchronization techniques, the RSSI data can be annotated by the precise
position estimations to obtain position annotated RSSI datasets.

The results show that with AR markers we can achieve very high precision
positioning (under 0.05 meters) in controlled indoor areas. This makes the
proposed AR based positioning method far more accurate than the state-of-art
WiFi or BLE based positioning. Of course, the aim is not to compare the
positioning systems. Instead, the proposed setup and method can serve the
indoor positioning domain as a benchmark tool to measure and verify other yet
less precise localization performances in the indoor positioning domain.
Moreover, we believe that there still exists some area to make the AR based
positioning system more precise, especially with better illumination
conditions, more rigorous measurements or larger and more markers.

In the literature, indoor positioning algorithms are evaluated mostly by
manually measuring the predefined reference points. The lack of high frequency
position annotations prevents the community from standardizing the evaluation
procedures. The position estimations of the tracking algorithms are assessed
using various distinct techniques. This need makes high frequency ground truth
data very valuable. Besides, we bring forth not only a method for fine grained
ground truth data collection, but also a method for collecting data
practically. Considering the lack of high precision data sets in the indoor
positioning community, we make the data publicly available for the researchers
who want to test their algorithms based on RSSI data.

We regard the AR based localization as a one-time prior assistive tool for the
calibration of a wireless indoor positioning system. It is intended for
collecting position annotated data which can later be used for parameter
estimation of the tracking algorithms, fingerprinting and more importantly
evaluation of the whole system in terms of position estimation accuracy.

Whereas, in this work, the data are shown to be processed and recorded by a
single node, the infrastructure and the dataset allow the researchers to
approach the wireless based indoor tracking problem through a distributed and
decentralized manner. As an example, as a future work, it can be considered to
handle the problem in a pure decentralized way, where each computer or sensor
node is to estimate the position of an emitter and share its belief with the
neighboring processors. This approach will make the setup easily scalable to
any area.

The main focus of this work is to propose a method for collecting position
annotated wireless data. Nevertheless the most prominent deficiency of the
output is that the dataset is not constructed to reflect the human body
influence on the wireless signals, even though there is a person that drives
the navigable platform. As wireless based positioning systems are usually
considered to be human or living things centered, the lack of this influence
in the dataset may be considered important. Moreover, the dataset cannot be
projected to the true nature of the crowded environments in terms of the
signal quality. However, for the purposes of tracking non-living things or
living things in sparse environments, the dataset still reflects the true
nature and is valid to become a benchmark testing of wireless based
positioning and tracking systems. Besides, the method itself and the system
architecture are not prone to human body influence as they are already
designed to be redundant and scalable.

Collecting wireless data using the current setup and constructing a new
dataset influenced by the human body are left as future work highly due to the
current pandemic situation. The future data set should include varying
versions of the same experiment where the object to be tracked is held,
pocketed or worn as an accessory. Moreover varying the movement and the number
of people is another issue that will present different natural conditions.

As more future work, with proper tunings, this method can be used for
automatic, online and fast fingerprinting for wireless signal based
localization, which is substantially a challenging task for the researchers
working in the domain. With the markers installed in the area, the camera-
beacon bundle may be navigated to collect data from any point, and these data
can be processed to form the fingerprints fast and easily. Moreover, a self
navigating robot can be used for real automatization of fingerprinting.
Furthermore, the technique also permits collecting data in three dimensions,
which will prove to be essential for indoor localization of robots,
specifically drones.

## Acknowledgment

This work is supported by the Turkish Directorate of Strategy and Budget under
the TAM Project number 2007K12-873. The authors would like to thank Sergio
Bromberg for his contribution of the back to back GoPro camera mount 3D model
to Thingiverse and PoiLabs for supplying POI Beacons for this work.

## References

## References

  * [1] I. Bisio, F. Lavagetto, M. Marchese, A. Sciarrone, Smart probabilistic fingerprinting for wifi-based indoor positioning with mobile devices, Pervasive and Mobile Computing 31 (2016) 107–123.
  * [2] R. Faragher, R. Harle, Location fingerprinting with bluetooth low energy beacons, IEEE Journal on Selected Areas in Communications 33 (11) (2015) 2418–2428.
  * [3] J. Torres-Sospedra, Óscar Belmonte-Fernández, G. M. Mendoza-Silva, R. Montoliu, A. Puertas-Cabedo, L. E. Rodríguez-Pupo, S. Trilles, A. Calia, M. Benedito-Bordonau, J. Huerta, 3 - lessons learned in generating ground truth for indoor positioning systems based on wi-fi fingerprinting, in: Geographical and Fingerprinting Data to Create Systems for Indoor Positioning and Indoor/Outdoor Navigation, Intelligent Data-Centric Systems, Academic Press, 2019, pp. 45 – 67.
  * [4] G. M. Mendoza-Silva, M. Matey-Sanz, J. Torres-Sospedra, J. Huerta, Ble rss measurements dataset for research on accurate indoor positioning, Data 4 (1) (2019) 12.
  * [5] P. Baronti, P. Barsocchi, S. Chessa, F. Mavilia, F. Palumbo, Indoor bluetooth low energy dataset for localization, tracking, occupancy, and social interaction, Sensors 18 (12) (2018) 4462.
  * [6] M. Girolami, F. Mavilia, F. Delmastro, A bluetooth low energy dataset for the analysis of social interactions with commercial devices, Data in Brief 32 (2020) 106102.
  * [7] E. Sansano-Sansano, F. J. Aranda, R. Montoliu, F. J. Álvarez, Ble-gspeed: A new ble-based dataset to estimate user gait speed, Data 5 (4) (2020) 115.
  * [8] O. S. Eyobu, A. Poulose, D. S. Han, An accuracy generalization benchmark for wireless indoor localization based on imu sensor data, in: 2018 IEEE 8th International Conference on Consumer Electronics - Berlin (ICCE-Berlin), 2018, pp. 1–3.
  * [9] C. M. de la Osa, G. G. Anagnostopoulos, M. Togneri, M. Deriaz, D. Konstantas, Positioning evaluation and ground truth definition for real life use cases, in: 2016 International Conference on Indoor Positioning and Indoor Navigation (IPIN), 2016, pp. 1–7.
  * [10] A. Poulose, D. S. Han, Hybrid deep learning model based indoor positioning using wi-fi rssi heat maps for autonomous applications, Electronics 10 (1) (2021) 2.
  * [11] H. Ai, K. Tang, W. Huang, S. Zhang, T. Li, Fast fingerprints construction via GPR of high spatial-temporal resolution with sparse RSS sampling in indoor localization, Computing 102 (3) (2020) 781–794.
  * [12] S. Adler, S. Schmitt, K. Wolter, M. Kyas, A survey of experimental evaluation in indoor localization research, in: 2015 International Conference on Indoor Positioning and Indoor Navigation (IPIN), 2015, pp. 1–10.
  * [13] J. Kim, H. Jun, Vision-based location positioning using augmented reality for indoor navigation, IEEE Transactions on Consumer Electronics 54 (3) (2008) 954–962.
  * [14] G. C. La Delfa, S. Monteleone, V. Catania, J. F. De Paz, J. Bajo, Performance analysis of visualmarkers for indoor navigation systems, Frontiers of Information Technology & Electronic Engineering 17 (8) (2016) 730–740.
  * [15] R. S. Xavier, B. M. F. da Silva, L. M. G. Goncalves, Accuracy analysis of augmented reality markers for visual mapping and localization, in: 2017 Workshop of Computer Vision (WVC), 2017, pp. 73–77.
  * [16] E. Ujkani, J. Dybedal, A. Aalerud, K. B. Kaldestad, G. Hovland, Visual marker guided point cloud registration in a large multi-sensor industrial robot cell, in: 2018 14th IEEE/ASME International Conference on Mechatronic and Embedded Systems and Applications (MESA), 2018, pp. 1–6.
  * [17] A. Alnabhan, B. Tomaszewski, Insar: Indoor navigation system using augmented reality, in: Proceedings of the Sixth ACM SIGSPATIAL International Workshop on Indoor Spatial Awareness, ISA ’14, Association for Computing Machinery, New York, NY, USA, 2014, p. 36–43.
  * [18] I. A. Koc, T. Serif, S. Gören, G. Ghinea, Indoor mapping and positioning using augmented reality, in: 2019 7th International Conference on Future Internet of Things and Cloud (FiCloud), 2019, pp. 335–342.
  * [19] D. Byrne, M. Kozlowski, R. Santos-Rodriguez, R. Piechocki, I. Craddock, Residential wearable rssi and accelerometer measurements with detailed location annotations, Scientific Data 5 (1) (2018) 180168.
  * [20] F. J. Romero-Ramirez, R. Muñoz-Salinas, R. Medina-Carnicer, Speeded up detection of squared fiducial markers, Image and Vision Computing 76 (2018) 38–47.
  * [21] S. Garrido-Jurado, R. Munoz-Salinas, F. Madrid-Cuevas, R. Medina-Carnicer, Generation of fiducial marker dictionaries using mixed integer linear programming, Pattern Recognition 51 (2016) 481–491.
  * [22] K. K. Liang, Efficient conversion from rotating matrix to rotation axis and angle by extending rodrigues’ formula (Oct. 2018).  arXiv:1810.02999.
  * [23] M. A. Fischler, R. C. Bolles, Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography, Commun. ACM 24 (6) (1981) 381–395.
  * [24] J. Wang, A. Hu, X. Li, Y. Wang, An improved pdr/magnetometer/floor map integration algorithm for ubiquitous positioning using the adaptive unscented kalman filter, ISPRS International Journal of Geo-Information 4 (4) (2015) 2638\.
  * [25] Z. Chen, Q. Zhu, Y. C. Soh, Smartphone inertial sensor-based indoor localization and tracking with ibeacon corrections, IEEE Transactions on Industrial Informatics 12 (4) (2016) 1540–1549.
  * [26] M. Quigley, B. Gerkey, K. Conley, J. Faust, T. Foote, J. Leibs, E. Berger, R. Wheeler, A. Ng, Ros: an open-source robot operating system, in: Proc. of the IEEE Intl. Conf. on Robotics and Automation (ICRA) Workshop on Open Source Robotics, Kobe, Japan, 2009.
  * [27] G. H. An, S. Lee, M.-W. Seo, K. Yun, W.-S. Cheong, S.-J. Kang, Charuco board-based omnidirectional camera calibration method, Electronics 7 (12) (2018) 421.
  * [28] S. Tomažič, I. Škrjanc, An automated indoor localization system for online bluetooth signal strength modeling using visual-inertial slam, Sensors 21 (2021) 2857.
  * [29] H. Obeidat, W. Shuaieb, O. Obeidat, R. Abd-Alhameed, A review of indoor localization techniques and wireless technologies, Wireless Personal Communications 119 (1) (2021) 289–327.
  * [30] F. S. Daniş, A. T. Cemgil, C. Ersoy, Adaptive sequential monte carlo filter for indoor positioning and tracking with bluetooth low energy beacons, IEEE Access 9 (2021) 37022–37038.


# Anomalous criticality coexists with giant cluster in the uniform forest
model

Hao Chen<EMAIL_ADDRESS>School of the Gifted Young, University
of Science and Technology of China, Hefei, Anhui 230026, P.R. China  Jesús
Salas<EMAIL_ADDRESS>Departamento de Matemáticas, Universidad Carlos III
de Madrid, Avenida de la Universidad 30, 28911 Leganés, Spain  Grupo de
Teorías de Campos y Física Estadística, Instituto Gregorio Millán Barbany,
Universidad Carlos III de Madrid, Unidad Asociada al Instituto de Estructura
de la Materia, CSIC, Spain  Youjin Deng<EMAIL_ADDRESS>Department of
Modern Physics,
University of Science and Technology of China, Hefei, Anhui 230026, P.R. China
Hefei National Laboratory,
University of Science and Technology of China, Hefei, Anhui 230088, P.R. China

(August 28, 2024)

###### Abstract

In percolation theory, the general scenario for the supercritical phase is
that all clusters, except the unique giant one, are small and the two-point
correlation exponentially decays to some constant. We show by extensive
simulations that the whole supercritical phase of the three-dimensional
uniform forest model simultaneously exhibits an infinite tree and a rich
variety of critical phenomena. Besides typical scalings like algebraically
decaying correlation, power-law distribution of cluster sizes, and divergent
correlation length, a number of anomalous behaviors emerge. The fractal
dimensions for off-giant trees take different values when being measured by
linear system size or gyration radius. The giant tree size displays two-length
scaling fluctuations, instead of following the central-limit theorem. In a
non-Gaussian fermionic field theory, these unusual properties are closely
related to the non-abelian continuous OSP$(1|2)$ supersymmetry in the
fermionic hyperbolic plane ${\mathbb{H}}^{0|2}$.

## Introduction.

Percolation studies connectivity in random geometric systems [1, 2]. In
statistical physics, percolation has been of immense theoretical interest,
providing a simple example that undergoes a non-trivial phase transition. The
celebrated Ising and Potts models [3] can be described as a correlated
percolation through the exact Fortuin-Kasteleyn transformation [4, 5, 6].
Thanks to the percolation approach, it was established [7, 8, 9, 10, 11] that
the Ising model in three dimensions (3D) has a sharp continuous phase
transition, and in 4D, it exhibits mean-field critical behavior, proving the
triviality of the 4D Euclidean scalar quantum field theory. Percolation has
also intensively been studied in mathematics, including a list of variations
like $k$-core and explosive percolation, etc [12, 13, 14, 15, 16, 17]. Also,
percolation has been applied in diverse branches of science and industry [18,
19].

In the basic bond percolation, one randomly occupies each lattice edge with
probability $p$ and constructs clusters of connected components. Clusters are
small for small $p$, and the probability (two-point correlation) that two
sites with distance $r$ are in the same cluster decays exponentially as
$g(r)\sim\exp(-r/\xi)$, and the correlation length $\xi$ diverges as
$\xi\sim(p_{c}-p)^{-\nu}$ as threshold $p_{c}$ is approached. At $p_{c}$, the
size $s$ of fractal clusters follows a universal power-law distribution as
$n(s)\sim s^{-\tau}$, and the correlation decays algebraically as $g(r)\sim
r^{2-d-\eta}$. For the supercritical phase ($p>p_{c}$), an infinite cluster of
size $C_{1}$ occupies a nonzero fraction of the lattice—i.e., $m\equiv
C_{1}\,L^{-d}$ converges to a constant as the linear size $L\to\infty$. In
percolation, $m$ plays a role as the order parameter and behaves as
$m\sim(p-p_{c})^{\beta}$ as $p\downarrow p_{c}$. Further, the second moment of
cluster sizes, $\chi\equiv L^{-d}\,\sum_{i}|C_{i}|^{2}$, acting as the
magnetic susceptibility, diverges as $\chi\sim|p-p_{c}|^{-\gamma}$. Among
these critical exponents $\nu,\beta,\gamma,\eta,\tau$, two are independent and
the others can be obtained from (hyper)scaling relationships [20].

For $p>p_{c}$, the infinite cluster can be even proved to be unique–i.e.,
there is one and only one giant cluster. All off-giant clusters are small with
finite correlation length $\xi^{\prime}$, and the correlation $g^{\prime}(r)$,
with the giant cluster being excluded, vanishes exponentially. In 2D, the
supercritical and subcritical ($p<p_{c}$) phases are dual to each other.
Further, the smallness of $\xi^{\prime}$ predicts that the size fluctuation of
the giant cluster would obey the central-limit theorem and follow a normal
(Gaussian) distribution.

In this work, we study the percolative properties of the supercritical phase
for the (weighted) uniform forest (UF) model in 3D [21, 22, 23, 24]. The UF
model (also called the arboreal gas) consists of a spanning forest of trees
(acyclic clusters), in which each tree is weighted by a factor $w$ per
occupied bond; the statistical weight of any configuration ${\cal{A}}$ can be
written as $\pi({\cal A})=w^{|\cal A|}\cdot\delta_{c({\cal A}),0}$. This is
similar to that for bond percolation with probability $p=w/(1+w)$, except the
$\delta$-function constraint on zero cyclomatic number $c({\cal A})=0$.
Further, as bond percolation, the 3D UF model undergoes a continuous
transition $w_{c}$ [Fig. 1(a)]. It was found that $w_{c}=0.43365(2)$,
$\nu=1.28(4)$, and $\beta/\nu=0.4160(6)$ [25]. Moreover, the supercritical
phase also has an infinite and unique tree [26], similar to percolation. As
shown in the inset of Fig. 1(a) for $w=0.9>w_{c}$, $m$ quickly saturates to
$m_{0}=0.685(6)$, and is clearly long-ranged.

## Results.

Despite these analogues to bond percolation, we find that the whole
supercritical phase of the 3D UF model exhibits the simultaneous emergence of
a surprisingly rich variety of critical behaviors and of a unique giant tree,
providing a counter example for the standard percolation scenario.

Figure 1:  Simultaneous emergence of critical behaviors and a giant cluster.
(a) The percolation transition and the giant tree of size $C_{1}$ in the
supercritical phase. The approximately common intersection of the ratios $Q$
for different sizes $L$ indicate the threshold at $w_{c}\approx 0.43365$. The
inset shows that, while the order parameter $m(w_{c})\equiv C_{1}\,L^{-3}$
algebraically vanishes, it quickly saturates to a constant, which is
$m_{0}=0.685(6)$ for $w_{s}=0.9$. (b) Emergent critical behaviors at $w=w_{s}$
demonstrated by the power-law divergence of the Fourier-transformed
susceptibility $\chi_{\bm{k}}\sim L^{1.99(2)}=L^{2}$, and the algebraic
decaying of the off-giant correlation $g^{\prime}(\bm{r})\sim|\bm{r}|^{2-d}$.

In the supercritical phase, the off-giant clusters are fractal and display
critical scaling behaviors that are generally expected from percolation theory
at criticality. (i) As shown in Fig. 1(b), the Fourier-transformed
susceptibility $\chi_{\bm{k}}$ (the definition will be given later) diverges
as $\chi_{\bm{k}}\sim L^{\gamma/\nu}=L^{1.99(2)}$. This result is in good
agreement with $\chi_{\bm{k}}\sim L^{2}$, following the standard finite-size
scaling (FSS) _Ansatz_. The off-giant correlation $g^{\prime}({\bm{r}})$
algebraically decays as $g^{\prime}(\bm{r})\sim|\bm{r}|^{2-d-\eta}$ with
$\eta=0$ [24]. Indeed, the scaling relation $2-\eta=\gamma/\nu$ is satisfied.
(ii) The size distribution of the clusters $n(s,L)$ has two terms: one
accounts for the contribution of the off-giant trees, while the other takes
care of the giant one. The former term contains a power law $s^{-\tau}$ times
the size distribution of the off-giant clusters $\tilde{n}^{\prime}(s,L)$,
which is governed by the second-largest cluster of size $C_{2}\sim
L^{d_{C_{2}}}$ with $d_{C_{2}}=2.29(2)$. The latter term has a prefactor
$L^{-d}$, and it is governed by the size of the largest cluster $C_{1}\sim
L^{d_{C_{1}}}=L^{3.000(2)}=L^{d}$. Putting all together, we have [see Fig.
2(a)]:

$n(s,L)\;=\;s^{-\tau}\,\tilde{n}^{\prime}(s\,L^{-d_{C_{2}}})+L^{-d}\,n_{1}(s,L)\,.$
(1)

The Fisher exponent is $\tau\equiv\tau_{2}=2.31(2)$, and the distribution
peak, arising from $C_{1}$, defines an effective Fisher exponent $\tau_{1}=2$.
The hyperscaling relations $\tau_{i}=1+d/d_{C_{i}}$ for $i=1,2$ are satisfied.

In addition, the supercritical phase exhibits a variety of unusual critical
behaviors. (iii) The standard FSS theory predicts that the critical
correlation length is $\xi\sim O(L)$, and, indeed, the gyration radius $R_{1}$
of the largest tree scales as $R_{1}\sim L$ for $w\geq w_{c}$. However, for
$w>w_{c}$, the off-giant correlation length, characterized by the gyration
radius $R_{2}$ of the second-largest cluster, diverges sublinearly as
$R_{2}\sim L^{\kappa_{2}}$ with $\kappa_{2}=0.76(2)$ [see the inset of Fig.
2(b)]. This indicates that the supercritical phase has two length scales—i.e.,
$L$ and $L^{0.76}\ll L$. Typically, the size $s$ of a fractal object depends
on its gyration radius $R$ as $s\sim R^{d_{f}}$, and the generic fractal
dimension $d_{f}$ of the giant cluster coincides with the finite-size fractal
dimension $d_{C_{1}}$ since $R_{1}\sim L$. (iv) For the off-giant clusters at
$w>w_{c}$, however, the generic and finite-size fractal dimensions take
different values—e.g., $d_{C_{2}}\neq d_{f_{2}}$ for the second-largest
cluster. Instead, one has $d_{f_{2}}=d_{C_{2}}/\kappa_{2}$, since $C_{2}\sim
L^{d_{C_{2}}}\sim R_{2}^{d_{f_{2}}}$ and $R_{2}\sim L^{\kappa_{2}}$. We obtain
$d_{f_{2}}=3.01(8)\approx 3$, which agrees well with $d$ and $d_{C_{1}}$; this
is further demonstrated for all the off-giant clusters in Fig. 2(b). Thus, the
off-giant clusters share the same generic fractal structure as the giant one.

Figure 2: Fractal structures and two-length scales at $w=w_{s}$. (a) Scaling
behavior of the cluster-size distribution (1). The inset displays $C_{2}\sim
L^{d_{C_{2}}}$ and with a line of slope $d_{C_{2}}=2.29$. (b) Power-law
dependence of cluster size $s$ on the gyration radius $R$. The curves seem to
collapse around two parallel lines of slope $d_{f}=3$. The inset shows the
radii of the two largest clusters, of which the fits give exponents
$\kappa_{1}=0.999(4)$, and $\kappa_{2}=0.76(2)$.

The giant cluster, occupying about $70\%$ of the lattice for $w=0.9$, exhibits
interesting critical behaviors. In a traditional supercritical phase with
finite correlation length, the central limit predicts that the giant-cluster
size should follow a normal (Gaussian) distribution and the normalized
fluctuation $F_{1}\equiv{\rm Var}({\cal C}_{1})\,L^{-d}$ should converge to
some constant. (v) As in Fig. 3(b), however, we find $F_{1}\sim L^{d_{F_{1}}}$
with $d_{F_{1}}=2.03(6)$. This is counter intuitive because one would naively
expect the central-limit theorem to apply because the ratio
$\xi^{\prime}/L\sim L^{\kappa_{2}-1}\to 0$ as $L\to\infty$ (here
$\xi^{\prime}$ is the off-giant correlation length).

We then consider the probability density function (PDF)
$f_{\mathcal{C}_{1}}(C_{1},L)$ of the random size $\mathcal{C}_{1}$ of the
giant cluster in a lattice of linear size $L$. In the standard FSS theory, by
rescaling
$\mathcal{X}=(\mathcal{C}_{1}-\langle\mathcal{C}_{1}\rangle)\,L^{-d_{C_{1}}}$,
and transforming
$f_{\mathcal{X}}(x)\,dx=f_{\mathcal{C}_{1}}(C_{1},L)\,dC_{1}$, one can obtain
a universal and $L$-independent function $f_{\mathcal{X}}(x)$. (vi) In the
supercritical phase of the UF model, however, we find that the
$f_{\mathcal{C}_{1}}(C_{1},L)$ data for different values of $L$ cannot be
collapsed onto a unique curve by any single exponent like $d_{C_{1}}=3$. Thus,
we consider the probability $f_{\mathcal{X}_{1}}(x_{1},L)\,dx_{1}$ for the
rescaled random deviation
$\mathcal{X}_{1}\equiv(\mathcal{C}_{1}-\langle\mathcal{C}_{1}\rangle)\,L^{-d_{C_{2}}}$
with $d_{C_{2}}=2.29$. Then, the $f_{\mathcal{X}_{1}}(x_{1},L)$ data
approximately collapses well near $x_{1}=0$ and for $x_{1}>0$ [Fig. 3(a)].
Nevertheless, $f_{\mathcal{X}_{1}}(x_{1},L)$ has a wide-range shoulder for
$x_{1}\ll 0$, for which an approximate data collapse can be achieved by
$L^{\delta}\,f_{\mathcal{X}_{1}^{\prime}}(x_{1}^{\prime})\,dx_{1}^{\prime}\equiv
f_{\mathcal{X}_{1}}(x_{1},L)\,dx_{1}$ with
$\mathcal{X}_{1}^{\prime}\equiv(\mathcal{C}_{1}-\langle\mathcal{C}_{1}\rangle)\,L^{-3}$
and $\delta=0.77$. This means that the whole configuration space is roughly
partitioned into two sectors: one takes up a finite configuration-space volume
while the other vanishes asymptotically as $L^{-\delta}$. In the dominant
sector, the critical fluctuation of ${\cal C}_{1}$ is governed by $d_{C_{2}}$
for off-giant clusters. In the vanishing sector, the variance ${\rm Var}({\cal
C}_{1})$ is $\sim O(L^{2d})$. Note that this exponent takes the largest
possible value. We further sample ${\rm Var}({\cal C}_{1})$ conditioned on
${\cal C}_{1}-C_{1}\geq 0$ for the dominant (dom) sector, and $L^{-d}({\cal
C}_{1}-C_{1})\leq a$ for the vanishing (van) sector, where we take $a=-0.1$.
We obtain $d_{F_{1}}(\mathrm{van})=2.96(4)$ and
$d_{F_{1}}(\mathrm{dom})=1.58(2)$ [see Fig. 3(b)], the latter of which gives
$d_{C_{2}}=2.29(1)$ from relation $2d_{C_{2}}-3=d_{F_{1}}$. Note that
$d_{F_{1}}({\rm total})=2.03(6)$ for the total configuration space is distinct
from $d_{F_{1}}(\mathrm{dom})$ or $d_{F_{1}}(\mathrm{van})$, indicating that
the crossover regime also plays an important role.

Figure 3: Two-length-scale fluctuation of $\mathcal{C}_{1}$ at $w=w_{s}$. (a)
PDF $f_{\mathcal{X}_{1}}(x_{1},L)$ of the rescaled random variable
$\mathcal{X}_{1}\equiv(\mathcal{C}_{1}-\langle\mathcal{C}_{1}\rangle)\,L^{-d_{C_{2}}}$
with $d_{C_{2}}=2.29$. A good collapse of data is observed in the peak region.
The inset shows the PDF
$L^{\delta}\,f_{\mathcal{X}^{\prime}_{1}}(x_{1}^{\prime})$ of the random
variable
$\mathcal{X}_{1}^{\prime}\equiv(\mathcal{C}_{1}-\langle\mathcal{C}_{1}\rangle)\,L^{-3}$
with $\delta=0.77$, showing an approximate collapse in the shoulder region.
(b) Normalized variance $F_{1}=\mathrm{Var}(\mathcal{C}_{1})\,L^{-3}$ in
different sectors. The fits give exponents $d_{F_{1}}(\text{total})=2.03(6)$,
$d_{F_{1}}(\text{dominant})=1.58(2)$, and
$d_{F_{1}}(\text{vanishing})=2.96(4)$.

## Theoretical Insights.

The critical behaviors of the 3D UF model can be partially understood from its
relation to statistical mechanical systems and from the perspective of quantum
field theory. The UF model corresponds to the $q\to 0$ limit of the $q$-state
Potts model [27, 28, 3, 29] in the Fourtuin–Kasteleyn random-cluster
representation [4, 5, 6]. Particularly, this limit should be taken such that
$v/q\equiv(e^{J}-1)/q=w$ is held fixed ($J$ is the reduced nearest-neighbor
coupling).

Caracciolo  _et al._ [30] mapped the UF model onto a non-Gaussian fermionic
theory with the non-abelian continuous OSP$(1|2)$ supersymmetry by
generalizing Kirchhoff’s matrix-tree theorem. They also showed how to map this
model in perturbation theory to all orders in $1/w$ onto a $N$-vector model
analytically continued to $N=-1$. They concluded that in 2D, the model is
asymptotically free, implying the absence of phase transition for any finite
$w>0$ [22, 23]; the criticality occurs at $w=+\infty$ [31, 28]. The relation
of the UF model and the supersphere non-linear sigma model with
${\mathbb{S}}^{0,2}$ was studied in [32, 33, 34].

Recently, the UF model was interpreted [23] as a non-linear sigma model with
the fermionic hyperbolic plane ${\mathbb{H}}^{0|2}$ as the target space. The
hyperbolic symmetry of spin models has many interesting properties [35], and
for the UF models, the relation leads [24] to the existence of a percolation
transition at finite $w_{c}>0$ for $d\geq 3$. It was also proven [26] that the
supercritical phase has, in the infinite-lattice limit, a unique infinite tree
for $d=3,4$. Moreover, it was established [24] that, at large enough $w>0$,
there are massless power-law correlations for $d\geq 3$, as if the model were
at criticality,

$g(\bm{r})\;=\;g_{0}+c\,|\bm{r}|^{2-d}+\cdots\,,$ (2)

where $g_{0}=m^{2}$ comes from the giant cluster, $c=c_{0}+O(1/w)$ is a
constant, and the dots stand for higher-order corrections. Nevertheless,
despite Eq. (2), the percolative properties in the supercritical phase remain
elusive.

## Algorithms and observables.

We use the Sweeny algorithm [36] and simulate the UF model on the simple-cubic
lattice with periodic boundary conditions, for $8\leq L\leq 128$. At every
step, one randomly picks up an edge $e_{ij}$ between sites $i$ and $j$. If
$e_{ij}$ is occupied, the bond is removed with probability $\min(1,1/w)$. If
$e_{ij}$ is empty, it is occupied with probability $\min(1,w)$ if $i$ and $j$
belong to different trees, and, otherwise, it is left unoccupied since a bond
on $e_{ij}$ would generate a cycle. The nontrivial operation is to detect the
connectivity of $i$ and $j$ in a dynamical setting. Using the link-cut tree
data structure [37], the connectivity query can be efficiently implemented in
$O(\log L)$ amortized time. Note that no critical slowing down occurs in 3D
[36, 25, 38].

For a random configuration, we denote the forest of trees as $\\{{\cal
C}_{k}\\}$, and specifically leave $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ for
the largest and second-largest clusters, respectively. For a tree
$\mathcal{C}_{k}$, an arbitrary site is chosen as the origin, and the
“unwrapped” coordinate $\bm{x}_{k}^{i}$ of each site $i$ is obtained by
growing the tree from the origin. This coordinate $\bm{x}_{k}^{i}$ is well
defined, since the path connecting any two sites in a tree is unique. The
mass-center coordinate $\overline{\bm{x}}_{k}=(1/{\cal
C}_{k})\sum_{i\in\mathcal{C}_{k}}\bm{x}_{k}^{i}$, and the squared gyration
radius ${\cal R}_{k}^{2}=(1/{\cal
C}_{k})\sum_{i\in\mathcal{C}_{k}}(\bm{x}^{i}_{k}-\overline{\bm{x}}_{k})^{2}$
are calculated.

By detecting the connectivity between sites $i$ and $j$ over configurations,
we measure the two-point correlation $g({\bm{r}}=\bm{r}_{i}-\bm{r}_{j})$, as
well as the off-giant correlation $g^{\prime}({\bm{r}})$, where $\bm{r}_{i}$
is the standard Euclidean coordinate of site $i$. For simplicity, we choose
$\bm{r}=(r,0,0)$ along the $x$-axis. Moreover, to study the algebraic decaying
behavior of the correlation, an auxiliary Ising spin $s_{i}\in\\{\pm 1\\}$ is
introduced for every site $i$: Independently for each tree, we assign the same
value $s_{i}=1$ or $-1$ with equal probability to all the sites in the tree.
By definition, $g(\bm{r}_{i}-\bm{r}_{j})=\langle s_{i}s_{j}\rangle$. The
magnetization ${\cal M}=\sum_{m}s_{m}$ and its Fourier transform ${\cal
M}(\bm{k})=\sum_{m}s_{m}\exp(i\bm{k}\cdot\bm{r}_{m})$ are sampled, where the
summation is over the whole lattice. The smallest nonzero momenta in the $x$
direction, $\bm{k}=(2\pi/L,0,0)$ is used for simplicity.

Statistical average and probability distribution are then taken over the
configurations generated in simulations—e.g., $C_{k}\equiv\langle{\cal
C}_{k}\rangle$ and $R_{k}\equiv\langle{\cal R}_{k}\rangle$ for $k=1,2$. Also,
we define the normalized fluctuation
$F_{1}\equiv\mathrm{Var}(\mathcal{C}_{1})\,L^{-3}$, the susceptibility
$\chi\equiv\langle{\cal M}^{2}\rangle\,L^{-3}$, the dimensionless ratio
$Q\equiv\langle{\cal M}^{2}\rangle^{2}/\langle{\cal M}^{4}\rangle$, and the
Fourier-transformed susceptibity $\chi_{\bm{k}}\equiv\langle{\cal
M}({\bm{k}}){\cal M}(-{\bm{k}})\rangle\,L^{-3}$.

| $d_{C_{1}}$ | $\kappa_{1}$ | $d_{C_{2}}$ | $\kappa_{2}$ | $d_{F_{1}}$ | $d_{\chi_{\bm{k}}}$
---|---|---|---|---|---|---
tot | 3.000(2) | 0.999(4) | 2.29(2) | 0.76(2) | 2.03(6) | 1.99(2)
dom | 3.002(3) | 1.001(5) | 2.28(2) | 0.78(2) | 1.58(2) | 1.63(3)
van | 2.997(4) | 0.997(6) | 3.00(2) | 1.01(2) | 2.96(4) | 2.83(5)
Table 1: Estimated critical exponents for the supercritical phase at
$w_{s}=0.9$, for the total (tot) configuration space, and for the dominant
(dom) and vanishing (van) sectors.

## Fits.

As a powerful quantity for locating a continuous phase transition, the
crossings of the $Q(w)$ curves for diffent sizes $L$ [Fig. 1(a)] clearly
support the previously determined percolation threshold $w_{c}=0.43365(2)$
[25]. We carry out extensive simulation at $w_{c}$ and $w_{s}=0.9$, deeply in
the supercritical phase. The critical behaviors, shown in Figs. 1–3, have been
qualitatively presented in Results.

We perform the least-squares fits to a power-law _Ansatz_ for any observable
$\mathcal{O}(L)$

$\mathcal{O}(L)\;=\;L^{d_{\mathcal{O}}}\,\left(a_{0}+a_{1}L^{-\omega_{1}}+a_{2}L^{-\omega_{2}}\right)+b_{0}\,.$
(3)

In most cases, we set $\omega_{1}=1$ and $\omega_{2}=2$, and $b_{0}=0$ is
fixed for observables that vanish for $L\to\infty$. As a precaution against
FSS corrections not included in the _Ansatz_ (3), we have performed each fit
by allowing only data with $L\geq L_{\text{min}}$. By studying how the
estimates of the parameters, as well as the $\chi^{2}$ per degree of freedom,
vary as a function of $L_{\text{min}}$, we determine our final estimates and
their error bars.

In particular, we consider the sizes and radii of the largest and second-
largest clusters, the normalized fluctuation $F_{1}$, and the Fourier-
transformed susceptibility $\chi_{\bm{k}}$, which scale as ($j=1,2$)

$\displaystyle C_{j}$ $\displaystyle\sim$ $\displaystyle
L^{d_{C_{j}}},\;\;R_{j}\sim L^{\kappa_{j}},\;\;F_{1}\sim
L^{d_{F_{1}}},\;\;\chi_{\bm{k}}\sim L^{d_{\chi_{\bm{k}}}}\,.\quad$ (4)

At the critical value $w_{c}$, the scaling behaviors follow the standard FSS
theory, which predicts $\kappa_{1}=\kappa_{2}=1$, $d_{C_{1}}=d_{C_{2}}=d_{f}$
($d_{f}$ is the generic fractal dimension), and $d_{F_{1}}=2d_{f}-d$. There is
only one non-trivial exponent $d_{C_{1}}$, which is determined to be
$d_{C_{1}}=2.5840(6)$.

In the supercritical phase with $w=w_{s}$, the final estimates are given in
Table 1. As expected, the largest cluster, occupying a finite fraction of the
lattice, has trivial exponents $d_{C_{1}}=3.000(2)=3$ and
$\kappa_{1}=0.999(4)=1$. The effective Fisher exponent $\tau_{1}$, governing
the decreasing of the distribution peak in Fig. 2(a), is also trivial
$\tau_{1}=1+d/d_{C_{1}}=2.000(1)$.

The finite-size fractal exponent of the second-largest cluster is
$d_{C_{2}}=2.29(2)$, which has not yet been reported to our knowledge. This
gives the Fisher exponent $\tau_{2}=1+d/d_{C_{2}}=2.31(2)$ in Eq. (1). The
gyration radius scales sublinearly versus $L$ with exponent
$\kappa_{2}=0.76(2)$, unexpected from the standard FSS theory. The generic
fractal dimension, $C_{2}\sim R_{2}^{d_{f_{2}}}$, is calculated as
$d_{f_{2}}=d_{C_{2}}/\kappa_{2}=3.01(8)=3$. Surprisingly, $d_{f_{2}}$ is just
the spatial dimension, and Fig. 2(b) further gives $d_{f_{2}}=3$ for all the
off-giant clusters.

By definition, the Fourier-transformed susceptibility is
$\chi({\bm{k}})=L^{-3}\,\sum_{m,n}\langle
s_{m}s_{n}\rangle\exp(i\bm{k}\cdot(\bm{r}_{m}-\bm{r}_{n}))$, where the
contribution from the background term $g_{0}$ in Eq. (2) is eliminated. Thus,
$\chi({\bm{k}})\sim L^{2}$ is expected, and this is strongly supported by the
estimated exponent $d_{\chi_{\bm{k}}}=1.99(2)$. Interestingly, the normalized
fluctuation of the largest tree is governed by exponent $d_{F_{1}}=2.03(6)=2$.

Despite the simplicity of scalings like $C_{1}\sim L^{3}$ and $F_{1}\sim
L^{2}$, the distribution of $\mathcal{C}_{1}$ is sophisticated [Fig. 3(a)].
Thus, we perform separate least-squares fits for the dominant and the
vanishing sectors, respectively conditioned on ${\cal C}_{1}-C_{1}\geq 0$ and
${\cal C}_{1}-C_{1}\leq-0.1\,L^{3}$. The results in Table 1 suggest that the
two sectors have dramatically different scaling behaviors. Particularly, in
the vanishing sector, we find that both the largest and the second-largest
clusters have $d_{C_{j}}=3$ and $\kappa_{j}=1$ for $j=1,2$. Further,
$d_{\chi_{\bm{k}}}=2.83(5)$ suggests that the $r$-dependent decaying of
correlation $g(r)|_{\mathrm{van}}$ is extremely slow (i.e.,
$g(r)|_{\mathrm{van}}\sim r^{-0.17}$). This behavior, together with
$d_{C_{j}}=3$, gives a strong hint for a logarithmic decay:
$g(r)|_{\mathrm{van}}\sim 1/\log(r)$.

## Summary.

While undergoing a typical continuous percolation transition, the 3D UF model
exhibits a variety of critical behaviors in the supercritical phase. The
simultaneous existence of anomalous criticality and of a unique giant cluster
is unexpected from the standard percolation theory. The critical scaling
behaviors not only appear in the off-giant clusters, but also in the
fluctuation of the giant tree. Unlike a conventional critical point, the whole
configuration space can be approximately divided into two configuration
sectors of distinct critical exponents. Further, the overall scaling behaviors
arise from some delicate interplay of the two sectors and of the crossover
regime in between.

Some insight can be borrowed from the fermionic field theory, but a complete
and deep understanding is still needed for this extremely rich critical
behavior. As a return, we believe that our work may also bring some insight
for critical phenomena in the Potts model and the nonlinear-sigma model, which
are two important classes of systems in statistical mechanics and condensed-
matter physics. For instance, for the XY model with long-range interaction,
which was recently found [39, 40] to exhibit critical behaviors in the low-
temperature phase, it would be desired to study such critical behaviors from
percolation perspective.

Some open quesions arise. For instance, what is the upper spatial
dimensionality $d_{u}$ for the supercritical-phase criticality for the UF
model? It is known that the zero-temperature UF model (the uniform tree model)
has $d_{u}=d_{c}=4$ and, at criticality, it has $d_{u}=d_{p}=6$. Thus,
$d_{c}=4$ or $d_{p}=6$ can equally serve as a candidate of $d_{u}$ for the
supercritical UF model. From recent studies for the Fortuin-Kasteleyn
representation of the Ising model [41, 42], we may have the third scenario
that the low-temperature UF model has simultaneously two upper dimensions at
both $d_{c}=4$ and $d_{p}=6$.

###### Acknowledgements.

This work was initiated by private communications with Tyler Helmuth, Roland
Bauerschmidt, and Nicholas Crawford, to whom we are indebted. H.C and Y.D.
have been supported by the National Natural Science Foundation of China (under
Grant No. 12275263), the Innovation Program for Quantum Science and Technology
(under grant No. 2021ZD0301900), Natural Science Foundation of Fujian province
of China (under Grant No. 2023J02032). J.S. was partially supported by Grant
No. PID2020-116567GB-C22 AEI/10.13039/501100011033, and by the Madrid
Government (Comunidad de Madrid-Spain) under the Multiannual Agreement with
UC3M in the line of Excellence of University Professors (EPUC3M23), and in the
context of the V PRICIT (Regional Programme of Research and Technological
Innovation).

## References

  * Stauffer and Aharony [2018] D. Stauffer and A. Aharony,  _Introduction to Percolation Theory_ (CRC press, 2018).
  * Duminil-Copin [2018] H. Duminil-Copin, Sixty years of percolation, in  _Proceedings of the International Congress of Mathematicians (ICM 2018)_ (World Scientific, 2018) pp. 2829–2856, arXiv:1712.04651 .
  * Wu [1982] F. Y. Wu, The Potts model, Rev. Mod. Phys. 54, 235 (1982).
  * Kasteleyn and Fortuin [1969] P. Kasteleyn and C. Fortuin, Phase transitions in lattice systems with random local properties, J. Phys. Soc. Jpn. Suppl. 26, 11 (1969).
  * Fortuin and Kasteleyn [1972] C. M. Fortuin and P. W. Kasteleyn, On the random-cluster model: I. Introduction and relation to other models, Physica 57, 536 (1972).
  * Grimmett [2006] G. Grimmett,  _The Random-Cluster Model_, Vol. 333 (Springer, 2006).
  * Duminil-Copin [2016] H. Duminil-Copin, Random currents expansion of the Ising model, in  _European Congress of Mathematics, Berlin, 18–22 July, 2016_ (EMS, 2016) pp. 869–889, arXiv:1607.06933 .
  * Duminil-Copin and Tassion [2016] H. Duminil-Copin and V. Tassion, A new proof of the sharpness of the phase transition for Bernoulli percolation and the Ising model, Comm. Math. Phys. 343, 725 (2016), arXiv:1502.03051 .
  * Aizenman and Duminil-Copin [2021] M. Aizenman and H. Duminil-Copin, Marginal triviality of the scaling limits of critical 4D Ising and $\lambda\phi_{4}^{4}$ models, Ann. Math. 194, 163 (2021), arXiv:1912.07973 .
  * Duminil-Copin [2022] H. Duminil-Copin, 100 Years of the (Critical) Ising Model on the Hypercubic Lattice, in  _Proc. Int. Cong. Math._, Vol. 1 (International Mathematical Union, 2022) arXiv.2208.00864 .
  * Grimmett [2022] G. R. Grimmett, The work of Hugo Duminil-Copin, (2022), arXiv:2207.02022 .
  * Achlioptas  _et al._ [2009] D. Achlioptas, R. M. D’Souza, and J. Spencer, Explosive percolation in random networks, Science 323, 1453 (2009).
  * Ziff [2009] R. M. Ziff, Explosive growth in biased dynamic percolation on two-dimensional regular lattice networks, Phys. Rev. Lett. 103, 045701 (2009), arXiv:0903.4678 .
  * Grassberger  _et al._ [2011] P. Grassberger, C. Christensen, G. Bizhani, S.-W. Son, and M. Paczuski, Explosive percolation is continuous, but with unusual finite size behavior, Phys. Rev. Lett. 106, 225701 (2011), arXiv:1103.3728 .
  * Riordan and Warnke [2011] O. Riordan and L. Warnke, Explosive percolation is continuous, Science 333, 322 (2011).
  * Li  _et al._ [2023] M. Li, J. Wang, and Y. Deng, Explosive percolation obeys standard finite-size scaling in an event-based ensemble, Phys. Rev. Lett. 130, 147101 (2023), arXiv:2301.09774 .
  * Ziff [2013] R. M. Ziff, Getting the jump on explosive percolation, Science 339, 1159 (2013).
  * Araújo  _et al._ [2014] N. Araújo, P. Grassberger, B. Kahng, K. Schrenk, and R. M. Ziff, Recent advances and open challenges in percolation, Eur. Phys. J. Special Topics 223, 2307 (2014), arXiv:1404.5325 .
  * Saberi [2015] A. A. Saberi, Recent advances in percolation theory and its applications, Phys. Rep. 578, 1 (2015), arXiv:1504.02898 .
  * Fernández  _et al._ [1992] R. Fernández, J. Frölich, and A. D. Sokal,  _Random Walks, Critical Phenomena, and Triviality in Quantum Field Theory_ (Springer, 1992).
  * Stephen [1976] M. Stephen, Percolation problems and the Potts model, Phys. Lett. A 56, 149 (1976).
  * Jacobsen  _et al._ [2005] J. L. Jacobsen, J. Salas, and A. D. Sokal, Spanning forests and the $q$-state Potts model in the limit $q\to 0$, J. Stat. Phys. 119, 1153 (2005), arXiv:cond-mat/0401026 .
  * Bauerschmidt  _et al._ [2021] R. Bauerschmidt, N. Crawford, T. Helmuth, and A. Swan, Random spanning forests and hyperbolic symmetry, Comm. Math. Phys. 381, 1223 (2021), arXiv:1912.04854 .
  * [24] R. Bauerschmidt, N. Crawford, and T. Helmuth, Percolation transition for random forests in $d\geq 3$,  arXiv:2107.01878 .
  * Deng  _et al._ [2007a] Y. Deng, T. M. Garoni, and A. D. Sokal, Ferromagnetic phase transition for the spanning-forest model ($q\to 0$ limit of the Potts model) in three or more dimensions, Phys. Rev. Lett. 98, 030602 (2007a), arXiv:cond-mat/0610193 .
  * Halberstam and Hutchcroft [2023] N. Halberstam and T. Hutchcroft, Uniqueness of the infinite tree in low-dimensional random forests, (2023), arXiv:2302.12224 .
  * Potts and Domb [1952] R. B. Potts and C. Domb, Some generalized order-disorder transformations, Proc. Cambridge Philos. Soc. 48, 106 (1952).
  * Baxter [2016] R. J. Baxter,  _Exactly Solved Models in Statistical Mechanics_ (Elsevier, 2016).
  * Wu [1984] F. Y. Wu, Potts model of magnetism (invited), J. Appl. Phys. 55, 2421 (1984).
  * Caracciolo  _et al._ [2004] S. Caracciolo, J. L. Jacobsen, H. Saleur, A. D. Sokal, and A. Sportiello, Fermionic field theory for trees and forests, Phys. Rev. Lett. 93, 080601 (2004), arXiv:cond-mat/0403271 .
  * Wu [1978] F. Y. Wu, Absence of phase transitions in treelike percolation in two dimensions, Phys. Rev. B 18, 516 (1978).
  * Jacobsen and Saleur [2005] J. L. Jacobsen and H. Saleur, The arboreal gas and the supersphere sigma model, Nucl. Phys. B 716, 439 (2005), arXiv:cond-mat/0502052 .
  * Bedini  _et al._ [2009] A. Bedini, S. Caracciolo, and A. Sportiello, Phase transition in the spanning-hyperforest model on complete hypergraphs, Nucl. Phys. B 822, 493 (2009), arXiv:0906.4503 .
  * Caracciolo  _et al._ [2017] S. Caracciolo, A. D. Sokal, and A. Sportiello, Spanning forests and OSP$(N|2M)$-invariant $\sigma$-models, J. Phys. A 50, 114001 (2017).
  * [35] R. Bauerschmidt and T. Helmuth, Spin systems with hyperbolic symmetry: a survey,  arXiv:2109.02566 .
  * Sweeny [1983] M. Sweeny, Monte Carlo study of weighted percolation clusters relevant to the Potts models, Phys. Rev. B 27, 4445 (1983).
  * Sleator and Endre Tarjan [1983] D. D. Sleator and R. Endre Tarjan, A data structure for dynamic trees, J. Comput. Syst. Sci. 26, 362 (1983).
  * Deng  _et al._ [2007b] Y. Deng, T. M. Garoni, and A. D. Sokal, Critical speeding-up in the local dynamics of the random-cluster model, Phys. Rev. Lett. 98, 230602 (2007b), arXiv:cond-mat/0701113 .
  * Giachetti  _et al._ [2021] G. Giachetti, N. Defenu, S. Ruffo, and A. Trombettoni, Berezinskii-Kosterlitz-Thouless Phase Transitions with Long-Range Couplings, Phys. Rev. Lett. 127, 156801 (2021), arXiv:2104.13217 .
  * Giachetti  _et al._ [2022] G. Giachetti, A. Trombettoni, S. Ruffo, and N. Defenu, Berezinskii-Kosterlitz-Thouless transitions in classical and quantum long-range systems, Phys. Rev. B 106, 014106 (2022), arXiv:2201.03650 .
  * Fang  _et al._ [2022] S. Fang, Z. Zhou, and Y. Deng, Geometric upper critical dimensions of the Ising model, Chin. Phys. Lett. 39, 080502 (2022), arXiv:2201.12954 .
  * Fang  _et al._ [2023] S. Fang, Z. Zhou, and Y. Deng, Geometric scaling behaviors of the Fortuin-Kasteleyn Ising model in high dimensions, Phys. Rev. E 107, 044103 (2023), arXiv:2212.08544 .