Triangulating a Polygon with Holes in Optimal (Deterministic) Time

We give the first optimal algorithm for triangulating a polygon with holes in $O(n+h\log h)$ time. The algorithm is deterministic.

Many applications follow [7, 19, 24, 29]. (For example, combining with Wang’s result [29] yields an optimal $O(n+h\log h)$-time algorithm for finding the Euclidean shortest path between two given points amidst $h$ disjoint polygonal obstacles with $n$ vertices total.)

Known approaches, as well as our new algorithm, work with the reformulation of the problem in terms of the trapezoidal decomposition of $h$ disjoint polygonal chains.

One could try to directly modify Chazelle’s linear-time algorithm for $h=1$ [5], but this appears difficult (considering how notoriously complicated his algorithm was). Bar-Yehuda and Chazelle’s approach was to directly reduce the general $h$ case to the $h=1$ case: join the multiple polygonal chains into one connected component by adding vertical line segments, and then apply Chazelle’s algorithm as a black box. To do so, the first observation is that polygonal chains intersecting a common vertical line can be joined together this way by Jordan sorting (the problem of sorting intersections of polygonal chains along a line), which is solvable by known techniques [18]. Roughly, their idea was then to build an interval tree [11, 26] for the $x$-spans of the chains, and solve the problem for each level of the interval tree (since each level consists of vertically separated instances in which chains intersect a common vertical line). The remaining task is merging the results of these $O(\log h)$ levels. Bar-Yehuda and Chazelle proposed a binary divide-and-conquer over the $O(\log h)$ levels, yielding $O((n+h\log h)\log \log h)$ time overall, and alternatively a divide-and-conquer with a larger fan-out, yielding $O(n+h{\log }^{1+\epsilon }h)$. This approach seems hard to improve further.

A different approach is randomized divide-and-conquer, as done by Clarkson et al. and Seidel [10, 27]. The trapezoidal decomposition of a random sample of $n/\log n$ edges allows us to divide the problem into subproblems of roughly logarithmic size. Assigning input elements to subproblems requires some form of “point location” and would normally cost $O(\log n)$ per element, but can be done faster knowing that the input elements comes from $h$ chains. Since the recursion has $O({\log }^{\ast }n)$ levels, the running time inherently increases by a ${\log }^{\ast }n$ factor, unfortunately. Even for $h=1$ (the case of a simple polygon), Clarkson et al. and Seidel’s methods do not achieve linear time (Amato, Goodrich, and Ramos [1] described a randomized linear-time algorithm which is simpler than Chazelle’s, but still relies on some of Chazelle’s ingredients and is still complicated). Furthermore, randomization seems essential: such random sampling methods where the sample size is large is generally harder to derandomize.

The difficulty of improving the random sampling approach turns out to be overstated. We observe that in combination with another known divide-and-conquer strategy via planar graph separators [22] or $r$-divisions [15, 16], the random sampling approach can indeed achieve optimal $O(n+h\log h)$ expected running time. The reason the ${\log }^{\ast }n$ factor is avoidable is that we can recurse for a constant number of levels (in fact, just once) and then switch to Chazelle’s algorithm (or more precisely Bar-Yehuda and Chazelle’s algorithm, which calls Chazelle’s) in the base case when the input size is sufficiently small – the $h\log h$ term is big enough to absorb the cost of the switch. The idea of combining sampling with $r$-divisions and invoking Chazelle’s in the base case is not original and can be viewed as a variant of an existing randomized algorithm by Eppstein, Goodrich, and Strash [13] for a related problem (see below). Still, it is interesting that this implication has not been noted earlier, and we include this optimal randomized solution in Appendix A.

Our main contribution is in obtaining an optimal deterministic algorithm. To do so, we combine ideas from both the random sampling approach and Bar-Yehuda and Chazelle’s deterministic approach, even though the two approaches appear different and incompatible. To derandomize the sampling approach, we need to compute a subset of the input edges that “behave” like a random sample, known as an $\epsilon$-net in the literature [17, 23] (with $\epsilon$ being $\frac{polylogn}{n}$ in our case). Constructing a net deterministically appears challenging, especially since we can only afford linear time. We observe a simple way to construct a net in linear time when the trapezoidal decomposition is given, by using planar graph separators or $r$-divisions (the net size is suboptimal, but is good enough for our purposes).

