Find a unit square containing most of the points

Question

Given points $p_1,\ldots,p_n$ in $\mathbb{R}^2$, the task is to find an axis-aligned unit square containing the maximum number of points. I came up with and $O(n^3)$ algorithm as follows.

Observation 1: For $n\geq 2$, there is an optimal solution containing at least two points on the boundary not lying on a line parallel to the coordinate axes.

Pf. Suppose not, say it contains only 1 or all are on one boundary. Then we may move the rectangle along the boundary until we hit some point.

We can take all pairs and for each pair $p_i = (x_1,y_1)$ and $p_j = (x_2,y_2)$, where $x_1\neq x_2$ and $y_1\neq y_2$. Suppose wlog that $y_2>y_1$ and $x_2>x_1$. Other cases are symettrical. we try to find the intersection $a_1$ of lines $x=x_1$ and $y=y_2$ and try to place the upper left corner to $a_1$ and then intersection $a_2$ of lines $x=x_2$ and $y=y_1$ and place bottom right corner to $a_2$.

In this way, we get $O(1)$ candidate squares for each of $O(n^2)$ pairs of points. We can calculate the number of points contained in each of them in $O(n)$ time, givin us $O(n^3)$ algorithm.

For a sanity check, is the above algorithm correct?

Now, is there a faster way to solve this?

Answer 1

"A window" will also mean "an axis-aligned unit square".

As observed in the question, there is an optimal window such that

• either the bottom left corner is a given point.
• or the left side contains a given point not on the bottom side and the bottom side contains a given point not on the left side.

Hence we need to check only the windows with bottom left corners at $(x_i, y_j)$. There are $n^2$ windows (with possible duplicates among them) to check.

---

Here is a simple $O(n^2)$-time algorithm. Basically, we apply the sliding interval algorithm to check all windows with their bottom sides on the same line.

The input is $p_i=(x_i,y_i)$, $1\le i\le n$.

1. Sort $p_i$'s by $x$-coordinate.
   Now $p_1, p_2, \cdots$ are sorted by $x$-coordinate.
2. For $i$ from $1$ to $n$ (which is meant for $y_i$), do the following.
   1. Initialize an empty queue that will be used to track the points in the sliding interval.
   2. Run the usual sliding interval algorithm on the array of numbers $x_1,x_2,\cdots, x_n$ using the dequeue and the sliding interval that will be $[x_1, x_1+1]$, $[x_2,x_2+1]$ and so on successively. However, when the interval is $[x_j, x_j+1]$, a number $x_k$ (with $k$ noted) in the dequeue will be counted if and only if $(x_k, y_k)$ is in the window with bottom left corner at $[x_j, y_i]$, i.e., $y_i\le y_k\le y_i+1$. Track the maximum number of counted numbers in the sliding interval.
3. Return the maximum number of all maximum numbers obtained in step 2.2.

---

It takes $O(n\log n)$ time to sort. Each run of the sliding interval algorithm costs $O(n)$ time. All $n$ runs cost $O(n^2)$ time. Other computations are insignificant in time-cost. Hence the time-complexity of the algorithm is $O(n^2)$.

Since the "sliding" of the intervals is the same for all runs, we could speed up the algorithm somewhat. However, that alone does not reduce the time-complexity to below $O(n^2)$ because of work needed to track the (number of) points in each sliding interval.

Answer 2

Yes, the algorithm is correct. It can be sped up by using implementing the 'count number of points within a box' using smarter data structures. In particular, a running time of $O(n^2 \log n)$ is achievable.

Let $N(x_\ell,y_\ell,x_u,y_u)$ denote the number of points within the box $[x_\ell,x_u] \times [y_\ell,y_u]$, and define $N'(x_u,y_u) := N(-\infty,-\infty,x_u,y_u)$. Note that

$$N(x_\ell,y_\ell,x_u,y_u) = N'(x_u,y_u) - N'(x_u,y_\ell) - N'(x_\ell,y_u) + N'(x_\ell,y_\ell).$$

I will describe below two ways to compute $N'(\cdot,\cdot)$ efficiently, which immediately leads to more efficient versions of your algorithm.

Pragmatic approach

Store all points in a 2D range tree, 2D interval tree, or 2D segment tree. I'm not sure, but I think these data structures allow you to efficiently compute $N'(\cdot,\cdot)$, though I'm not sure if there are any worse-case running time bounds.

Worst-case complexity

