Emergency Landing Spot Detection Algorithm for Unmanned Aerial Vehicles

: The use and research of Unmanned Aerial Vehicle (UAV) have been increasing over the years due to the applicability in several operations such as search and rescue, delivery, surveillance, and others. Considering the increased presence of these vehicles in the airspace, it becomes necessary to reﬂect on the safety issues or failures that the UAVs may have and the appropriate action. Moreover, in many missions, the vehicle will not return to its original location. If it fails to arrive at the landing spot, it needs to have the onboard capability to estimate the best area to safely land. This paper addresses the scenario of detecting a safe landing spot during operation. The algorithm classiﬁes the incoming Light Detection and Ranging (LiDAR) data and store the location of suitable areas. The developed method analyses geometric features on point cloud data and detects potential right spots. The algorithm uses the Principal Component Analysis (PCA) to ﬁnd planes in point cloud clusters. The areas that have a slope less than a threshold are considered potential landing spots. These spots are evaluated regarding ground and vehicle conditions such as the distance to the UAV, the presence of obstacles, the area’s roughness, and the spot’s slope. Finally, the output of the algorithm is the optimum spot to land and can vary during operation. The proposed approach evaluates the algorithm in simulated scenarios and an experimental dataset presenting suitability to be applied in real-time operations.


Introduction
Presently, Unmanned Aerial Vehicles (UAVs) are gaining more interest in the scientific and industrial research community due to their autonomy, maneuverability, and payload capacity, making them suitable robots to perform real-world tasks in different scenarios. There have been several research topics related to hardware development, human-system interaction, obstacle detection, and collision avoidance [1]. The UAVs are designed to be remotely controlled by a human operator or execute a mission autonomously [2]. In the latter case, the degree of autonomy and the mission they can achieve depends on the sensors used. Concerning the UAV classification, there are several ways in which UAVs can be defined: aerodynamics, landing, weight, and range [3]. Considering the state of the art of UAVs, commercial or research, there is a large spectrum of applications, such as search-and-rescue operations [4,5], delivery, surveillance, inspection, and interaction with the environment [6,7].
Given the application scenarios, there are missions where the UAV must fly in civilian airspace; i.e., they must fly over populated areas. However, they are susceptible to external disturbance or electromechanical malfunction. Different failure scenarios could impact the operation, thus leading to an emergency landing: • Global Positioning System (GPS) failures. In general, UAVs use GPS messages to navigate. Although several sensors could aid the navigation, the vehicle may need to land in the occurrence of loss of GPS signal. • Loss of communication. In the event of loss of communication between the UAV and a base station, one possible action is to perform an emergency landing. • Battery failure. If a battery failure is detected, the vehicle may not be able to continue its operation, and as a result, an emergency landing is necessary. • Software and hardware errors. A UAV can experience a mechanical fault during a mission, like a broken propeller or a fail in the motor, or even a software issue that could require that the vehicle must perform an emergency landing. • Environment factors. Bad weather conditions such as strong winds and rain make the vehicle unable to carry out the mission, forcing it to land.
In these situation, UAVs must safely land to reduce damage to themselves and avoid causing any injury to humans. Consequently, detecting a reliable landing spot is essential to safe operation. In order to estimate a landing spot, a set of conditions must be evaluated when analyzing the sensor data. These conditions are typically restraints on the landing surface. The reliability of a safe landing site depends on several factors, such as the aircraft's distance to the landing site and the ground conditions, such as (a) the slope of the plane: the landing spot must have a slope smaller than a threshold so that the vehicle does not slide when landing; (b) the roughness of the area: there are cases when the vehicle will land on a surface with vegetation or other small obstacles-in this case, these obstacles' height must not be larger than a maximum value predefined before the mission, avoiding contact with the propellers; (c) the size of the spot: the landing spot must be large enough for the UAV; (d) the presence of obstacles: the approximation to landing spot depends on the area's obstacles, namely building, vegetation, or humans. Therefore, the presence of obstacles must be taken into account when evaluating the spot. Another rule also imposed on the algorithm is the distance to the landing spot. The UAV must be able to reach the desired spot with the remaining battery power. Since one factor is the aircraft's distance to the landing site, the landing spot's suitability varies during the mission. Figure 1 illustrates the conceptual approach for emergency landing spot detection with a UAV in real-time. Considering the robustness problems that a UAV faces during an autonomous mission, this paper addresses developing an algorithm for detecting, storing, and selecting real-time emergency landing spots based on LiDAR data. The main focus is to develop an algorithm that will continuously evaluate the terrain based on LiDAR data and establish a list of possible landing spots. Additionally, a continuous evaluation of the detected landing spots is required given that the best spot can vary during the UAV operation. The developed algorithm accumulates the 3D point cloud provided by the sensor and then generates downsampled point clouds. Following that, a plane is detected by applying the PCA algorithm in a spherical neighborhood. The potential spots are subsequently evaluated regarding different parameters, and finally, the resulting score is assigned to the spot. This paper is an extended version of work published in [8,9]. Compared to the previous papers [8,9], this work presents a more extensive review of the related works. Furthermore, the proposed approach was evaluated with field tests during a 3D reconstruction survey mission that took place in the Monastery of Tibães, Braga, Portugal. The results were also improved with the evaluation of spots' score for each dataset. This work intends to contribute by applying a Voxel Grid Filter to the point clouds provided by a LiDAR. The application of this filter speeds up the plane detection step by reducing the size of the point cloud and removing redundant points. Moreover, it is shown that the algorithm can obtain a list of safe landing spots without segmenting the entire point cloud by analyzing spherical regions within it.
The paper outline is as follows: Section 2 presents a preliminary study of the related works regarding emergency landing and landing detection given LiDAR and camera data. Then, Section 3 shows the high-level hardware and software architecture, detailing their components. Section 4 introduces the developed algorithm and its requirements. In addition, the algorithm sequence of data processing blocks is displayed and explained. In Section 5, the results obtained are presented, as well as the analyses and performance of the algorithm. Finally, in Section 7, we discuss the conclusions of the method for detecting safe landing spots. Here we also present suggestions or directions for further work.

