Computing Largest Subsets of Points Whose Convex Hulls Have Bounded Area and Diameter

We describe the algorithm of Eppstein, Overmars, Rote and Woeginger [5].

Take any point $p\in P$ and choose it to be the bottommost vertex of a set of candidate solutions. In other words, we search for the optimal solution among all convex polygons with $k$ points inside that have $p$ as the bottommost vertex. By iterating over all points $p\in P$, we will find the global optimum with a linear factor extra in the time bound.

Given the bottommost point $p$ in the solution, we can use a recursive structure when setting up the dynamic program by considering a bottom-vertex triangulation. This is a triangulation where all vertices on the boundary connect to $p$. A convex polygon $C$ is now partitioned into a fan of triangles around $p$ and above it.

Assume the vertices are listed in clockwise order around $C$ as $p,p_1,\mathrm{\dots },p_m$. Then the last (most clockwise) triangle in the fan around $p$ is $\mathrm{△}p{p}_{m-1}{p}_m$.

If $C$ is an optimal solution with $k$ points from $P$, then the convex polygon $p,p_1,\mathrm{\dots },p_{m-1}$ is an optimal solution with $k-K(p,p_{m-1},p_m)-1$ points among the solutions that have $p_{m-1}$ as the counterclockwise neighbor of $p$ and where $p_{m-2},p_{m-1},p_m$ makes a right turn (is convex at $p_{m-1}$ in $C$). This suggests a 3-dimensional table $T[q,r,k]$, defined as: the minimum area of a convex polygon with $p$ as the bottommost vertex, $q$ as its ccw (counterclockwise) neighbor, and $r$ as the ccw neighbor of $q$. See Figure 2. This leads to a standard dynamic program with the recurrence:

where $s$ is chosen over the points above the horizontal line through $p$, left of the directed line from $p$ to $r$, and that leads to convexity at $r$. Since $K(p,q,r)$ counts only points strictly inside $\mathrm{△}pqr$, we subtract one more for point $q$. As the base case we let $T[q,r,k]=A(p,q,r)$ if triangle $\mathrm{△}pqr$ has exactly $k$ points, $p$, $q$, and $r$ included.

It is easy to see that filling an entry takes $O(n)$ time, and the whole dynamic program requires $O(n^{3}k)$ time. Together with the fact that we try all $n$ points as the bottommost point $p$, the whole algorithm takes $O(n^{4}k)$ time.

In [5] a technique is described that reduces the running time by a linear factor. A similar table is used, assuming a given point $p$ as bottommost point, but with a subtly different definition: $T[q,r,k]$ is the smallest area convex polygon with $k$ points inside where $p$ is the bottommost vertex, $q$ is its ccw neighbor, and the convex polygon lies on the same side of the line through $q$ and $r$ as $p$ does. Note that $r$ might not be part of the convex polygon itself. With a rotational sweep around $q$, entries are filled in only $O(1)$ time per point $r$. We consider all points $r$ above a horizontal line through $p$, not just the ones left of the directed line through $p$ and $q$. We consider the full lines through $r$ and $q$ and order them by angle around $q$, starting clockwise from the line through $p$ and $q$.

The rotational sweep maintains a value $T_q^{\ast }$ that is the minimum area encountered so far with $k$ points. The sweep has two cases: $(i)$ Point $r$ lies left of the directed line through $p$ and $q$ (in that order), in which case we may use $T[q,r,k-K(p,q,r)-1]$ to potentially find a smaller $T_q^{\ast }$, in which case we update $T_q^{\ast }$ and fill $T[q,r,k]$; $(ii)$ Point $r$ lies to the right of the directed line, in which case we fill $T[q,r,k]$ with $T_q^{\ast }$. This leads to an $O(n^{3}k)$ time algorithm. Details can be found in [5].

In this section we show how to compute a minimum area convex polygon with bounded diameter and $k$ points inside. We do this by trying all pairs of points $p,q,\in P$, assuming they form the diameter, and then computing an optimal convex polygon under this assumption. The algorithm is an extension of the basic, $O(n^{4}k)$ time algorithm from [5]. At the end of this section we show how to use this algorithm to find a convex polygon with bounded diameter and area and the maximum number of points inside. The full version [8] contains omitted proofs.

