{'fraction_non_alphanumeric': 0.02691480314172795, 'fraction_numerical': 0.009505227875331432, 'mean_word_length': 4.6822057987492895, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Drone-based Volume Estimation in Indoor Environments

Volume estimation in large indoor spaces is an important challenge in robotic inspection of industrial warehouses. We propose an approach for volume estimation for autonomous systems using visual features for indoor localization and surface reconstruction from 2D-LiDAR measurements. A Gaussian Process-based model incorporates information collected from measurements given statistical prior information about the terrain, from which the volume estimate is computed. Our algorithm finds feasible trajectories which minimize the uncertainty of the volume estimate. We show results in simulation for the surface reconstruction and volume estimate of topographic data.


INTRODUCTION
Autonomous robotic platforms are increasingly used for data collection, for example in structural inspection Almadhoun et al. (2016), agriculture surveillance, search and rescue, and industrial environments. In industrial warehouses it is common to store raw materials as stockpiles, and determining the current amount of material in stock is of paramount importance for logistics. However, to the best of our knowledge, the task of estimating the volume with an automated robotic platform has not been addressed. On the contrary, current business practice is to take differential measurements of the volume added or removed, which is prone to drift over time, and on periodic inspections from experts, which are costly and inaccurate. A quadcopter equipped with an adequate suite of sensors could be used for this purpose, since it can fly above the pile of bulk material taking advantage of its maneuverability to take measurements from poses otherwise unreachable while avoiding obstacles in cluttered environments. Estimating the volume using autonomous quadcoptor in indoor environment imposes several requirements, e.g.: (1) indoor localization, where the location uncertainty depends on the drone state, (2) accurate measurements in environments with uneven light and dust, (3) an efficient surface reconstruction method, able to cope with large amounts of data, and (4) path planning, due to the limited time budget to fly the drone and the presence of obstacles.
Localizing a mobile platform in a GPS-deprived environment with a known map can be achieved using information This research has been supported by NCCR Automation, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 180545), and by the Swiss Innovation Agency (Innosuisse) via grant number 52481. obtained from dead-reckoning, infrared, radio or soundbased distance measurements, visual information using a motion capture system, or from an onboard camera. In this work localization is inferred from an onboard camera, detecting a set of fixed, previously mapped features, given the lower overall system cost, flexibility and taken into account the accuracy requirements.
Surface reconstruction can be performed from images, with photogrammetry methods such as structure from motion (SfM) (Newcombe et al. (2011)). This methods are however problematic for objects with homogeneous surfaces or improper lighting. An alternative way is via active laser scanners (LiDAR) which project a laser beam and measure the time of flight of the reflected light. This sensors are precise and are not affected by the effect of scale uncertainty present in vision based measurements. Furthermore 2D LiDAR systems are lighter than their 3D counterparts, allowing the use of more maneuverable quadrotors. There are commercial examples of drone-based solutions that use this method in outdoor environments (van Rees (2015)). The reconstruction quality depends to a large degree on the availability and quality of measurements. Classic approaches for quality-driven and automated 3D scanning use volumetric (Khalfaoui et al. (2013)) and Poisson meshbased metric. Algorithms also been proposed for multiview stereo reconstruction (Hepp et al. (2018)), defining heuristics to decide on the utility of the next measurements and optimize set viewpoints based on initial scans. Alternatively, 2D laser scanner has been used for 3D mapping in Park et al. (2018); Singh (2014, 2017), often mounting the 2D scanner on a rotating motor to emulate 3D LiDAR properties.
In this paper, we consider the problem of estimating the volume of material within a given domain using a drone mounted 2D LiDAR unit operating in an indoor environment, leveraging information about the surface. This setting has challenges unique to GPS denied environments, notably the uncertainty in the localization depends on the position of the drone, which must follow trajectories that keep enough features in view to maintain localization accuracy. In Hollinger et al. (2012) the authors consider visual inspection of ship hulls using underwater vehicles. In ? the authors consider inspection of 3D object, and Popović et al. (2020) deals with a similar problem setting, but localization uncertainty is not taken into account. This work differs since it deals with indoor environments and carefully considers the uncertainty of the measurements. Our contributions are threefold: (1) we derive measurement models and uncertainty estimates for the camera-based localization scheme and the LIDAR system used to measure the surface; (2) we propose a scalable methodology for estimating the volume of material based on LIDAR measurements and qualifying the uncertainty of our estimate; (3) we propose a preliminary informative path planning method that greedily minimizes the uncertainty in the volume estimate.