Related Work
The main problem to be addressed is the integration in a UAV to detect safe landing zones. Thus, it is necessary to process different types of data from several sensors so that it is possible to obtain information about the surface where we intended to land. Each sensor onboard the UAV provides specific information regarding terrain and obstacles and therefore has advantages and disadvantages that must be evaluated. Consequently, there are multiple literature approaches regarding safe landing, either for emergencies or for other purposes.

LiDAR Based Detection Systems
Among the primary sensors used to solve the problem of detecting emergency landing spots is LiDAR. Generally, the point clouds are processed, determining terrain features such as slope and roughness. In addition, LiDAR data present usefulness in other tasks, for instance, hazard detection and avoidance. Typically, a geometric approach is a common approach concerning LiDAR data.
In 2002, in order to perceive obstacles and land spacecraft, a 2.5D grid map technique was used by Johnson et al. [10]. Their work applied a Least Mean Square (LMS) algorithm to estimate a plane in the grid map, computing the incidence angle and the roughness to generate a landing cost map. Another geometric approach was presented by Whalley et al. [11,12] in 2009. In the author's strategy, each LiDAR scan is mapped to a grid in which a sliding window moves along to calculate the terrain restraints. Finally, the optimal spot is chosen after the entire grid is covered.
Regarding unknown environments, Chamberlain et al. [13], in 2011, presented a proposal that permits a full-scale unmanned helicopter to land without human control or input. Their methodology utilizes a 3D scanning LiDAR in two modes. A forward scan is performed to detect obstacles, and a downward scan is subsequently carried out to map the environment. A 3D virtual model of the helicopter is placed on each cell of the resulting map. An assessment of the area is executed considering the skid contact, wind direction, and presence of obstacles. This work is expanded in [14][15][16] by presenting experimental results in urban and natural environments. In the evaluation step, a 2D Delaunay triangulation [17] of the potential landing sites data is generated in order to establish the intersections with the landing skids and the roll and pitch of the helicopter. Ultimately, the volume between the triangulation and the 3D model is used to estimate bad contact.
One of the drawbacks of the geometric approach is the presence of low vegetation in the landing spot, as it makes the spot extremely rough. In order to overcome this hindrance, Maturana and Scherer [18] proposed a 3D Convolutional Neural Network [19] algorithm. The strategy is to generate a volumetric density map given the globally registered LiDAR point clouds. The volumetric map is then divided into sub-volumes that are the input of the neural network. The algorithm was able to detect small obstacles using synthetic datasets and simulation.
In 2017, Lorenzo et al. [20] proposed a landing site detection algorithm on many core systems by using parallel processing. An octree is built in order to increase process speed during a neighborhood search step. Then, the normal of all points are computed in parallel, and geometric restraints are used to assess the spot.
Recently, the authors of [21] proposed an approach to select landing zones based on a region growing algorithm. The initial strategy is to find flat regions from LiDAR point clouds and then apply a Progressive Sample Consensus (PROSAC) algorithm to fit a plane. The detected planes are assessed conforming slope, roughness and maximum difference of height. However, the proposed method was only evaluated in simulation experiments and the PROSAC algorithm can be time consuming.

