Revisiting Boustrophedon Coverage Path Planning as a Generalized Traveling Salesman Problem

TL;DR

This paper introduces an obstacle-aware coverage path planner for MAVs that optimizes sweep paths in polygonal environments, improving efficiency and applicability in real-world, obstructed terrains.

Contribution

It extends boustrophedon coverage planning by incorporating obstacle considerations and sweep optimization, providing an open source framework for terrain coverage in complex environments.

Findings

Achieves 14% lower path costs than conventional planners.

Successfully tested on 320 synthetic maps and real terrain.

Operates efficiently enough for field deployment.

Abstract

In this paper, we present a path planner for low-altitude terrain coverage in known environments with unmanned rotary-wing micro aerial vehicles (MAVs). Airborne systems can assist humanitarian demining by surveying suspected hazardous areas (SHAs) with cameras, ground-penetrating synthetic aperture radar (GPSAR), and metal detectors. Most available coverage planner implementations for MAVs do not consider obstacles and thus cannot be deployed in obstructed environments. We describe an open source framework to perform coverage planning in polygon flight corridors with obstacles. Our planner extends boustrophedon coverage planning by optimizing over different sweep combinations to find the optimal sweep path, and considers obstacles during transition flights between cells. We evaluate the path planner on 320 synthetic maps and show that it is able to solve realistic planning instances fast enough to run in the field. The planner achieves 14% lower path costs than a conventional coverage planner. We validate the planner on a real platform where we show low-altitude coverage over a sloped terrain with trees.

Tables

Table 1. Table 1: Detailed computation times for the instance in Figure 2 with 15 holes ( 52 52 52 hole vertices). Creating all edges is the greatest computational effort. The reduced problem thus has the smallest computation time. Our planner results in the lowest path cost.

Graph elements

Computation time (s)

Path cost (s)

Cells

Nodes

Edges

Cells

Sweeps

Nodes

Pruning

Edges

Solve

our_bcd

$52$

$254$

$\mathrm{63\hspace{0.17em}056}$

$\mathrm{5.680\hspace{0.17em}28}$

$\mathrm{1.100\hspace{0.17em}210}$

$\mathrm{2.984\hspace{0.17em}410}$

$\mathrm{2.392\hspace{0.17em}840}$

$22.2329$

$\mathrm{1.195\hspace{0.17em}960}$

$2391.21$

our_tcd

$87$

$512$

$\mathrm{258\hspace{0.17em}704}$

$\mathrm{6.920\hspace{0.17em}08}$

$\mathrm{1.107\hspace{0.17em}450}$

$\mathrm{4.940\hspace{0.17em}300}$

$\mathrm{3.679\hspace{0.17em}660}$

$84.6797$

$\mathrm{6.555\hspace{0.17em}060}$

$2879.15$

one_dir_gk

$52$

$52$

$2652$

$\mathrm{5.535\hspace{0.17em}08}$

$\mathrm{0.089\hspace{0.17em}057}$

$\mathrm{0.246\hspace{0.17em}038}$

$\mathrm{0.026\hspace{0.17em}695}$

$0.8067$

$\mathrm{0.513\hspace{0.17em}497}$

$2664.98$

Table 2. Table 2: Flight experiments parameters.

Area	Flight Time	\ac AGL	Sweep Offset	$v_{\max }$	$a_{\max }$
$1950m^2$	$1050s$	$4m$	$1.5m$	$3.0m{s}^{-1}$	$0.5m{s}^{-2}$

Equations

$$w$$

$$c_{ij}=t_{ij}+t_j$$

$$t$$

Full text

1 Introduction

\Acp

MAV such as the DJI M600 shown in Figure 1 present ideal platforms for supporting demining efforts and other critical remote sensing tasks. These commercially available platforms are fast to deploy and can collect useful data that are often inaccessible to other modes of sensing pajares2015overview . For demining applications, this means that \iacMAV carrying \iacSAR system can be flown over a minefield to collect in-situ, high precision measurements for locating buried landmines schartel2018uav .