After choosing a pair $p,q$, we first select all points of $P$ that are at most $|pq|$ from both $p$ and from $q$. They are the points in a lune formed by the intersection of two disks with radius $|pq|$, centered at $p$ and at $q$, see Figure 3. For ease of description and without loss of generality, we let $pq$ be vertical with $p$ below $q$. We now need to form the right chain of a convex polygon, clockwise (cw) from $q$ to $p$, and the left chain, clockwise from $p$ to $q$. While doing so, we must make sure that we do not consider solutions that use a point $\mathrm{\ell }$ on the left chain that is too far from a point $r$ on the right chain. Note that two points on the right chain can never be too far from each other, and the same is true for the left chain (since we restricted our points to the ones inside the lune).

To solve the problem for a given $p,q$ we combine the dynamic programming solution from [5] with a rotating calipers algorithm to make sure that all antipodal pairs (one from the left chain and one from the right chain) have distance no more than $|pq|$.

Let $R[r,r^{\prime },\mathrm{\ell },k]$ be the minimum area of a convex polygon with $pq$ as a diameter, where $r$ is the ccw neighbor of $p$, $r^{\prime }$ is the ccw neighbor of $r$, $\mathrm{\ell }$ is the ccw neighbor of $q$, that has exactly $k$ points inside, and where $\mathrm{\ell }$ is antipodal to $r$ and to $r^{\prime }$. Analogously, let $L[\mathrm{\ell },{\mathrm{\ell }}^{\prime },r,k]$ be the minimum area of a convex polygon with $pq$ as a diameter, where $\mathrm{\ell }$ is the ccw neighbor of $q$, ${\mathrm{\ell }}^{\prime }$ is the ccw neighbor of $\mathrm{\ell }$, $r$ is the ccw neighbor of $p$, that has exactly $k$ points inside, and where $r$ is antipodal to $\mathrm{\ell }$ and to ${\mathrm{\ell }}^{\prime }$. We allow $r^{\prime }$ to be $q$ and ${\mathrm{\ell }}^{\prime }$ to be $p$. Note that $p$ and $q$ are antipodal: there is a pair of horizontal lines through $p$ and through $q$ that have any convex polygon with $p,q$ as diameter in between. We study antipodality more closely in the following lemma.

From the lemma above, it follows that $\mathrm{\ell }$ is not antipodal to $r^{\prime }$ or $r$ is not antipodal to ${\mathrm{\ell }}^{\prime }$. Now consider an optimal solution $C_{\text{opt}}$ with a given diameter $pq$, and let $\mathrm{\ell }$, $r$, ${\mathrm{\ell }}^{\prime }$ and $r^{\prime }$ be as in the proof. If $\mathrm{\ell }$ is not antipodal to $r^{\prime }$, then we can see $C_{\text{opt}}$ as being composed of a triangle $\mathrm{△}\mathrm{\ell }{\mathrm{\ell }}^{\prime }q$ and an optimal subsolution $C_{\text{opt}}^{\prime }$ where ${\mathrm{\ell }}^{\prime }$ is the ccw neighbor of $q$. Since $\mathrm{\ell }$ is antipodal to no more than $r$ and possibly $p$, by induction we can assume that $C_{\text{opt}}^{\prime }$ has all antipodal pairs no larger than the diameter $|pq|$, we need to test only the distance $|\mathrm{\ell }r|$ to conclude that $C_{\text{opt}}$ also has all antipodal pairs close enough. If $r$ is not antipodal to ${\mathrm{\ell }}^{\prime }$, the argument is symmetric.

If $|r\mathrm{\ell }|>|pq|$ then $R[r,.,\mathrm{\ell },.]=\mathrm{\infty }$. Similarly, if $|r^{\prime }\mathrm{\ell }|>|pq|$ then $R[.,r^{\prime },\mathrm{\ell },.]=\mathrm{\infty }$.

The recurrence optimizes over two options: use a subproblem of the $R$-type or use a subproblem of the $L$-type. The optimization of the $R$-type using a step on the right side is:

