A Linear Time Algorithm for the Maximum Overlap of Two Convex Polygons Under Translation

The computational geometry literature is replete with various types of problems about polygons [36], e.g., polygon containment [14] (which has seen exciting development even in the convex polygon case, with translations and rotations, as recently as last year’s SoCG [13]). The maximum overlap problem may be viewed as an outgrowth of polygon containment problems.

The maximum overlap problem may also be viewed as a formulation of shape matching [7], which has numerous applications and has long been a popular topic in computational geometry. Many notions of distance between shapes have been used in the past (e.g., Hausdorff distance, Fréchet distance, etc.), and the overlap area is another very natural measure (the larger it is, the closer the two shapes are) and has been considered in numerous papers for different classes of objects (e.g., two convex polytopes in higher dimensions [3, 5, 4], two arbitrary polygons in the plane [18, 25], two unions of balls [10], one convex polygon vs. a discrete point set¹¹1In the discrete case, we are maximizing the number of points inside the translated polygon. [9, 1], etc.) as well as different types of motions allowed (translations, rotations, and scaling – with translations-only being the most often studied). Problem 1 is perhaps the simplest and most basic version, but is already quite fascinating from the theoretical or technical perspective, as we will see.

Due to the concavity of (the square root of) the objective function [21], it is not difficult to obtain an $O((n+m)polylog(n+m))$-time algorithm for the problem by applying the well-known parametric search technique [30, 2] (or more precisely, a multidimensional version of parametric search [38]). De Berg, Cheong, Deviller, van Kreveld, and Teillaud [21] took more care in eliminating extra logarithmic factors to obtain their $O((n+m)\log (n+m))$-time algorithm, by avoiding parametric search and instead using more elementary binary search and matrix search techniques [24]. But any form of binary search would seem to generate at least one logarithmic factor, and so it is tempting to believe that their time bound might be the best possible.

On the other hand, in the geometric optimization literature, there is another well-known technique that can give rise to linear-time algorithms for various problems: namely, prune-and-search (pioneered by Megiddo and Dyer [31, 22, 32]). For prune-and-search to be applicable, one needs a way to prune a fraction of the input elements at each iteration, but our problem is sufficiently complex that it is unclear how one could throw away any vertex from the input polygons and preserve the answer.

Perhaps because of these reasons, no improved algorithms have been reported since de Berg et al.’s 1998 work. This is not because of lack of interest in finding faster algorithms. For example, a (lesser known) paper by Kim, Choi, and Ahn [28] considered an unbalanced case of Problem 1 and gave an algorithm with $O(n+m^2{\log }^{3}n)$ running time, which is linear when $m\ll \sqrt{n/{\log }^{3}n}$. Ahn et al. [6] investigated approximation algorithms for the problem and obtain logarithmic running time for approximation factor arbitrarily close to 1, while Har-Peled and Roy [25] more generally gave a linear-time algorithm to approximate the entire objective function, which was then used to obtain a linear-time approximation algorithm for an extension of the problem to polygons that are decomposable into $O(1)$ convex pieces.

We obtain a randomized (exact) algorithm for Problem 1, running in linear (i.e., $O(n+m)$) expected time! This is the first improvement to de Berg et al.’s result [21] in 26 years, and is clearly optimal, and so fully settles the complexity of Problem 1.

Besides our result itself, the way we obtain linear time is also interesting. As mentioned, if one does not care about extra $\log (n+m)$ factors, one can use (multidimensional) parametric search to solve the problem. Specifically, parametric search allows us to reduce the optimization problem to the implementation of a line oracle (determining which side of a given line the optimal point $t^{\ast }$ lies in), which in turn can be reduced to a point oracle (evaluating the area of $(P+t)\cap Q$ for a given point $t$). A point oracle can be implemented in linear time by computing the intersection $(P+t)\cap Q$ explicitly. (In de Berg et al.’s work [21], they directly implemented the line oracle in linear time.)