At present, the target coverage area is $1\text{\,}\mathrm{ha}$ per battery charge. Furthermore, since the \acMAV must fly at relatively low altitudes to emit enough energy into the ground, structures in the environment such as trees or buildings can force \acpNFZ across the space. These obstacles can be especially problematic since they are often omitted from maps and must be dealt with at deployment time. This results in a need for an autonomous coverage planning and execution solution that can rapidly replan trajectories in the field.

Many existing commercial mission planners, e.g., Ardupilot Mission Planner111http://ardupilot.org/planner/, Pix4DCapture222https://www.pix4d.com/product/pix4dcapture, Drone Harmony333https://droneharmony.com, and DJIFlightPlanner444https://www.djiflightplanner.com/ provide implementations of 2D coverage planners which allow specifying polygon flight corridors and generating lawnmower patterns. However, these mission planners aim at high-altitude, open space, top-down photography and thus neither allow specifying obstacles within the polygon nor consider the polygon edges as strict boundaries. Furthermore, due to safety limitations of commercial platforms, most commercial planners enforce a minimum altitude which is above the required $3\text{\,}\mathrm{m}$ of our \acSAR system.

While solutions for generating coverage paths in general polygon environments exist, the implementations are either not publicly available, computationally unsuitable for complex large environments or lack the full system integration that would allow a user to freely designate the flight zone.

In this work we develop a complete pipeline for rapidly generating 2D coverage plans that can account for \acpNFZ. Our proposed solution takes as input a general polygon (potentially containing \acpNFZ) and performs an exact cell decomposition over the region. Initial sweep patterns are computed for each cell from which feasible flight trajectories are computed. Finally, the \acE-GTSP across the complete cell adjacency graph is solved to minimize the total path time.

Results comparing our proposed coverage planner to exhaustively searching a coverage adjacency graph demonstrate an order of magnitude speedup in the computation time for each trajectory. Specifically, we show that our solution’s computation time grows reasonably with increasing map complexity, while the exact solution grows unbounded. Furthermore, the flight path computed by our method is $14\text{\,}\mathrm{\char 37\relax}$ shorter than those found by solving a regular \acTSP while providing the same level of coverage.

The contributions of this work are:

• •

\Iac

E-GTSP based fast 2D coverage planning algorithm that allows specifying polygonal flight zones and obstacles.

• •

Benchmarks against existing solutions.

• •

An open source \acROS implementation.555https://github.com/ethz-asl/polygon_coverage_planning

We organize the remainder of the paper as follows. Section 2 presents related work. In section 3 we present the background in computational geometry to generate the sweep permutations. In section 4 we present our coverage planner formulation to solve for the optimal coverage pattern. Finally, we validate our method in section 5 before we close with concluding remarks in section 6.

2 Related Work

Coverage planning in partially known environments is an omnipresent problem in robotics and includes applications such as agriculture, photogrammetry and search. It is closely related to the \acAGP and \acTSP, which are NP-hard o1983some ; papadimitriou1977euclidean , and for which heuristic solutions have been developed to cope with high dimensional problems. The coverage planning problem can be stated as: given a specified region and sensor model, generate a plan for a mobile robot that provides complete coverage, respects motion and spatial constraints, and minimizes some cost metric (often path length or time). The underlying algorithms to solve sweep planning in a 2D polygon map are usually either based on approximate or exact cellular decomposition. The state of the art in robotic coverage planning is summarized in galceran2013survey ; cabreira2019survey .

