The structure of the remainder of the paper is as follows. Section 2 provides a succinct summary of the authors' previous work with CPP for fixed‐wing UAVs, alongside the details of new work which has been used to develop the cost function in this paper. It also includes experimental evidence to show that the assumption of the uniform steady wind is valid over the area and time‐scale of a typical survey. Section 3 explains the proposed recursive‐greedy decomposition algorithm in detail. Section 4 presents and compares the paths generated from the proposed method to that of previous methods in simulation. A number of examples real‐world fields are used to demonstrate the algorithms functionality, and a Monte‐Carlo simulation using 30 randomly generated polygons is used to make a numerical comparison of the methods. Section 5 details flight testing performed, to assess time prediction accuracy and further comparison to previous techniques are made. Section 7 and 7.1 provide concluding remarks and suggest extensions to this study for the future. Appendix A shows image coverage results and orthophoto outputs of one of the flight tests, demonstrating that total image overage has been achieved in the real‐world.

Note that this study is in parallel to the work in Coombes, Fletcher, et al. ( 2018 ). Whereas the focus of the previous paper was to produce a fully optimal path given large computational resources, this paper develops a suboptimal algorithm that trades a small loss of flight time efficiency for a large decrease in calculation time, enabling real‐world use in the field.

An example output is shown in Figure 1 where the UAV path is defined for two decomposed cells. Results will be provided in the form of simulations and extensive flight testing to compare its effectiveness to previous methods, and to show the validity of the flight time in wind model and the steady uniform wind assumption.

This paper aims to develop a CPP algorithm for complex and irregular field shapes that is fast enough to run when on‐site conducting an aerial survey. The algorithm should take into account real‐world considerations like wind, NFZ, and launch points and output a set of way‐points that can be uploaded to most commercial fixed‐wing UAV platforms. This enables operators to conduct a custom survey for user‐defined aircraft, imaging systems, and survey requirements. The specific novel contributions from this paper are given as follows:

In literature, research into complex large scale surveys by a fixed‐wing UAV is currently relatively sparse. For complex‐shaped fields; work is limited to extensions of ground robot decomposition techniques (Huang, 2001 ; Li et al., 2011 ; Xu, Viriyasuthee, & Rekleitis, 2011 ), and does not account for the additional complexities previously mentioned. In contrast, UAV path planning in the wind and the associated modeling is well understood (Ceccarelli, Enright, Frazzoli, Rasmussen, & Schumacher, 2007 2007 ; Jennings, Ordonez, & Ceccarelli, 2008 ; Schopferer & Pfeifer, 2015 ), however, these techniques are not suited to CPP for aerial surveys. Only a single paper outside of our own work has been found that combines these two areas. Richards ( 2018 ) plots a full boustrophedon path in the same orientation across the whole concave field polygon, then calculates a time cost in steady uniform wind for each sweep and uses a modified traveling salesman solver to find the fastest order and direction to traverse all the sweeps without forcing the aircraft to perform maneuvers that it is incapable of. It is mentioned that this technique and others are dependent on an accurate flight time model in the real world. Therefore, comprehensive flight testing should be performed to verify such flight time models. While there are a few papers that do perform flight testing on fixed‐wing UAV CPP (Avellar, Pereira, Pimenta, & Iscold, 2015 ; Paull et al., 2014 ; Xu et al., 2011 ) they did not attempt to predict flight time, rather just showing that complete coverage can be achieved. In addition, this paper assume a steady uniform wind, which is a standard assumption for many large scale planning problems. However, there are some examples for use with UAV gliders that use dynamic‐wind maps to maximize soaring times, where actual flight tests were conducted successfully (Depenbusch, Bird, & Langelaan, 2018a , 2018b).

Much of the work in this area assumes a negligible effect from the wind. It is shown by Li et al. ( 2011 ) that minimizing the number of turns (NT) in the survey, that is, aligning the angle of a boustrophedon path with the long axis of each decomposed convex polygon, results in the most efficient survey path. This works by minimizing the time wasted in the turns between sweeps, where no images can be taken due to the aircraft roll angle. Huang ( 2001 ) tries to minimize flight time in a survey by using the similar heuristic cost function, the sum of long axis lengths (minimal sum of altitudes [MSA]) from each decomposed convex polygon. In addition to effectively trying to limit the number of turns, it also has the effect to limit the total number of decomposed cells which serves to decrease the extra wasted transit time between them. However, the wind strength and the flight path angle to the wind will have a massive impact on the flight time of a survey (Coombes, Chen, & Liu, 2017 ; Hovenburg, de Alcantara Andrade, Rodin, Johansen, & Storvold, 2018 ; Richards, 2018 ) due to high wind velocity relative to airspeed (typically 20%–50%). Flight times can be significantly improved by flying perpendicular to the wind direction (Coombes, Fletcher, Chen, & Liu, 2018 ). This means that in wind, finding the best path is a complex balancing act between flying perpendicular to the wind, minimizing the number turns and minimizing transit times between decomposed cells. Furthermore, CPP can be complicated further by being close to areas in which overflight is not allowed (e.g., built‐up areas) or not advised for flight safety (e.g., power lines) known as NFZs.

