Path Planning Method for UAVs Based on Constrained Polygonal Space and an Extremely Sparse Waypoint Graph

: Finding an optimal/quasi-optimal path for Unmanned Aerial Vehicles (UAVs) utilizing full map information yields time performance degradation in large and complex three-dimensional (3D) urban environments populated by various obstacles. A major portion of the computing time is usually wasted on modeling and exploration of spaces that have a very low possibility of providing optimal/sub-optimal paths. However, computing time can be signiﬁcantly reduced by searching for paths solely in the spaces that have the highest priority of providing an optimal/sub-optimal path. Many Path Planning (PP) techniques have been proposed, but a majority of the existing techniques equally evaluate many spaces of the maps, including unlikely ones, thereby creating time performance issues. Ignoring high-probability spaces and instead exploring too many spaces on maps while searching for a path yields extensive computing-time overhead. This paper presents a new PP method that ﬁnds optimal/quasi-optimal and safe (e.g., collision-free) working paths for UAVs in a 3D urban environment encompassing substantial obstacles. By using Constrained Polygonal Space (CPS) and an Extremely Sparse Waypoint Graph (ESWG) while searching for a path, the proposed PP method signiﬁcantly lowers pathﬁnding time complexity without degrading the length of the path by much. We suggest an intelligent method exploiting obstacle geometry information to constrain the search space in a 3D polygon form from which a quasi-optimal ﬂyable path can be found quickly. Furthermore, we perform task modeling with an ESWG using as few nodes and edges from the CPS as possible, and we ﬁnd an abstract path that is subsequently improved. The results achieved from extensive experiments, and comparison with prior methods certify the efﬁcacy of the proposed method and verify the above assertions.


Introduction
Unmanned aerial vehicles (UAVs) are highly useful for executing diverse missions not only in urban environments but also in hazardous areas that are not easily reachable, such as forests, deserts, and hilly areas. UAVs (being lightweight, low-cost, and with the abilities to fly at lower altitudes) are now extensively used for a wide range of both military and civilian tasks. Owing to military and civilian investments in UAV technology, this field continuously advances with the passage of time. Based on a forecast by the Teal Group, the market for UAVs is constantly growing globally, and yearly spending on this technology is expected to be higher than US $12 billion by 2024 [1]. Advancements in the technology, such as improved computation capacity, low-cost sensors, artificial intelligence-based algorithms, and fuzzy logic-based decision-making abilities, enable UAVs to easily perform many practical applications in complex environments that otherwise would take a long time and require significantly high costs. The economic and potential applications of UAVs in the real world are most lucrative, including distribution of vaccines [2], tourism security and safety [3], vegetable inspection [4], document delivery for libraries [5], industrial applications [6], forest and urban firefighting [7], sensing of large reduce exploration and modeling of low-probability spaces to overcome these computing issues. Various space modification methods have been designed to increase the pathfinding speed, such as abstractions in a hierarchical form [40], symmetry breaking [41], subgoal graphs [42], jump-point searches [43], accurate heuristics [44], compressed path databases [45], pruning dominant states [46], swamp hierarchies [47], influence-aware pathfinding [48], and constraint-aware methods for navigation [49]. Besides the validity of these latest developments, in most cases, either many low-priority locations of a map are searched uselessly or path quality significantly degrades. Recently, a number of studies considered reducing computation times by dealing with pertinent obstacles that are only crossed along a straight axis in the pathfinding process [50,51]. However, these mechanisms have higher computational complexity and yield non-taut paths if obstacle density is high. Hence, these methods are vulnerable to either returning longer paths or demanding more computing power in determining a path. To address the above limitations, this study presents a new PP method that significantly lowers pathfinding computing time without impacting path lengths by leveraging a Constrained Polygonal Space (CPS) and an Extremely Sparse Waypoint Graph (ESWG) while finding a working path from a 3D urban environment.
The rest of this paper is structured as follows: Section 2 presents the background and related work on renowned PP algorithms. Section 3 illustrates the proposed PP method and describes its main steps. Section 4 explains the results obtained from the simulations. Finally, the conclusions and future avenues for research are discussed in Section 5.

Background and Related Work
In this section, we briefly discuss the UAV operating environment's modeling techniques, the pathfinding algorithms, and geometric-and sampling-based PP methods. The initial step of the global PP is to model the real environment with correct geometric shapes. It is closely linked to the choice of search algorithm because most search algorithms yield good performance when they are collectively employed with a particular environment's illustration. A comprehensive discussion about the performance impacts of distinct environment modeling techniques collectively tested with their respective search methods was given by Sariff et al. [52]. Many UAV operating environment methods have been discussed in the published studies. These modeling methods are categorized as RoadMap (RM), Cell Decomposition (CD), and Potential Fields (PF). Hyungil et al. [53] presented a comprehensive survey on environment modeling techniques used in PP. Each modeling method differs in terms of the scale of space/time complexity, the modeling method's accuracy, and the path quality. For example, when the cell sizes are relatively small, CD-based methods yield poor path quality. In contrast, if the cells are too wide, they are vulnerable to very high time and space complexity. The PF-based methods are prone to getting trapped in local minima, and, thereby, solution quality can be degraded. After modeling the environment with a visibility/waypoint graph, a search algorithm is utilized for the graph's exploration in order to find a path.
Most of the existing search algorithms explore and model whole maps during the PP that can lead to various overheads, such as resource-hogging, needless exploration of many parts of a map, and latency issues during pathfinding. Generally, they do not take advantage of the available useful knowledge related to obstacles' geometries from underlying environments in order to lower the complications in path computing. While finding a path from a provided graph, they mostly hold all edges that are visible in the memory, thereby memory requirements of these algorithms are high. Current bio-inspired search algorithms are vulnerable to pre-mature convergence by relying solely on the specified parameters that can lead to poor path quality. In addition, they were mostly tested in semi-urban environments, and their completeness property may yield infeasible results in realistic-urban environments. To address these technical problems, we proposed a new PP method for computing low-cost paths in order to facilitate UAV's aerial missions in urban environments.