If you care about provable worse-case complexity, here is an approach that will enable you to solve your problem in $O(n^2 \log n)$ time, using a sweep line algorithm and persistent data structures (a standard technique for algorithms on 2D points).

Build a self-balancing binary tree with one leaf per unique $y$-coordinate of the points (thus at most $n$ leaves). We sweep a vertical line from left to right. We'll have one tree per $x$-coordinate of the points (with repetition). Each node of the tree corresponds to a range of $y$-coordinates (correspondings to the leaves that are descendants of that node). We'll augment the tree for $x$-coordinate $x_u$ so that the node corresponding to range $[y_\ell,y_u]$ stores the count $N(-\infty,y_\ell,x_u,y_u)$.

Start with the leftmost $x$-coordinate, and build the tree for this $x$-coordinate. Then, as you sweep the line to the right by one position, you pick up one additional point $(x_u,y_u)$, so you need to increment the count for the leaf corresponding to $(x_u,y_u)$, as well as ancestors of that leaf. Instead of modifying the tree in place, use a persistent data structure to implement this modification, so we retain the tree before this modification and the tree afterwards; these trees differ in about $\lg n$ nodes and the remaining $n-\lg n$ nodes are common to both trees. Continue sweeping the sweep line to the right, modifying the tree as you encounter each point.

Once you have built all the trees, you can compute $N'(x_u,y_u)$ by looking up the tree for $x$-coordinate $x_u$, then finding the set of internal nodes whose ranges form a partition of $[-\infty,y_u]$. This a set of $O(\log n)$ nodes. Summing up the counts in those nodes gives us the value of $N'(x_u,y_u)$.

What is the running time? It takes $O(n \log n)$ time to build the persistent data structure with all the trees. Also, to compute $N'(x_u,y_u)$, you can find the tree for a particular point, find the partition of nodes, and sum up all of them in $O(\log n)$ time. Therefore, computing $N(x_\ell,y_\ell,x_u,y_u)$ can be done in $O(\log n)$.

Therefore, the total running time of the algorithm becomes $O(n^2 \log n)$, as you have to compute $N(\cdot)$ for $n^2$ pairs of points.

I suspect perhaps there is a clever algorithm to achieve $O(n^2)$ running time, but I can't see how at the moment.

Answer 3

For every input there exists a unit square with maximal number of points within such that there is a point on the bottom and on the left sides of the square (it can be the same point which coincides with bottom left corner of the square). So let's consider only squares with this property.

We can sort all x-coordinates of given points and move from least to greatest coordinate. For coordinate $x_i$ let's consider a vertical line $l$ given by equation $x = x_i$. Let's also consider a vertical strip $Strip_l$ which is a part of plane between vertical lines $x = x_i$ (aka $l$) and $x = x_i + 1$. Let's search for unit squares which left side lies on line $l$ and which bottom side lies on line $y = y_k$, where $y_k$ is a coordinate of some of the given points which lie within current $Strip_l$. We will such square $sq_{y_k}^l$.

Now let's assume that we have a data structure which can maintain a set of y-coordinates of points, which currently lie within $Strip_l$, and for each such coordinate $y_k$ store a number $num$ of points within a unit square $sq_{y_k}^l$. That is, we should be able:

1. To put a pair $(y, num)$ to this data structure. And to update values associated with y-coordinates, which are already in the data structure, appropriately.
2. To remove a coordinate $y$ from this data structure.

Moreover, we should be able query a maximum value which is currently stored in the data structure. Note, that a question about whether a given point from $Strip_l$ with coordinate $y$ lies within $sq_{y_k}^l$ for a fixed $Strip_l$ is equivalent to a question whether $y$ lies within a one-dimensional segment $[y_k, y_k+1]$.

But now we got exactly the same problem which is described in this question.

So we can sort all given points by their coordinates. And then move $Strip_l$ from left to right along a plane, maintaining a data structure that we've just described. Each point will be accessed no more than 2 times: put and remove from a data structure. And processing of a single point will take $O(\log n)$ time. So the total time complexity of this algorithm will be $O(n\log n)$.

Note, that in a comparison-based model this solution is optimal. Indeed, let $a_1, \ldots, a_n$ be an instance of an element distinctness problem. Then convert it to an instance $(2a_1, 0), \ldots, (2a_n, 0)$ of a problem from the current question. Then the answer for this problem equals 1 iff the answer for original element distinctness problem is "all elements are distinct".