Vision-Based Detection Systems
The vision-based system has been the most popular approach to detect and evaluate a landing spot. There are multiple strategies to land a UAV using vision-based algorithms. Cameras have the advantage of being inexpensive and lighter compared to LiDARs. Bosch et al. [22] developed an algorithm to autonomously detect landing sites with monocular images. First, a cluster of potential points in a plane is selected by applying a homography estimation process. Then, the algorithm distinguishes planar from non-planar surfaces. This information is stored in the stochastic 2D grid cell. However, this proposal stores the probability of the region being flat and hence cannot be used for other applications, such as obstacle avoidance.
Several works have used machine vision to detect landing spots [23][24][25]. The authors proposed a technique that generates two binary images by applying a Canny Edge Detector and a line-expansion algorithm. Given the fusion of the two images, the preliminary map is built to label safe and unsafe areas. Then, the surface is classified, and fuzzy logic is used to select the safest spot. The work was later expanded in 2015 by Warren et al. [26]. A 3D reconstruction using Structure-from-Motion is carried out to analyze the surface.
Eendebak et al. [27] presented a technique for emergency landing in real-time operations given camera movement. Their work is based on a Background Estimation on stabilized video in which the camera frames are compared to distinguish moving objects and structures. Then, a distance map to all detected obstacles is generated. The maximum distance on the map is calculated and chosen as the best landing area.
Forster et al. [28] proposed an elevation map-based technique to landing spot detection for micro-aerial vehicles using a monocular camera. A depth map is generated given the camera images and the vehicle's pose. Furthermore, an elevation map is also built and consequently updated according to the depth maps. During the algorithm selection stage, a flat surface with a radius depending on the vehicle size is considered a landing spot.
Another proposal to detect landing spots in unknown environments in real-time is presented by Hinzmann et al. [29]. A segmentation step is applied to the camera images, classifying the regions as "grass" or "not grass." Next, 3D reconstruction and elevation maps are achieved given the results of the segmentation stage. Finally, the area is assessed considering the ground features, such as slope and roughness.
In 2019, Kaljahi et al. [30] proposed a system for detecting landing sites using Gabor Transform and Markov chain codes. The authors proposed a method that obtains flat surfaces from images using the Gabor Transform. Then, histogram operations are applied in order to determine the pixels that contribute to the highest peak. The next step consists of using Markov Chain Codes to estimate and group potential pixels as candidate regions. For each candidate region, Chi square distance is computed and compared to a reference to define the similarity between the regions, resulting in the final safe landing zone.
In 2020, Bektash et al. [31] presented a machine-vision-based algorithm to recognize landing areas. The model chosen for the network was the convolutional neural network that receives split images provided by a camera as an input. However, manual vision inspection of the camera frames is used as criteria for classifying the landing site.

Other Approaches
In addition to techniques that use one primary sensor to solve landing detection, other approaches combine multiple sources of information. Serrano [32] presented a probabilistic framework using Radio Detection Furthermore, Ranging (RADAR), LiDAR, and camera data to improve the robustness of the selection of landing sites stage in space operation. By applying a Least Median of Squares regression to fit a plane using the RADAR and LiDAR data, the authors calculated the terrain slope and roughness. Furthermore, an edge detection algorithm is applied in the camera images to distinguish craters and rocks. Finally, a Bayesian Network [33] is used to evaluate the area safety considering the information acquired in the previous steps.
Using a similar set of sensors, Howard and Seraji [34] built three hazard maps from RADAR, LiDAR, and camera data. These maps were then used to obtain measurements and the landing site features in which each map is labeled with a confidence variable. Finally, the maps were fused using fuzzy logic presenting the landing site safety.

