Line Segment Visibility in Simple Polygons: Exact, Robust, Scalable Computation and Applications

In the following, $P$ denotes a simple polygon and $n$ its number of vertices. Two points $x$ and $y$ are visible to each other if the line segment $xy\subset P$. The visibility polygon of a point $s\in P$, denoted by $V_P(s)$, is the set of all points in $P$ visible to $s$. Given a segment $\mathrm{\ell }=pq$ inside of $P$, the weak visibility polygon $V_P(\mathrm{\ell })$ is the set of all points in $P$ that are visible from some point on $\mathrm{\ell }$, i.e., $V_P(\mathrm{\ell })≔\{p\in P\mid \exists s\in \mathrm{\ell }:p\in V_P(s)\}$; see Figure 1.

For two points $p,q\in P$, the link distance $d_L(p,q)$ between them is the minimum number of edges in any polygonal path within $P$ between $x$ and $y$ (see Figure 1, Left). Consider a point $s\in P$ and its visibility polygon $V_P(s)$. The vertices of $V_P(s)$ are either vertices of $P$ or the shadows of reflex vertices of $P$ seen by $s$. The edges of $V_P(s)$ between reflex vertices and their shadows are called windows. A window $w$ is a chord of $P$ and divides $P$ into two parts; define $P(w,s)$ to be part of $P$ which does not contain $s$. The window partition ${𝒲}_P(s)$ is the partition of $P$ returned by Algorithm 1 described by Suri [50] given input $s$. We drop the subscript $P$ from ${𝒲}_P(s)$ when it is clear from context which polygon we are referring to.

Let $R(w)$ denote the region generated by a window $w$. The rooted dual tree $T_s$ of $𝒲(s)$ has as vertices the windows defining $𝒲(s)$ and edges between two vertices if their regions share a common boundary. The root of $T_s$ is $s$. We refer to $T_s$ as the window tree of $s$. See Figure 1, Right. Let $g(u,v)$ denote the minimum distance of two nodes $u,v$ in $T_s$ and consider a point $p\in R(w)$. The link distance $d_L(s,p)$ is $g(s,w)+1$ [50].

Constructing a triangulation $𝒯$ of a polygon $P$ takes linear time [11]. Thus, computing $V_P(\mathrm{\ell })$ for a segment $\mathrm{\ell }$ in $P$ takes linear time – first triangulate $P$, then use the algorithm from [24]. As we are unaware of any implementation of [11], we use a constrained Delaunay triangulation of $P$ provided by CGAL [5]. This has a worst-case complexity of $𝒪(n^2)$ [8] but works extremely well in practice, avoids degenerate triangles and thus simplifies our workflow. Our implementation runs in linear time given $𝒯$ as in [24]. To achieve this runtime, we implemented a specialized data structure called a finger tree [30, 29] (see Section 2.1). To the best of our knowledge, all previous implementations of finger trees are in functional programming languages like Haskell or Scala [29, 22]; we used C++ for their iterative versions. We implemented an improvement to [24] suggested in [50]; this allows us to compute $V_P(\mathrm{\ell })$ in time linear in the number of triangles in $𝒯$ it intersects rather than $n$. Apart from allowing us to reliably compute $V_P(\mathrm{\ell })$ even in extremely large input polygons, this is crucial to compute window trees in $𝒪(n)$. All our implementations use exact number types and handle degeneracies. In Section 2, we review the algorithms introduced in [24, 50] and discuss relevant implementation details. We recommend readers familiar with these results to skim through our implementation notes in Section 2 before skipping to Section 3.

Computing visibility polygons is central to computational geometry. We focus on weak visibility polygons of segments inside simple polygons; for a more general overview of visibility problems, see [56, 46, 48, 21, 15]. Given a triangulation of a polygon, Guibas, Hershberger, Leven, Sharir and Tarjan [24] presented an $𝒪(n)$ algorithm that computes the visibility polygon. A number of algorithms have been proposed for data structures that can efficiently compute these for a query segment $\mathrm{\ell }$, each with its own trade-offs in terms of preprocessing time, query time, and space complexity. We summarize the most relevant ones in Table 1.