There are three main categories of techniques developed in previous literature to perform sensor CPP tasks; cell decomposition (Choset, 2000 ), grid‐based methods (Arney, 2007 ), and informative path planning (IPP) (Galceran & Carreras, 2013 ). IPP guides the vehicle around the ROI in real‐time, aiming to maximize information gain while ensuring total coverage (Paull, Thibault, Nagaty, Seto, & Li, 2014 ). These methods are not suited to fixed‐wing UAVs because they tend to have quite long erratic paths and require a gimbal to keep the imaging system orthogonal to the ground during the frequent turns (Acar, Choset, Rizzi, Atkar, & Hull, 2002 ). Alternatively, grid‐based methods discretize the whole ROI polygon into a grid (Elfes, 1987 ) and plan a path the covers each unoccupied cell. This only approximates the region, also generating paths that would lead to a large number of turns, making grid‐based methods unsuited for fixed‐wing aircraft surveys for the same reasons as IPP. Grid‐based methods also tend to require maneuvers that a nonholonomic vehicle, those with a constant forward motion, could not feasibly perform. Cell decomposition methods, on the other hand, decompose the ROI into a number of smaller cells which are individually covered with a simple and efficient path. These methods lend themselves to aerial survey CPP, because the long, thin decomposed cells will require fewer turns and have the aircraft in straight level flight for longer periods of time where high quality stable images can be obtained (Li, Chen, Er, & Wang, 2011 ).

Current commercial mission planning software packages are often simplistic and are not suited for large highly complex field shapes (Coombes, Chen, & Liu, 2018 ). They tend to take no account for the shape of the field, wind or any other operational factors. While there have been a number of CPP techniques developed in the literature to generate paths for a total coverage of an area, the vast majority have been primarily designed for ground robots (Latombe, 2012 ). Contrary to ground robots, aircraft are not constrained to operate within the bounds of the polygon and are subject to additional complexities such as no‐fly zones (NFZs). Therefore, although similar, algorithms developed for ground robots are not directly applicable for UAVs.

Farms tend to be large areas, and due to historic land boundaries and local features such as roads, rivers, forests, and towns, are often complex irregular shapes (especially in Europe). A study of agricultural fields in Finland found that only 13% of farms were simple convex shapes (Oksanen & Visala, 2009 ). Until recently, large scale unmanned aerial surveys have not generally been performed because beyond visual line‐of‐sight flight for UAVs is legislatively difficult. However this is about to change; the Civil Aviation Authority (CAA) in the UK (Department for Transport (UK), 2016 ) and the Federal Aviation Administration (FAA) in the United States (Alta Devices, 2017 ) are pushing to make this more straightforward. This will enable farmers and operators, for the first time, to legally perform large‐scale complex surveys over a number of fields in a single flight. The next limit to the size of the surveys that are possible is the limited endurance of UAVs. While UAV flights are limited by endurance, the coverage path planning (CPP) mission plan can make a huge difference to how long the flight can take, thus the total possible coverage area for given flight time. This can make the difference between performing the survey over single or multiple flights, which is to be avoided as it would add operational time, complexity, and cost.

The process of remote sensing from a small fixed‐wing UAV begins by planning a path that enables large sets of images to be taken orthogonal to the ground to cover the entire region of interest (ROI), which is usually represented by a user‐defined polygon (Colomina & Molina, 2014 ). For a simple convex field, the most common technique is a simple back and forth sweep motion known as a boustrophedon path. These images are then stitched together using structure‐from‐motion algorithms to generate a single high‐resolution image, also known as an “orthophoto.” A digital elevation model (DEM) and/or nonvisible spectral images of the ROI could also be produced using one of a number of commercially available photogrammetry software packages (Gini et al., 2013 ). These can then be used to extract actionable information on the surveyed field.

UAV remote sensing can be performed by either a rotary or fixed‐wing UAV, equipped with an imaging system which can be fixed directly to the airframe, or mounted on a gimbal (Anderson, 2014 ). There is a wide range of imaging systems that can be used depending on the application, including visual, multispectral, or hyperspectral systems which sample over a wide range of electromagnetic frequencies. Rotary wing UAVs are ideal for low altitude, slow flights, providing high‐resolution images but have very low endurance. Depending on the platform, payload, and batteries, flight times of around 15–20 min are typical. This means that only a small area can be imaged. Alternatively, fixed‐wing UAVs tend to fly higher and faster, resulting in lower resolution imagery, but they have much higher endurance and can survey a much larger area in a single flight. This paper will deal only with fixed‐wing remote sensing due to the focus of this paper on large aerial surveys.