However, there is an obvious circularity issue: we are not given the trapezoidal decomposition, since that is what we want to compute! Our key idea is to go back to Bar-Yehuda and Chazelle’s interval-tree approach, which gives a partition of the input into $O(\log n)$ subsets where we already know how to compute the trapezoidal decomposition of each subset. Although merging their trapezoidal decompositions is not easy, merging their nets is: the union of the nets of the subsets is a net for the whole input (the net size is worsen by an $O(\log n)$ factor, but is still good enough for our purposes). By building the trapezoidal decomposition of the net in combination with an $r$-division, we can then divide the problem into subproblems of polylogarithmic size. The detailed description of the whole algorithm is given in Section 3, and is just about 2 pages long.

In the last part of the paper, we address an extension of the problem: computing the trapezoidal decomposition of polygonal chains that may not be simple, i.e., may self-intersect, and may intersect each other. The decomposition is defined in the same way, except that we also add vertical segments from each intersection point to the edge immediately above/below. The crux of the problem lies in the computation of the intersections. (Once all the intersection points are known, we could actually reduce to the non-intersecting case to obtain the final trapezoidal decomposition.)

The problem of computing all self-intersections of one chain (the $h=1$ case) is already challenging. Chazelle [5] posed the (still open) question of whether there is an $O(n+X)$-time algorithm, where $X$ denotes the number of intersections. Known line-segment intersection algorithms [9, 25, 6, 2] can solve the problem in $O(n\log n+X)$ time, which is optimal when $X\gg n\log n$ but not when $X$ is sublinear; they do not exploit the fact that the input line segments form a connected chain. Clarkson, Cole, Tarjan [8] showed that the random sampling approach yields an $O(n{\log }^{\ast }n+X)$ expected time bound for this problem. On the other hand, Eppstein, Goodrich, and Strash [13] used random sampling in combination with $r$-divisions to obtain an alternative $O(n+X{\log }^{(c)}n)$ expected time bound for an arbitrarily large constant $c$, where ${\log }^{(c)}$ denotes the logarithm function iterated $c$ times. Thus, their algorithm runs in optimal linear expected time when $X\ll n/{\log }^{(c)}n$. We are not aware of any deterministic result better than $O(n\log n+X)$, however.

We present a deterministic algorithm for computing the trapezoidal decomposition of a non-simple polygonal chain with $X$ intersections in $O(n+X{n}^{\epsilon })$ time. Thus, our algorithm runs in optimal linear time when $X\ll n^{1-\epsilon }$. (This case is already sufficient for various applications; for example, Eppstein and Goodrich [12] provided empirical evidence that road networks tend to have $X\approx \sqrt{n}$ intersections.)

More generally, when there are multiple chains ($h>1$), Clarkson, Cole, and Tarjan [8] obtained an $O(n{\log }^{\ast }n+h\log h+X)$ expected time bound. Our deterministic algorithm, presented in Section 4, takes $O(n+h\log h+X{n}^{\epsilon })$ time. (In Appendix A, we also note an extension of Eppstein, Goodrich, and Strash’s randomized algorithm with $O(n+h\log h+X{\log }^{(c)}n)$ expected time bound.)

For a set $S$ of $n$ line segments in the plane, the trapezoidal decomposition, denoted $\text{TD}(S)$, is the subdivision of the plane formed by $S$ together with vertical line segments drawn from each vertex $v$ of the arrangement (i.e., each endpoint and each intersection point) to the segment immediately above and below $v$. The faces in this subdivision are indeed trapezoids (or triangles) where the left and right sides (if exist) are vertical and the top and bottom sides (if exist) are parts of the input segments of $S$. If $S$ is disjoint, $\text{TD}(S)$ has $O(n)$ trapezoids and can be constructed in $O(n\log n)$ time [11, 26]. If $S$ has $X$ intersections, $\text{TD}(S)$ has $O(n+X)$ trapezoids and can be constructed in $O((n+X)\log n)$ time by the classic Bentley–Ottmann sweep [4, 26] or in $O(n\log n+X)$ time by the fastest known algorithms [9, 25, 6].