Our results are given in bold. We remark that implementing the more recent results with improved output sensitive query time [3, 6, 9, 13] with exact predicates seems extremely challenging. It is also unclear if these will lead to improvements in runtime in practice (see Section 4 for experiments).

Computing segment visibility is a crucial subroutine in several important geometric data structures; perhaps the most well-known one being window partitions and window trees [50]. Other than computing minimum link paths, window partition is also used in convex decomposition algorithms [7, 38] as well as in pursuit evasion problems [52].

Minimum link paths have numerous applications across domains: in motion planning, where turns are expensive and their number should be minimized [47]; in communication systems design under line-of-sight constraints, where a min-link path minimizes the number of necessary repeaters; in placement of telescoping manipulators with flexible telescopic links [34]; in curve simplification by approximating curves with minimum link paths [25, 31, 43]; and in VLSI layouts, where wires typically run along orthogonal axes, motivating the study of link path problems in rectilinear settings [54, 14]. For a comprehensive overview, we refer the reader to [39, Chapter 12] and [40]. There also exist several works on implementing minimum link paths: [53] is restricted to rectilinear paths and is the only work with a publicly available implementation we could find, [20, 57] also consider rectilinear paths, and [17, 58] are restricted to rectilinear environments. Baum, Bläsius, Gemsa, Rutter, and Wegner use minimum link paths in [4] to simplify the boundary of a polygon in a corridor; this is the only other experimental paper on (non-rectilinear) minimum link paths we are aware of. While the authors provided parts of their non-public implementation to us and the paper provides valuable insights into optimizing minimum link path computations for a specific source-target pair, their implementation is tailored to their specific use case and does not offer a method to compute optimal paths between two arbitrary points in a simple polygon.

Let $P$ be a simple polygon with $n$ vertices and let $s$ be a given source vertex of $P$. For two points $x,y\in P$, let $\pi (x,y)$ denote the Euclidean shortest path between $x$ and $y$ inside $P$. Lee and Preparata [37] showed that for all vertices $v$ of $P$, $\pi (s,v)$ is a polygonal path whose corners are vertices of $P$ and that ${\bigcup }_v\pi (s,v)$ over all vertices $v$ forms a tree rooted at $s$. We use $Q_P(s)$ to denote this shortest-path tree rooted at $s$ – the vertices of $Q_P(s)$ are vertices of $P$ and $(u,v)$ is an edge of the tree if it is an edge of either $\pi (s,u)$ or $\pi (s,v)$.

Let $𝒯$ be a triangulation of $P$ and $G$ the dual graph of $𝒯$. Let $d=uw$ be a diagonal or and edge of $P$ and let $a$ be the deepest common ancestor of $u$ and $w$ in $Q_P(s)$, see Figure 2 Left. Then the path $\pi (a,u)$ and $\pi (a,w)$ are two outward convex chains [37], i.e., the convex hull of each subpath lies outside of the region bounded by $\pi (a,u)$, $\pi (a,w)$, and $uw$. Following [37], we call $F=F_{uw}=\pi (a,u)\cup uw\cup \pi (a,w)$ the funnel of $d$ with cusp $a$.

We compute $Q_P(s)$ in linear time by maintaining funnels of diagonals in $𝒯$ as we process its vertices. We maintain the current funnel $F_{uw}$ as a sorted list $[u_l,u_{l-1},\mathrm{\dots },u_1,a,w_1,\mathrm{\dots },w_k]$ where $a=u_0=w_0$ is the cusp of $F$, $\pi (a,u)=[u_0,u_1,\mathrm{\dots },u_l]$ and $\pi (a,w)=[w_0,w_1,\mathrm{\dots },w_k]$ with $u_l=u$ and $w_k=w$. We call Algorithm 2 [24] with $s$ and an adjacent vertex $v_1$ as the initial funnel $F$; for an interior source point $s$, we start by locating the triangle containing $s$ and create three funnels, one for each edge, with cusp $s$.