Early work in robot coverage planning developed approximate cellular decomposition methods. In zelinsky1993planning , the target region is first decomposed into connected uniform grid cells such that coverage is achieved by visiting all cells. The authors presented a coverage wave front algorithm that obeyed the given starting and goal cells. Unfortunately, grid based methods grow linearly in memory with the area to be covered and exponentially in finding an optimal path, making them unsuitable for large areas with small sensor footprints. Furthermore, guiding the wave front algorithm to pick a route that meets all mission requirements involves designing and tuning a specific cost function for each application. Spanning-trees are an alternative method to solve this problem gabriely2001spanning . However, by default those result in a rather large number of turns which is undesirable for \acpMAV.

A simple approach stems from the idea that a polygonal region can be covered using sequential parallel back-and-forth motions, i.e. lawnmower or boustrophedon paths. This approach for 2D top-down lawnmower patterns in polygonal maps dominates commercial aerial remote sensing. The user defines a polygonal target region and a specified sweep direction, and the planner generates a path consisting of a sequence of equally-spaced sweep lines connected by turns. The generated coverage paths are time-efficient, complete, and intuitive to the user. However, these planners cannot account for \acpNFZ nor do they satisfy our targeted flight altitude requirements. Furthermore, except for Ardupilot Mission Planner, these commercial systems do not allow modification as they are closed source.

Exact cellular decomposition is a geometric approach that can handle \acpNFZ by dividing a configuration space into simpler component cells, whose union is the complete free space. A full coverage path can then be found by creating a boustrophedon path in each cell and connecting all cells. The individual cell coverage plans are connected in an adjacency graph and solving the resulting \acTSP solves the coverage problem. Compared to grid-based methods, exact decomposition generally results in significantly fewer cells in simple large environments and thus a smaller \acTSP. Choset and Pignon choset1998coverage present the standard solution in 2D environments, which has been adapted by many planners for robot coverage galceran2015coverage ; bahnemann2017decentralized . Several improvements have also been proposed to minimize the number of cell traversals and the number of turns along the coverage path. Namely, methods for optimizing over the sweep line direction rather than using a fixed direction for all cells have been proposed li2011coverage ; torres2016coverage . However, these come at the cost of increased complexity when solving for the optimal path through the cell-connectivity graph.

In our work we revisit the exact cellular decomposition method and also generate multiple possible sweep directions per cell, however we formulate the resulting search as \iacE-GTSP, allowing us to use a state of the art genetic solver that handles significantly larger problem sizes. The \acE-GTSP, also known as \acTSP with neighborhoods, is a generalization of the classical \acTSP mitchell2000geometric . The goal is to find a shortest tour that visits exactly one node in each of a set of $k$ neighborhoods. As first proposed by Waanders waanders2011coverage , we cluster all possible sweep patterns of a cell in a neighborhood and search for the shortest path that includes exactly one sweep pattern per cell.

Bochkarev and Lewis bochkarev2016minimizing ; lewis2017semi took this \acE-GTSP formulation even further. Instead of using per-cell predefined sweep patterns, they precompute a set of globally good sweeping directions for each monotone cell and build a graph over all individual straight segments where a neighborhood is defined by traversing a segment in either direction. Compared to our approach their approach can give better solutions where the robot traverses from one cell to another cell without covering it completely first at the expense of a more complex \acE-GTSP and a predefined sweep direction.

To solve the \acE-GTSP, several options exist. Exact solvers like fischetti1997branch only work reliably for small problem sizes. Converting the problem into a directed graph using a product graph rice2012exact and solving it with an optimal graph solver, e.g. Dijkstra, basically falls back to solving it exhaustively as in choset1998coverage and torres2016coverage , which remains intractable for larger problem sizes. Thus a practical solution for our problem sizes is to use heuristic solvers. Helsgaun helsgaun2015solving transforms the \acE-GTSP into a \acTSP and uses an approximate \acTSP solver. The \acGK, on the other hand, is a memetic algorithm that approximately solves the \acE-GTSP directly and is faster but with reduced performance for large problems gutin2010memetic .

3 Geometric Path Generation