Geometric Path Planning Methods
A geometric PP method that assists in determining a good-quality path in a 3D environment of relatively higher complexity was given in [54]. The environment is represented by using a height reduction strategy to solve the trade-off between path-finding efficiency and accuracy in environment modeling. Unfortunately, it does not reduce the searches in the left-over parts of the area, and, thereby, computing time can be higher in most cases. An incremental PP algorithm considering both local and global constraints for good quality pathfinding was designed by Hu et al. [55]. It is fast, and it reduces the set of good-quality path candidates to only four, with minimal computing time. However, the study ignores space reduction to efficiently find a candidate solution. A geometric PP method considering the minimum turn radius of a UAV for optimal paths in a 3D space was designed by Sikha et al. [56]. The suggested concept is reliable and assists in determining a path with the least complexity. An enhanced heuristic-based PP method to find a good-quality path efficiently by considering UAV flight limits was designed by Kun et al. [57]. A new over-segmentation-based method to determine the free-space overlay of a connected region set was suggested by Plaku et al. [58]. This method quickly finds a safe and good-quality path. However, the approach does not take into account information about sharp turns, narrow passages, and other environmental constraints, which may degrade the suggested method's utility. Furthermore, to augment both efficiency and accuracy, it is extremely important to find irrelevant areas that can be discarded if they cannot help to find an optimal/quasi-optimal path in an environment [59]. Several studies have designed closely related PP methods with undoubtedly reduced time cost, such as the Approximation with Visibility Line (ApVL) method [60]. The ApVL PP method [60] is an improvement of the Base Line-Oriented Visibility Line (BLOVL) algorithm [50], and it is regarded as a highly suitable algorithm for finding an approximate shortest path in 3D urban environments. It reduces the obstacle count significantly (e.g., it processes obstacles that are on a straight line only), and constructs visibility graphs from the chosen obstacles' corners only to incrementally find a low-cost path. Meanwhile, it has relatively higher time complexity. In addition, it either yields longer paths or requires more processing to find a working path. In some scenarios, it is even unable to find a flyable path owing to connected obstacles with straight-line obstacles.

Sampling-Based Path Planning Methods
Sampling-based methods include the Probabilistic Road Map (PRM) [61], Rapidly Exploring Random Trees (RRTs) [62], and their refined versions. These methods have demonstrated effectiveness at quickly generating near-optimal/optimal global solutions. Their algorithmic simplicity makes sampling-based methods applicable to solving both real-time and single-query PP problems. The RRT PP method and its subsequent versions such as informed RRT* [63], Transition-aware RRT (T-RRT) [64], RRT-connect [65], and AnyTime-RRT (AT-RRT) [66] are all complete probabilistically. Most RRT-based methods yield slow convergence rates in complex environments, and they mostly fail to resolve the trade-off between length and time while finding reasonable-quality paths. Sertac et al. [67] designed a better version of the original RRT, named RRT*. This method has a fast rate of convergence compared to RRT, and it has an ability to find a quasi-optimal path with minor post-processing. However, computing issues such as pre-mature convergence, high space complexity, path searching from a whole map, and discarding beneficial samples while converging into a solution pre-maturely make it unreliable for solving practical missions. Jauwairia et al. [68] designed a new variant of RRT* named RRT*-Smart. It has a faster convergence rate, compared to the traditional RRT* algorithm, by using smart-sampling and optimization techniques. However, the main limitations of RRT*-Smart are higher sensitivity to the operating environment, too many iterations, and extensive memory consumption. Yanjie et al. [69] suggested a sampling-based PP method with improved convergence rate. Iram et al. [70] presented a concept relatively closer to our PP method, called RRT*-AB (Adjustable Bounds), to determine low-cost paths. It shows better results than the traditional RRT* method. However, it yields time performance issues due to the near-neighbor search and extensive rewiring operations while optimizing the path lengths. Hence, the shortest paths determined by the existing methods have higher time complexity. Accordingly, a constrained space complexity analysis with in-depth complexity parameters and an obstacle's geometry information has not been simultaneously explored to find a good-quality path with the least time cost.

The Proposed Method
A constrained polygonal space and an extremely sparse waypoint graph-based PP method are imperative for addressing the time complexity issues that emerge due to unnecessary path exploration of low-probability spaces on an obstacles-rich map. The proposed PP method limits path exploration to only the constrained spaces that have a higher probability of containing optimal/quasi-optimal paths, and it safely removes the unlikely spaces in order to hasten the pathfinding computations. It removes the only spaces from a map that likely cannot assist in finding a low-cost solution with high probability. It effectively resolves the two competing goals of efficiency and path length while finding paths for UAVs in urban environments. This section provides a brief overview of our proposed PP method and outlines its workings. In Figure 1, we demonstrate the proposed PP method's conceptual overview. To find a path, P, between source s and target location t for a UAV, while safely bypassing obstacles in UAV flying environment W, the following six key conceptualizations are introduced: (i) operating-environment modeling by using data from a real environment map; (ii) generation of a constrained polygonal space by exploiting obstacle geometry information; (iii) determining and analyzing the complexity of the constrained polygonal space using a multiple criteria-based method leveraging six different complexity parame-ters; (iv) providing need-based extension of the constrained polygonal space to the next level; (v) task-modeling with an extremely sparse waypoint graph that has as few nodes and edges as possible by utilizing the concepts of far distance reachability and direction guidance; and (vi) abstract pathfinding with the A* algorithm from the CPS using an ESWG, and enhancing the abstract path quality by generating additional nodes and edges in the vicinity of abstract path nodes. Concise descriptions of each main component, with relevant equations/procedures, are summarized below.

