Maximum disc containing no more than k points

Hello. I hope this is the right place to ask (and if not could
somebody kindly direct me to where it is?)

I have the following problems: given a (large, roughly 9M) set of
points, I need to find for each point the largest disc centered at
that point that contains no more than "k" other points.

If there is no other way I can live with "the largest square" instead
of "the largest disc".

This problem can be solved off-line and not necessarily in real-time/
on-line.

Another variation on the problem that I have: each point has an
integer positive weight, and I need to find the largest disc/square
that contains no more than "k" sum-weights (that is: the sum of
weights of points in the inscribed region is no more than k).

tia!

Hi,

It is.


Equivalent to: find the k:th nearest neighbor under the Euclidean norm.


Same problem with maximum norm.


Essentially the same problem (above the weights can be thought to be 1).

In low to moderately high dimension (say 32) you'll do very fine by 
constructing a kd-tree and searching in that. For higher dimensions a 
good answer is brute force.

-- 
 http://www.**--****.com/ 

The simple method requires N*N ops (brute force). There is a way to
solve this problem by a smaller number (roughly proportional N).
Assume we have a 32 bpp bitmap. Put all of your points onto it.
Initially this bitmap is filled by zeroes. When you put a point number
i onto this bitmap, zero at corresponding point is to be replaced by
number k (number starts from one). Then you need to check each point
and find a circle around it, which contains no more than k another
points. It can be done quickly using some simple approaches (check
circle with radius r1, if it is too small multiply r1 by 1.5, if r1 is
too large divide it by 1.5 and so on. When you have two radiuses -r1
and r2 then check an average radius). There are two problems:
coordinates will be rounded and two (or more ) points will be placed
just onto one pixel. The first problem is rather simple- after you
have got set of points with rounded coordinates, check their precise
values. The second problem requires another approach. Simple (and not
elegant) is to create two, three bitmaps and if a pixel at bitmap
number one is already occupied, place a point to the bitmap number two
and so on. Much better approach is to place onto bitmap not numbers,
but pointers to buffers with numbers. This method probably will be
rally fast. You needn't have a large bitmap- ropugh approximation of
actual coordinates will be more than enough. It is an interesting task
for a programmer- to quickly allocate all beforenamed buffers.

A grid generalizes this idea: pick an axis-aligned box to bound the 
points. Then divide it into a grid of equal-sized voxels (where each 
voxel is an aligned box). In each voxel keep up a list of points that 
are contained in it. No rounding is done.

Grids are good for point data that is uniformly distributed: no voxel 
should contain many points, and there should be not many empty voxels. 
Due to exponentially growing storage needs, grids are not suitable for 
dimensions > 3.


Using a grid, the k nearest neighbors are found by examining voxels in 
order of distance from the query point until it is known that all the 
rest of the voxels can not contain closer candidates.

-- 
 http://www.**--****.com/ 

I've just noticed, that what I'm talking about solves not the original
problem, but similar- how many points lays within a circle with
predefined radius R. Let's see, which approach will be the best for
the original problem (to find radius of a circle, what contents k
points) in 1D event. The quick sort algorithm, of course. Once we have
an ordered set of values, all what we need is to check 2*k neighbors
of some point (k to the left side of its position and k to the right
one). Similar method can be used for 2D event too. Let's make two
copies of our set of points. One will be ordered accordingly to x-
coordinate, the second to y- coordinate. This method gives us
significant economy too. Assume all points lays inside a square. Using
this approach, we need to check not the whole square, but a cross
within it. It means, that total number of ops will be proportional not
to N*N, but N*sqrt(N). It seems to me, that this method will be better
than grid based- at least it not depends on the way in what points are
distributed.

In 1D, the problem can be solved in O(n log n), by sorting as you say. 

In 2D, the problem can be solved in O(n log n) too, by using Voronoi 
graphs.

But even in 2D, I'd go for a kd-tree. While there are no non-trivial 
worst-case analysis of kd-trees (as far as I know), they are a very 
practical tool for approximity queries. Expected time analyses show 
O(n log n) under some distribution assumptions. Of course, as 
dimensionality increases, the gap between the performance of kd-tree and 
brute-force gets smaller.

There are a lot of open questions on what can be assumed and what not. A 
practical approach allows for different distance measures, scales to 
higher dimensions, allows for dynamic updates on the point set, and 
works efficiently.

-- 
 http://www.**--****.com/ 