In order to solve the coverage path problem we follow the route of exact cellular decomposition techniques outlined in Figure 2. We decompose a \acPWH into cells, in a way that guarantees that each cell can be fully covered by simple boustrophedon paths. Our algorithm creates a permutation of sweeping directions for each cell and finds a shortest route that connects and covers every cell to define the coverage path.

3.1 Sweep Pattern Permutation

A sweep pattern describes the combination of straight segments and transition segments that covers a single polygon cell. A continuous parallel sweep pattern can be generated perpendicular to any monotone direction of a simple polygon. Without loss of generality we can consider the monotone direction the $y$-direction of the polygon. A polygon $P$ is considered $y$-monotone, if any line in the $x$-direction intersects $P$ at a single point, a segment or not at all. In other words the line intersects $P$ at most twice de1997computational .

Figure 3 (a) shows the straight segment generation in a $y$-monotone polygon. We initialize the first straight segment at the bottommost vertex parallel to the $x$-axis, which we refer to as sweep direction. In general, we restrict the sweep directions to be collinear to one of the edges of our polygon, as these directions have been proven to lead to a minimum number of straight segments to cover the polygon huang2001optimal . The individual straight segments are generated by alternating between intersecting a line in the $x$-direction with the polygon and offsetting the line from the bottommost towards the topmost vertex. The sensor footprint hereby defines the offset distance.

Based on this construction criterion, the straight segments can start clockwise or counter-clockwise at the bottom vertex and go from the bottom to the top or from top to bottom. Thus we generate four possible sweep patterns per sweepable direction as shown in Figure 4.

To connect two straight segments (and later to connect two sweep patterns) we calculate the Euclidean shortest path that avoids collision with the \acNFZ. The Euclidean shortest path in a \acPWH is computed along the reduced visibility graph latombe2012robot . Figure 3 (b) shows an example solution (red) using A*∗* to search the graph hart1968formal . The graph node set (circles) consists of all non-convex hull vertices and convex hole vertices.

3.2 Polygon Decomposition

To cover a general \acPWH, we partition it into monotone non-intercepting cells whose union is again the original \acPWH. In general any monotone decomposition would be feasible but we only consider the \acTCD and the \acfBCD choset1998coverage ; de1997computational . Both decompositions have different advantages as shown in Figure 5. The \acTCD provides a partitioning that adjusts well in rectangular scenes with multiple directions of extension. The \acBCD also adjusts well to rectangular scenes and usually leads to fewer cells and thus fewer redundant sweeps and traversal segments choset1998coverage . On the downside it can lead to degenerate sweeping behaviour in cells with narrow protrusions.

Since both decompositions result from scan line algorithms, the scan line direction determines the set of cells. To find a good decomposition direction we calculate the decomposition for every individual edge direction. A potentially good decomposition is the decomposition with the smallest altitude sum $w$, where altitude refers to the monotone extension of a cell. The minimum altitude sum corresponds to a minimum number of sweeps in the case of a fixed sweep direction huang2001optimal .

[TABLE]

where $m$ is the number of monotone cells, $y_{\mathrm{max,}i}$ is the $y$-coordinate of the uppermost vertex and $y_{\mathrm{min,}i}$ is the $y$-coordinate of the lowermost vertex in a $y$-monotone polygon cell.

4 Coverage Path Planning

After decomposing the input \acPWH into simple cells and generating a set of sweep patterns for each cell the planner has to find a shortest sequence of sweep patterns such that every cell, and thus the whole \acPWH, is covered. To solve this problem efficiently we formulate it as an \aclE-GTSP.

Figure 6 sketches the adjacency graph $G=(N,A)$, where $N=\{n_{1}\ldots n_{n}\}$ is the node set and $A=\{(n_{i},n_{j}):n_{i},n_{j}\in N,i\neq j\}$ is the set of directed arcs. The node set $N$ is divided into $m$ mutually exclusive and exhaustive clusters $N_{1}\ldots N_{m}$, i.e., $N=N_{1}\cup\ldots\cup N_{m}$ with $N_{i}\cap N_{j}=\emptyset~{}\forall i,j\in\{1\ldots m\},i\neq j$.

