OpenGL: fitting a circle boundary to set of 2D points

OpenGL: fitting a circle boundary to set of 2D points

OpenGL programmers,

Here is a nice small program that lets the user interactively
specify a set of 2D points and then automatically determines
the circle boundary that best fits those points using a least
squares fit algorithm implemented by Dave Eberly.

You'll need OpenGL (or Mesa) and GLUT to compile the program.
To compile, do something such as:

cc -o circlefit circlefit.c -lglut -lGLU -lGL -lXmu -lXext -lX11 -lm

If you need the OpenGL Utility Toolkit (GLUT), you can get
it from:

http://reality.sgi.com/mjk_asd/glut3/glut3.html

I hope this is interesting.

- Mark

/* Copyright (c) Mark J. Kilgard, 1997. */

/* This program is freely distributable without licensing fees
and is provided without guarantee or warrantee expressed or
implied. This program is -not- in the public domain. */

/* This is a small interactive demo of Dave Eberly's algorithm
that fits a circle boundary to a set of 2D points. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>

/* Some <math.h> files do not define M_PI... */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

typedef struct {
double x, y;

} Point2;

/****************************************************************************
Least squares fit of circle to set of points.
by Dave Eberly (ebe...@cs.unc.edu or ebe...@ndl.com)
ftp://ftp.cs.unc.edu/pub/users/eberly/magic/circfit.c
---------------------------------------------------------------------------
Input:  (x_i,y_i), 1 <= i <= N, where N >= 3 and not all points
are collinear
Output:  circle center (a,b) and radius r

Energy function to be minimized is

E(a,b,r) = sum_{i=1}^N (L_i-r)^2

where L_i = |(x_i-a,y_i-b)|, the length of the specified vector.
Taking partial derivatives and setting equal to zero yield the
three nonlinear equations

E_r = 0:  r = Average(L_i)
E_a = 0:  a = Average(x_i) + r * Average(dL_i/da)
E_b = 0:  b = Average(y_i) + r * Average(dL_i/db)

Replacing r in the last two equations yields

a = Average(x_i) + Average(L_i) * Average(dL_i/da) = F(a,b)
b = Average(y_i) + Average(L_i) * Average(dL_i/db) = G(a,b)

which can possibly be solved by fixed point iteration as

a_{n+1} = F(a_n,b_n),  b_{n+a} = G(a_n,b_n)

with initial guess a_0 = Average(x_i) and b_0 = Average(y_i).
Derivative calculations show that

dL_i/da = (a-x_i)/L_i,  dL_i/db = (b-y_i)/L_i.

---------------------------------------------------------------------------
WARNING.  I have not analyzed the convergence properties of the fixed
point iteration scheme.  In a few experiments it seems to converge
just fine, but I do not guarantee convergence in all cases.
****************************************************************************/

int
CircleFit(int N, Point2 * P, double *pa, double *pb, double *pr)
{
/* user-selected parameters */
const int maxIterations = 256;
const double tolerance = 1e-06;

double a, b, r;

/* compute the average of the data points */
int i, j;
double xAvr = 0.0;
double yAvr = 0.0;

for (i = 0; i < N; i++) {
xAvr += P[i].x;
yAvr += P[i].y;
}
xAvr /= N;
yAvr /= N;

/* initial guess */
a = xAvr;
b = yAvr;

for (j = 0; j < maxIterations; j++) {
/* update the iterates */
double a0 = a;
double b0 = b;

/* compute average L, dL/da, dL/db */
double LAvr = 0.0;
double LaAvr = 0.0;
double LbAvr = 0.0;

for (i = 0; i < N; i++) {
double dx = P[i].x - a;
double dy = P[i].y - b;
double L = sqrt(dx * dx + dy * dy);
if (fabs(L) > tolerance) {
LAvr += L;
LaAvr -= dx / L;
LbAvr -= dy / L;
}
}
LAvr /= N;
LaAvr /= N;
LbAvr /= N;

a = xAvr + LAvr * LaAvr;
b = yAvr + LAvr * LbAvr;
r = LAvr;

if (fabs(a - a0) <= tolerance && fabs(b - b0) <= tolerance)
break;
}

*pa = a;
*pb = b;
*pr = r;

return (j < maxIterations ? j : -1);

}

enum {
M_SHOW_CIRCLE, M_CIRCLE_INFO, M_RESET_POINTS, M_QUIT

};

#define MAX_POINTS 100

int num = 0;
Point2 list[MAX_POINTS];
int circleFitNeedsRecalc = 0;
int showCircle = 1;
int circleInfo = 0;
int windowHeight;
double a, b, r = 0.0;   /* X, Y, and radius of best fit circle.
*/

void
drawCircle(float x, float y, float r)
{
double angle;

glPushMatrix();
glTranslatef(x, y, 0);
glBegin(GL_TRIANGLE_FAN);
glVertex2f(0, 0);
for (angle = 0.0; angle <= 2 * M_PI; angle += M_PI / 24) {
glVertex2f(r * cos(angle), r * sin(angle));
}
glEnd();
glPopMatrix();

}

void
display(void)
{
int i;

if (circleFitNeedsRecalc) {
int rc;

rc = CircleFit(num, list, &a, &b, &r);
if (rc == -1) {
fprintf(stderr, "circlefit: Problem fitting points to a circle
encountered.\n");
} else {
if (circleInfo) {
printf("%g @ (%g,%g)\n", r, a, b);
}
}
circleFitNeedsRecalc = 0;
}
glClear(GL_COLOR_BUFFER_BIT);

if (showCircle && r > 0.0) {
glColor3ub(0xbb, 0xbb, 0xdd);
drawCircle(a, b, r);
}
glColor3ub(0, 100, 0);
glBegin(GL_POINTS);
for (i = 0; i < num; i++) {
glVertex2d(list[i].x, list[i].y);
}
glEnd();
glutSwapBuffers();

}

void
reshape(int w, int h)
{
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
gluOrtho2D(0, w, 0, h);
glMatrixMode(GL_MODELVIEW);
windowHeight = h;

}

void
{
if (num + 1 >= MAX_POINTS) {
fprintf(stderr, "circlefit: limited to only %d points\n", MAX_POINTS);
return;
}
list[num].x = x;
list[num].y = y;
num++;
circleFitNeedsRecalc = 1;
glutPostRedisplay();

}

void
mouse(int button, int state, int x, int y)
{
if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
}

}

void
{
switch (value) {
case M_SHOW_CIRCLE:
showCircle = !showCircle;
break;
case M_CIRCLE_INFO:
circleInfo = !circleInfo;
break;
case M_RESET_POINTS:
num = 0;
r = 0.0;
break;
case M_QUIT:
exit(0);
}
glutPostRedisplay();

}

int
main(int argc, char **argv)
{
glutInitWindowSize(400, 400);
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
glutCreateWindow("Least squares fit of circle to set of points");

printf("\n");
printf("Least squares fit of circle to set of points\n");
printf("--------------------------------------------\n");
printf("Click left mouse button to position points.  The\n");
printf("program then shows the circle whose boundary best\n");
printf("fits the set of points specified.  Try clicking\n");
printf("points in a near circle.\n");
printf("\n");

glClearColor(125.0 / 256.0, 158.0 / 256.0, 192.0 / 256.0, 1.0);
glPointSize(3.0);
glutReshapeFunc(reshape);
glutMouseFunc(mouse);
glutDisplayFunc(display);