Precision agriculture is a technique based on observing, measuring, and responding to inter‐ and intra‐field variability in crops, with the aim of maximizing returns and/or minimizing costs. The types of data that are most often gathered for precision agriculture include weed & disease distribution, soil pH, crop yields, and various crop health indices. These data are traditionally gathered by hand, by agronomists, dividing a field into sections and sampling. This actionable agricultural data can be used to help increase yields, while lowering fertilizer, herbicide, and water use. However, farms can be vast (the average US farm is 175 ha; United States Department of Agriculture, 2016 ) so gathering any data by hand can be prohibitively time‐consuming and expensive (Rew & Cousens, 2001 ). This necessitates a more effective way of agricultural surveying, such as using UAVs. An example would be site‐specific weed management (SSWM) (Zhang & Kovacs, 2012 ); where weed maps are generated from remote sensing aerial images and herbicides are more efficiently applied to only the affected areas of the field. Other examples include disease mapping (Fornace, Drakeley, William, Espino, & Cox, 2014 ) and yield prediction (Bendig et al., 2014 ).

Both commercial and hobbyist unmanned aerial vehicle (UAV) markets have grown significantly over the last few years. The worldwide market for agricultural drones was $494 million in 2016 and is anticipated to reach $3.69 billion by 2022, suggesting that agricultural drones will likely dominate all other drone sectors (WinterGreen Research Inc., 2016 ). They have reduced in price and have matured in both safety and reliability. UAV‐based remote sensing has been shown to be effective across a huge range of applications with growing popularity, for example, forestry and agriculture (Grenzdörffer, Engel, & Teichert, 2008 ) and coastal and environmental remote sensing (Klemas, 2015 ). This is especially true in the field of precision agriculture (Anderson, 2014 ) due to their low cost, greater flexibility, and better spatial and temporal resolution compared with more traditional methods such as manned aircraft or satellite remote sensing (Whitehead & Hugenholtz, 2014 ). Therefore, efforts to improve the safety, functionality, and efficiency of agricultural drones, and their algorithms, are essential.

Two alternative cost functions have been chosen from the literature for comparison. These are the NT (Equation 9 ; Li et al.,) and the MSA (Equation 10 ; Huang,) algorithms previously mentioned in Section 1.1 whereis the total number of sweeps.

Example boustrophedon path over a decomposed cell, showing the formulation of EFZ from the turn zone and corner‐cutting zone. As a part of the EFZ intersects with the NFZs this is an unacceptable decomposition. EFZ, extended flight zone; NFZ, no‐fly zone; ROI, region of interest [Color figure can be viewed at wileyonlinelibrary.com]

Shown in Figure 7 is an example of EFZ of a single convex survey cell from an already decomposed ROI. Here the NFZ is shown to be power pylons that come very close to the boundary of the ROI. Shown in the exploded view is how each edge is shifted outward from the original polygon. This creates an additional green area which is where the turns are performed. The areas in yellow are portions of the flight that the aircraft will overfly areas outside of the ROI due to corner cutting. In the case of this example, the pylon NFZ intersects with the EFZ so the cost J Ω will be infinite and another solution will need to be used, as shown in Equation 7 .

To calculate if a NFZ will be violated, an area called the extended flight zone (EFZ) is created. The EFZ is a polygon that represents the total area that the aircraft will overfly for each decomposed cell, inclusive of the modified Dubins paths for turning and traversal. Each EFZ is represented by polygon P E , k with the same index k as its parent P C , k . Each edge, E k j , of each P E , k is shifted outside P C , k by the turn radius of the aircraft where R t = V ψ ˙ at the survey angle ψ s . To approximate the area that the aircraft will turn in, we assume zero wind for the turn radius.

If it is assumed that there are no NFZs within the operator‐defined ROI, then there are three instances where the survey could violate this rule: (a) overflying an NFZ during cell traversal, (b) overflying NFZ while turning between sweeps, and (c) overflying an NFZ while corner cutting. Corner cutting is a technique used in the decomposition phase to find simpler faster decompositions, this will be detailed in Section 3 . If any of these cases occur then a hard constraint is imposed by adding an infinite cost to the cost function.

Finally, while the overall aim is to minimize the flight time, there may be other operational factors that need to be considered. To demonstrate a method for their inclusion, NFZs have been included. These are areas that the survey is prohibited to overfly or perform turns over. This may be due to structures or personnel not under control by the operator: In this case, the aircraft is not to come within 50 m of these areas (Civil Aviation Authority, 2012 ).

An example concave polygon decomposition with boustrophedon paths on each cell. Each boustrophedon path has four possible start and endpoints, this figure shows the fastest route to traverse each cell to and from the aircraft start location.(Coombes, Fletcher, et al.,). ROI, region of interest [Color figure can be viewed at wileyonlinelibrary.com]

The traversal path is complicated by the fact that each cell ( c k ) has four corners that the aircraft can start its survey from c n k , where n ∈ [ 1 , 2 , 3 , 4 ] (Note that the entry corner will also define the exit corner depending on whether there is an even or odd number of sweeps). The algorithm must, therefore, decide on how best to traverse a given set of decomposed cells. Shown in Figure 6 is an example ROI that has been decomposed into three convex polygonal cells. In this example, there are 3 ! x 4 3 different possible routes to take between the start and endpoint. To find the fastest one, a graph search is performed using Dijkstra's algorithm. Also shown in Figure 6 is the fastest path which has ensured that the transversal path is minimized.