Overall Discussion
The decision between LiDARs or cameras to serve as the primary sensor to detect emergency landing spots depends on the vehicle's characteristics and the mission in which it will be applied. Generally, the feasibility of vision-based techniques depends on visibility restrictions. However, this drawback can be easily overcome by LiDARs sensors. Nevertheless, the processing power, weight, and price of LiDARs are more demanding for UAVs than cameras.
The use of LiDAR is justified due to the greater robustness to situations of luminosity variations [14,15,21], already obtaining depth information at the instant of the SCAN, not requiring a stereo baseline or depth estimation by sequence of images, and being robust in flights nocturnal.
In terms of processing time of LiDAR-based systems, the current works indicate the importance of knowing the sensor data structure to reduce the number of operations and, thereby, to improve data access efficiency and execution time of the algorithm. Hence, the point clouds are generally spatially structured. It is worth noting that the approach chosen for structuring can cause loss of spatial information, such as 2.5D grid and plane-projected images. Consequently, additional statistical data must be stored in these scenarios.
Regarding the selection of landing spots, the exposed works revealed the importance of classifying the area in a probabilistic manner. These probabilistic values should consider the possibility of landing, presence of obstructions, battery power, ground restraints. Consequently, the levels vary during the operation.

System Design
This section describes the hardware and software architectures that were considered during the implementation of the proposed algorithm. The landing spot detection algorithm was developed in Robotic Operating System (ROS) [35], due to providing a modular structure able to integrate available ROS packages straightforwardly with other software modules already available in the UAV.
The proposed framework processes the 3D point clouds provided by a LiDAR to detect landing zones in real time. Furthermore, another layer of software is responsible for analyzing the UAV status, such as remaining battery power, the status of hardware and software components. Therefore, in a failure situation in which a forced landing is necessary, this layer will search for the current best landing spot detected by the proposed algorithm.

Hardware
The landing spot detection algorithm is designed to be executed simultaneously with other UAV tasks during operation. The VLP-16 LiDAR provides 3D point clouds that are processed using an onboard computer. Considering that modern LiDARs produce hundreds of thousands of points in each scan, a high transmission rate between the sensor and the computer is necessary (Figure 2 in blue).
The point clouds are given in the LiDAR reference frame. Consequently, knowing the UAV pose is essential to transform the data to a fixed frame. There are several approaches to estimate the vehicle pose. In this case, the method is to fuse the Global Navigation Satellite System (GNSS) and Inertial Measurement Unit (IMU) measurements. Another critical issue is to correlate in real-time the LiDAR scan with all the UAV onboard sensors. Therefore, to ensure synchronization, the timestamp data from the GNSS is provided to both systems, LiDAR and PC, with Chrony [36] service. E th e r n e t Figure 2. High-level hardware architecture (adapted from Refs. [37,38]). Figure 3 presents the high-level software architecture to detect a potential landing spot given a point cloud generated by a LiDAR sensor. All the pipeline elements were developed within the ROS framework, except for the sensor's input data. The software is divided into several data processing-blocks: Generally, the algorithms consist of estimating the parameters from the plane equation in a limited region of the original point cloud. Moreover, the algorithms have to fit a plane in the presence of outliers, i.e., points that do not fit the plane model. Besides the landing spot detection for aerial vehicles, estimating the ground conditions is also useful for detecting clear paths and making further processing less complex. • Spot Evaluation: the segmented point cloud containing the detected plane is then evaluated to classify the spot's reliability. The roughness of the landing zone can be assessed by computing the standard deviation in the z-axis. The higher the standard deviation, the rougher the area. In addition, a high value can also indicate the presence of obstacles, such as trees or buildings. The spots are also evaluated regarding their slope as the landing zone cannot be steep enough to destabilize the vehicle when landed. Furthermore, the size of the area and the distance to the UAV can also be used to evaluate the spot. Each evaluation factor has different importance according to the environment. Thus, it becomes interesting to assign different weights to these factors. • Spot Register: the assessed spots are then registered as a landing spot. However, the suitability of a landing point to be the optimal choice varies during operation. In this way, the points registered have to be periodically reassessed.

Emergency Landing Spot Detection Algorithm
The developed algorithm procedure is divided into several steps: the frame transformation of the input data, the downsampling of the point cloud, the spatial structuring of the data, detection of planes, filtering of potential candidates, and classification of detected spots. In addition, there are some rules created to obtain the desired objective. Figure 4 shows the flowchart of the developed software.