We explain Lines 4 and 11 of Algorithm 2 in more detail. As $\pi (a,u)$ and $\pi (a,w)$ are outward convex, the slopes of the segments defining them are monotonic. Thus, we can do a binary search on these slopes to find the point of tangency for $x$ to $F$ in Line 4; see Figure 2, Left. We do not store shortest paths explicitly as suggested by Line 11, we simply store back pointers from $x$ to $v$; $\pi (s,x)$ is computed by following these pointers until we reach $s$.

Although using binary search as above is straightforward with a binary search tree, it requires $\log n$ time, as $\pi (s,u)$ and $\pi (s,v)$ (and thus our list representing $F$) might contain $𝒪(n)$ elements. For an overall linear runtime, we need to find $v$ in (amortized) constant time. This is achieved by storing the list $F$ as a so-called finger tree [26]. Our finger trees are modified 2-3 search trees which include pointers (called fingers) to the last and first element of $F$. They allow for access to $v$, as well as the splitting of $F$ at $v$ in $𝒪(\log \delta )$ time, where $\delta$ is the distance from $v$ to the nearest finger. See Guibas, McReight, Plass and Roberts [26] for details on properties of finger trees. Gibiansky [22] also provides valuable insights.

Given a polygon $P$ and a segment $\mathrm{\ell }=ab$ inside of $P$, recall that the weak visibility polygon $V_P(\mathrm{\ell })$ contains all points in $P$ that are visible from a point on $\mathrm{\ell }$, i.e., $V_P(\mathrm{\ell })=\{s\in P\mid \exists t\in \mathrm{\ell }:s\text{ is visible from }t\}$. The following structural lemma is crucial.

Consider such an edge $e^{\prime }=cd$. The structure formed by $\mathrm{\ell }$, $e^{\prime }$, $\pi (a,c)$ and $\pi (b,d)$ is called an hourglass, see Figure 2, Right. When running the shortest path tree algorithm with $a$ as the source, the funnel $F_{e^{\prime }}$ has $x$ as its cusp and $y$ as an adjacent vertex. Similarly, when running the algorithm with $b$ as the source, the funnel $F_{e^{\prime }}$ has $w$ as its cusp and $z$ as an adjacent vertex. These two tangents define the visible section of $e^{\prime }$.

During the computation of the $Q_P(s)$, we check, for each vertex $v$, in constant time if $\pi (s,v)$ is outward convex or not. To extract the visibility polygon, we iterate over the vertices of the polygon and check whether an edge is (partially) visible as explained above. Then we connect continuous visible portions of the boundary, skipping over invisible edges, by connecting visible vertices with their shadow vertices (defined by the tangents). While the algorithm defined in [24] is not output sensitive (taking $𝒪(n)$ time regardless of the size of $V_P(\mathrm{\ell })$), we explain in Section 2.3 modifications to this algorithm to have this characteristic.

A naive implementation of Algorithm 1 using [24] would, for each window, compute a complete shortest path tree of $P$, resulting in a worst case runtime of $\mathrm{\Omega }(n^2)$. Suri [50] suggested modifications to [24] to make the computation of segment visibility output sensitive. The main idea is to synchronize the shortest path tree computation of both endpoints of a segment $\mathrm{\ell }$ and abort the computation of the visibility polygon if an edge is no longer visible from $\mathrm{\ell }$.

