# Find a unit square containing most of the points

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?

• Sort all your points by pairs of their coordinates in a lexicographical order. Then create a sliding window of size 1x1. Move this window from top to bottom and from left to right so that the lines through its left and top sides contain points from a given set. If you move it in a monotone (in some appropriate sense) manner then you can recount number of points within a window so that every point will be viewed no more than twice. Commented May 6, 2023 at 21:23
• @VladislavBezhentsev, would you care to write a full answer? I don't think I understand what you mean by a 1x1 window, by moving in a monotone manner, or by recounting number of points within a window.
– D.W.
Commented May 6, 2023 at 21:33
• Well, if I managed to sort all the squares in a way such the set of squares each point belongs to forms a segment in this sorted sequence, then this can be solved using segment trees. However, I don't see, how to sort it. I think this somehow corresponds to the idea of the sliding window. Commented May 6, 2023 at 21:39
• By sweepline: sort the points by increasing $y$ then move a sliding window of height $1$; it contains at most $n$ points and is moved at most $n$ times. For every position, maintain the points by increasing $x$. The total work is $O(n\log n)$. And finally, for a given position, scan horizontally with a sliding window of width $1$. This approach takes time $O(n^2\log n)$ at worse.
– user16034
Commented May 7, 2023 at 7:10
• @MichalDvořák Nit-picking: the observation is missing the case when the optimal solution contains only $1$ point or all given points in each optimal solution are on the same horizontal line or on the same vertical line. Commented May 7, 2023 at 11:32

"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.

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.

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".