may be, method, what I described, requres N*log(N) too (it is an
average number of ops, of course- not for the worst case scenario). .
N*sqrt(N) will be if one needs to check the whole cross. Real
efficiency depends on many things, some of them hardly can be known
without experiments. There is a gap between pure theory and practical
realizations. I met not once situation, when there is no adequate
theory for some events. I wrote a subroutine for quicksort, using
different approaches (recursive call and just-in-place sorting). My
input set was relatively large- millions of input points, and
recursive call was not the best- it requires too large stack. Using
assembler, I rewroute this procedure without chains of calls. This one
was better than written using C- assembler variant worked three times
faster.

Thank you all for the answers so far. Kaba's original reply was right
and I should have understood it right away -- what I need a simple
modification on KNN.
I went back and brushed up on my kd-tree theoretical stuff, and I
understand that it is O(n*sqrt(n)).

After a little bit of reading I found out about cover trees, which
perform in n*log(n). I am now looking to generalize NN to kNN on cover
trees, which I think I have.
The only(?) question remaining is whether there is an all-kNN
algorithm for cover trees, which is what I really need, and whether it
can be easily modified for my case (where I require that
sum(weight(P'))<k, not that |P'|<k). (Reminder: each point has a
weight, and I want the largest radius that contains points whose total
weight < k, I am not looking for the radius where num points < k).

Thanks again!

That's worst case (for 2D). The worst case almost never applies, it is 
very conservative. You should expect O(n log n) behaviour, particularly 
in low dimensions such as 2. In practice, you will have hard time 
finding a faster algorithm for this problem. Of course, you need to be 
careful on the implementation. 

Something that works really good:
* A kd-tree, where the leaf nodes contain the points, intermediate nodes 
do not.
* The sliding midpoint rule is used to split nodes (on the axis on which 
the node is longest).

-- 
 http://www.**--****.com/ 

I've just made a small experiment. Input data contents 190000 points.
the first attempt was to find for each point a circle with k=20
another points within it. It took 60 seconds, and average number of
points what was checked was 793. The second attempt was with k=40. It
took 90 second, and averale number was 1135. My comp has a single-core
Athlon 1300 GHz and less than 500 M of a memory. Beat me, boys, if you
can.

I have Intel Q6600 2.40Ghz quad-core. I'm using only a single core. The 
point sets are generated from a gaussian distribution with mean 0 and 
standard deviation 1, the hardest test set I know. The norm is the 
Euclidean norm. I am searching for 20 nearest neighbors.

Using a kd-tree I described:

2D:

190000 points: 1.984s.
380000 points: 4.047s.
760000 points: 8.297s.

I am doubling the number of points in each step. It can be seen that the 
performance scales as O(n log n).

4D:

190000 points: 5.703s.
380000 points: 11.66s.
760000 points: 23.66s.

Performance scales as O(n log n).

6D:

190000 points: 21.89s.
380000 points: 46.13s.
760000 points: 98.17s.

As the dimension gets larger, the performance deteriorates from the 
O(n log n) towards O(n^2). It has already happened here, showing 
behaviour of O(n^(1.13)).

To test the 8D, I will now switch on all of the four cores to save time.

8D:

190000 points: 36.19s.
380000 points: 91.06s.
760000 points: 213.2s.

This shows behaviour similar to O(n^(1.23)). This is still far from the 
worst case O(n^(1.88)).

In summary, kd-trees a practical tool for nn searching. Explaining their 
surprising performance theoretically is still largely an open problem.

-- 
 http://www.**--****.com/ 

I don't understand what you are asking. Maybe how the searching 
algorithm on a kd-tree works? Here's how:

0. Set the culling distance to infinity (or to some number for fixed 
distance searches), and the current node to the root node.

1. From the current node, traverse the tree to the leaf node closest to 
the query point. On the way down, store the nodes you do _not_ visit in 
a priority queue based on distance.

2. Examine the points in the leaf node. Keep up accumulating an ordered 
list of nearest candidates. Once you have more than k of them, throw the 
excess out, and set the culling distance to the current k:th nearest 
neighbor.

3. If there are no more nodes in the priority queue, you are done. 
Otherwise, retrieve the closest node from the priority queue. If it is 
beyond the current culling distance, do step 3 again. Otherwise, go to 
step 1.

The trick is that the hierarchical space subdivision can throw out whole 
regions of points without ever considering them. 


For a space subdivision structure, clusters are easy because they are 
separated by empty space: empty space can be skipped trivially, and 
makes culling of large portions of points possible.


For the 2D point set, 61 on average.

-- 
 http://www.**--****.com/