Here the $r^{\prime \prime }$ must be such that $r^{\prime \prime }{r}^{\prime }r$ is a right turn at $r^{\prime }$, vertex $r^{\prime \prime }$ is right of the line through $r^{\prime }$ and $q$, and $|r^{\prime \prime }\mathrm{\ell }|\le |pq|$. In such a case, $r^{\prime }$ is antipodal to only $\mathrm{\ell }$. We can choose $r^{\prime \prime }$ to be $q$. If $r^{\prime }=q$, then this option cannot be chosen because the right boundary is finished. Finally, $\mathrm{\ell }$ can be $p$ to capture the case that the left side is empty.

The optimization of the $R$-type using a step on the left side is:

Here, ${\mathrm{\ell }}^{\prime }$ must be such that $q\mathrm{\ell }{\mathrm{\ell }}^{\prime }$ is a left turn at $\mathrm{\ell }$, vertex ${\mathrm{\ell }}^{\prime }$ is to the left of the line through $\mathrm{\ell }$ and $p$ and right of the line through $\mathrm{\ell }$ that is parallel to the line through $r$ and $r^{\prime }$ (showing that ${\mathrm{\ell }}^{\prime }$ and $r$ are not antipodal), and $|{\mathrm{\ell }}^{\prime }{r}^{\prime }|\le |pq|$. We can choose ${\mathrm{\ell }}^{\prime }$ to be $p$, which would complete the left boundary. Also here, $r^{\prime }$ might be $q$.

As base case of the recursion, we take all triangles including the points $p$ and $q$. In particular, we define $R[r_1,q,p,k]=A(p,r_1,q)$ if $k=3+K(p,r_1,q)$, and $\mathrm{\infty }$ otherwise.

The definition of $L[\mathrm{\ell },{\mathrm{\ell }}^{\prime },r,k]$ is symmetric, having analogous options. Note that since $q$ is the unique highest point and $p$ the unique lowest point, “left” and “right” of any line through $p$ (or $q$) and another point is well-defined.

Intuitively, the algorithm constructs an optimal solution as follows, see Figure 5. Consider a triangulation of the solution $C$ for a given diameter $p,q$ where $p$ has chords to all vertices on the right and $q$ has chords to all vertices on the left. Let $p=r_0,r_1,\mathrm{\dots },r_j=q$ be the vertices incident to edges on the right boundary. Then the edges $r_1{r}_2,\mathrm{\dots },r_{j-1}{r}_j$ are the sequences of edges of the triangles in the right half that are opposite to $p$ in these triangles. Let $q={\mathrm{\ell }}_0,\mathrm{\dots },{\mathrm{\ell }}_h=p$ be defined analogously for the left half.

Consider the sorted sequence (by positive turn angle from the positive $x$-axis vector) of edges ${\overrightarrow{e}}_1,\mathrm{\dots },{\overrightarrow{e}}_{j+h-2}$ from $\overrightarrow{r_1{r}_2},\mathrm{\dots },\overrightarrow{r_{j-1}{r}_j}$ and $\overrightarrow{{\mathrm{\ell }}_2{\mathrm{\ell }}_1},\overrightarrow{{\mathrm{\ell }}_3{\mathrm{\ell }}_2},\mathrm{\dots },\overrightarrow{{\mathrm{\ell }}_h{\mathrm{\ell }}_{h-1}}$. This sorted sequence defines an order on the triangles in the left and right halves of $C$ that is a merge of the two sequences. The solution $C$ can be deconstructed by peeling off the triangles according to that ordered sequence. If triangle $\mathrm{△}p{r}_1{r}_2$ is first in the sequence, then $C$ is composed of this triangle and the convex polygon $C$ where that triangle is removed. As we argue next, this solution is found when computing $R[r_1,r_2,{\mathrm{\ell }}_1,k]$, where $k$ is the number of points of $P$ in $C$.

As a consequence of this lemma, $r_1$ is not antipodal to ${\mathrm{\ell }}_2$, so removing $r_1$ from consideration the recursive problem does not remove any antipodal pairs in $C$ that must be checked.

We use three types of data for the experiments:

• $\mathrm{■}$

  Uniformly distributed point sets inside a square of size $20\times 20$ mm² of increasing cardinalities $100,\mathrm{\hspace{0.17em}110},\mathrm{\dots },200$.