Let $P$ be a simple polygon and $𝒯$ be a triangulation of the interior of $P$. For a given segment $\mathrm{\ell }=ab$, we start by computing the shortest path tree rooted in $a$ and $b$ along the triangles intersected by $\mathrm{\ell }$. When processing the next point $x$ in a triangle $\mathrm{△}uwx$, we only continue with a diagonal $ux$ if either $\pi (a,u)$ and $\pi (b,x)$, or $\pi (a,x)$ and $\pi (b,u)$ are both outward convex. Symmetrically, we only continue with diagonal $wx$ if either $\pi (a,w)$ and $\pi (b,x)$, or $\pi (a,x)$ and $\pi (b,w)$ are both outward convex. The resulting runtime to compute a visibility polygon $V_P(\mathrm{\ell })$, is $𝒪(k_{\mathrm{\ell }})$, where $k_{\mathrm{\ell }}$ is the number of triangles intersected by $V_P(\mathrm{\ell })$. To ensure that the implementation achieves this runtime, we carefully adjust the traversal of the polygon boundary to only consider vertices that are part of the considered triangles.

Note that $k_{\mathrm{\ell }}$ can be arbitrarily large compared to the size of the visibility polygon $V_P(\mathrm{\ell })$, e.g., for a triangular visibility polygon that intersects many triangles of $𝒯$. However, exploiting the fact that each triangle is intersected by at most three partitions [50], we obtain a linear-time algorithm to compute the window partition $𝒲(s)$ of a point $s$ in $P$. Note that to start the shortest-path tree computation, we have to locate a triangle intersecting $\mathrm{\ell }$. We do this using CGAL’s built-in locate function, which in the worst case has a linear runtime.

To compute the window partition $𝒲(s)$, we combine Algorithm 1 with our output sensitive visibility polygon computation. The algorithm initially computes a single triangulation $𝒯$ of $P$, which we use to compute the non-regularized visibility polygon of $s$, giving us the initial set of windows. Then, for each window $w$, we can compute the visibility polygon $R(w):=V_{P(w,s)}(w)$ on the remaining polygon $P(w,s)$. We reuse $𝒯$ by checking if the next vertex is in $P(w,s)$. This window partition forms the basis for our implementation of the minimum link path algorithms in Section 3.

To construct exact minimum link paths from a source point $s$ to any target point $t$ in a simple polygon $P$, we use the window partition $𝒲(s)$ rooted at $s$. First, we locate the region $R(w)$ that contains $t$ and then traverse the unique path $(w,w_1\mathrm{\dots }s)$ from $w$ to $s$ in $T_s$. Our construction of $𝒲(s)$ informs us of a minimum link path from $s$ to $t$ in $P$ that has a vertex on each of the windows $w,w_1,\mathrm{\dots },s$ in that order. Constructing $T_s$ takes $𝒪(n)$ time and space [50] and locating the target node takes $𝒪(n)$. Point location could be reduced to $𝒪(\log n)$ time with preprocessing, but the preprocessing takes super-linear time [15, Chapter 6]. In Section 3.3, we explain how we reconstruct the minimum link path in $P$ in linear time.

To compute an additive ($\text{OPT}+4$) approximative minimum link path from $p$ to $q$ in a simple polygon $P$, we can use the window partition $𝒲(s)$ of an arbitrary point $s\in P$. For two vertices $u$ and $v$ of $T_s$, let $g(u,v)$ denote the distance in $T_s$ between them.

The idea is to first identify the two regions $R(w_p)$ and $R(w_q)$ that contain $p$ and $q$ in $𝒲(x)$ respectively. Then, we identify the lowest common ancestor $w_0$ of $w_p$ and $w_q$ in $T_s$. This allows us to construct two paths: one from $w_p$ to $w_0$ and one from $w_q$ to $w_0$. Finally, we have to connect the two paths with an additional link on $w_0$.