The trapezoidal decomposition of a collection $P$ of (not necessarily simple) polygonal chains, denoted $\text{TD}(P)$, is defined as the trapezoidal decomposition of its edge set.

We use the following weak definition of $r$-divisions [15]: For a planar graph $G$ with $n$ vertices, an $r$-division is a subset $B$ of $O(n/\sqrt{r})$ vertices such that every connected component in $G-B$ has size at most $r$. An $r$-division can be viewed as an extension of planar graph separators, and can be constructed in $O(n)$ time by known algorithms [16, 21] for any given $r$.

In our application, we consider an $r$-division of the dual²²2In the dual graph, the vertices are the trapezoids, and two trapezoids are adjacent if their boundaries (vertical sides or parts of input line segments) overlap. of a trapezoidal decomposition. Here, $B$ is a subset of trapezoids. Construct the subdivision of the plane formed by the sides of the trapezoids in $B$. The cells of this subdivision will be referred to as the regions of the $r$-division in this paper. Each region is a polygon possibly with holes, which is the union of at most $r$ trapezoids (a region could also just be a single trapezoid, for example, when it is from $B$). The total number of sides in all the regions³³3We use sides instead of edges to avoid confusion between geometric and graph terminologies. is $O(n/\sqrt{r})$. (With a stronger version of $r$-divisions [15, 21], it is possible to ensure that there are $O(n/r)$ regions, each having at most $O(\sqrt{r})$ sides and at most $O(1)$ holes, though this will not be necessary for our application.)

For a set $ℐ$ of $h$ intervals, the interval tree [11, 26] is a binary tree defined recursively as follows. If $ℐ=\mathrm{\varnothing }$, the tree is just a single node. Otherwise, let $a$, $b$, $m$ be the minimum, maximum, and median endpoint. At the root, store $[a,b]$, which we call the range, and $m$, which we call the pivot value. At the root, also store all intervals of $ℐ$ containing $m$. Recursively build the interval tree for the subset of all intervals of $ℐ$ completely left of $m$ (this is the left subtree) and the interval tree for the subset of all intervals of $ℐ$ completely right of $m$ (this is the right subtree). The resulting tree has $O(\log h)$ height, and can be constructed in $O(h\log h)$ time. Note that each interval of $ℐ$ is stored in exactly one node.

Let $s$ and $t$ be parameters to be set later. Our algorithm is described in a series of steps below. The purpose of Steps 0–2 is to compute a subset $R$ satisfying a certain property. Afterwards, Steps 3–5 use its trapezoidal decomposition $\text{TD}(R)$ to divide the problem into smaller subproblems and then solve these subproblems.

1. 0.

   For each chain $p\in P$, let $I_p$ be the interval determined by the smallest and largest $x$-coordinate in $p$; these can be generated in $O(n)$ time. Build the interval tree for these $O(h)$ intervals in $O(h\log h)$ time. For each node $v$ of the tree, let $P_v$ be the list of all chains $p\in P$ whose intervals $I_p$ are stored at $v$. Note that each chain appears in exactly one list.