Various linear-time prune-and-search algorithms (e.g., for low-dimensional linear programming [31, 22, 32], ham-sandwich cuts [33], centerpoints [26], Euclidean 1-center [23], rectilinear 1-median [39, 35], etc.) also exploit line or point oracles; using (in modern parlance) cuttings [16], such oracle calls let them refine the search space and eliminate a fraction of the input at each iteration. Unfortunately, for our problem, we can’t technically eliminate any part of the input polygons $P$ and $Q$. Instead, our new idea is to reduce the cost of the oracles iteratively as the search space is refined, which effectively is as good as reducing the input complexity itself. As the algorithm progresses, oracles get cheaper and cheaper, and we will bound the total running time by a geometric series.

To execute this strategy, one should think of oracles as data structures with sublinear query time (i.e., sublinear in the original input size $n+m$). Our key idea to designing such data structures is to divide the input polygons into blocks of a certain size based on angles/slopes. This idea is inspired by the recent work of Chan and Hair [13] (they studied the convex polygon containment problem with both translations and rotations, and although our problem without rotations might seem different, their divide-and-conquer algorithms also crucially relied on division of the input convex polygons by angles/slopes). We show that at each iteration, as we refine the search space, we can use the current “block data structure” to obtain a better “block data structure” with a larger block size. (The refinement process from one block structure to another block structure shares some general similarities with the idea of fractional cascading [17].)

We are not aware of any prior algorithms in the computational geometry literature that work in the same way (which makes our techniques interesting). A linear-time algorithm by Chan [12] for a matrix searching problem (reporting all elements in a Monge matrix less than given value) has a recurrence similar to the one we obtain here, but this appears to be a coincidence. An algorithm for computing so-called “$\epsilon$-approximations” in range spaces by Matoušek [29] also achieved linear time (for constant $\epsilon$) by iteratively increasing block sizes, but the algorithm still operates in the realm of traditional prune-and-search (reducing actual input size). Chazelle’s infamous linear-time algorithm for triangulating simple polygons [15] also iteratively refines a “granularity” parameter analogous to our block size, and as one of many steps, also requires the implementation of oracles (for ray shooting, in his case) whose cost similarly varies as a function of the granularity (at least from a superficial reading of his paper) – fortunately, our algorithm here will not be as complicated!

Recently, Jung, Kang, and Ahn [27] have announced a linear-time algorithm for a related problem: given two convex polygons $P$ and $Q$ in the plane, find a vector $t\in {ℝ}^2$ that minimizes the area of the convex hull $(P+t)\cup Q$. Although the problem may look similar to our problem on the surface, it is in many ways different: in the minimum convex hull problem, the optimal point $t^{\ast }$ is known to lie on one of $O(n+m)$ candidate lines (formed by an edge of one polygon and a corresponding extreme vertex of the other polygon), whereas in the maximum overlap problem (Problem 1), $t^{\ast }$ lies on one of $O({(n+m)}^2)$ candidate lines (defined by pairs of edges of the two polygons). Thus a more traditional prune-and-search approach can be used to solve the minimum convex hull problem, but not the maximum overlap problem. Indeed, Jung et al.’s paper [27] considered the maximum overlap problem, but they do not give an improvement over de Berg et al.’s $O((n+m)\log (n+m))$ time algorithm for the two-dimensional case. Instead, they give some new results for the problem in higher dimensions $d$, but their time bound is $O(n^{2\lfloor d/2\rfloor })$, which is much larger than linear.

We observe that, in the specific case of maximum overlap for convex polytopes in dimension $d=3$, one can obtain a significantly faster $O((n+m)polylog(n+m))$-time algorithm by directly using multi-dimensional parametric search [38], since a point oracle (computing the area of the intersection of two convex polytopes in ${ℝ}^3$) can be done in $O(n+m)$ time by explicitly computing the intersection [15, 11], and this is parallelizable with $O((n+m)polylog(n+m))$ work and $O(polylog(n+m))$ steps by known parallel algorithms for intersection of halfspaces in ${ℝ}^3$, i.e., convex hull in ${ℝ}^3$ in the dual (e.g., see [8]). This observation seems to have not appeared in literature before (including [40, 27]).