• $\mathrm{■}$

  Normally distributed point sets of 100 points centered inside a square of size $20\times 20$ mm² of increasing standard deviations ranging from $\sigma =0.5,1.0,1.5,\mathrm{\dots },6.5$ mm. Points outside the square are discarded and re-generated.

• $\mathrm{■}$

  Medical data used in [9]. We selected ten data sets with up to $700$ points, in particular the ten most populated sets.

We will always search for an area of at most $4$ mm².¹¹1While [9] lists 2 mm² as the desired area of the region, personal communication with the authors of the paper showed that they use 4 mm² currently in the application. For the AD algorithm, we examine maximum diameters of $2$, $3$, $4$, $5$, and $6$ mm. Note that for a diameter of $2$ mm, the maximum area possible is $\pi$ mm², so in this setting, we cannot “use” the full allowed area of $4$ mm².

Since the first two types of data are randomly generated, we generated $100$ point sets in each setting and report averages over these $100$ runs.

The real-world data is obtained by preprocessing the data set provided by the Department of Digital Pathology at UMC Utrecht, containing $1982$ point sets representing the mitotic figures identified in Whole-Slide-Images (WSIs) from actual patients [9]. The identification was carried out by the classification model introduced in [11], which produced $2$-dimensional points with detection accuracies, since the points have a different likeliness of being a mitotic cell. We apply a preprocessing step to filter out any point with detection accuracy below $0.86$ and convert the coordinates from pixels to millimeters. We consider the top $10$ most populated point sets for our experiments. Table 1 provides an overview of this data. The average nearest neighbor distances provide an indication of how clustered the points are.

Both the synthetic generation and preprocessing of real-world data are fully automated in a Python script which we provide in a repository.

We implemented the A algorithm from [5] and the AD algorithm in C++²²2Our code is available at https://github.com/gianmarcopicarella/area-diameter-algorithm.. We implemented all other necessary data structures and algorithms that are not already in the std ourselves. The implementations do not use parallelism and run in a single thread. The experiments were run on a MacBook Pro using an Apple M1 Pro (3228 MHz) and 16GB (6400 MT/s LPDDR5) RAM. It runs macOS Sonoma 14.2.1, and we compiled it for 64 bits using clang 19.1.6 with the -O3 option. We used the Google Benchmark library to set up our experiments and perform data collection. The data is then post-processed in Python with several scripts that we provide in our repository. Regarding the Area Selection method (AS), we used the Python implementation made available by [9]. We use this implementation mostly to compare output, not to compare resources, although we report times.

The AS algorithm from [9] first partitions the data in patches of $3\times 3$ mm, allowing overlap so that solutions across the boundary of a patch are not missed. For each patch, if it contains at most $25$ points, a brute-force method determines the optimum. Otherwise, a heuristic is used that takes the convex hull of a patch and incrementally discards points from the convex hull, until the criteria are met. More precisely, it considers each convex hull point $p$, analyzes the decrease in area and cardinality if it were discarded, but also if it were replaced by any point inside the triangle formed by $p$ and its two neighbors (details in [9]).

The A and AD algorithms use a few optimizations. First, we preprocess for quickly determining the point count in triangles, so that the dynamic program can avoid spending more than constant time. We use only $O(n^2)$ time for preprocessing and storage by computing the number of points below each of the $\left(\genfrac{}{}{0pt}{}{n}{2}\right)$ possible edges. To represent all table entries of the AD algorithm, we used hashing, since many table entries are not used. In particular, we use an array with $k$ hash tables for $R$ and for $L$: to retrieve $R[r,r^{\prime },\mathrm{\ell },k^{\prime }]$ we inspect the $k^{\prime }$-th hash table and use the indices of $r$, $r^{\prime }$ and $\mathrm{\ell }$ as the key, and similar for $L$. In the A algorithm, most table entries are used, so the A algorithm does not benefit as much from hashing. Hence it is not used. Finally, both the A and AD algorithms use an early exit strategy. For the AD algorithm point pairs are treated in decreasing order of number of points in the lune. If a solution of cardinality $k$ is known, then all point pairs whose lune contains fewer points can be skipped. For the A algorithm, the early exit strategy is much less strong, but we can still skip processing points $p$ as bottommost point if there are not enough points above it to improve the current best solution. Any preprocessing time is included in the timings. The implementations handle degenerate cases like collinear points properly.