PROBLEM STATEMENT
We consider the problem of estimating the volume of a pile of bulk material inside a region of interest using a quadrotor-based mobile sensor. Let h(x) = h(x 1 , x 2 ) be the true surface function of the height of the pile, defined and the dynamics of the quadcopter are given bẏ where χ full (t) = (p, θ, v, ω, b) : [0, T ] → R 13 is the drone state, consisting of its position p, orientation θ, velocity v and angular velocity ω, as well as the battery state of charge b, u(t) : [0, T ] → R 4 are the rotor voltages, g is a feedback controller that stabilizes the system to a commanded position and yaw setpoint r, andχ full (t) is an estimate of the current state of the drone. In this work we assume that the low-level controller can follow the reference closely, such that if the time derivatives of r(t) up to order four (velocity, acceleration, jerk and snap) are within specified bounds, derived from actuator limits U Bouabdallah and Siegwart (2005).
The quadcopter is equipped with a 2D LiDAR, measuring the radial distance d l from the sensor to the surface of the pile, that is used to estimate the volume of the pile based on a reconstruction of the surface. The drone is also equipped with a camera and computer vision algorithms that allow it to detect and identify features placed in the environment at known locations. These features are used to localize the drone which must keep a minimum number of features within its field of view at all times.
We have now all the ingredients needed to formulate our volume estimation problem. We aim to find a reference trajectory r(t) is the trajectory time and C free is the set of 3D positions and yaw where the drone is allowed to fly and is able to detect enough features to localize itself, that minimizes the volume estimate uncertainty, where V is the variance operator. This is a hard problem as the uncertainty in the volume estimate is a consequence of uncertainty in the LiDAR measurements which in turn depends on the uncertainty in the camera-based localization system. In the subsequent sections we derive detailed measurement models for the localization system (Section 3) and the LIDAR, propose a method for fusing them to estimate the volume (Section 4) and propose a greedy algorithm to approximate a solution of (3) in Section 5.

LOCALIZATION
We are primarily interested in operations in indoor GPS denied environments. Our drone is equipped with an IMU and a front facing camera, which is used to localize the drone based on markers placed in known locations in the environment. Both markers and natural features are detected using computer vision algorithms (??) which also have been demonstrated by Tinamu Labs (Jiang (2022)). We require that enough previously mapped features are visible at any given point of the flight to assure that the state can always be determined, independently from other sensory information. set of N features detected and identified in the image plane of the onboard camera, as exemplified in Figure 2.
We can compute the coordinates in the image frame P im from the global frame P g as a function of the position and orientation of the drone where T g c (χ) is the geometric transform from the global to the camera frame, and C c im the transform from camera frame to image plane given by the camera model. We can solve a nonlinear least squares problem to find an estimate of χχ Linearizing the nonlinear least squares problem around the solution, for normally distributed measurement errors, the estimateχ is also normally distributed with variance where σ is the standard deviation of a measurement in the image plane, and J = ∇P im (χ)| χ=χ is the Jacobian. Assumingχ ≈ χ, we can compute the variance of the state estimate Σ χ as a function of the drone coordinates Σ χ (χ) = V(χ)|χ =χ .
Let q pos (Σ χ (χ)) = 1/ T r(Σ χ ) be a function mapping the covariance matrix to a scalar quality of position fix metric. This metric is used to define the feasible flying domain C free as C free = {χ ∈ C : q pos (Σ χ (χ)) > τ | W} (7) where χ is the drone state, C is the configuration space, τ is the minimum quality of fix threshold, and W the available information about the visual features positions.
The IMU and position measurements derived from the vision system are fused using an extended Kalman filter Einicke and White (1999) to produce an estimate of the full state of the drone. The full estimate is used by the inner loop controller g while the planning algorithm only use the position and yaw localization.