Representation of the Environment Where the UAV Operates
The initial step in the PP process is to represent the UAV's moving/flying environment from a real environment map with the help of relevant geometrical shapes. Generally, it is a process of dividing W into obstacle-free regions (ξ f ree ) and the obstacle regions (ξ obstacles ). An example of ξ f ree and ξ obstacles is presented in Figure 2, in which black and yellow regions represent the ξ f ree and the ξ obstacles regions, respectively. The obstacles in W can be modeled with sufficient accuracy by leveraging geometric shapes, such as cubes, rectangles, cylinders, circles, polygons, prisms, etc. In our work, obstacle modeling from a raw environment map is carried out with the help of 3D convex polyhedrons having six faces each. The minimum height (e.g., z min ) of each obstacle is 0, and the three other dimensions are random. To generate a convex-obstacles set, we extract the digital map's elevation readings in accordance with digital environment elevation data standards, and we find a convex-hull to accomplish the obstacle modeling task. After the real map's processing, we create a map with fixed convex-polyhedron 3D obstacles. The beginning location of the mission is denoted with s (i.e., a 3D point), where s = (x s , y s , z s ). The target location is denoted with t (also a 3D point), where t = (x t , y t , z t ). Taking into account W's representations, the UAV profile information, and path searching from a 3D environment, the objective of the proposed PP method is to find good-quality paths with the least computing complexity. The proposed PP method fulfills the stated assertions by finding a path from high-probability regions of the map, and the UAV is considered a single 3D point, just like s or t.

Generation of the First Constrained Polygonal Space
Searching for a path by leveraging full map information can be very time consuming, and it may result in serious computing overhead in complex and large urban environments. To address these issues, we convert the full 3D map into a CPS that guarantees P from s to t, and at the same time makes pathfinding computations faster. The space can be reduced using five main steps: (i) stipulating both s and t for a mission, (ii) sketching a straightline l o between s and t, (iii) extracting pertinent obstacles (i.e., only those obstacles where edges/vertices cross with the l o ), (iv) analyzing the geometry of the extracted obstacles, and (v) drawing a minimum-span ζ min polygonal space S 3 1 from the base vertices of pertinent obstacles in such a way that a cross section of the reduced space (e.g., S 3 1 ) is not completely blocked by all of the pertinent obstacles' cross sections. More specifically, the ζ is selected in such a way that the path can be found in each scenario from the constrained space. Obstacles that are on the l o between s and t are called pertinent obstacles. After receiving the s and t locations for the flight, we sketch a l o from s to t. After sketching the l o , four outputs can be derived, which are: (I) neither penetration nor collision, (II) not penetration, but the possibility of a collision exists, (III) penetration, but all obstacles have a lower height, and the UAV can go over them safely, and (IV) both penetration and collision. All four results are demonstrated in Figure 3. If the l o does not collide with any obstacle, and the collision possibility is zero, then path P = l o , which is an optimal path (e.g., a straight line) as given in Figure 3a. Moreover, if the l o does not penetrate any obstacles, but there remain some obstacles in close proximity to the l o with which the UAV could collide with a higher probability, then we deal with these obstacles to generate a safe P, as demonstrated in Figure 3b. In the third case, as given in Figure 3c, few obstacles are penetrated by the l o , all obstacles have a lower height, and the UAV can go over them safely. In the fourth case, as given in Figure 3d, some obstacles are penetrated by the l o , and we need to bypass them at a lower cost to find a P between s and t while fulfilling the stated objectives. Such obstacles are extracted from the map and utilized for pathfinding from the CPS. The complete pseudo-code utilized for extracting pertinent obstacles is illustrated in Algorithm 1. In Algorithm 1, map W encompassing N distinct obstacles, s, and t, are given as input. The set E, where (E ⊆ N) of the pertinent obstacles is retrieved as output. Line 2 can do l o sketching between s and t. Lines 3-7 can do obstacles' extraction that are crossed by the l o . At the end, set E of pertinent-obstacles is collected. Moreover, if no intersection occurs between the obstacles and the l o , and, then, E = ∅ will be the output.
After getting set E of pertinent obstacles, we enlarge the obstacles by a safe distance (d sa f e ), and apply the flying minimum and maximum limits, denoted as h min and h max , respectively. Subsequently, we analyze the pertinent obstacle cross sections that are on the l o between s and t. Then, we draw a bottom boundary, Ω b , around the bottom vertices of the pertinent obstacles, and a top boundary, Ω t , from minimal height h with path guarantees keeping s and t as two endpoints. This forms a 3D constrained space where the shape can resemble a 3D polygon, and we call this constrained region the CPS, represented with S 3 1 . This CPS converts the more difficult problem of UAV pathfinding into the relatively easier problem of pathfinding for a 3D point. A visual overview of S 3 1 is demonstrated in Figure 4. The S 3 1 can be simply defined as full map W partitioned into a small space/region where the outline is the same as a 3D polygon, with s and t as two endpoints. It is obtained by drawing a boundary around the outermost vertices of the pertinent obstacles in such a way that a path can be found from it regardless of its quality. The process of transforming a full map into S 3 1 is given in Figure 4a-c. In Figure 4a, environment map W is shown, which will be converted into S 3 1 . Figure 4b shows the l o drawing, and identifies the corresponding pertinent obstacles (e.g., yellow obstacles) that were crossed by it; these obstacles will be used to subsequently generate S 3 1 . Figure 4c shows the outline of S 3 1 in 2D form that is obtained by drawing a boundary around obstacles identified in Figure 4b as a consequence of l o penetration.
End for 8

