Find the most points enclosed in a fixed size circle

Question

This came up when a friend talked about a programming competition, and we wondered what the best approach was:

Given a list of points, find the centre of a circle of predetermined size that covers the most points. If there are several such circles, its only important to find one of them.

Example input: 1000 points, in a 500x500 space, and a circle of 60 diameter.

Answer 1

Very quick idea, not necessarily right one:

• for each point P you calculate a "candidate covering area" - a continuum of points where the center of its covering circle might be. Naturally it is also a circle of diameter D with the center in P.
• for each point P you intersect its candidate covering area with the corresponding areas of other points. Some of the candidate covering areas may intersect with the P's and with each other. For each intersection you count the number of intersected areas. A figure which is intersected by the most of candidate areas is a candidate area for the center of a covering circle which covers P and as many of other points as possible.
• find a candidate area with the highest number of intersections

Seems to be N^2 complexity, provided that calculating intersections of circle shaped areas is easy

Answer 2

Hi

How about using a clustering algorithm to identify the cluster of points. Then ascertain the cluster with the maximum number of points. Take the mean point of the cluster having the maximum points as the center of your circle and then draw the circle.

MATLAB supports implementation of k-means algorithm and it gives back a 2-d array (a matrix to be precise) of cluster means and corresponding cluster ids.

One well known flip side of k-means is deciding on k(number of clusters) before hand. This can be resolved though - one can learn the value of k from the data points. Please check this paper.

I hope this helps.

cheers

Answer 3

Unless I've missed something obvious I think there is a simple answer.

For a rectangular area MxN, number of points P, radius R:

• Initialise a map (e.g. 2D array of int) of your MxN area to all zeroes
• For each of your P points
  • increment all map points within radius R by 1
• Find map element with maximum value - this will be the centre of the circle you are looking for

This is O(P), assuming P is the variable of interest.

Answer 4

My best approach so far is:

Every circle containing points must have a left-most point. So it makes a list of all the points to the right of a point that are potentially within the bounds of a circle. It sorts the points by x first, to make the sweep sane.

It then sorts them again, this time by the number of neighbours to the right that they have, so that the point with the most neighbours get examined first.

It then examines each point, and for each point to the right, it computes a circle where this pair of points is on the left perimeter. It then counts the points within such a circle.

Because the points have been sorted by potential, it can early-out once it's considered all the nodes that might potentially lead to a better solution.

import random, math, time
from Tkinter import * # our UI

def sqr(x):
    return x*x

class Point:
    def __init__(self,x,y):
        self.x = float(x)
        self.y = float(y)
        self.left = 0
        self.right = []
    def __repr__(self):
        return "("+str(self.x)+","+str(self.y)+")"
    def distance(self,other):
        return math.sqrt(sqr(self.x-other.x)+sqr(self.y-other.y))

def equidist(left,right,dist):
    u = (right.x-left.x)
    v = (right.y-left.y)
    if 0 != u:
        r = math.sqrt(sqr(dist)-((sqr(u)+sqr(v))/4.))
        theta = math.atan(v/u)
        x = left.x+(u/2)-(r*math.sin(theta))
        if x < left.x:
            x = left.x+(u/2)+(r*math.sin(theta))
            y = left.y+(v/2)-(r*math.cos(theta))
        else:
            y = left.y+(v/2)+(r*math.cos(theta))
    else:
        theta = math.asin(v/(2*dist))
        x = left.x-(dist*math.cos(theta))
        y = left.y + (v/2)
    return Point(x,y)