VOLUME ESTIMATION
Now that we know the position and orientation of the drone, as well as their uncertainties we consider the LiDAR model. The sensor rotates in a plane with constant angular velocity, collecting approximately 10K distance samples per second, at discrete angles as illustrated in Figures 1  and 3. We can compute the hit point coordinates of the LiDAR measurements with the geometric transformation z m (χ, α) = T l g (χ, α)d l , where T l g (χ, α) is the geometric transformation converting points in the LiDAR scan line frame to global coordinates, α is the angle of the sensor, d l is the measured distance, and z m ∈ R 3 are the coordinates of the hit point. Propagating the uncertainty of each of the variables we obtain an estimate of the covariance of the measurement z m , where again we use the Jacobian of the geometric transformation where Σ m , Σ χ , Σ α and Σ d l are the covariance matrices of a measurement, the drone position and orientation, the LiDAR angle, and measured distance, respectively. Note that due the intrinsic properties of the sensor, only obstacles in a certain distance range can be detected d l ∈ [d min , d max ]. Figure 3 shows how the LiDAR scan obtains data about the surface. We condense the position uncertainty Σ m , propagating the errors-in-variables to errors in height, given known statistical information about the slope of the terrain where σ t is the standard deviation of a Gaussian distribution fitted to an histogram of the slopes of this terrain, as shown in Figure 4.

Surface model
A common way to store the information about a 3D object is to use voxels, a generalization of the concept of pixel in 3D. Variations include sparse octree representations Laine and Karras (2010). However, if we assume that the surface can be described with a function from x, y to h we can simplify the representation. One option is to use kriging-based methods (Krige (1951)) such as Gaussian process (GP) (Williams and Rasmussen (2006)). This is the method used in Popović et al. (2020) where an informative path planning framework is proposed.
We parameterize the surface using a grid of heights. Each point of the grid is described by a univariate normal distribution where N and M are the dimensions of the grid. Between the grid points, we represent the surface distribution (the surface height is a random variable due to uncertainty) using a Gaussian Process model where the height grid provides inducing points.
We use a specific kernel tailored to the physics of our surface. Due to their physical properties, the materials piled up in the region of interest have a volumetric organization which can be described statistically. The Matérn Kernel is a good choice to estimate the correlation of the heights of two points separated by a distance s where l is the lengthscale, ν is a positive parameter controlling the smoothness of the function, and K ν and Γ ν are Bessel and Gamma functions, respectively. We describe the surface as a sparse Gaussian Process with fixed inducing points X which are obtained from the height grid, f (x) ∼ GP(0, k(X, X )). At an arbitrary set of points X * ∈ D, we can predict the expected value and variance of the height of the surface using the equations where E is the expected value, V is the variance, K X * X is the covariance between the points X and X * computed with the kernel (12), where s is the pairwise distance between the points (Williams and Rasmussen (2006)), Z is the vector of the z coordinates of the inducing points, and Σ the vector of their uncertainties. For a certain kernel lengthscale, the vector Θ . = [Z, Σ] contains all the information needed to make predictions. For computational reasons, the kernel is made sparse by setting to zero the correlation between points outside of a ball with radius γl. With this model for the surface, we can readily compute the volume, its expected value and variance using the equations where A is the area of the base of each square cuboid.

Update
When a new measurement is obtained we first propagate it to the grid points. This is accomplished using the kernel Σ m (z) + σ 2 t (e s/l − 1), where σ t is the standard deviation of the normal distribution that approximates the slope distribution, s is the distance between the measurement point and the grid coordinates, and l is the lenghtscale used in (12). We then use a Kalman filter to update the parameters of the height model (11) µ where µ k−1 , µ k , σ k−1 , and σ k are the mean and standard deviation of each point of the grid, before and after the update step, z k and Σ m (z) k are the measurement value and covariance from (9). K k is the Kalman gain.