A node $n_{i}$ represents an individual sweep pattern that covers a single monotone polygon cell as shown in Figure 4 or the start or goal point. Every monotone cell (and the start and goal point) represents an individual cluster $N_{i}$. The arcs are the shortest path connecting the end of one sweep pattern with the start of the next sweep pattern in another cluster. The start node has outgoing arcs to all nodes, and the goal node has incoming arcs from all nodes.

Every arc $(n_{i},n_{j})$ has a non-negative cost $c_{ij}$ which is the sum of the shortest path cost $t_{ij}$ from $n_{i}$ to $n_{j}$ and the cost $t_{j}$ to execute sweep pattern $n_{j}$.

[TABLE]

Because the start and the end of a sweep pattern do not coincide, the cost matrix $\boldsymbol{C}=(c_{ij})$ is asymmetric, i.e., $c_{ij}\neq c_{ji}$.

Since \acpMAV can hover and turn on the spot, we can plan rest-to-rest segments between waypoints which obey the straight-line assumption of our coverage planner. The trajectories are modeled with velocity ramp profiles with instantaneous acceleration and deceleration and a maximum velocity. The cost of a trajectory is the sum of all segment times. The segment time $t$ between two waypoints is a function of the distance $d$, the maximum velocity $v_{max}$, and the maximum acceleration $a_{max}$,

[TABLE]

$t_{a}$ is the time to accelerate to maximum velocity and $d_{a}$ is the distance travelled while accelerating. Figure 7 shows that by tuning $v_{max}$ and $a_{max}$ with respect to our platform and sensor constraints the optimization finds a good compromise between minimizing distance and minimizing the number of turns.

Algorithm 1 summarizes the process of setting up the \acE-GTSP problem. First we decompose the polygon into monotone cells. For each cell we compute the sweep permutations and make each sweep pattern a node in the graph where the neighborhood is defined by the cell id. Finally, we create the edges between all sweeps patterns of different cells using the precomputed reduced visibility graph. Once the cost matrix is fully defined, we solve the \acE-GTSP using \acGK666https://csee.essex.ac.uk/staff/dkarap/?page=publications&key=Gutin2009a as an off-the-shelf open source solver gutin2010memetic .

4.1 Pruning

While our problem size (tens of clusters, hundreds of nodes, hundreds of thousands of edges) is no problem for the \acGK, generating the edges is the bottleneck of the algorithm because every sweep pattern needs to be connected to almost any other sweep pattern through a Euclidean shortest path. The total number of arcs grows quadratically with the number of nodes.

Fortunately, the optimization problem is modular, i.e., any sweep pattern combination that visits every cell will achieve full coverage and the path cost is the only optimization criterion. We can safely prune a node $n_{i}\in N_{k}$ if it is cheaper to first traverse from $n_{i}$’s start point to the start point of $n_{j}\in N_{k}$, perform coverage $n_{j}$ and return to $n_{i}$’s goal point.

5 Results

The algorithm has been implemented in C++ using \acCGAL cgal:eb-18b . \acCGAL provides efficient, reliable, and exact geometric algorithms which we extended to generate the sweep lines and the \acBCD. In order to find an exact solution we implement rice2012exact to convert the \acE-GTSP into a directed graph for an exhaustive exact solution using Dijkstra.

5.1 Simulation Benchmark

We setup a simulation benchmark to evaluate the performance of our algorithm (our) against the classical sweep planner with one sweep direction (one_dir) choset1998coverage , to compare \acBCD and \acTCD, and to demonstrate the superiority of a designated \acE-GTSP solver over exact brute-force search. Our test instances, e.g. Figure 2, are automatically generated from the EPFL aerial rooftop dataset sun2014free to provide polygon maps with realistic obstacles and dimensions. Our synthetic dataset consists of $320$ rectangular worlds with an area of $1\text{\,}\mathrm{ha}$ and [math] to $15$ rooftop obstacles. The benchmark was executed on an Intel®Core™i7-7820HQ CPU @ $2.90\text{\,}\mathrm{GHz}$.