We want to find a translation maximizing

De Berg, Cheong, Devillers, van Kreveld, and Teillaud showed several interesting properties of the area as a function of $t$ [21]. For our purposes, we only need to import the following.

The restriction on $t$ is important, as $A(t)$ suddenly transitions from a convex function to a constant function $A(t)=0$ when $t$ is outside of the specified region. For our algorithm, however, this technicality will have little impact.

First, we give a definition of the type of cuttings we use.

Endowed with the ability to sample random hyperplanes, we can actually produce such a cutting with very high probability in sublinear time, as per Clarkson and Shor [19, 20] (roughly, we can take logarithmically many random hyperplanes, and return a canonical triangulation of their arrangement).

Let $v_1,\mathrm{\dots },v_n$ be the vertices of $P$ listed in counterclockwise order, and let $v_1^{\prime },\mathrm{\dots },v_m^{\prime }$ be the vertices of $Q$ listed in counterclockwise order. Let $P[v_x:v_y]$ denote the convex hull of the set of points containing: (i) the origin, (ii) the vertices $v_x$ and $v_y$, and (iii) all vertices encountered when traversing the boundary of $P$ in counterclockwise order from $v_x$ to $v_y$. We define $Q[v_x^{\prime }:v_y^{\prime }]$ analogously.

By a standard application of prefix sums, we can produce a data structure to report $\text{Area}(P[v_x:v_y])$ and $\text{Area}(Q[v_x^{\prime }:v_y^{\prime }])$ with $O(N)$ preprocessing time and $O(1)$ query time. This is since (taking $\text{Area}(P[v_x:v_y])$ as an example):

Throughout, we refer to the above data structure as the area prefix sums for $P$ and $Q$.

Consider two convex polygons $A$ and $B$, or more generally, two polygonal chains $A$ and $B$. Let $e$ be an edge of $A$ and let $e^{\prime }$ be an edge of $B$. These two edges define a parallelogram $\pi (e,e^{\prime })$ in configuration space such that, for any vector $t\in {ℝ}^2$, we have that $e+t$ intersects $e^{\prime }$ iff $t\in \pi (e,e^{\prime })$. Let $\partial \pi (e,e^{\prime })$ denote the set of four line segments defining $\pi (e,e^{\prime })$. We define the set of all such line segments as

We define $\overleftrightarrow{\mathrm{\Pi }}(A,B)$ to be $\mathrm{\Pi }(A,B)$ where each line segment is extended to a line.

The angle of an edge $e$ on the boundary of an arbitrary convex polygon $A$, denoted ${\theta }_A(e)$, is the angle (measured counterclockwise) starting from a horizontal rightward ray and ending with a ray coincident with $e$ that has $A$ on its left side. We always have that ${\theta }_A(e)\in [0,2\pi )$. For any subset $X$ of the edges of $A$, let ${\mathrm{\Lambda }}_A(X)$ denote the minimal interval in $[0,2\pi )$ such that ${\theta }_A(e)\in {\mathrm{\Lambda }}_A(X)$ for all $e\in X$. We call ${\mathrm{\Lambda }}_A(X)$ the angle range for the subset $X$ with respect to $A$.

By definition of the block structure, we have a partition of $S=[N/b]$ into two subsets $S_{\text{Good}}$ and $S_{\text{Bad}}$ such that $|S_{\text{Bad}}|\le \mu$, and for each $i\in S_{\text{Good}}$ we have a constant complexity quadratic function $f_i(t)=\text{Area}((P_i+t)\cap Q_i)$ for all $t\in 𝒯$.

We can thus compute