return E
The CPS can enclose the pertinent obstacles-just part of them or as a whole. In some cases, due to complex 3D environments, S 3 1 can enclose obstacles that do not belong to set E but that are part of S 3 1 , either partially or completely, and we include such obstacles in E and utilize them in pathfinding. The essence of S 3 1 is that it guarantees a flyable path between s and t. However, S 3 1 may or may not be an ideal choice for a good-quality path (e.g., optimal or quasi-optimal) due to several complexity parameters (as illustrated in Figure 1) about obstacles. By considering such potential complexity parameters and a probabilistic analysis of optimal paths, we conducted a CPS complexity analysis leveraging six complexity parameters (also known as complexity constraints), prevailing in S 3 1 that relate to the obstacles' geometries, and a low-cost path tends to lie outside S 3 1 , in most cases, with a significantly higher probability.

Determining and Analyzing the Constrained Polygonal Space Complexity Using a Multicriteria-Based Method
In order to check whether S 3 1 is good enough for optimal/quasi-optimal pathfinding or not, we performed a multi-criteria-based complexity analysis of S 3 1 using multiple complexity parameters before task modeling and pathfinding. We computed the complexity, χ, of S 3 1 by leveraging detailed information regarding the obstacles' geometries. We employed six complexity parameters: the proportion of free spaces, the obstacle occupancy in distinct regions of the CPS, the complexity of obstacle-avoidance options, the proportion of connected obstacles, the length deviations from the optimal path, and obstacles' tendency in the CPS that hinders the solution quality. Through extensive simulations and analysis, it is found that there exists a very firm relationship between the complexity parameters of the CPS and path quality. The total χ is the weighted sum of six parameters cited above. Brief overviews of those six complexity parameters, with their procedures and equations, are described below.

Free Spaces' Ratio
To compute a feasible, safe, and smooth P, it is highly enticing that the amount of free spaces must be high in the CPS. To measure the obstacle-free spaces in the CPS, we first determine the size (ξ) of the S 3 1 , blocked-spaces (ξ obstacles ), and free-spaces (ξ f ree ). The overall size (ξ) of the CPS that is in the form of a polygon can be obtained by the Gauss determinant using Equation (1): where x, y are the coordinate values of the CPS boundary, |g| denotes total vertices of the CPS boundary, and h is the height of the CPS. Out of the ξ-sized CPS, we find the amount occupied by the obstacles (ξ obstacles ) using Equation (2): where n denotes total obstacles' count present in a S 3 1 and ωO i denotes the obstacles' occupancy. The occupancy of an O i obstacles can be determined using Equation (3): where O W , O H , O L denote the width, height, and length of an obstacle, respectively. The free space (X f ree ) amount, where UAV can fly safely, can be determined using Equation (4): The ratio (r f ) of the free spaces can be calculated using Equation (5): The value of r f ranges between 1 and 0. We represent this ratio as (1 − r f ) to compute the occupied spaces value in overall CPS S 3 1 complexity computation.

Deviation in Length from an Optimal Path
The proposed method is a global PP approach, and all information about the obstacles geometries is known in advance. By utilizing the obstacles geometry information, we can estimate the path length without calculating the actual path. We estimate length of an optimal path L o , where L o = l o as an optimal path and estimate the deviation in it that can occur due to obstacles that are on l o . We call it deviation D from the optimal path and formalization used to estimate is explained below.

•
Computing the optimal path L o that is a straight line between s and t using Equation (6).
• Calculating the deviation D in optimal path L o due to the presence of obstacles on the l o between s and t locations. The D in the paths length due to an obstacle (e.g., O i ) is calculated using Equation (7).
• Estimating the total deviation D that can likely occur due to the presence of the obstacles between s and t in the CPS using Equation (8).
• Calculating the length of the estimated paths (L 1 ) avoiding all obstacles that are on the l 0 using Equation (9).
where L 1 denotes the lengths of the paths avoiding obstacles in the selected space, L 0 denotes the Euclidean distance between s and t locations, n represents obstacles' strength in the selected space, and D i denotes the degradation in path length due to obstacles. • Computing the complexity C p of the estimated path that can assist in analyzing the S 3 1 complexity using Equation (10): The value of the C p ranges between 1 and 0. We use this value to represent the estimated path complexity in terms of length in overall CPS S 3 1 complexity evaluation in Equation (20).

Complexity of the Obstacles' Avoidance Options
There are usually four options in total to bypass any obstacle present in a W such as right, left, up, and down (in the case of hanged obstacles or flying obstacles). Meanwhile, after the space reduction, the number of options to avoid obstacle will likely be reduced, and there can be an increase in the complexity of remaining options. Because of this, a path may be taut and path length can be prolonged. Hence, while determining the CPS complexity, we take into account the complexity of options needed to bypass obstacles. To calculate the complexity (C AO ) of each option, the entropy concept is employed. Entropy is acknowledged as the most effective and accurate measurement for similar tasks in numerous fields. In this work, we consider the urban environment; therefore, the P cannot go beneath the obstacles since bottom height of each obstacle is zero, and, hence, there are only three options in total to avoid any obstacle. The strategy below is employed to calculate the C AO .
1. Find the proportion (p i ) value of every avoidance option (i.e., AO l , AO r , AO t ) category using Equation (11).
where b denotes the total AO, and its value can be determined using Equation (12).
2. The complexity C AO of avoidance options can be calculated using Equation (13): The normalized value of the C AO lies between 0 and 1, denoted as C AO ∈ [0, 1]. The C AO with 0 value means that avoidance complexity is low (e.g., all obstacles can be avoided from the same side). In contrast, the C AO value 1 means that enough variations exist in options to avoid all obstacles, and the path can contain many turns. In the CPS analysis, we take into account the C AO values.