Note that, even if $w_0=w_p$ or $w_0=w_q$, i.e., the path is directed in $T_s$, we may still require an additional link using the above algorithm: let $w_0=w_p$ and $w_0,w_1,\mathrm{\dots },w_q$ be a directed path in $T_s$. As in Section 3.1, we construct a path from $q$ to $p$ with one vertex $v_i$ on each window $w_i$ for $0\le i\le t$. Unlike the exact algorithm, we cannot guarantee that $p$ is visible from $v_0$, we can only guarantee the existence of a point $v_0^{\prime }$ on $w_0$ that is visible from $p$. Therefore, we may need to add an additional link from $v_0$ to $v_0^{\prime }$ and then to $p$.

We precompute window partitions rooted at a set of seed vertices and answer queries by probing a small subset of them. At query time, a small set of candidate seed vertices is chosen, e.g., the $k$ nearest to the endpoints. Additionally, we post-process the candidates by removing any introduced loops. The algorithm returns a path with a minimum number of links among all candidates and runs in time proportional to the number of candidates and the cost of locating nodes and reconstructing paths, i.e., $𝒪(k\cdot n)$ in the worst case for $k$ candidates. Note that Lemma 2 also provides a lower bound on the link distance between $p$ and $q$ that can be used to evaluate the quality of the approximated path.

Now we reconstruct the path from $s$ to $t$ for the exact algorithm detailed in the prior sections. Specifically, we look to construct $(s,p_0,p_1,\mathrm{\dots },p_j,t)$, where $p_i$ belongs to the window $w_i$. This is done backwards from $t$. For the step from $w_0$ to $s$, any point on $w_0$ works because $s$ sees all points on $w_0$. For going from a window $w_{i+1}$ to $w_i$, this can easily be done in $O(1)$ by extending $w_{i+1}$ to a line and determining its intersection with $w_i$. However, especially if $t$ is not fixed when constructing the window partition, e.g., because we are answering multiple queries with the same source $s$, we need to find a point on the last window $w_j$ that sees $t$. There is a straightforward algorithm that can do this in linear time: compute $V_{R(w_j)}(t)$ and intersect this with $w_j$ to obtain the maximal connected subsegment of $w_j$ that sees $t$. In practice, this approach is not preferable due to the lack of efficient implementations of visibility algorithms for non-regularized polygons that often appear in the generated window partitions, see Figure 3. The linear-time algorithm available in CGAL [28] is not directly applicable, as it requires the polygon to be simple. Additionally, the triangular expansion approach implemented in CGAL, which is often preferred in practical applications, requires a triangulation of each window, which cannot be directly derived from our triangulation of $P$, as windows are not necessarily edges of the triangulation; even if we had such a triangulation, the interface of CGAL’s visibility package does not allow specifying a triangulation as input. Even without these issues, triangulating each region separately would result in a (practical) superlinear worst-case runtime.

To avoid these problems, we compute the maximal connected subsegment $w_j^{\prime }$ of $w_j$ that is visible from $t$ in linear (in the complexity of $R(w_j)$) time as follows. Consider the ray between $t$ and the left endpoint of $w_j$; this possibly intersects $\partial R(w_j)$. Consecutive hitting points define chains of edges on one side of the ray that occlude parts of $w$. March through the vertices in these chains and choose the vertex defining the rightmost ray from $t$ through it. We find the supporting ray for the right endpoint symmetrically. If these supporting rays cross, $w_j$ is not visible from $t$ (an impossibility). Otherwise, they define the maximal subsegment on $w_j$ visible from $t$; see Figure 4 for an illustration. Because each other step of path reconstruction takes time $O(1)$, this allows path reconstruction for source $s$ in time $O(|R(w_j)|+j)$ once $𝒲(s)$ is computed and the window containing $t$ is identified.

For the approximation algorithm, we apply the above procedure for both the windows containing $s$ and $t$; apart from this and an added segment on the window defining their deepest common ancestor, the situation is analogous.

For evaluating our implementation, we consider the following research questions.

RQ1

How does our implementation scale to large inputs, and how is the total runtime distributed onto the individual subroutines?

RQ2

Is there a measurable difference between using finger trees and CGAL’s multiset? Which of the two approaches is faster?