2. 1.

   For each node $v$ of the interval tree with pivot value $m_v$:

   1. (a)

      Sort the intersections of $P_v$ with the vertical line $x=m_v$ (note that all chains in $P_v$ intersect this line). This is an instance of the Jordan sorting problem for multiple disjoint polygonal chains. Bar-Yehuda and Chazelle [3] showed that this can be done in $O(n_v+h_v\log h_v)$ time (by adapting a known Jordan sorting algorithm for a single polygonal chain [18]), where $n_v$ and $h_v$ are the number of vertices and chains in $P_v$ respectively. Note that ${\sum }_v{n}_v=n$ and ${\sum }_v{h}_v=h$.

   2. (b)

      Join all chains in $P_v$ into a single simple, connected component $P_v^{\prime }$ by adding vertical line segments between all pairs of consecutive intersection points along $x=m_v$, by using the output of step 1(a). (Not all such vertical segments are needed to connect the chains, as shown in Figure 2 (left), but this would not matter.) Apply Chazelle’s algorithm [5] to compute $\text{TD}(P_v^{\prime })$, and thus $\text{TD}(P_v)$, in $O(n_v)$ time.

   3. (c)

      Compute an $s$-division of the dual of $\text{TD}(P_v)$ in $O(n_v)$ time. Preprocess the regions of the $s$-division for point location [26] in linear time. Let $R_v$ be the subset of all $O(n_v/\sqrt{s})$ edges of $P_v$ that participate in the region boundaries in this $s$-division (this excludes vertical sides). Then $R_v$ obeys the following property:⁴⁴4This property basically says that $R_v$ is an $O(s/n_v)$-net [17, 23] of the edge set of $P_v$, in the set system corresponding to edges intersecting vertical line segments. The size of our net ($O(n_v/\sqrt{s})$) isn’t the smallest possible (which would be $O(n_v/s)$) but is good enough for our purposes.

      ($\ast$) Given any vertical line segment $e$ not intersecting $R_v$, the segment $e$ can intersect at most $O(s)$ edges of $P_v$; furthermore, these edges can be found in $O(\log n+s)$ time.

      This property is easy to see: any such segment $e$ must lie completely inside a region in the $s$-division, which can be found by point location [26] in $O(\log n)$ time, and there are at most $O(s)$ edges of $P_v$ participating in the region.

3. 2.

   Let $R$ be the union of all the $R_v$’s, together with $O(h)$ extra line segments of zero length at the start vertices of the chains of $P$. Note that $|R|=O(n/\sqrt{s}+h)$, and $R$ satisfies the following property:⁵⁵5This property basically says that $R$ is an $O((s\log h)/n)$-net of the edge set of $P$ (ignoring the zero-length segments).

   ($†$) Given any vertical line segment $e$ not intersecting $R$, the segment $e$ can intersect at most $O(s\log h)$ edges of $P$; furthermore, these edges can be found in $O((\log n+s)\log h)$ time.

   This follows immediately from property ($\ast$), since the $x$-coordinate of a vertical segment $e$ can lie in the ranges of only $O(\log h)$ nodes $v$ of the interval tree.

4. 3.

   Compute $\text{TD}(R)$ in $O(|R|\log |R|)$ time, and compute a $t$-division of the dual of $\text{TD}(R)$ in $O(|R|)$ time.

5. 4.

   For each region $\gamma$ in the $t$-division and each vertical side $e$ of the boundary $\partial \gamma$, compute all intersections of $e$ with $P$. Since the interior of $e$ does not intersect $R$, we know that $e$ has at most $O(s\log h)$ intersections with $P$ by property ($†$), and they can be found in $O((\log n+s)\log h)$ time. (See Figure 2 (right).) By summing over all $O(|R|/\sqrt{t})$ choices of $e$, there are at most $O((|R|/\sqrt{t})\cdot s\log h)$ such intersections in total and they can be found in $O((|R|/\sqrt{t})\cdot (\log n+s)\log h)$ time.

   Cut the chains of $P$ at all these intersection points to obtain a new collection $P^{\prime }$ of polygonal chains. In order to cut, we may need to sort the intersections along each edge of $P$, but this costs at most a logarithmic factor per intersection, i.e., $O((|R|/\sqrt{t})\cdot s\log h\cdot \log n)$ in total. The number of vertices in $P^{\prime }$ is $n^{\prime }\le n+O((|R|/\sqrt{t})\cdot s\log h)$, and the number of chains in $P^{\prime }$ is $h^{\prime }\le h+O((|R|/\sqrt{t})\cdot s\log h)$.