Occupancy of Obstacles at Distinct Regions of the CPS
Besides the other complexity parameters described earlier, another important parameter that can lead to genuine performance concerns while pathfinding is the obstacles' occupancy at distinct regions of the CPS. If obstacles in large numbers are clustered at one location (e.g., obstacles' placement in the CPS is uneven), then the path quality likely degrades. The obstacle occupancy at one place introduces cycles/sharp-turns in a P because the P revolves around boundaries of many obstacles before approaching t. To calculate occupancy Π o of obstacles, we partition S 3 1 into n sub-spaces {s 1 , s 2 , s 3 , . . . , s n } and find the obstacles occupancy in each subspace s i . The obstacles' occupancy in a s i subspace can be determined by taking the ratio of the obstacles' occupancy ξ obstacles in the s i divided by overall obstacles' occupancy in the CPS. To calculate occupancy, the S 3 1 is partitioned in five equal-size sub-spaces. The occupancy Π o i of the s i can be mathematically expressed as where ξ obstacles denotes occupancy of all obstacles in total from the CPS as given in Equation (2) and s i obstacles represent the ith subspace's obstacles occupancy, and its value can be computed using Equation (15): where ωO i represents an obstacle's volume, and O denotes all obstacles count in the CPS. After computing the occupancy of five sub-spaces, we determine the overall occupancy Π o of the S 3 1 using the following equation: where Π o is the obstacles' occupancy in the S 3 1 . The rationale to choose maximum values is to effectively deal with the worst cases. The occupancy analysis assists with finding the smooth paths by giving considerable attention to the regions of high occupancy.

Ratio of Obstacles' Tendency in the CPS
In some scenarios, the CPS can enclose more obstacles compared to the W (e.g., the tendency of obstacles on l o is high compared to the whole W). Hence, it is viable to assess the impact of obstacles' tendency to yield a good quality path. To analyze the obstacles' tendency T o , we find the number of obstacles in a CPS, and take a ratio with the obstacles' count present in a W. We denote the number of obstacles present in the CPS with n and number of obstacles present in a full map with N, respectively. We determine the value of T o using Equation (17): The value of T o can lie between 0 and 1, T o ∈ [0, 1]. The higher value of T o means that more obstacles are present in the CPS. In our work, we take into account the T o value while determining and analyzing the complexity of the S 3 1 .

Ratio of the Connected Obstacles
In some cases, some obstacles exist that are not directly penetrated with l 0 , but they have connections with the pertinent obstacles (e.g., obstacles directly crossed with l 0 ). These obstacles can escalate the time of path computation and yield unnecessary turns in the path. Hence, while analyzing the CPS complexity, it is paramount to take into account the connected obstacles' effect along with other five complexity parameters. The ratio of the connected obstacles can be found by counting the connected obstacles in the CPS divided by the obstacles' count in the S 3 1 . The count of connected obstacles can be determined using Equation (18): where O j is the pertinent obstacle, and O CON denotes the connected obstacle with O j (e.g., pertinent obstacles). The overall ratio of the connected obstacles (r co ) can be computed using Equation (19): where n represents connected obstacles' strength, and n shows the number of obstacles in E. When all six complexity parameters' values have been calculated, the total complexity χ of the S 3 1 can be quantified using Equation (20): In Equation (20), 1 − r f denotes the ratio of spaces occupied by obstacles, C p denotes the deviation in path length from an optimal path, C AO means the complexity of options while avoiding obstacles, Π o is the occupancy of the obstacles at the distinct region of the CPS, T o denotes the tendency of obstacles in the CPS in relation to a full map, and r co denotes the ratio of connected obstacles. The drawing of all six complexity parameters described in prior subsections (e.g., Sections 3.3.1-3.3.6) is given in Figure 5.
For calculation simplicity, we used complexity parameters values in normalized form; therefore, the χ(S 3 1 ) ranges between 0 and 1. In Equation (20), w i , where i = 1, 2, 3, 4, 5, 6 represents each complexity parameter's weight, and they fulfill two conditions, (i) w 1 + w 2 + w 3 + w 4 + w 5 + w 6 = 1, and (ii) w i > 0. We adjust each parameter weight by taking into account the significance and influence of every parameter in the CPS complexity analysis. The probability σ of a path P to be found from the S 3 1 with good quality is given as follows: where χ(S 3 1 ) represents the complexity of the S 3 1 , and T denotes a threshold. If σ(P) = 1, no extension of space is required because S 3 1 is appropriate for low-cost pathfinding. Meanwhile, if σ(P) = 0, then an additional space will be needed to find good quality paths. The threshold T value relies on numerous global factors of the W, and local constraints (e.g., related to the UAV). In simulations, we set the threshold value to 0.7 to make a decision about the space expansion. We did substantial experiments to validate the T value using P's length as a main criteria. However, T's value can be tuned flexibly based on the UAV's workspace and resources.