in time $O(N/b)$ by simply evaluating each $f_i(t^{\ast })$. Since $|S_{\text{Bad}}|\le \mu$, we can compute

by evaluating the area for each $(P_i+t^{\ast })\cap Q_i$ directly. The time required is $O(\mu b)$ using known algorithms for intersecting two convex polygons whose edges are in sorted order. $◀$

We now turn our attention to the summation ${\sum }_{i\ne j}\text{Area}((P_i+t^{\ast })\cap Q_j)$. Because every term amounts to computing the intersection for two blocks with different angle ranges, we could use binary searches to evaluate each term in $O(polylog(b))$ time. However, this is too slow for our purposes, as there are a near quadratic number of terms. Instead, we take advantage of additional structure to decompose the summation into $O(N/b)$ intersections of unions of blocks, such that each intersection can still be computed via binary search.

As a tool, we need a way to efficiently intersect convex chains that meet certain conditions.

Because the angle ranges for $X$ and $Y$ are disjoint, they act as pseudo-disks, and hence can intersect at most twice. We can find the intersection point(s), if they exist, using a standard nested binary search in $O(\log |X|\log |Y|)$ time (or a more clever binary search [37] in $O(\log |X|+\log |Y|)$ time, though we don’t need this improvement). $◀$

Now we can show how to compute ${\sum }_{i\ne j}\text{Area}((P_i+t^{\ast })\cap Q_j)$.

We give a recursive algorithm for the problem. As a subroutine, we consider the problem of computing, given two integers $x,y\in [N/b]$ with $x>y$,

Let $T(y-x,b)$ be the time complexity of this subroutine. Since the original summation ${\sum }_{i\ne j}\text{Area}((P_i+t^{\ast })\cap Q_j)$ is equivalent to ${\sum }_{(i,j)\in W(1,\lceil N/b\rceil )}\text{Area}((P_i+t)\cap Q_j)$, the lemma amounts to showing that $T(\lceil N/b\rceil -1,b)=O(N\cdot \frac{{\log }^{2}b}{b})$, or more generally, $T(z,b)=O(z{\log }^{2}b)$.

Consider an instance of this subroutine with any parameters $y>x$. Let $m:=\lceil \frac{x+y}{2}\rceil$ and $z:=y-x$. We can decompose the target summation as follows:

The first two summations can be evaluated in time at most $2T(\lceil \frac{y-x}{2}\rceil ,b)$ via recursion (unless $x-y=1$, in which case they are automatically zero). The last two summations can actually be evaluated directly. We consider the case of the very last summation; the remaining one can be evaluated via a symmetric process. First note that

Let and $B={\bigcup }_{m\le j\le y}{Q}_j$. Observe that $A$ and $B$ are both convex polygons, since each is the union of consecutive blocks. By construction, the edges that $A$ shares with $P$ and the edges that $B$ shares with $Q$ form two convex polygonal chains that have disjoint angle ranges, i.e., ${\mathrm{\Lambda }}_P(E_P(A))\cap {\mathrm{\Lambda }}_Q(E_Q(B))=\mathrm{\varnothing }$. These chains account for all but two edges of $A$ and two edges of $B$. Thus we can decompose $A$ and $B$ into $O(1)$ convex polygonal chains satisfying the preconditions of Lemma 10. This implies that $A$ and $B$ intersect in at most a constant number of locations, and furthermore we can find the intersection points in time $O({\log }^2(bz))$.

By decomposing $A$ and $B$ into $O(1)$ pieces based on the intersection locations, we can use the area prefix sums⁸⁸8The area prefix sums for $P$ and $Q$ also act as ares prefix sums for $A$ and $B$., along with a direct computation on the $O(1)$ intersecting pairs of edges, to find the area of intersection in $O(1)$ additional time. So the recurrence is $T(z,b)=2T(\lceil z/2\rceil ,b)+O({\log }^2(bz))$ when $z>1$ and $T(1,b)=O({\log }^{2}b)$. This solves to $T(z,b)=O(z{\log }^{2}b)$. $◀$