Frames Transformation
The pose of a robot is generally given in the local navigation coordinates system. Consequently, the landing spot must be in the same frame. However, LiDARs report data in their own reference frame. Then, the first step of our proposal is to transform the point cloud from the LiDAR reference frame to the local navigation reference frame. This procedure can be fulfilled by using transformation matrices. Therefore, the relation between the frames is defined by Equation (1): Considering a 3D point p L (t) in the LiDAR reference frame, at instant t, the transformation matrix given by T b L transforms the point to the body reference frame. At this step, it is considered that the LiDAR is fixed to the robot, hence the T b L matrix not varying over time. Finally, the T n b (t) matrix multiplied by R b L p L (t) results in the data expressed in the local navigation reference frame.
The algorithm does not apply the landing zone detection step for each LiDAR scan. This is due to the following reasons: the vehicle may not have traveled a sufficient distance to need a new spot other than the take off point; in addition, the sensor may not return sufficient points of the environment. Therefore, the point cloud registered in the local frame is stored until it reaches the conditions to start the next step. The defined conditions are: • Traveled distance: if the robot has traveled a long distance in relation to the last landing spot, it becomes necessary to find a new location, since, in an emergency scenario, the vehicle may not be able to reach the landing zone. At the moment, this step is done considering only the vehicle pose. However, in future projects, it is worth considering the vehicle velocity.

Point Cloud Downsampling
Accumulating the data until one of the necessary conditions are reached will increase the processing effort. In order to obtain better performance in terms of execution time and memory consumption, it becomes necessary to perform the downsampling of the point cloud. Therefore, a Voxel Grid filter [39] is applied. The filter takes a spatial average of the points in the cloud. A set of 3D volumetric pixels (voxel) grid with size v f ilter is generated over the cloud, and the points are approximated with their centroid. Figure 5 shows the result of the downsampling. At first analysis, the point clouds displayed in Figure 5a,b are almost identical. However, the downsample technique created a point cloud three times smaller in size. This method allows decreasing the number of points of the cloud and gives a point cloud with approximately constant density. Finally, the algorithm resets the original point cloud in order to free the memory.

Data Structuring and Neighbor Search
The next step is to determine the neighborhood of a point to perform the plane identification method. The process of detecting a plane in a point cloud is the most timeconsuming process in the algorithm. Due to the large number of data, this step needs to be optimized in order to obtain a real-time analysis. For this reason, the point cloud is spatially structured using octree [40]. Each internal node of the octree is subdivided into eight octants. By using a tree structure like an octree, the execution time of a search algorithm is considerably reduced [41]. There are several approaches to identify the neighborhood of a point. In this case, the algorithm finds the spherical neighborhood of a randomly chosen point. Considering a point p(x p , y p , z p ) in the point cloud, the neighborhood N(p) of this point with radius r is determined by: where q(x q , y q , z q ) is any point in the cloud. The minimum radius (r min ) of a plane is defined by the user as depends on sensor and environment characteristics. By applying Equation (2), all the points placed inside the sphere are considered as part of the neighborhood.

Plane Detection
After the previous step, it is possible to calculate a planar surface given the neighborhood points. This is done satisfactorily by using the PCA algorithm. PCA applies an orthogonal transformation to map the data to a set of values called principal components that correspond to the eigenvectors of the covariance matrix.
Considering the plane equation given by: The normal vector is represented by: The eigenvectors determined with PCA serve as the three axes of the plane while the eigenvalues indicate the square sum of points deviations along the corresponding axis. Therefore, the eigenvector with the smallest eigenvalue represents the normal vector given by Equation (4), and the points are bounded by the other two axes. In this perspective, the slope of the plane is examined using the normal vector. The slope is the angle between the normal vector and the vertical vectorẑ = 0 0 1 T and is computed using Equation (5): Hence, a first evaluation can be done using the plane slope. If the resulting value is greater than the maximum slope (θ max ) permitted for the robot, the plane is rejected. Otherwise, the neighborhood radius used in the previous section is increased and the process is repeated until the maximum radius (r max ) is reached. By doing this procedure, the method tries to find regions with different sizes that can be considered a landing spot. Algorithm 1 describes the neighborhood and plane identification steps. The idea of the method is to detect planes for n points random search points with increasing radius. The PCA step is only performed if the number of points within the neighborhood sphere is greater than n min . The slope angle is immediately computed after detecting the plane and, if it is greater than θ max , the plane is rejected and the algorithm restarts the process for a new search point. In addition, when we increase the radius of the sphere and a previously detected plane has a slope greater than the threshold, only the point cloud corresponding to the previous slope and the respective radius of the sphere is sent to the evaluation step.