Expansion of the First Constrained Polygonal Space
Although S 3 1 always finds path P, it does not guarantee the quality of P in each scenario due to higher complexity in the obstacles. To circumvent this issue and ensure consistent quality for P in each scenario, the scenarios that require a bigger space are identified carefully through the CPS complexity analysis utilizing the six different parameters. With the assistance of the complexity analysis of S 3 1 , for a good-quality P, we can accurately identify the cases that require relatively more space than already in S 3 1 . Having sufficient information about the obstacles connected with the boundary of S 3 1 enables us to flexibly expand the space to the next level. We adopted this method to expand the space, since it yields less computing overhead and significantly enhances path quality. Hence, by processing obstacles that are penetrated by the boundary of S 3 1 , and by marking a polygonal boundary in an analogous way, the first CPS formation emanates into a second CPS of a relatively bigger size, compared to the S 3 1 as visually depicted in Figure 6b. We denote this expanded space with S 3 2 , and it encompasses S 3 1 fully. The S 3 2 includes pertinent obstacles fully, and it provides greater opportunities for P to be determined from S 3 2 solely. The utility of the S 3 2 is that it is highly desirable space for producing P of shortest lengths. The space can be extended to nested levels in identical manner. Meanwhile, we expand the spaces only up to 2-levels because an optimal/quasi-optimal P tends to lie in S 3 1 and S 3 2 with acceptable probability. After the selection of appropriate highest priority space, an ESWG is constructed for pathfinding.

Extremely Sparse Waypoint Graph Generation from the Selected CPS
Waypoint graph (WG) is one of the approaches for task modeling and pathfinding, respectively. The WG constructs an indirect and compact graph connecting s with t by catching the connectivity of the ξ f ree to form a multiple paths' network. However, generating a WG is very expensive in terms of computation. The overall time complexity of generating a WG with n nodes is O(n 3 ). Many studies that have focused on lowering the complexity of the WG generation have been reported by joining adjacent obstacles, altering obstacles' shape, and ignoring small obstacles. More recent research [60] highlights that WG's time complexity in the 3D environments can be decreased to the O(n 2 ) by only considering the obstacles crossed by the l o . To expedite the time complexity reduction of the WG, this paper suggests a new concept of an ESVG construction method which does not compose a dense WG. This method forms an ESVG from the CPS like a roadmap with connectivity between s and t through intermediate nodes. An ESVG is a double-edge type graph G of reachable and mutually-visible locations, mathematically expressed as: To construct a G from the 3D CPS, two steps are generally applied: making a nodes' set X and generating an edge set Y. The initial step is about creating nodes set X. We utilized three vertices of obstacles, bottom, top, and mid to make an ESWG for the first time. Both bottom and top vertices' geometry values are known, and mid vertices can be found leveraging the midpoint formula on top and bottom axis values. Every obstacle has total eight vertices (e.g., four bottom and four top). An ith obstacle's vertices and their respective values can be expressed mathematically in a matrix as demonstrated in Equation (22). The height of bottom four vertices of each obstacles are transformed to the h min and top vertices of the obstacles have the same height as of the CPS height (e.g., h). In below metrics, the value of z min is zero but after adjustment becomes h min = z min : x min y min z max ; x min y min z min x min y max z max ; x min y max z min x max y min z max ; x max y min z min x max y max z max ; x max y max z min For example, an O whose original z max is greater than the CPS height, and the bottom and top vertices are fully known, the pair of mid-points' two side faces denoted with f 1 and f 2 can be determined using Equations (23) and (24), respectively: A similar procedure can be utilized to find the pair of vertices around all obstacles. After computing set X from the pertinent obstacles, we add both s and t in a set X and generate set Y of edges through two novel strategies and visibility-checks. Any two nodes u and v in X are inter-visible if a uv segment of line connecting them is collision free. A function named-line-of-sight (LOS) determines the visible segments among pairs of nodes through visibility analysis. The time complexity of this mechanism heavily relies on the function of LOS checking, and number of nodes. Meanwhile, in our work, we incorporate two additional strategies of far-reachability (FR) and direction-guidance (DG), thereby time complexity is significantly reduced. In addition, it only adds the edge between a vertices' pair that are as far as possible from each other and that guide to the t's direction. We set the visibility to off/false using coordinate values for those pairs of vertices that are on the same obstacle but do not favor the direction of the t. Hence, the visibility checking function has less time complexity in making a G. The time complexity of an ESWG formation is the O((n f l) 2 ) time, where l denotes the number of levels, n denotes obstacles' count, and f represents the counts of obstacles' facets. However, the upperbound of the l is constant (e.g., l = 3 ); therefore, ESWG's time complexity is O(( f n) 2 ). With the help of X and Y, a G is obtained that has reliable connectivity between s and t, and it encompasses all characteristic of a roadmap.

Path Finding from an ESWG and Enhancing Obtained Path Quality
Once an ESWG is modeled, a path searching algorithm is employed to search a P from it. In this paper, we used A * algorithm for computing a P between s and t from an ESWG. The A* is reliable algorithms for extracting a P of low-cost. The evaluation function utilized by this algorithm is expressed in Equation (25): In Equation (25), the f (n) denotes the estimated path cost in total between s and t via a node n, the g(n) denotes the actual distance to reach node n, and the h(n) is a heuristic function that computes the distance from node n to t. This algorithm was selected to make P's computing process fast. By exploring an ESWG using this algorithm, an abstract P is found. We consider both length and time of the obtained P for evaluating its quality. Meanwhile, in some cases, the P cannot be of minimal length, which needs post-processing to shorten it. We present the working of the A * algorithm while finding a path between s and t in Figure 7.
In many UAV practical applications, the length of the P is paramount, and to preserve UAV resources, it should be minimum. To address this issue, we shorten P's length by including more nodes in close proximity of it, and refine the sharp turns. The path quality improvements is mainly carried out by determining the adjacent P's neighbor nodes, find the proximity between the adjacent neighbor nodes and P nodes, and in the close proximity of the P, we introduce new nodes with relatively denser resolution and add smooth edges. The reason to add more nodes closely is to retain visibility to improve P quality. After injecting additional nodes, a P of good quality is obtained by jointly using the newly added and the P nodes. This path-refining method has the potential to improve path quality significantly with reduced computing cost.