To relate our solution to the complexity of the maps, we plot coverage path cost and algorithm runtime against the number of hole vertices, as these are the events of the decomposition scan line algorithms. Figure 8 (a) shows the trajectory cost of the different planner configurations. Because our planner takes advantage of different sweep directions and start points, it gives better results than the classic solution with only one fixed sweep pattern per cell. Furthermore, \acBCD leads to shorter trajectories than \acTCD because it generates fewer cells and thus fewer redundant sweeps at the adjacent edges and transition segments between the cells. On a side note, whenever the exact solution is available it coincides with the approximate solution from \acGK.

The relative differences in path cost between our planner and the other configurations is shown in Figure 8 (b). We observe up to $14\text{\,}\mathrm{\char 37\relax}$ improvement over the optimal fixed direction and $29\text{\,}\mathrm{\char 37\relax}$ improvement of \acBCD over \acTCD. Figure 8 (c) shows the computation time of the planners. Exact solutions fail to solve within $200\text{\,}\mathrm{s}$ for any of our scenarios with more than twenty hole vertices. Generating the product graph from the adjacency graph $G$ and Boolean lattice of possible cluster sequences grows exponentially in the number of clusters rice2012exact . The reduced problem one_dir_gk achieves the best computation time. At most it takes $12\text{\,}\mathrm{s}$ for the most complex map. The full \acE-GTSP with \acBCD solves it within $65\text{\,}\mathrm{s}$ which is still reasonable for employing the planner in field experiments. Table 1 reveals that generating the dense \acE-GTSP graph is the greatest computational burden.

5.2 Experiment

We validate our planner in a real flight on a DJI M600 Pro. Our drone covers a non-convex open area with sparse trees, depicted in Figure 1, in a low altitude mapping scenario. The flight corridor is selected in a georeferenced map, an optimal sweep path is calculated with our planner and executed under GNSS control. The drone follows a velocity ramp profile as described in section 4 and additionally turns at every waypoint in the flight direction. The slope of the terrain is not considered during planning. To regulate the \acAGL, we fuse Lidar and Radar altimeter data into a consistent altitude estimate. Table 2 shows the general experiment setup. The sweep distance is chosen based on image overlap. The velocity and acceleration are chosen to meet controller constraints and avoid motion blur. To validate coverage with a nadir configuration sensor we record QXGA top-down imagery and generate the \acDTM shown in Figure 9 using the Pix4D mapping tool. The \acDTM shows good coverage of the designated area but imperfect reconstruction due to instantaneous acceleration movements. Collision-free trajectory smoothing may improve the turning maneuvers both in speed and smoothness and consequentially increase overall performance of the coverage planner.

6 Conclusion

In this work we presented a boustrophedon coverage path planner based on an \acE-GTSP formulation. We showed in comprehensive benchmarks on realistic synthetic polygon maps that our planner reliably solves complex coverage tasks in reasonable computation time, making it suitable for field deployment. Furthermore, we showed that our planner outperforms the classic boustrophedon coverage planner in terms of path cost. We validated in a field experiment the usability of our coverage algorithm on a real \acMAV and show that we can cover a $1950\text{\,}{\mathrm{m}}^{2}$ area with obstacles at low altitude. Future work includes optimizing coverage for a side looking \acSAR configuration, collision-free trajectory smoothing to improve turning times, and integration into an airborne mine detection system.

Acknowledgment

This work is part of the FindMine project and was supported by the Urs Endress Foundation. The authors would like to thank their student Lucia Liu for her initial work on the \acBCD, Florian Braun and Michael Riner-Kuhn for their hardware support, and the reviewers for their constructive advice.