PLANNING
To estimate the volume of the surface, the robot needs to collect information of the whole surface, moving on a path suitable for this purpose. Two main approaches exist in path planning, that achieve this goal. Coverage path planning Galceran and Carreras (2013) is a method to design paths that visit all points of interest while avoiding obstacles. Typical cost functions minimize the length of the path, whereas with informative path planning Stache et al. (2021) the objective is to maximize the amount of information in a feasible path. An overview of path planning algorithms can be found in LaValle (2006).
For the purposes of planning we (1) assume a constant nominal surface height h 0 , which allows us to compute d l ; (2) decouple the motion of the drone and the rotating LiDAR scan, since the latter is at least one order of magnitude faster; (3) neglect the possible effect of shadows, as defined in Figure 3; (4) we find a discrete set of waypoints instead of a continuous time trajectory. We use a greedy algorithm for planning, where the next point to be picked is the one with that minimizes the volume estimate uncertainty where σ V is computed with (14c), r is the reference, R is the radius of a ball where the next step can lie. The simulations of the next section consider only the  x and y components of r, and fix z and the yaw. The greedy algorithm is suboptimal but fast to evaluate, and allows us to test the surface reconstruction method. Future work will focus on implementing more advanced planning algorithms.

SIMULATIONS
We demonstrate our approach for volume estimation on a topographic map of the Alps mountain range, rich on features, which we scale down by a factor of 1000x. We show simulation results for two trajectories, where the first one is a fixed, manually created square wave pattern, and the second is the trajectory resulting from applying the path planning algorithm from Section 5, in the variables r x and r y , while keeping r z and r yaw fixed. The disposition of visually identifiable features in the example used in this simulations, shown in Figure 1, creates an uneven uncertainty of position map, as shown in the right plot of Figure 5, where darker colors represent lower uncertainty. Figure 6 shows the path of the drone and the altitude uncertainty map, and Figure 7 the surface reconstruction map for the fixed square wave pattern. Traversing the whole region results in reduced uncertainty of surface reconstruction. It is in general not trivial to manually design paths that avoid constraints or regions with insuficient localization quality. As for the greedy algorithm and its resulting trajectory, we show in Figure 8 the path of the drone and the altitude uncertainty map, and in Figure  9 the surface reconstruction map. The simulation results show that a feasible reference trajectory is found that visits most of the region of interest through regions with a high quality of position fix and drives the uncertainty of the volume down to σ V /µ V = 2.26% and a relative error of 2.53% for the greedy algorithm, comparable with σ V /µ V = 2.42% with a relative error of 2.30% for the square wave pattern. We speculate that the two methods perform similarly in this example due to the greedy nature of the planning algorithm used, which optimizes only for the next step ahead, thus finding only an approximate solution for the optimization problem defined in (3). Also note that in this preliminary result the yaw and altitude are kept fixed, so the planner has less degrees of freedom to explore. These factors taken together make it difficult to improve on the benchmark square wave pattern. Figure  10 shows the evolution of the volume estimate as well as the uncertainty as a function of the number of samples collected. The two paths have the same length and number of samples, and the evolution of the volume and its uncertainty is similar for both paths. Furthermore, the final reconstructions from Figures 7 and 9 approximate the ground truth shown in Figure 5 adequately. In summary the simple greedy planner finds a path that is by all metrics Ground truth Estimate µ Estimate µ ± 2σ Estimate µ ± 1σ Fig. 10. Evolution of the volume estimate and its uncertainty along the trajectory. The left plot shows the result for the square wave pattern and the right plot for the greedy algorithm.
similar to the manually designed square wave pattern, allowing the automation of the task of designing paths for experimental volume estimation campaigns.

CONCLUSION
We develop and implement a framework for volume and uncertainty estimation for an autonomous robotic platform collecting experimental data with a LiDAR scanner under variable position uncertainty. The surface reconstruction and volume estimation are validated in simulation with a feature rich surface, for two trajectories generated manually and with a simple greedy algorithm. Future work will focus on the implementation of an improved planning algorithm employing a longer prediction horizon, validation in more complex maps with intricate constraints as well as experimental validation of the method.