This follows from known observations by de Berg et al. [21]. $◀$

Before describing how to perform cascading, we need a tool to help classify the new blocks based on their relationship with ${𝒯}^{\prime }$.

Let the set of vertices for ${𝒯}^{\prime }$ be $V({𝒯}^{\prime })$. For each $t\in V({𝒯}^{\prime })$, we know that, by Lemma 10, $X+t$ and $Y$ can intersect in at most two points, and we can find the two pairs of edges that contain these intersection(s) in time $O(\log |X|\log |Y|)$.¹⁰¹⁰10There is an edge case where at least one of the intersections occurs at a vertex of $X$ or a vertex of $Y$, but this can be handled with slightly more work without increasing the asymptotic running time. Let the set of all $O(1)$ edge pairs $(e_X,e_Y)$ for all vertices in $V({𝒯}^{\prime })$ be $L$. We can determine if ${𝒯}^{\prime }\subseteq \pi (e_X,e_Y)$ for all $(e_X,e_Y)\in L$ via brute force in time $O(1)$. If this is not the case, then some line segment of $\mathrm{\Pi }(X,Y)$ intersects the interior of ${𝒯}^{\prime }$, and we are done. If this is the case, then because ${𝒯}^{\prime }$ is convex and each $\pi (e_X,e_Y)$ is convex, we have that each translation $t\in {𝒯}^{\prime }$ realizes a set of intersecting pairs that is a (possibly strict) superset of $L$.

It remains to determine if any translations $t\in {𝒯}^{\prime }$ realize an intersecting pair not in $L$, because this case (and only this case) would mean that a line segment in $\mathrm{\Pi }(X,Y)$ intersects the interior of ${𝒯}^{\prime }$. Let $X^{-}$ be $X$ with its edges participating in $L$ deleted, and let $Y^{-}$ be $Y$ with its edges participating in $L$ deleted. In time $O(\log |X|\log |Y|)$ we can use techniques similar to the ones given above to determine if the interior of the (implicitly generated) Minkowski sum of ${𝒯}^{\prime }$ and $X^{-}$ has any intersection with $Y$, and then if the interior of the (implicitly generated) Minkowski sum of ${𝒯}^{\prime }$ and $X$ has any intersection with $Y^{-}$. $◀$

We are now ready to prove Lemma 6. The main observation is that the new quadratic functions stored in this data structure only need to apply when $t\in {𝒯}^{\prime }$, not the more general case that $t\in 𝒯$, so many of the previously difficult blocks can now be handled. We will find a quadratic function for almost all indices $k\in [N/b^{\prime }]$ in two steps:

1. 1.

   Scan through the information provided in the input $(\mu ,b,𝒯)$-block structure, and use this to recover part of the function $\text{Area}((P_k^{\prime }+t)\cap Q_k^{\prime })$.

2. 2.

   Recover the missing information for $\text{Area}((P_k^{\prime }+t)\cap Q_k^{\prime })$ by intersecting $O(\frac{b^{\prime }}{b})$ pairs of convex polygons with disjoint angle ranges. We evaluate these using nested binary searches (Lemma 10 or 13) and area prefix sums.

For all $k\in [N/b^{\prime }]$, note that because $b$ divides $b^{\prime }$, we have a set of indices $C_k=\{(k-1)\cdot \frac{b^{\prime }}{b}+1,\mathrm{\dots },k\cdot \frac{b^{\prime }}{b}\}$ such that $P_k^{\prime }={\bigcup }_{i\in C_k}{P}_i$ and $Q_k^{\prime }={\bigcup }_{j\in C_k}{Q}_j$.¹¹¹¹11The set may be smaller for the very last value of $k$, but this does not impact the proof. If for all $(i,j)\in C_k\times C_k$ and for all $t\in {𝒯}^{\prime }$ we can summarize $\text{Area}((P_i+t)\cap Q_j)$ with a constant-complexity quadratic function, then we can simply add these to get a new function $f_k(t)=\text{Area}((P_k^{\prime }+t)\cap Q_k^{\prime })$.

