Drexel dragonThe Math ForumDonate to the Math Forum

Ask Dr. Math - Questions and Answers from our Archives
_____________________________________________
Associated Topics || Dr. Math Home || Search Dr. Math
_____________________________________________

Computing Ellipse Parameters


Date: 9/11/96 at 6:5:25
From: symah
Subject: Computing Ellipse Parameters

Hello Dr. Math,

I want to find the five parameters of an ellipse (center, axes,
orientation) given five or more points. What is the best and fastest 
method?

I want to find the best fitting ellipse to an array of points (five 
parameters to find, all points are random points of the ellipse). The 
only information I have is the coordinates of these points. Least mean 
squares leads to very complicated equations. Are there other 
interesting methods?

Thanks a lot.

  Nicolas


Date: 9/11/96 at 19:6:7
From: Doctor Tom
Subject: Re: Computing Ellipse Parameters

This isn't an easy problem.

If you have more than 5 points, I simply don't know any technique
other than a least-squares method to get a good fit.

I do know a trick that's very numerically stable for getting the
equation of an ellipse given exactly 5 points.  (I can even supply
C code that does it!)

After you get the equation of the conic, you've got to check that
it is, in fact, an ellipse and not a circle, hyperbola, or some
degenerate case, and then use standard methods to get the center,
axes, and orientation.

The method is based on Pascal's theorem, which states that if you take
any hexagon inscribed in a conic and find the intersections of the
opposite edges, those three intersections must lie on a line!

Since you've got 5 points, call the sixth point (x, y, w) (in
homogeneous coordinates).  Now write down, based on Pascal's theorem,
the equations of the opposite lines, find their intersections, and
then write down the condition that the three intersections meet on
a line.  Multiply this whole mess out, and you'll get an equation
like this:

  Ax^2 + Bxy + Cy^2 + Dxw + Eyw + Fw^2 = 0

Since you're not interested in the projective case with points at
infinity, just set w = 1, and you get the standard equation of a 
conic.  The A, B, C, D, E, and F will just be numbers that depend on 
the coordinates of the 5 fixed points.

To save you a pile of time, I'll include the code I wrote to do
this.  Needless to say, it took a couple of hours to get it right.

***********************************************************************

 /* Code to find the equation of a conic */
 /*               Tom Davis              */
 /*             April 12, 1996           */


#include <stdio.h>
#define FLOAT double

#define cross(a, b, ab) ab[0] = a[1]*b[2] - a[2]*b[1]; \
			ab[1] = a[2]*b[0] - a[0]*b[2]; \
			ab[2] = a[0]*b[1] - a[1]*b[0];

#define abs(x) ((x > 0.0) ? x : (-x))

/* toconic takes five points in homogeneous coordinates, and returns the
 * coefficients of a general conic equation in a, b, c, ..., f:
 * 
 * a*x*x + b*x*y + c*y*y + d*x + e*y + f = 0.
 * 
 * The routine returns 1 on success; 0 otherwise.  (It can fail, for
 * example, if there are duplicate points.
 * 
 * Typically, the points will be finite, in which case the third (w)
 * coordinate for all the input vectors will be 1, although the code
 * deals cleanly with points at infinity.
 * 
 * For example, to find the equation of the conic passing through (5, 0), 
 * (-5, 0), (3, 2), (3, -2), and (-3, 2), set:
 * 
 * p0[0] =  5, p0[1] =  0, p0[2] = 1, 
 * p1[0] = -5, p1[1] =  0, p1[2] = 1, 
 * p2[0] =  3, p2[1] =  2, p2[2] = 1, 
 * p3[0] =  3, p3[1] = -2, p3[2] = 1, 
 * p4[0] = -3, p4[1] =  2, p4[2] = 1.
 *
 * But if you want the equation of the hyperbola that is tangent to the
 * line 2x=y at infinity,  simply make one of the points be the point at
 * infinity along that line, for example:
 *
 * p0[0] = 1, p0[1] = 2, p0[2] = 0.
 */