Results and Discussion
This section explains the simulation experiments and corresponding results. The performance of the proposed PP method was analyzed using two criteria: time complexity and path length compared to prior studies. To make the proposed PP method a benchmark, we compared the simulation results with randomized motion planning and visibility graph-based algorithms. The simulation tests were performed and compared using Matlab v. 9.8.0.1451342 (R2020a) on a computer running Windows 10 with 8 GB of RAM and a 2.6 GHz CPU. In the tests, we assumed a 25 kg fixed-wing UAV similar to ones used in existing studies. We took into account both global and local constraints in the simulations. The parameters of the local constraints (i.e., on the UAV) were a 1 m wingspan and a maximum turning angle at a radius of π/6 . The minimum and maximum flying altitudes were h min = 23 m, h max = 155 m. The global constraints belonged to the geometry of obstacles in W. We consider six complexity parameters that can significantly hinder the quality of a P and the UAV's safety while selecting a space size. We assumed that the UAV had enough battery power to complete the task in one flight. The safe distance to avoid collisions with obstacles was 10m (d sa f e = 10 m). We assumed that wind was negligible during the flight. The proposed method finds P using an ESWG that respects both global and local constraints. The weights of space complexity parameters were w 1 = 0.2, w 2 = 0.1, w 3 = 0.2, w 4 = 0.2, w 5 = 0.1, and w 6 = 0.2. We assigned values to these weights by considering the significance and influence of each parameter on the accurate space selection for an optimal/quasi-optimal pathfinding. We tested our method with diverse combinations of the weight values, and analyzed the accuracy of space selection for optimal/quasi-optimal paths. Subsequently, we determined the best combination of these weight indexes' values that make accurate space selection consistently. Furthermore, these weight values were validated via numerical analysis by computing optimal paths from numerous maps using whole W, and analyzing the number of times optimal paths tend to lie in the selected space. The locations for s and t were chosen randomly during experiments. We compared our PP method's performance with two existing algorithms: the ApVL algorithm proposed by Guillermo et al. [60] and the RRT*-AB algorithm proposed by Noreen et al. [70]. Both comparison algorithms are state-of-the-art for PP. We tested them on our maps to compare the performance of our method with them. We show a sketch of the 3D maps employed in the experiments and three exemplary paths from each method in Figure 8. The path produced by the proposed PP method was more smooth and shorter than the paths from the other methods. To analyze and compare the proposed PP method's results, we designed three distinct scenarios with sufficient obstacles in the 3D environment maps. Each scenario was tested with all three methods, and the results were analyzed. All obstacles had a random width, depth, and height with a rectangular-shaped base. Comprehensive details on maps counts, map sizes, s and t locations, numbers of obstacles and their geometric information, etc. are given in each scenario description below.