Algorithm 1 Algorithm for neighborhood and plane identification steps.
Input: pointcloud downsampled Output: point cloud cluster 1: for j = 1 to n points do 2: Select random point in cloud as the search point Declare radius and vectors of cloud indices 3: float r min , vector radiusInx While the plane is accepted, do radius search 4: while isPlane = true do 5: Start radius search; 6: if (radiusInx ≥ n min ) then 7: Get points inside the neighborhood; 8: Start PCA; 9: Compute plane parameters; 10: Compute plane slope; 11: if slope ≤ θ max then 12: Increase radius; In summary, only the planes that have a slope less than or equal to θ max are sent to the register stage. Table 1 describes the parameters used in the algorithm.

Parameters
Description n points Number of random points chosen from the cloud cluster n min Minimum points used to fit a plane r max Maximum radius considered for a plane r min Minimum radius considered for a plane θ max Maximum slope accepted for the drone

Registration and Classification
Given the algorithm presented in the previous section, the next step consists in evaluating the detected planes. In this context, many factors are considered to decide the best landing spot. Initially, the point cloud centroid is computed and is considered the spot center. After this, each factor is computed individually. A score (g n ) from 0 to 20 is computed for each parameter. It is worth noting that the degree of importance for each parameter depends on the operation. Consequently, each grade has a different weight previously defined by the user. Finally, each spot is classified regarding the following equation: Table 2 shows the parameters computed to classify the plane and their respective weights. In this stage, the planes that have standard deviation in the z-axis greater than the maximum standard deviation σ max are immediately rejected. Distance from the spot to the vehicle Using these parameters, the algorithm evaluates the landing spot in terms of terrain roughness, vehicle stability when landed, obstacle clearance of a location and distance to the UAV. Ultimately, the spots are stored and sorted from highest rate to smallest. In general, a spot quality depends on the vehicle trajectory. In this perspective, the stored spots are re-evaluated periodically.

Algorithm Output
The algorithm returns the current best landing spot in the local navigation frame and the score obtained in the classification step. Concurrently, the algorithm keeps updating the spots score as the vehicle continues the mission by computing the new distance to the UAV and recalculating the score. Therefore, the landing spot detection and selection algorithm is always being performed for each scan provided by the LiDAR.