6. 5.

   For each region $\gamma$ in the $t$-division, let $P_{\gamma }^{\prime }$ be the chains of $P^{\prime }$ inside $\gamma$, and compute $\text{TD}(P_{\gamma }^{\prime })$ by applying Bar-Yehuda and Chazelle’s previous algorithms [3] in $O(n_{\gamma }^{\prime }+h_{\gamma }^{\prime }{\log }^{O(1)}{h}_{\gamma }^{\prime })$ time, where $n_{\gamma }^{\prime }$ and $h_{\gamma }^{\prime }$ are the number of vertices and chains in $P_{\gamma }^{\prime }$ respectively.⁶⁶6We could avoid Bar-Yehuda and Chazelle’s algorithm and use some slower algorithm running in $O(n_{\gamma }^{\prime }+{(h_{\gamma }^{\prime })}^{O(1)})$ time, if we bootstrap one more time. Note that ${\sum }_{\gamma }{n}_{\gamma }^{\prime }=n^{\prime }$ and ${\sum }_{\gamma }{h}_{\gamma }^{\prime }=h^{\prime }$. Furthermore, ${\max }_{\gamma }{h}_{\gamma }^{\prime }\le O(t\cdot s\log h)$ (since there are at most $O(t)$ sides in $\gamma$, and at most $O(t)$ start vertices of chains inside $\gamma$ due to those extra segments of zero length added to $R$).

   Afterwards, we can piece together $\text{TD}(P_{\gamma }^{\prime })$ (after clipping to $\gamma$) over all regions $\gamma$, to obtain $\text{TD}(P)$ in $O(n^{\prime })$ additional time.

The total running time for Steps 0–1 is $O(n+h\log h)$. The running time for Step 3 is $O(|R|\log |R|)=O((n/\sqrt{s}+h)\log n)$ time, as $|R|=O(n/\sqrt{s}+h)$. The total running time for Steps 4–5 is at most $O(n+((|R|/\sqrt{t})\cdot s\log h\cdot \log n)+(h+(|R|/\sqrt{t})\cdot s\log h){\log }^{O(1)}(t\cdot s\log h))$.

Setting $s={\log }^{2}n$ and $t={\log }^{8}n$ gives an overall time bound of $O(n+h\log n+h{(\log \log n)}^{O(1)})=O(n+h\log n)$, which is the same as $O(n+h\log h)$ (since the second term dominates only if $h\ge n/\log n$, but then $\log n=O(\log h)$). $◀$

As a corollary, we can test whether a collection $P$ of $h$ polygonal chains has any intersection in $O(n+h\log h)$ deterministic time: simply run the above algorithm to attempt to compute $\text{TD}(P)$, and run Chazelle’s linear-time algorithm [5] to check simplicity of the resulting (connected) diagram.

Let $\mathrm{\ell }$ be the maximum number of vertices per chain in $P$. Let $s$ and $t$ be parameters to be set later. We describe an algorithm below to compute all intersections of $P$. Steps 1–2 compute a subset $R$ satisfying a certain property. Afterwards, Steps 3–5 use its trapezoidal decomposition $\text{TD}(R)$ to divide the problem into smaller subproblems and then solve these subproblems.

1. 1.

   For each $i\in \{1,\mathrm{\dots },b\}$:

   1. (a)

      Apply Theorem 1 to compute $\text{TD}(P_i)$ in $O(n_i+h_i\log h_i)$ time, where $n_i$ and $h_i$ are the number of vertices and chains in $P_i$ respectively.

   2. (b)

      Compute an $s$-division of $\text{TD}(P_i)$ in $O(n_i)$ time. Preprocess the regions of the $s$-division for point location [26] in linear time. Let $R_i$ be the subset of all $O(n_i/\sqrt{s})$ edges of $P_i$ that participate in the boundaries of regions in this $s$-division.

2. 2.

   Let $R$ be the union of all the $R_i$’s, together with $O(h)$ extra line segments of zero length at the start vertices of the chains of $P$. Note that $|R|=O(n/\sqrt{s}+h)$, and $R$ satisfies the following property, similar to before:

   ($†$) Given any vertical line segment $e$ not intersecting $R$, the segment $e$ can intersect at most $O(bs)$ edges of $P$; furthermore, these edges can be found in $O(b(\log n+s))$ time.

