Find a Circle ((x,y,r)) that has maximum number of points 'on' it; given a set of points(x,y) in a 2D plane

Question

Given N points on a plane (of form (x, y)), find the circle that have maximum number of points on it ? P.S. : The point should lie on the circumference of the circle. What is the most efficient algorithm to solve this problem and how does it work? Which Data structures would you use to solve this problem. This was asked in one of the FANG coding interviews.

Answer 1

As a starting point, the simple O(N³) solution is to find the circle corresponding to each unique triple of points, while counting the number of occurrences of each circle you find.

If a circle has N points on it, then you will find it N-choose-3 times, so the circle you find most often is the one with the most points on it.

There are complications in any actual implementation, but they are different complications depending on how your points are represented and whether you want exact or approximate answers.

Answer 2

Case 1 : Hough Transform

In computer vision problem solving we often search for circles amongst edge information. This problem is characterised by having many data points, possibly originating from many different circles in the presence of a great deal of noise.

The usual approach to solving this problem is the Hough Transform https://en.wikipedia.org/wiki/Circle_Hough_Transform. The basic idea is to sum the evidence for the circles which can pass through each point (x, y).

We create an integer array called Hough [a, b, r] which parametrises all the possible circles that can pass through your point (x,y). This is equivalent to drawing a circle of radius 1 in the r=1 plane centred on (x,y); a circle of radius 2 in the r=2 plane centred on (x,y) etc.

Each time a circle is drawn through a point in [a, b, r] we add 1 to the corresponding value. Some points accumulate a lot of evidence. These points correspond to the circles of interest.

Image from cis.rit.edu illustrates what happens in one of the r-planes.

Doing this for each point (x,y) will generate evidence towards each of the points in [a,b,r] corresponding to the circle you seek. So just scan this array to find the point with the maximum evidence. That is your circle.

Example of Hough Transform

Knowing the radius of the circle reduces this from an O(n^3) problem to a O(n^2) problem as only one plane needs to be constructed and scanned. I have also had good results plotting the radii in a log space to give a (less accurate) O(n^2 log n) algorithm.

Case 2 : Circle Fitting

If the points are known to lie near the boundary of a single circle, and/or if the points are not very numerous and/or else we are very sure that there is very little noise then the Hough transform is a poor solution as it is computationally intensive, memory hungry and because of the raster nature of the accumulator array possibly not very accurate.

In this case we may want to fit a circle by analogy to line fitting techniques that use linear regression. A discussion of circle fitting techniques can be found in https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf.

A (rather simple minded) implementation of this algorithm is presented below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct { 
    float x;
    float y;
    } point;

/*
 * Function using a modified least squares approach to circle fitting.
 *
 * Reference :
 *
 * Umbach, D. and Jones, K. N. "A few methods for fitting circles to data", 
 * IEEE Trans Instrumentation and Measurement
 * vol XX(Y) 2000
 *
 * https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf
 *
 * NOTES 
 *
 * The code below has not been checked for numerical stability or conditioning.
 */

int circle_fit_MLS (point P[], int n, double *x_pos, double *y_pos, double *radius)
{
 int i;
 double sum_x=0.0, sum_y=0.0, sum_xx=0.0, sum_xy=0.0, sum_yy=0.0, sum_xyy=0.0, sum_yxx=0.0, sum_xxx=0.0, sum_yyy=0.0;
 double A, B, C, D, E;
 double x2, y2, xy, F, xdif, ydif;


 for (i=0; i<n; i++)
 {
    sum_x += P[i].x;
    sum_y += P[i].y;
    x2 = P[i].x * P[i].x;
    y2 = P[i].y * P[i].y;
    xy = P[i].x * P[i].y;
    sum_xx += x2;
    sum_yy += y2;
    sum_xy += xy;
    sum_xyy += xy*P[i].y;
    sum_yxx += P[i].y*x2;
    sum_xxx += x2*P[i].x;
    sum_yyy += y2*P[i].y;
 }

 A = n * sum_xx - sum_x * sum_x;
 B = n * sum_xy - sum_x * sum_y;
 C = n * sum_yy - sum_y * sum_y;
 D = 0.5 * ( n * (sum_xyy + sum_xxx) - sum_x * sum_yy - sum_x * sum_xx);
 E = 0.5 * ( n * (sum_yxx + sum_yyy) - sum_y * sum_xx - sum_y * sum_yy);

 F = A*C - B*B;
 *x_pos = (D*C - B*E) / F;
 *y_pos = (A*E - B*D) / F;

 *radius = 0;
 for (i=0; i<n; i++)
 {
     xdif = P[i].x - *x_pos;
     ydif = P[i].y - *y_pos;
     *radius += sqrt(xdif * xdif + ydif * ydif);
 }
 *radius /= n;

 return 0;
}

The main program below can be used for testing the code. Please post back any results/observations/suggestions for improvement in the comments.

int main()
{
 point *P;
 int    n, i;
 double xpos, ypos, radius;

 printf ("Please enter the number of points \n> ");
 scanf ("%d", &n);
 P = malloc (n * sizeof(point));

 for (i=0; i<n; i++)
 {
     printf ("x%d = ", i);
     scanf ("%f", &P[i].x);
     printf ("y%d = ", i);
     scanf ("%f", &P[i].y);
 }

 circle_fit_MLS (P, n, &xpos, &ypos, &radius);

 printf (" a = %f\n", xpos);
 printf (" b = %f\n", ypos);
 printf (" r = %f\n", radius);
}