Simulation Setup
Before proceeding to field tests, one requirement was to evaluate the algorithm's robustness and quality under a simulation environment. Therefore, the simulation environment that was considered was the Modular Open Robots Simulation Engine (MORSE) [42][43][44]. MORSE is an open-source simulator with several features, such as ROS support and virtual sensors, including a generic 3D LiDAR that performs a 180-degree scan and publishes as a point cloud message. Moreover, it is developed in Python and uses the Blender Game Engine to model and render the simulation.
The simulation was performed on a computer with an Intel Core i7-6500U CPU @ 2.50 GHz with 4 cores and 8 GB RAM, with a Linux 4.15.0-50-generic Ubuntu. Moreover, it used the ROS Kinetic Kame distribution (http://wiki.ros.org/kinetic accessed on 14 May 2021).

Environment I
The first scenario is a standard MORSE environment, as shown in Figure 6. This environment is relatively flat with only a few obstructions, namely the buildings. Thus, this simulation's objective was to analyze the algorithm's behavior in terms of time consumption, specifically in the downsampling and plane detection steps. The algorithm waits for the vehicle to travel the minimum distance of tv d = 10 m or the accumulated point cloud reaches its maximum size of c size = 10 6 points to start the downsampling. Each voxel of the Voxel Grid filter was defined to v f ilter = 5 cm. The environment was simulated for several values of maximum sphere radius and search points in the plane detection stage. Finally, potential spots slope greater than θ max = 15 • or standard deviation along the z-axis greater than σ max = 0.20 m were rejected. Figure 7 shows a comparison of the execution time between the downsampling step and the PCA stage for several simulations with increased number of search points. The downsampling technique is relatively fast, taking about 20 ms. Moreover, increasing the number of search points implies more iterations for the PCA algorithm. Considering that the PCA is applied in a downsampled point cloud, downsampling the point cloud is justified because it limits the sample size.  Regarding the number of detected spots, it is not possible to obtain a precise analysis of the performance of the algorithm in this simulation. This is due to the fact that the environment has many regions for landing. Table 3 summarizes the results for the simulation. It also presents the simulation results with fixed search points but increasing maximum radius. The PCA run-time increases rapidly to a small increase in radius. This implies that the value chosen for the neighborhood radius has greater weight in this step. A second environment (displayed in Figure 8) was developed in Blender with the intention of evaluating the performance of the algorithm in a bad scenario with few landing spots. The set consists of six rectangular surfaces that represent the desired landing spots surrounded by a mountain like structures. Considering the environment and the dimension of each surface, the simulations were realized for four different search points. In this simulation, the standard deviation was decreased to σ max = 0.15 m.   Figure 9a,b, the run-time of the PCA algorithm was similar to the downsampling stage. The reason for this is that the planes are immediately rejected. Since the environment has a high variance in z, the angle between the normal vector calculated by the PCA and the vertical vector (ẑ = [0, 0, 1] T ) is generally greater than the maximum allowed. As a result, the algorithm does not perform many iterations for the same search point.
Considering the detected spots, Figure 10 shows the results obtained from the simulation. The landing surfaces are the red rectangles, and the black crosses are the detected spots. However, the algorithm was not able to find a landing spot in all planes. In this scenario, one drawback of the algorithm is that if the random search point falls near a high deviation zone like the surface edges, and the algorithm may discard the plane or detect a spot outside the surface, such as a result for 100 search points (Figure 10d).  For a better understanding of how the best landing spot varies during the simulation, Figure 11 displays the graphs of the best spot score and the score of each detected area for the case with 50 search points. Furthermore, when a spot is detected or is considered the best spot, it is marked on the graph. As the vehicle moves away from the spot, the grade associated with it decreases. If another spot has a higher score, then it becomes the new best spot. Table 4 summarizes the results for the second simulation. The PCA mean run-time is smaller than the downsampling mean time in the first test. The difference between rejected and accepted plans increased considerably as expected.
After this procedure, new simulations with different values for the voxel size were carried out to comprehend the influence of the voxel size on the other steps of the algorithm. Table 5 presents the results obtained. It was chosen as the number of search points n points = 20. As the value of the voxel increases, the size of the new point cloud decreases considerably. This result shows that it is necessary to adjust the variable correctly, as it can lead to a great loss of information about the environment. Furtherore, after 1 m of voxel size, the algorithm stopped detecting landing points. In terms of the PCA performance, the increase in size did not reflect a major change in the execution time. However, the number of detected planes increased, indicating that the downsampled point cloud did not present as many obstructions as it should.       The multirotor UAV STORK (shown in Figure 12) is an autonomous aerial vehicle with ss designed to achieve real-time data acquisition and processing efficiently. It has been used for several applications, including power lines inspection and mapping. The UAV was built, allowing a modular payload assembly, i.e., the sensor can be replaced by others that supply different types of data and using the same frame [45]. The UAV can operate in both manual and autonomous modes. Currently, the UAV Stork's low-level control is accomplished by a customized autopilot (INESC TEC Autopilot), while an onboard computer ODROID-XU4 performs the high-level control. In terms of navigation, the UAV possesses two IMUs and two GNSS receivers. In addition to the low-cost sensors, STORK also has the high-performance IMU STIM300 and the single-band GNSS receiver ComNav K501G that supports Real-Time Kinematic (RTK). Generally, this setup of sensors succeeds in satisfying the requirements for many applications.
Regarding the perception assembly, the sensors to extract information about the surrounding area, the UAV STORK, has two visual cameras (Teledyne Dalsa G3-GC10-C2050 and FLIR Point-Grey CM3-U3-13S2C-CS) and a 3D spinning LiDAR sensor (Velodyne VLP-16). These sensors provide data used as input for processing algorithms in several modules, such as navigation and 3D reconstruction. For instance, the Velodyne VLP-16 data is the input for the emergency landing spot detection algorithm.

Dataset
The experimental dataset corresponds to the staircase of the Monastery of Tibães, Braga, Portugal. The stairway consists of several water fountains that are intertwined by stairs. This area has numerous trees that hinder the navigation and localization of the UAV and are considered obstacles for detecting landing points.
Two experiments were performed to assemble this dataset. First, the vehicle started by mapping a field of vegetation near the staircase (see Figure 13a). After that, the UAV traveled through the different levels of the stairway. The second experiment was done at the monastery building, specifically on the fence of the monastery (Figure 13b). This region consists mainly of buildings and vegetation. It is an open area with more landing spots.  In this experiment, the landing spot detection algorithm was evaluated with more restrictive parameters. Planes with standard deviation along the z-axis greater than σ max = 5 cm were rejected. In addition, the maximum allowed radius for a spot was r max = 4 m, and the maximum point cloud size for starting the downsampling algorithm was 30,000 points. In addition, it was chosen as the number of search points n points = 5.

Dataset Results
For a better understanding of the algorithm procedure, in Figures 14 and 15, we identified the detected spots along with the camera image in the current timestamp of the operation. The red lines represent the center of the spot found during the plane detection step. The small red lines are detected spots that were not chosen as the best landing zones. It is worth noting that some of the detected landing spots are part of a larger spot. Besides the fact that the search points are chosen at random, another reason that justifies this behavior is the need for the maximum radius not to be too large so that the plan detection algorithm does not consume much time, as previously demonstrated by the simulations.
Moreover, the best spot grade has been stored and shown in Figure 16. The result obtained illustrates the effect of the weights in the spot classification stage. For instance, since the maximum allowed value for the standard deviation is very restrictive, relatively safe zones can be chosen as the best spot even with a low score. Table 6 presents some statistics obtained for the experiment. In order to obtain the influence of the algorithm on the computer's energy consumption, the consumption rate was measured before initializing the algorithm and afterwards. An average increase of 2.96 W was observed. Additionally, the detected spots had an average radius of 3.24 m.   Similarly to the simulation of the second environment, the algorithm's performance was evaluated for several sizes of voxels. Table 7 displays the results. The results obtained are similar to the simulation in relation to the loss of information with the voxel increase. However, a considerable reduction in the execution time of the plane detection step was observed. One explanation is that the dataset has more landing zones compared to the second simulated environment. Therefore, many iterations are performed in the PCA step, unlike the simulation, where there are only a few planar surfaces. In addition, the number of detected spots reduced with the increasing voxel size.

Discussion
Some considerations concerning the acquired results have been already discussed in Section 5. The algorithm was able to detect safe landing spots in simulation and also during the experimental dataset. Downsampling the accumulated point cloud proved to be very useful. The plane detection algorithm is the most demanding in terms of processing time and memory consumption. This is expected due to the loop nature of this part of the algorithm. Therefore, in the downsampling step, reducing the number of points while maintaining a faithful representation of the environment becomes essential to decrease the execution time in the following steps.
Regarding the experimental field test, the point cloud was generated based on LiDAR data at an approximate rate of 10 Hz. Consequently, the algorithm can run in real-time as the algorithm processes the data and outputs a spot in less than 100 ms. However, this depends on the parameters chosen by the user, namely the size of the accumulated point cloud and the number of search points in the plane-detection loop. The plane detection algorithm had an average execution time of 37.87 ms. Therefore, we can see the good computational efficiency of the algorithm. In terms of physical memory, the average usage was 14.90 MB for a total of 8GB of the computer. It was also possible to conclude the algorithm's energy only added an average of 2.96 W to the consumption rate regarding energy consumption.
One limitation of the algorithm is that it does not guarantee that a landing point will be found due to the random approach to choosing the search point from the point cloud. However, the results showed that, by gradually increasing the neighborhood sphere's radius, the algorithm could obtain a satisfactory region to detect the planes.

Conclusions and Future Work
This paper addressed the development of a real-time emergency landing detection method based on LiDAR for UAVs. Methods like the one presented in this paper will allow us to have more safety during UAV out-of-sight missions. The proposed algorithm was evaluated in a simulation environment to characterize the robustness under different conditions and also in an experimental dataset during a 3D reconstruction mission in the staircase of the Monastery of Tibães, Braga, Portugal. The overall results showed the algorithm's ability in both cases, i.e., simulation and experimental datasets, to detect the landing spot in real-time and continuously recalculate new landing spots during a UAV mission based on the parameters defined in Table 2. To the best of the authors' knowledge, this paper contributes to the global system solution to provide emergency landing spots and the implementation of Voxel Grid Filter for point cloud compression. This step allows the algorithm to considerably decrease the number of points while maintaining an adequate representation of the environment.
As future work, we intend to validate the proposed algorithm under different surfaces with water and snow. These regions can reflect the light beams of the LiDAR, and the algorithm would consider areas with slopes equal to zero. Therefore, it would be chosen as a landing point. Furthermore, it is also interesting to benchmark the algorithm performance with other data-structuring techniques regarding computational time and spot detection. Another research work will integrate more weighting factors (wind speed, battery status, etc.) that could be considered during the spot evaluation layer. Funding: This work is financed by National Funds through the Portuguese funding agency, FCT -Fundação para a Ciência e a Tecnologia, within project UIDB/50014/2020.