3. 3.

   Compute $\text{TD}(R)$ in $O(|R|\log |R|+X)$ time, and compute a $t$-division of $\text{TD}(R)$ in $O(|R|+X)$ time.

4. 4.

   For each region $\gamma$ in the $t$-division and each vertical side $e$ of the boundary $\partial \gamma$, compute all intersections of $e$ with $P$. Since the interior of $e$ does not intersect $R$, there are at most $O(bs)$ such intersections and they can be found in $O(b(\log n+s))$ time by property ($†$). By summing over all $O((|R|+X)/\sqrt{t})$ choices of $e$, there are at most $O(((|R|+X)/\sqrt{t})\cdot bs)$ such intersections in total and they can be found in $O(((|R|+X)/\sqrt{t})\cdot b(\log n+s))$ time.

   Cut the chains of $P$ at all these intersection points to obtain a new collection $P^{\prime }$ of polygonal chains. In order to cut, we may need to sort the intersections along each edge of $P$, but this costs at most a logarithmic factor per intersection, i.e., $O(((|R|+X)/\sqrt{t})\cdot bs\cdot \log n)$ in total. The number of vertices in $P^{\prime }$ is $n^{\prime }\le n+O(((|R|+X)/\sqrt{t})\cdot bs)$, and the number of chains in $P^{\prime }$ is $h^{\prime }\le h+O(((|R|+X)/\sqrt{t})\cdot bs)$.

5. 5.

   For each region $\gamma$ in the $t$-division, let $P_{\gamma }^{\prime }$ be the set of all chains of $P^{\prime }$ with start vertices in $\gamma$, let $Q_{\gamma }$ be the set of all non-vertical sides of $\partial \gamma$, and apply Theorem 1 to test whether $P_{\gamma }^{\prime }\cup Q_{\gamma }$ has any intersection in $O(n_{\gamma }^{\prime }+(h_{\gamma }^{\prime }+|Q_{\gamma }|)\log n)$ time, where $n_{\gamma }^{\prime }$ and $h_{\gamma }^{\prime }$ are the number of vertices and chains in $P_{\gamma }^{\prime }$ respectively. Note that ${\sum }_{\gamma }{n}_{\gamma }^{\prime }=n^{\prime }$ and ${\sum }_{\gamma }{h}_{\gamma }^{\prime }=h^{\prime }$ and ${\sum }_{\gamma }|Q_{\gamma }|=O((|R|+X)/\sqrt{t})$. If the answer is no, mark the region $\gamma$ good and mark all the chains in $P_{\gamma }^{\prime }$ good; otherwise, mark the region $\gamma$ bad and mark all the chains in $P_{\gamma }^{\prime }$ bad. If $\gamma$ is good, all chains in $P_{\gamma }^{\prime }$ do not intersect each other and stay inside $\gamma$ (since chains in $P^{\prime }$ do not intersect the vertical sides of $\partial \gamma$). Thus, good chains do not intersect each other. On the other hand, the number of bad regions is at most $O(X)$, since we can charge each bad region to a unique intersection (note that each chain is in just one $P_{\gamma }^{\prime }$); so, the number of bad chains is at most $O(bstX)$, since each $P_{\gamma }^{\prime }$ has at most $O(bs\cdot t)$ chains.

6. 6.

   For each bad chain $p^{\prime }$ in $P^{\prime }$, find all regions $\gamma$ in the $t$-division that are intersected by $p^{\prime }$. To do so, begin with the region $\gamma$ containing the start vertex of $p^{\prime }$. Find the next line segment $e\in Q_{\gamma }$ intersected by $p^{\prime }$, if exists, by brute force in $O(\mathrm{\ell }t)$ time, since $p^{\prime }$ has size at most $O(\mathrm{\ell })$ and $Q_{\gamma }$ has size at most $O(t)$. Then reset $\gamma$ to the neighboring region defined by $e$, and repeat. The total time is at most $O(\mathrm{\ell }t\cdot bstX)$, since every region $\gamma$ visited by $p^{\prime }$ after the first can be charged to a unique intersection, and there are at most $O(bstX)$ bad chains $p^{\prime }$.