RQ3

How does the performance of our implementation depend on the size of the output?

To address these questions, we randomly generate segments for simple polygons from the Salzburg Database of Geometric Inputs²²2We used all valid simple polygons from https://zenodo.org/records/3784789 not labeled as contrived. [19], the largest of which has $\mathrm{18\hspace{0.17em}556\hspace{0.17em}888}$ vertices. These polygons are of various types and include polygons with many collinear points or other degeneracies. To generate segments, we sample a source point uniformly at random within the given polygon, compute the visibility polygon of that point, and then either sample a target point in the visibility polygon or pick a random vertex of it. The former approach corresponds to sampling segments inside polygons uniformly at random, whereas the latter focuses on placing one of the endpoints on the boundary of the polygon, exercising edge cases. Using each polygon and each approach twice results in a total of $\mathrm{24\hspace{0.17em}096}$ instances. For each such instance, and each of the two shortest path tree data structures, we compute the exact visibility polygon five times, measuring the runtime (both total and of various subroutines). Figure 6 depicts the resulting runtimes.

We can see that the overall runtime is dominated by computing a triangulation of the input polygon, which we do via CGAL’s constrained Delaunay triangulation (CDT). Even CGAL’s mark_in_domain routine, which simply flags the CDT’s triangles that constitute the polygon, is typically more expensive than the actual visibility computation. The fact that the runtime for the visibility computation seems to be almost unaffected by the size of the input polygon can be explained by the fact that it – except for a few small operations such as finding the segment endpoints in the triangulation – mostly depends on the size of the output polygon; see Figure 7, Left. On our inputs, the implementation using finger trees is slightly faster than the one using CGAL’s Multiset, see Figure 7 Right. In a preliminary microbenchmark based on building funnels by left insertion, our finger trees became significantly faster already at around $15$ funnel vertices; this makes a difference even for the typically small funnels we encounter in our algorithm. We therefore used our finger trees in our other experiments.

Recall that computing a window partition is a basic routine that is useful in different contexts. Furthermore, we only need one window partition from a source point $s$ to compute shortest link-distance paths to any target $t$, and only a constant number to approximate shortest link-distance paths for any source-target pair inside our polygon. We therefore study the performance of our window partition on its own, as well as in the context of minimum link paths, addressing the following research questions.

RQ4

How does our window partition implementation scale to large polygons?

RQ5

How does it compare to the other parts of our exact minimum link path implementation?

RQ6

How does the approximate minimum link path algorithm perform w.r.t. runtime and solution quality?

To answer these questions, we randomly selected a sample of $1000$ simple polygons with at most $\mathrm{2\hspace{0.17em}000\hspace{0.17em}000}$ vertices from [19]. For each such polygon $P$, we generate $10$ source points $s\in P$ uniformly at random and computed a window partition for each source. Furthermore, for each source point $s$, we generate a random target point $t\in P$ outside of the visibility polygon of $s$ and computed an exact minimum link $s$-$t$-path, as well as an approximate minimum link path using $k=7$ randomly chosen source points for the window partitions in our approximation. The resulting runtimes are depicted in Figure 8.

For each path computed using the approximation algorithm, Figure 9 shows the absolute, additive error depending on the number of window partitions $k$; the quality of the approximation notably improves with higher $k$. One should, however, note that increasing $k$ makes both preprocessing and path reconstruction more expensive, both scaling linearly in $k$.

Overall, we find that both our window partition and minimum link path implementation scale well to large polygons; in fact, in our experiments, time spent inputting, verifying and outputting polygons was typically orders of magnitude higher than pure computation time. If we want to query minimum link paths from a single source $s$, each subsequent query helps amortizing the cost of computing a window partition; similarly, for doing many queries for different sources and targets in the same polygon, the approximate algorithm can be used to amortize the cost of computing a constant number of window partitions. Depending on the relative importance of runtime versus solution quality, the approximation may be preferable.