We first argue that each quadratic function $f_k(t)$ can be found efficiently, if it exists. Consider just a single value of $k$. We use the decomposition

For the first summation on the right hand side, we scan through the list $S_{\text{Good}}$ from the original $(\mu ,b,𝒯)$-block structure. Each function $\text{Area}((P_i+t)\cap Q_i)$ with $i\in S_{\text{Good}}$ already has a quadratic function, and reading these takes time $O(\frac{b^{\prime }}{b})$ for each $k$, or $O(N/b)$ time for all $k$. Globally, there are at most $\mu$ indices $i$ not in $S_{\text{Good}}$. For each of these, we can check whether ${𝒯}^{\prime }$ lies in a single cell of $\mathrm{\Pi }(P_i,Q_i)$ via directly examining each line segment in $\mathrm{\Pi }(P_i,Q_i)$ in time $O(b^2)$. If this is the case, then by Lemma 12, $\text{Area}((P_i+t)\cap Q_i)$ can be summarized by a quadratic function for all $t\in {𝒯}^{\prime }$, and we find the function via directly examining each parallelogram in $\mathrm{\Pi }(P_i,Q_i)$ in time $O(b^2)$; the total additional time is $O(\mu b^2)$. (If this fails for some $i\in C_k$, we will add the index $k$ to $S_{\text{Bad}}^{\prime }$.)

Now consider the second summation on the right hand side. We can rewrite this as

where $Q_{i,k}^{-}$ is the union of all $Q_j$’s in the original sum with , and $Q_{i,k}^{+}$ is the union of all $Q_j$’s in the original sum with $j>i$. We evaluate each of the terms individually.

Consider the case of $\text{Area}((P_i+t)\cap Q_{i,k}^{-})$, as the case of $\text{Area}((P_i+t)\cap Q_{i,k}^{+})$ is symmetric. Now, $Q_{i,k}^{-}$ is a convex polygon, because it is the union of consecutive blocks. By construction, ${\mathrm{\Lambda }}_P(E_P(P_i))\cap {\mathrm{\Lambda }}_Q(E_Q(Q_{i,k}^{-}))=\mathrm{\varnothing }$. Therefore, we can break $P_i+t$ and $Q_{i,k}^{-}$ into $O(1)$ pieces that satisfy the preconditions of Lemma 10, showing that $P_i+t$ and $Q_{i,k}^{-}$ can only intersect at a constant number of locations for any translation $t\in {ℝ}^2$. These same pieces allow us to apply Lemma 13 to determine if ${𝒯}^{\prime }$ lies in a single cell of $\mathrm{\Pi }(P_i,Q_{i,k}^{-})$, which takes time $O(\log |P_i|\log |Q_{i,k}^{-}|)=O(\log b\log b^{\prime })$. If this is the case, then by Lemma 12 there exists a quadratic function equal to $\text{Area}((P_i+t)\cap Q_{i,k}^{-})$ for all $t\in {𝒯}^{\prime }$. The nonconstant components of this quadratic function only depend on the $O(1)$ intersecting pairs of edges, and we have already found these via Lemma 13. Thus we can break $P_i$ and $Q_{i,k}^{-}$ into $O(1)$ pieces based on the intersection points, and use the area prefix sums from the $(\mu ,b,𝒯)$-block structure (which also serve as area prefix sums for $P_i$ and $Q_{i,k}^{-}$), to find the quadratic function equal to $\text{Area}((P_i+t)\cap Q_{i,k}^{-})$ in $O(1)$ additional time. (If this fails for some $i\in C_k$, we will add the index $k$ to $S_{\text{Bad}}^{\prime }$.)