class Vis:
    def __init__(self):
        self.frame = Frame(root)
        self.canvas = Canvas(self.frame,bg="white",width=width,height=height)
        self.canvas.pack()
        self.frame.pack()
        self.run()
    def run(self):
        self.count_calc0 = 0
        self.count_calc1 = 0
        self.count_calc2 = 0
        self.count_calc3 = 0
        self.count_calc4 = 0
        self.count_calc5 = 0
        self.prev_x = 0
        self.best = -1
        self.best_centre = []
        for self.sweep in xrange(0,len(points)):
            self.count_calc0 += 1
            if len(points[self.sweep].right) <= self.best:
                break
            self.calc(points[self.sweep])
        self.sweep = len(points) # so that draw() stops highlighting it
        print "BEST",self.best+1, self.best_centre # count left-most point too
        print "counts",self.count_calc0, self.count_calc1,self.count_calc2,self.count_calc3,self.count_calc4,self.count_calc5
        self.draw()
    def calc(self,p):
        for self.right in p.right:
            self.count_calc1 += 1
            if (self.right.left + len(self.right.right)) < self.best:
                # this can never help us
                continue
            self.count_calc2 += 1
            self.centre = equidist(p,self.right,radius)
            assert abs(self.centre.distance(p)-self.centre.distance(self.right)) < 1
            count = 0
            for p2 in p.right:
                self.count_calc3 += 1
                if self.centre.distance(p2) <= radius:
                    count += 1
            if self.best < count:
                self.count_calc4 += 4
                self.best = count
                self.best_centre = [self.centre]
            elif self.best == count:
                self.count_calc5 += 5
                self.best_centre.append(self.centre)
            self.draw()
            self.frame.update()
            time.sleep(0.1)
    def draw(self):
        self.canvas.delete(ALL)
        # draw best circle
        for best in self.best_centre:
            self.canvas.create_oval(best.x-radius,best.y-radius,\
                best.x+radius+1,best.y+radius+1,fill="red",\
                outline="red")
        # draw current circle
        if self.sweep < len(points):
            self.canvas.create_oval(self.centre.x-radius,self.centre.y-radius,\
                self.centre.x+radius+1,self.centre.y+radius+1,fill="pink",\
                outline="pink")
        # draw all the connections
        for p in points:
            for p2 in p.right:
                self.canvas.create_line(p.x,p.y,p2.x,p2.y,fill="lightGray")
        # plot visited points
        for i in xrange(0,self.sweep):
            p = points[i]
            self.canvas.create_line(p.x-2,p.y,p.x+3,p.y,fill="blue")
            self.canvas.create_line(p.x,p.y-2,p.x,p.y+3,fill="blue")
        # plot current point
        if self.sweep < len(points):
            p = points[self.sweep]
            self.canvas.create_line(p.x-2,p.y,p.x+3,p.y,fill="red")
            self.canvas.create_line(p.x,p.y-2,p.x,p.y+3,fill="red")
            self.canvas.create_line(p.x,p.y,self.right.x,self.right.y,fill="red")
            self.canvas.create_line(p.x,p.y,self.centre.x,self.centre.y,fill="cyan")
            self.canvas.create_line(self.right.x,self.right.y,self.centre.x,self.centre.y,fill="cyan")
        # plot unvisited points
        for i in xrange(self.sweep+1,len(points)):
            p = points[i]
            self.canvas.create_line(p.x-2,p.y,p.x+3,p.y,fill="green")
            self.canvas.create_line(p.x,p.y-2,p.x,p.y+3,fill="green")

radius = 60
diameter = radius*2
width = 800
height = 600

points = []

# make some points
for i in xrange(0,100):
    points.append(Point(random.randrange(width),random.randrange(height)))

# sort points for find-the-right sweep
points.sort(lambda a, b: int(a.x)-int(b.x))

# work out those points to the right of each point
for i in xrange(0,len(points)):
    p = points[i]
    for j in xrange(i+1,len(points)):
        p2 = points[j]
        if p2.x > (p.x+diameter):
            break
        if (abs(p.y-p2.y) <= diameter) and \
            p.distance(p2) < diameter:
            p.right.append(p2)
            p2.left += 1

# sort points in potential order for sweep, point with most right first
points.sort(lambda a, b: len(b.right)-len(a.right))

# debug
for p in points:
    print p, p.left, p.right

# show it
root = Tk()
vis = Vis()
root.mainloop()