7. 7.

   For each bad chain $p^{\prime }$ in $P^{\prime }$ and each good region $\gamma$ intersected by $p^{\prime }$ (all such regions have been identified from step 6), compute all intersections of $p^{\prime }$ with the good chains of $P_{\gamma }^{\prime }$ by brute force in $(\mathrm{\ell }\cdot \mathrm{\ell }bst)$ time (since $P_{\gamma }^{\prime }$ has at most $O(\mathrm{\ell }\cdot bst)$ edges). The total time is at most $O(\mathrm{\ell }\cdot \mathrm{\ell }bst\cdot bstX)$.

8. 8.

   Lastly, compute all intersections among the bad chains, which have $n^{\prime \prime }=O(\mathrm{\ell }\cdot bstX)$ edges, e.g., by the Bentley–Ottmann sweep [4] in $O((n^{\prime \prime }+X)\log n^{\prime \prime })$ time.

The total running time for Step 1 is $O(n+h\log h)$. The running time for Step 3 is $O(|R|\log |R|+X)=O((n/\sqrt{s}+h)\log n+X)$ time. The total running time for Steps 4–5 is at most $O(n+((|R|+X)/\sqrt{t})\cdot bs\cdot \log n+h\log n)$. The total running time for Steps 6–8 is $O({(b\mathrm{\ell }st)}^{O(1)}X)$.

Setting $s={\log }^{2}n$ and $t=b^2{\log }^{6}n$ gives an overall time bound of $O(n+h\log h+{(b\mathrm{\ell })}^{O(1)}X{\log }^{O(1)}n)$.

By initially subdividing the input chains into subchains of length $\log n$, we can lower $\mathrm{\ell }$ to $\log n$ while increasing $h$ to $h+O(n/\log n)$; the $O(h\log h)$ term would then be increased to $O(n+h\log h)$.

After computing all intersections, we can turn $P$ into a collection of at most $h$ disjoint simple chains with $O(n+X)$ vertices, by making tiny cuts to eliminate cycles and taking an Euler tour of each connected component (to traverse the components, we may need to sort intersection points along each edge, but this costs at most $O(X\log n)$ time). Afterwards, we can apply Theorem 1 to compute $\text{TD}(P)$ in $O(n+X+h\log h)$ additional time. $◀$

Finally, we combine the above lemma with an $n^{\mathrm{\Theta }(\epsilon )}$-way divide-and-conquer to solve the general problem:

Divide $P$ into $b$ subcollections $P_1,\mathrm{\dots },P_b$, each with $n/b$ vertices, possibly by splitting at most $b$ of the input chains. Suppose $P_i$ has $h_i$ chains and $X_i$ intersections. Note that ${\sum }_i{h}_i\le h+b$ and ${\sum }_i{X}_i\le X$. Recursively compute $\text{TD}(P_i)$, whose output gives us the $X_i$ intersections of each $P_i$. Subsequently, $P_i$ can be turned into a non-intersecting collection of at most $h_i$ simple chains with $O(n/b+X_i)$ vertices. Now, apply Lemma 2 to compute $\text{TD}(P)$ in $O(n+X+h\log h+b^{O(1)}X{\log }^{O(1)}n)$ time.

The total running time satisfies the recurrence

which solves to $T(n,h,X)=O((n+h\log h+b^{O(1)}X{\log }^{O(1)}n)\cdot {\log }_{b}n)$. Setting $b=n^{\epsilon /c}$ for a sufficiently large constant $c$ yields the result. $◀$

The main open question is whether the running time of the above theorem could be improved to $O(n+h\log h+X)$. A more modest goal is to find a deterministic algorithm with running time $O(n+h\log h+X{\log }^{O(1)}n)$.