Previous work by the authors (Coombes, Fletcher, et al., 2018 ) assumed that the sum of the flight times for the individual cells could be used as an estimate of the overall flight time. This was done to meet the dynamic programming requirement for optimal substructure. However, it has since been found that in some circumstances this algorithm tends to slightly over‐decompose the ROI due to the lack of a penalty for doing so. As a result, one of the developments for this study has been to include the transfer paths from the launch site, between each decomposed cell, and back to the landing site into the cost function. These are generated as modified Dubins paths in the same way as the turn portion of the boustrophedon paths, as described in Section 2.1.2 .

The other three (#2–#4) wind survey flights gave similarly consistent results as shown in Table 1 . Weaker wind strength leads to slightly more variation in wind direction, however as the wind is weaker, this will have less effect on the aircraft. Using the same simulation techniques described in Section 4 , it has been calculated that the variation in wind observed will result in less than 3% variation in the flight time for a typical mission.

The wind field is measured at the same survey site where the accuracy assessment flights in Section 4 are conducted, however, on four different days. The average measurements across the first wind field measurement flights were 7.3 m/s and 230.2°. Shown as (#1) in Table 1 . This closely matched the wind forecast information (225° at 7 m/s) which was taken from a meteorological station located 5 miles north of the test site. The wind strength and direction remained quite consistent across the whole survey area and across the 14 min flight time, with standard deviations of 0.627 m/s and 10.72°. The wind field with the S1000 flight path is shown on Figure 5 .

To achieve this, a sweep pattern was flown over the survey area, stopping and hovering at 16 separate way‐points and sampling wind for 15 s at each. This was repeated on four different days. The aircraft used was a DJI S1000 octorotor, with an FT205EV ultrasonic wind sensor mounted above the vehicle out of the rotor wash, as shown in Figure 4 . This allows direct measurement of the wind at any altitude. It has a 15 min flight time using a 520 Wh battery, and a gross weight of 5 kg. Using a Raspberry Pi 2 running robot operating system (ROS), it interfaces with the autopilot and wind sensor to log wind speed and direction, position, heading, and altitude data.

The flight‐time calculations for both the sweep and turn portions of flight assume a steady uniform wind over the relatively short time and area of the survey. If this were to be incorrect, this could significantly affect the time prediction accuracy. To show that the assumption does hold, the wind field (wind at multiple locations) was measured at a typical survey altitude of 100 m by a multirotor UAV. The measured wind field was assessed for spacial and temporal variations. Note that this step would not be necessary in a real‐world survey, it has been included here to test the assumption of wind prediction.

The turns joining two sweeps are represented by Dubins paths modified to account for wind. This algorithm is laid out in detail in Coombes et al. (), and as such, will only have a brief overview here. A circular coordinated turn in a steady uniform moving body of air causes the path flown by aircraft relative to the ground to be trochoidal. Each turn is made up of three continuous connected paths; two trochoidal turns paths linked by a straight path tangential to both trochoids. The time taken for this full turn maneuverconsists of the time taken to fly each trochoidal turn, and the straight tangential path between them, shown in Figure 3 and Equation 7 . For further detail about the derivations of these terms, please refer to Coombes et al. () and Techy and Woolsey ().

The time to fly along a single sweep line () is calculated as the length of the sweep () divided by the aircraft's ground speed (). If a steady uniform wind field across the survey is assumed, the aircraft's ground‐speed can be calculated by wind triangle vector subtraction; this is standard book work but is illustrated in Richards (). Ground‐speed is calculated for each sweep in Equation 5 , whererepresents the vehicle's airspeed,represents the wind speed and the direction of the flight path () is either equal to the sweep direction for outward legs, or opposite to it for return legs, see Equation 6

The complete single‐cell path consists of two different states of flight; the straight sweep paths where the images are captured, and the trochoidal turn maneuvers used to transition between sweeps. The flight time for a single cell can be calculated as the sum of the time spent during each of these states, see Equation 4

This section describes the method used to calculate the flight time in wind cost () for such a set of decomposed cells, which is to be used for the recursive optimization process in Section 3 . The cost function is made up of two parts; the flight time () and an additional cost () representing a hard constraint to prevent overflying of NFZs. The total flight time for a set of cells made up of the time to fly; (a) the launch (ascent) path (), (b) the sum of the times to fly each individual cell using boustrophedon paths, (c) the total time required to traverse between the decomposed cells, and finally (d) the landing (descent) path (). See Equation 3 , whereis the number of decomposed cells.

Although it may be possible to overfly the complete convex hull of the ROI, this may not be efficient, or desirable. In this case, the ROI can be decomposed into a number of convex polygonal cells that can be surveyed individually and, when combined, provide total coverage of the ROI. Figure 2 shows an example of complex polygonal decomposition into three convex polygonal cells, with a boustrophedon path displayed for one of these.

The GSD represents the physical size of a pixel in an image and is measured in cm/pixel, different applications have different GSD requirements. So the height above ground of the survey () can be set based camera parameters and the desired GSD as shown in Equation 2

The most efficient path to ensure complete image coverage of a single convex polygon by a fixed‐wing UAV is a boustrophedon path (Li et al.,), which consists of back and forth sweeps across the polygonal ROI. To minimize the number of turns, the survey angle () is fixed based on the angle of the long axis of the polygon. The distance between sweeps () is set as shown in Equation 1 . It is based on the horizontal angular field of view (FOV) of the sensor (), the number of pixels the sensor has in the horizontal direction (), the required ground sample distance (GSD) () and the lateral overlap ratio (or sidelap,) is defined by the user based on the application. A higherof around 0.6–0.8 is needed for DEM or a lower value of around 0.2–0.4 for orthophotos.

Algorithm 2 explains the polygon rotation angle optimization using a brute force search. A set of angles ( ψ r ) is defined (line 1) and the optimal cost ( J λ ) is initialized as infinity (line 2). Using a loop (line 3), the initial polygon ( P i ) is rotated by the angle to be checked ( ψ r ) (line 4) and the decomposition algorithm is called (line 5) to calculate the optimal cost ( J λ ′ ) and set of decomposed convex cells ( P λ ′ ) for the given angle of rotation. The decomposed cells are rotated back to the original frame of reference (line 6). Finally, if the optimal cost at the current angle is lower than the current overall optimal cost ( J λ ), then J λ and the corresponding set of optimized cells ( P ψ λ ) are updated.

For each polygon rotation angle ( ψ r ), there will be an optimal set of polygons ( P λ ′ ) which will have an associated cost ( J λ ′ ) selected by the recursive greedy function. The rotation angle which produces the lowest overall cost ( J λ ), is output as a set of polygons ( P ψ λ ) after re‐rotation back to the original angle. Note that only angles between 0 ≤ ψ r < π are required because π ≤ ψ r < 2 π will produce identical results. Also note that the wind direction must also be rotated for use in the cost function to keep it in the same relative angle.

As previously mentioned, the algorithm attempts polygon splits at each concave vertex. It achieves this by using a simplified version of the trapezoidal decomposition algorithm found in Seidel ( 1991 ). As explained in Coombes, Fletcher, et al. ( 2018 ), the long thin trapezoids often produced, lend themselves very well to efficient fixed‐wing surveys. The simplification of the algorithm used here is, that instead of splitting the polygon at every vertex to make n v − 2 (where n v is number of vertices) trapezoids, it only makes the split at a single concave vertex at each step (line 18 of Algorithm 1). As trapezoidal decomposition splits along vertical sweep lines, this often leads to decomposed cells that have similar ψ s angles. If the vertical sweep line is aligned with the wind, this could lead to longer flight times, as the UAV will be flying directly into and away from the wind. However, if the ROI is rotated, the dominant ψ s across multiple cells can be changed. Due to the efficiency of the recursive greedy search (RGS) algorithm, this is used to enlarge the search space by running the recursive function for a range of initial polygon rotations ( ψ r ), the number of which is dictated by a user‐selected positive integer number n r . This can be run using parallel computing techniques, even on a portable computer.

The order of the function can be approximated using Equation 11 , whererepresents the number of concave vertices. This represents a significant improvement over the dynamic programming (DP) solution from our previous work (Coombes, Fletcher, et al.,) which was approximated as, whererepresented the number of decomposed cells, a much larger number for a given ROI.

The algorithm works using a “top‐down” architecture; calculating the cost of the convex hull and iteratively splitting concave regions, one concave vertex at a time until no concave cells remain. The pseudo‐code for this function can be found in Algorithm 1. The steps of the algorithm can be summarized as follows;

Steps 3–5 are repeated until each of the concave vertices have been examined (line 21). The function outputs the final value of the current best cost ( J λ ′ ) and its associated set of decomposed convex polygons ( P λ ′ , line 37).

The sum of the costs of the subregions ( J v ) identified in Step 3 is calculated (lines 24 and 28). If this is less than the current best cost value ( J λ ′ ), then the current best cost ( J λ ′ , line 32) and optimal decomposition set ( P λ ′ , line 33) are updated.

The cost of each subregion ( P s k ) is calculated directly using the relevant cost function if the subregion is convex (line 23). If the subregion is concave, the function recurses using the concave subregion as its input, which will output the optimal cost of that concave subregion (line 26).

Trapezoidal decomposition is used to split the original polygon at the first of the remaining concave vertices ( V c j ) . This will result in two or more polygons (aka “subregions,” P s ) being generated (line 18).

First, the convex hull ( P h ) of the ROI polygon ( P ψ ) is determined and the relevant cost function (Equation 3-9 , or 10 ) is used to calculate the associated cost (line 3). This cost is used to initialize the current best cost value ( J λ ′ , line 5) and the optimal set of decomposed polygons ( P λ ′ , line 4).

Example recursive‐greedy decomposition search graph for a polygon with two concave vertices. The decomposition algorithm splits the cell using a vertical line that intersects with each concave vertex in turn. Note that steps ③ and ⑤ result in three decomposed cells (rather than two) due to the angle and shape of the original polygon [Color figure can be viewed at wileyonlinelibrary.com]

The proposed method uses trapezoidal decomposition to break a concave irregular polygon into two (or more) subregions, known as “cells.” At the same time, the convex hull of the original polygon is calculated and the optimizer will select the option which produces the smallest cost. Because each of the two cells may or may not be convex, the algorithm works recursively until every cell included in the cost calculation is convex. Figure 8 demonstrates how the decomposition is performed. The algorithm is “greedy” in that if, at any time, the calculated flight time of the convex hull of any cell is found to be more efficient than the decomposition, then this option is selected. Although this does not guarantee an optimal solution, it ensures that a solution is found quickly, even for highly irregular polygons with a large number of concave vertices. This is an important operational consideration because this algorithm is designed to be used in the field using portable computers with limited computational power.

To use boustrophedon paths to survey an ROI described by an irregular concave polygon, it must first be converted into one or a number of convex polygons. Depending on the shape of the ROI, it may be more efficient to either split the area into multiple convex cells and survey these separately, or to simply calculate the convex hull of the concave polygon. This section describes the optimization algorithm used to make this decision so that one of the cost functions (Equation 3-9 , or 10 ) laid out in the previous section can be minimized.

4 SIMULATION RESULTS AND DISCUSSION

This section presents simulated mission plans for a number of complex real‐world example fields and randomly generated field‐style polygon shapes to demonstrate the improvements that the proposed algorithm offers over previous cost functions. Three complex, concave, agricultural fields will be used as examples, which are each shown in Figure 9. The launch/landing locations are from feasible accessible locations and are also marked in Figure 9. All of these field polygons have a large recess which is a notoriously hard feature to plan an efficient mission around. This is because, when decomposed, they tend to generate a high number of convex polygons resulting in a large search space for path optimization (Coombes, Fletcher, et al., 2018). The example fields have a relatively high complexity (11–16 vertices) and a wide range of areas (35–129 ha) for this type of problem.

Figure 9 Open in figure viewer PowerPoint Example complex concave field survey polygons used to explore the effectiveness of the proposed algorithm in simulation (a) 11 vertex arable field polygon 34.89 ha, with adjacent NFZ in red; (b) 15 vertex arable field polygon 47.86 ha; (c) 16 vertex arable field polygon 129 ha. NFZ, no‐fly zone [Color figure can be viewed at wileyonlinelibrary.com]

The aircraft simulated is a senseFly eBee; a commercially available and widely used survey UAV. The eBee is a relatively low‐velocity fixed‐wing aircraft. The imaging system is a MicaSense RedEdge multispectral camera. Only an orthophoto is required, therefore, the requirements for image overlap/sidelap have been set to 30%. The RedEdge has a narrow FOV of 47.2°, so to get the required GSD of 4.3 cm/pixel, the survey is conducted at an altitude of 70 m. From Equation 1, this makes the distance between the sweeps d x to be 46.2 m. To show how the proposed method handles the presence of wind, a steady uniform southerly wind ( ψ w = 0 ∘ ) of 5 m/s is also applied to the simulation. The aircraft and survey parameters are summarized below in Table 2.

Table 2. Simulated senseFly eBee and survey parameters Parameter Value Wind 0 ∘ /3 m/s V 7 m/s ψ ˙ 0.7 rad/s w s , w o 0.3, 0.3 h 70 m GSD 4.3 cm/pixel D x 46.2 m

The polygon decomposition and mission planner algorithms laid out in Section 3 are applied to each of the example polygons. To compare methods, three missions per polygon will be planned, one for each of the cost function options; FTIW (Equation 3—proposed algorithm), MSA, and NT (Equation 9 and 10, both previously discussed in Sections 1 and 2.5). As previously discussed, the FTIW cost function has two main advantages over the cost functions found in the literature. First, the FTIW cost takes into account the effect of wind on the flight time. This is expected to cause the suggested sweep direction to be aligned more perpendicular to the wind, which has already been found to be beneficial (Coombes, Fletcher, et al., 2018) for real‐word survey efficiency. Second, the proposed algorithm may choose to survey the convex hull of a concave polygon rather than decomposing further (“corner cutting”). This is expected to result in fewer decomposed cells and hence a reduction in the transition time between them and will be selected by the optimizer only when it is beneficial overall.

The flight‐times for each field, and for each cost function are shown in Table 3. It can be seen that the FTIW algorithm outperforms both cost functions from the literature by between 5% and 25% in all three examples. This is to be expected because these examples have been selected to showcase the benefits of the proposed algorithm. A fairer numerical comparison is presented later in the paper using Monte‐Carlo simulation of randomly generated polygons.

Table 3. Calculated time comparisons between FTIW, NT, and MSA cost functions for three example fields Simulated flight times FTIW % improvement vs. # Polygon FTIW (s) NT (s) MSA (s) NT (%) MSA (%) 1 11 Vertices 1,552 1,631 1,672 4.84 7.18 1 11 Vertices with NFZ 1,603 – – – – 2 15 Vertices 1,121 1,182 1,436 5.16 21.9 2 15 Vertices 9 0 ∘ Wind 1,077 1,234 1,424 12.7 24.37 3 16 Vertices 2,953 3,277 3,342 9.89 11.64

Examining the results in more detail, Figure 10 shows the optimized mission plan for each cost function for the 11‐vertex field from Figure 9a. It can be seen that the MSA and NT paths use six and eight decomposed cells, respectively, but the FTIW output uses just four. This is because the FTIW algorithm is able to take advantage of the ability to cut corners rather than relying on decomposition. As a result, the aircraft will spend less time transitioning between cells where no surveying is taking place and is the primary reason for the 4.8% and 7.18% improvements over the NT and MSA cost functions, respectively. While the proposed algorithm will take advantage of prior knowledge of the wind, in this case, the proposed method cannot find an efficient low‐cell‐count decomposition perpendicular to the wind. This is simply due to most of the polygons' edges being aligned close to 4 5 ∘ to the wind direction.

Figure 10 Open in figure viewer PowerPoint Field polygon from Figure 9 a decomposed and flight paths generated using FTIW, MSA, and NT methods. Predicted flight times are shown. FTIW, flight time in wind; MSA, minimal sum of altitude; NT, number of turns [Color figure can be viewed at wileyonlinelibrary.com]

The second example from Figure 9b is more complex with 15 vertices, five of which are concave. For this field, mission plans have been generated for both southerly and westerly winds ( ψ w = 0 , 90 ). Their solutions are shown in Figure 11. For the southerly wind, the FTIW solution shows lower flight times than both NT and MSA algorithms, although the advantage over the MSA (21.9%) is much more pronounced than that over the NT algorithm (5.2%). This is because many of the polygon long axes are already perpendicular to the wind for the NT solution. The MSA algorithm, on the other hand, not only generates more cells but also flies largely parallel to the wind direction.

Figure 11 Open in figure viewer PowerPoint Field polygon from Figure 9 b decomposed and flight paths generated using FTIW, MSA, and NT methods. This is done for both southerly and westerly winds. Predicted flight times are shown. FTIW, flight time in wind; MSA, minimal sum of altitude; NT, number of turns [Color figure can be viewed at wileyonlinelibrary.com]

For comparison, Figure 11 also shows the flight paths for the three algorithms for a westerly wind, in which case the advantage of accounting for the wind becomes very obvious. It can be seen that the flight paths for the NT and MSA cost functions are identical because the do not take into account the wind in the planning, but the flight times have changed. Note that the flight time for the NT algorithm has increased due to the fact that the dominant sweep direction has gone from perpendicular to parallel to the wind direction. In contrast, the flight time for the MSA algorithm has reduced for the opposite reason. The FTIW algorithm, however, has accounted for the wind and changed the dominant sweep direction, and in this case, the convex hull of the entire ROI is found to be faster than decomposition resulting in a substantial 24.4% improvement over MSA and a 12% improvement over NT.

The final example, shown in Figure 9c, has similar complexity to the previous example (16 vertices and five concaves) but is a much larger field with an area of 129 ha. Shown in Figure 12 are the decomposition solutions and flight paths for each cost function. Using the FTIW algorithm gives the fastest flight time at 2,953 s, approximately 10% faster than the other algorithms. Contrary to the other two examples, the advantage of the FTIW algorithm comes primarily from flying perpendicular to the wind. In this case, cutting corners is less beneficial because it requires unnecessarily surveying much more land area. The magnitudes of the times are also significant in this case because the maximum flight time of the senseFly eBee, as configured, is approximately 50 min (3,000 s). Therefore, the 10% flight time reduction from around 3,300 s to just below 3,000 s could make the difference as to whether the survey could be performed in a single flight.

Figure 12 Open in figure viewer PowerPoint Field polygon from Figure 9 c decomposed and flight paths generated using FTIW, MSA, and NT methods. Predicted flight times shown. FTIW, flight time in wind; MSA, minimal sum of altitude; NT, number of turns [Color figure can be viewed at wileyonlinelibrary.com]

4.1 NFZ example Section 2.4 showed how violating NFZs can be avoided by adding a hard constraint to any decomposition which results in overflying an NFZ. The 11‐vertex field from Figure 9 is used to illustrate this, where the red area in Figure 9a is the NFZ containing uncontrolled buildings. Previously, the NFZ was not taken into account, and the NFZ was violated by a number of turns at the north‐east end of cell 3. Figure 13 shows a decomposition where the NFZ is no longer violated. This was achieved by flying the entire north‐east edge at an angle ( ψ r ) of 13 5 ∘ . As seen in Table 3, this solution is around 50 s slower than the previous solution but is still faster than both algorithms from the literature that do not account for the NFZ. Figure 13 Open in figure viewer PowerPoint Field polygon from Figure 9 a decomposed and flight paths generated using the FTIW and NFZ cost functions to plan a fast route in wind while avoiding NFZs. FTIW, flight time in wind; NFZ, no‐fly zone [Color figure can be viewed at wileyonlinelibrary.com]

4.2 Computational complexity One of the major advantages of this suboptimal method is that it runs much faster than the optimal method from Coombes, Fletcher, et al. (2018), which uses a full initial trapezoidal decomposition in its “bottom‐up” dynamic programming (DP) approach. As the number of initial trapezoidal decomposed cells rises, the computational complexity for DP grows exponentially; making it practically infeasible to find a mission plan for a field polygon with more than 18 initial cells. However, as the proposed method searches the solution space very differently to the DP method, it is hard to compare their computational time complexity in a rigorous way. As DP tries all combinations of each initial decomposed cell and attempts to merge them, the computation time increases with the initial number of cells, which is based on a number of polygon vertices N − 2 . However, the proposed method recursively slices cells apart at concave vertices, therefore, its time complexity increases with a greater number of concave vertices. To show how computational complexity increases with polygon complexity, one of 20 random polygons has its polygon vertex number reduced to just four vertices (a single trapezoid), then its optimization runtime is recorded for both RGS and DP methods. Then vertices are added back one at a time and the runtime recording is repeated until the mission for the entire polygon has been made. This is then repeated for each of the 20 polygons and an average runtime is calculated for each vertex count. The results of this computation time comparison is shown in Figure 14. In the simpler runs with a vertex count <17, the DP method is faster due to the less computationally intensive cost function. However, as the vertex count crosses 17 cells, the proposed technique becomes much faster; quickly becoming orders of magnitude faster than the DP method. For example, at 20 vertices, the DP method took 1,720 times longer than RGS and it was not feasible to run the DP algorithm on polygons with a vertex count greater than this due to its computational time. Figure 14 Open in figure viewer PowerPoint Average computation time for both DP and RGS methods against a number of polygon vertices in the field to plan a mission for. DP, dynamic programming; RGS, recursive greedy search; ROI, region of interest [Color figure can be viewed at wileyonlinelibrary.com]

4.3 Monte‐Carlo simulation To present a fair numerical comparison between the proposed method and those from the literature, 30 polygons have been randomly generated between 11 and 16 vertices. The improvement offered by using FTIW over NT and MSA is assessed. In this case, the platform used is the popular Skywalker X8, a flying‐wing aircraft that flies at a higher airspeed of 16 m/s (leading to higher turn radii). This has been chosen because it matches the aircraft which is used for flight testing in Section 5. The GSD requirements have been relaxed for this flight to 8.2 cm/pixel; meaning that the survey can be conducted at a higher flight altitude of 120 m for logistical reasons. The full aircraft and survey parameters are listed in Table 4. Table 4. Simulated X8 and survey parameters Parameter Value Wind 27 0 ∘ /5 m/s I 15.5 m/s ψ ˙ 0.7 rad/s w s , w o 0.3, 0.3 h 120 m GSD 8.2 cm/pixel D x 73 m Using identical sets of 30 random polygons, the simulation is run with five different wind strengths that range from 0 to 13 m/s in an easterly direction ( ψ w = 27 0 ∘ ). This is a large variation, which represents between 0% and 81% of the aircraft's airspeed. The results of these Monte‐Carlo simulations are shown in Table 5, and a few example polygons and the generated paths are shown in Figure 15. It can be seen that under all wind conditions, using the proposed method gives a significant flight time advantage. When the wind speed is increased, the advantage of the FTIW cost function becomes even greater. Table 5. Monte‐Carlo simulation results—30 random polygons runs—FTIW comparison to MSA and NT FTIW improvement over NT FTIW improvement over MSA # Description (m/s) Av time (s) Av percentage (%) Av time (s) Av percentage (%) 1 V w = 0 192 9.3 172 8.6 2 V w = 3 201 9.5 177 8.5 3 V w = 5 210 9.6 184 8.4 4 V w = 10 399 13.6 346 11.2 5 V w = 13 846 19 749.6 15.4 Figure 15 Open in figure viewer PowerPoint Three example random polygons from the 30 Monte‐Carlo simulations with flight paths generated for RGS with the FTIW cost function. FTIW, flight time in wind; RGS, recursive greedy search [Color figure can be viewed at wileyonlinelibrary.com] The first run has zero wind speed, which means that the advantage provided by corner cutting is isolated. It can be seen that this has a significant affect on the results, providing a 9.3% reduction in flight time compared to the NT algorithm and an 8.6% reduction compared to the MSA cost function. As the wind speed increases, the flight time advantage increases further. At low wind speeds (5 m/s, 30% of airspeed), the FTIW advantage increases marginally up to 9.6% compared to the NT cost function and reduces marginally down to 8.4% compared to the MSA algorithm. However, for high wind speeds (13 m/s, 81% of airspeed), the advantage of the prior wind knowledge is much more significant, increasing to 19% and 15.4% compared to NT and MSA, respectively. A weak correlation has also been found between the complexity of the polygon and the FTIW advantage as shown in Figure 16. With an increase in the number of vertices and especially concave vertices, the polygon becomes more complex. Because the alternative methods do not allow for corner cutting, their pathways tend to be based on a higher number of decomposed cells, increasing the path complexity and the transit time between cells. Therefore, it is concluded that with an increase in either polygon complexity or an increase in wind speed, the advantages of using the proposed method is maximized. Figure 16 Open in figure viewer PowerPoint Correlation between the number of polygon vertices and flight time reduction using FTIW cost function with corner cutting, against using NT cost function without corner cutting. FTIW, flight time in wind; NT, number of turns [Color figure can be viewed at wileyonlinelibrary.com]