Comparison with the Existing Approaches by Varying Map Sizes and Obstacle Counts
This scenario is defined with W at sizes ranging between 100 m × 100 m × 300 m-1000 m × 1000 m × 400 m. It encompassed 50 maps with distinct obstacle counts (e.g., . For the sake of simplicity and rational comparisons, we categorized all maps into 10 distinct groups considering both map size and obstacle strengths, as given in Table 1. The locations for s and t were marked in alternate places for every test/map. Furthermore, the obstacle density in W varied on each map. The ApVL algorithm [60] processed only obstacles that were on the l o and generated a dense graph to find P. However, the ApVL algorithm has higher complexity, and it produced a non-taut P in most scenarios due to the connected obstacles in urban environments, as shown in Figure 8b. The RRT*-AB algorithm restricts the space, but it explores many locations while finding P. The P generated by this algorithm had a longer length, and computing time was immense. In contrast, the proposed PP method processed fewer obstacles and employed an ESWG to find a P with good quality. The complete information about maps utilized in this particular scenario and the mean computing time for task modeling by our PP method and its comparisons with the two previous algorithms are given in Table 1. The computing time for modeling the UAV environment shown in Table 1 is the total time needed to constrain the space, an ESWG construction, and an ESWG's expansion for path-quality improvements. From Table 1, we can see that time surged with an increase in the map size and obstacle counts. Through extensive comparisons with prior PP algorithms, our method lowered computing time for task modeling by 15.05% on average. The pathfinding results (i.e., time needed and path length) and their comparisons with the two existing algorithms are depicted in Figure 9. Both computing time and path length are the mean of five maps in each map's group (given in Table 1) with arbitrary obstacles' placement. The simulation results emphasize that, for each method, there is a surge in the computation time with the increase in W's complexity. Moreover, the proposed PP method shows 15.4% and 33.6% curtailment in mean computing time compared to the ApVL algorithm and the RRT*-AB algorithm, respectively. In path lengths, the proposed method shows 5.34% improvements compared to the ApVL algorithm. Moreover, average improvements in the path lengths compared to the RRT*-AB algorithm are 6.34%. These results highlight that the proposed PP method is superior in terms of both computing time and path length over prior algorithms.

Comparison with the Existing Approaches by Varying the Number of Obstacles
This scenario is comprised of 10 maps with varying numbers of obstacles in a W of 1 km 2 to analyze the impact on our PP method's performance from increasing the number of obstacles. The obstacles were clustered between s and t in such a way that all methods would avoid them during pathfinding. Figure 10 presents the results of our proposed PP method and a comparison with the other algorithms from varying the number of obstacles. When W enclosed more obstacles, the suggested method could quickly determine a goodquality and safe P from W. It was better than the ApVL and RRT*-AB algorithms, based on the metrics, even when varying the number of obstacles in the CPS. The P determined by the proposed method was the shortest and smoother than the other two algorithms. Through simulation results and their comparison with previous methods on 10 obstacles' counts-based maps, the proposed method decreased the computing time for pathfinding by 27.06%, on average. For path lengths, the paths generated by our PP method, on average, were 4.6% shorter (i.e., produced at a lower cost) than the previous methods.

Comparison with the Existing Approaches by Varying Source and Target Locations
This scenario was tested using a W at 150 m × 150 m × 400-1000 m × 1000 m × 400 m. It encompassed five maps with obstacle counts of up to 25. We analyzed our method's performance through seven runs on every map with different coordinates for s and t in each run. By varying the positions of s and t, the number of obstacles to be modeled in each run/test can be distinct, and, accordingly, comparison metrics can vary with W complexity. The proposed method's averages, obtained from the seven runs on each map, are given in Table 2. The results indicate that the proposed method yielded comparable performance in all tests.
The proposed method yielded an average computing time of 0.71 s, compared to the ApVL and RRT*-AB algorithms, which had mean computing times of 1.01 s and 1.37 s, respectively. In addition, the proposed method lowered the path length, compared to both prior algorithms, by 5.05%. Although the proposed method gave better results, a relatively higher number of initial P nodes can degrade its performance in complex environments. The worst-case complexity with our method was O(n 3 ). However, the test results revealed that time complexity did not accelerate like O(n 3 ) in all test scenarios for finding good quality paths. The results obtained from all these scenarios showed that our method performed consistently better than the ApVL and the RRT*-AB algorithms. Aside from the path lengths and computing times, we analyzed and compared its performance against prior algorithms with respect to graph/tree nodes and path node counts. In Figure 11, the results of the proposed method in terms of average path nodes and graph/tree nodes for the above experiments are presented. As shown in Figure 11, both the graph nodes and path points of our ESWG method were lessened, compared to the prior methods.  Analysis of the memory requirements: As shown in Figure 11a, the proposed PP method generates a WG with fewer vertices. In addition, it does not register visibility of the edges that contribute minimally in an optimal/quasi-optimal path due to less coverage in terms of distance or they are not in the same direction as the target location while generating an ESWG. Therefore, the memory requirements of proposed method are not high compared to the existing methods and a complete graph. However, most global PP methods keep all visible edges in the memory that significantly increase the memory requirements. Furthermore, the visibility check function is called a substantial number of times in visibility graph-based PP methods, thereby space complexity drastically increases. The proposed method resolves these space complexity related issues through reduction in search space, modeling tasks with an ESWG, producing far lower but relevant edges and vertices, and reducing the visibility checks between vertices by incorporating farreachability and direction-guidance concepts while making an ESWG.
The proposed PP method is complete, and it can be used for many UAV practical applications in urban environments. The proposed method gives good performance due to two main concepts: (i) a new space reduction concept is proposed, which not only assists in lowering the time complexity by restricting the path exploration in the space of highest priority, but also assists in finding low-cost paths in most cases; and (ii) an ESWG, which models the tasks with far lower edges and vertices that curtails the computing time of path searching significantly by making a direction-guided search of the target location. It effectively resolves the trade-off between optimality and efficiency in pathfinding from an urban environment populated by various obstacles.

Conclusions and Future Directions
This article proposed a new PP method based on CPS and an ESWG to enable a UAV's safe navigation in 3D urban environments. The main objectives of the proposed PP method are to lower the time complexity in both task modeling and pathfinding without degrading the path quality for UAVs operating at lower elevations in urban environments with fixed, convex obstacles. The main contributions of this article are listed as follows: • We propose a new PP method based on CPS and an ESWG that has the potential to find an optimal/quasi-optimal path with considerably reduced time complexity. • We propose a new space reduction method that abstracts the full map into a 3D constrained polygonal space that guarantees a path for the UAV's mission. • We analyze the effectiveness of CPS for low-cost paths considering six complexity parameters, including the ratio of free space, obstacle density in distinct regions of the CPS, the complexity in the options for avoiding obstacles, deviation in the length from the optimal path, the ratio of connected obstacles to pertinent obstacles, and obstacle tendencies in the CPS. • The proposed method enlarges the CPS to the next level/space by including obstacles that are in close proximity to the first CPS if the first CPS fails to provide an opportunity for low-cost solutions owing to a higher complexity from obstacles in it. • The proposed method generates an ESWG from the CPS, leveraging the principles of maximum distance reachability, having only a few nodes and edges, and direction guidance, and it computes an abstract path that is further improved with the assistance of more nodes and edges around it. • This initial work makes use of obstacle information from an underlying environment in order to lower the computing overhead for pathfinding without compromising path quality in 3D urban environments.
The proposed method performance is substantiated through extensive tests, and, in most cases, it performs consistently better than prior PP methods. It lowers pathfinding time complexity considerably by restricting path exploration solely in the highest priority CPS that has the greatest chance of providing an optimal/sub-optimal path. While conducting the tests, we considered numerous parameters related to the underlying operating environment's complexity and UAV's safety. Meanwhile, during the PP at lower elevations in urban environments, we may need to consider hanged/thin obstacles (e.g., electrical wires and poles in the streets). Another group of evaluation parameters can be the wind/gust and wind/crosswind (e.g., wind's direction and speed), especially when passing through buildings. Hence, further testing with these parameters is yet to be investigated in future work. Finally, we intend to analyze the fidelity of our proposed PP method with other task modeling methods (e.g., Voronoi diagrams, grids, and navigation meshes, etc.).