int toconic(FLOAT p0[3], FLOAT p1[3], FLOAT p2[3], FLOAT p3[3], FLOAT p4[3], 
	FLOAT *a, FLOAT *b, FLOAT *c, FLOAT *d, FLOAT *e, FLOAT *f)
{
    FLOAT L0[3], L1[3], L2[3], L3[3];
    FLOAT A, B, C, Q[3];
    FLOAT a1, a2, b1, b2, c1, c2;
    FLOAT x0, x4, y0, y4, w0, w4;
    FLOAT aa, bb, cc, dd, ee, ff;
    FLOAT y4w0, w4y0, w4w0, y4y0, x4w0, w4x0, x4x0, y4x0, x4y0;
    FLOAT a1a2, a1b2, a1c2, b1a2, b1b2, b1c2, c1a2, c1b2, c1c2;
    
    cross(p0, p1, L0)
    cross(p1, p2, L1)
    cross(p2, p3, L2)
    cross(p3, p4, L3)
    cross(L0, L3, Q)
    A = Q[0]; B = Q[1]; C = Q[2];
    a1 = L1[0]; b1 = L1[1]; c1 = L1[2];
    a2 = L2[0]; b2 = L2[1]; c2 = L2[2];
    x0 = p0[0]; y0 = p0[1]; w0 = p0[2];
    x4 = p4[0]; y4 = p4[1]; w4 = p4[2];
    
    y4w0 = y4*w0;
    w4y0 = w4*y0;
    w4w0 = w4*w0;
    y4y0 = y4*y0;
    x4w0 = x4*w0;
    w4x0 = w4*x0;
    x4x0 = x4*x0;
    y4x0 = y4*x0;
    x4y0 = x4*y0;
    a1a2 = a1*a2;
    a1b2 = a1*b2;
    a1c2 = a1*c2;
    b1a2 = b1*a2;
    b1b2 = b1*b2;
    b1c2 = b1*c2;
    c1a2 = c1*a2;
    c1b2 = c1*b2;
    c1c2 = c1*c2;

    aa = -A*a1a2*y4w0
	 +A*a1a2*w4y0 
	 -B*b1a2*y4w0
	 -B*c1a2*w4w0
	 +B*a1b2*w4y0
	 +B*a1c2*w4w0
	 +C*b1a2*y4y0
	 +C*c1a2*w4y0
	 -C*a1b2*y4y0
	 -C*a1c2*y4w0;

    cc =  A*c1b2*w4w0
         +A*a1b2*x4w0
	 -A*b1c2*w4w0
	 -A*b1a2*w4x0
	 +B*b1b2*x4w0
	 -B*b1b2*w4x0
	 +C*b1c2*x4w0
	 +C*b1a2*x4x0
	 -C*c1b2*w4x0
	 -C*a1b2*x4x0;

    ff =  A*c1a2*y4x0
	 +A*c1b2*y4y0
	 -A*a1c2*x4y0
	 -A*b1c2*y4y0
	 -B*c1a2*x4x0
	 -B*c1b2*x4y0
	 +B*a1c2*x4x0
	 +B*b1c2*y4x0
	 -C*c1c2*x4y0
	 +C*c1c2*y4x0;

    bb =  A*c1a2*w4w0
	 +A*a1a2*x4w0
	 -A*a1b2*y4w0
	 -A*a1c2*w4w0
	 -A*a1a2*w4x0
	 +A*b1a2*w4y0
	 +B*b1a2*x4w0
	 -B*b1b2*y4w0
	 -B*c1b2*w4w0
	 -B*a1b2*w4x0
	 +B*b1b2*w4y0
	 +B*b1c2*w4w0
	 -C*b1c2*y4w0
	 -C*b1a2*x4y0
	 -C*b1a2*y4x0
	 -C*c1a2*w4x0
	 +C*c1b2*w4y0
	 +C*a1b2*x4y0
	 +C*a1b2*y4x0
	 +C*a1c2*x4w0;

    dd = -A*c1a2*y4w0
	 +A*a1a2*y4x0
	 +A*a1b2*y4y0
	 +A*a1c2*w4y0
	 -A*a1a2*x4y0
	 -A*b1a2*y4y0
	 +B*b1a2*y4x0
	 +B*c1a2*w4x0
	 +B*c1a2*x4w0
	 +B*c1b2*w4y0
	 -B*a1b2*x4y0
	 -B*a1c2*w4x0
	 -B*a1c2*x4w0
	 -B*b1c2*y4w0
	 +C*b1c2*y4y0
	 +C*c1c2*w4y0
	 -C*c1a2*x4y0
	 -C*c1b2*y4y0
	 -C*c1c2*y4w0
	 +C*a1c2*y4x0;

    ee = -A*c1a2*w4x0
	 -A*c1b2*y4w0
	 -A*c1b2*w4y0
	 -A*a1b2*x4y0
	 +A*a1c2*x4w0
	 +A*b1c2*y4w0
	 +A*b1c2*w4y0
	 +A*b1a2*y4x0
	 -B*b1a2*x4x0
	 -B*b1b2*x4y0
	 +B*c1b2*x4w0
	 +B*a1b2*x4x0
	 +B*b1b2*y4x0
	 -B*b1c2*w4x0
	 -C*b1c2*x4y0
	 +C*c1c2*x4w0
	 +C*c1a2*x4x0
	 +C*c1b2*y4x0
	 -C*c1c2*w4x0
	 -C*a1c2*x4x0;

    if (aa != 0.0) {
	bb /= aa; cc /= aa; dd /= aa; ee /= aa; ff /= aa; aa = 1.0;
    } else if (bb != 0.0) {
	cc /= bb; dd /= bb; ee /= bb; ff /= bb; bb = 1.0;
    } else if (cc != 0.0) {
	dd /= cc; ee /= cc; ff /= cc; cc = 1.0;
    } else if (dd != 0.0) {
	ee /= dd; ff /= dd; dd = 1.0;
    } else if (ee != 0.0) {
	ff /= ee; ee = 1.0;
    } else {
	return 0;
    }
    *a = aa;
    *b = bb;
    *c = cc;
    *d = dd;
    *e = ee;
    *f = ff;
    return 1;
}

FLOAT p0[3] = {0, 0, 1};
FLOAT p1[3] = {1, 1, 1};
FLOAT p2[3] = {-1, -1, 1};
FLOAT p3[3] = {2, 2, 1};
FLOAT p4[3] = {3, 3, 1};

main()
{
    FLOAT a, b, c, d, e, f, s0, s1, s2, s3, s4;
    FLOAT x, y, w;
    int i, ret;

/*
    ret = toconic(p0, p1, p2, p3, p4, &a, &b, &c, &d, &e, &f);
    if (ret == 1) {
	printf("success\n");
	printf("%g %g %g %g %g %g\n", a, b, c, d, e, f);
	x = p0[0]; y = p0[1]; w = p0[2];
	printf("%g ", a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w);
	x = p1[0]; y = p1[1]; w = p1[2];
	printf("%g ", a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w);
	x = p2[0]; y = p2[1]; w = p2[2];
	printf("%g ", a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w);
	x = p3[0]; y = p3[1]; w = p3[2];
	printf("%g ", a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w);
	x = p4[0]; y = p4[1]; w = p4[2];
	printf("%g\n", a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w);
    } else {
	printf("toconic failed\n");
    }
*/

    for (i = 0; i < 100000; i++) {
	p0[0] = (FLOAT) (rand()%30);
	p0[1] = (FLOAT) (rand()%30);
	p0[2] = (FLOAT) (rand()%30);
	p1[0] = (FLOAT) (rand()%30);
	p1[1] = (FLOAT) (rand()%30);
	p1[2] = (FLOAT) (rand()%30);
	p2[0] = (FLOAT) (rand()%30);
	p2[1] = (FLOAT) (rand()%30);
	p2[2] = (FLOAT) (rand()%30);
	p3[0] = (FLOAT) (rand()%30);
	p3[1] = (FLOAT) (rand()%30);
	p3[2] = (FLOAT) (rand()%30);
	p4[0] = (FLOAT) (rand()%30);
	p4[1] = (FLOAT) (rand()%30);
	p4[2] = (FLOAT) (rand()%30);
	if (toconic(p0, p1, p2, p3, p4, &a, &b, &c, &d, &e, &f)) {
	    x = p0[0]; y = p0[1]; w = p0[2];
	    s0=a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w;
	    x = p1[0]; y = p1[1]; w = p1[2];
	    s1=a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w;
	    x = p2[0]; y = p2[1]; w = p2[2];
	    s2=a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w;
	    x = p3[0]; y = p3[1]; w = p3[2];
	    s3=a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w;
	    x = p4[0]; y = p4[1]; w = p4[2];
	    s4=a*x*x+b*x*y+c*y*y+d*x*w+e*y*w+f*w*w;
	    if (abs(s0) > .00001 || abs(s1) > .00001 || abs(s2) > .00001
		|| abs(s3) > .00001 || abs(s4) > .00001) {
		    printf("%g %g %g %g %g\n", s0, s1, s2, s3, s4);
		    printf("\t%g %g %g\n", p0[0], p0[1], p0[2]);
		    printf("\t%g %g %g\n", p1[0], p1[1], p1[2]);
		    printf("\t%g %g %g\n", p2[0], p2[1], p2[2]);
		    printf("\t%g %g %g\n", p3[0], p3[1], p3[2]);
		    printf("\t%g %g %g\n", p4[0], p4[1], p4[2]);
		    printf("%g %g %g %g %g %g\n", a, b, c, d, e, f);
		}
	}
    }
}


-Doctor Tom,  The Math Forum
 Check out our web site!  http://mathforum.org/dr.math/   
    
Associated Topics:
College Conic Sections/Circles

Search the Dr. Math Library:


Find items containing (put spaces between keywords):
 
Click only once for faster results:

[ Choose "whole words" when searching for a word like age.]

all keywords, in any order at least one, that exact phrase
parts of words whole words

Submit your own question to Dr. Math

[Privacy Policy] [Terms of Use]

_____________________________________
Math Forum Home || Math Library || Quick Reference || Math Forum Search
_____________________________________

Ask Dr. MathTM
© 1994-2013 The Math Forum
http://mathforum.org/dr.math/