How to tell if up to round-off error in floating point, a collection of 2-d double precision floating point pairs might lie on some ellipse?

Question

So arbitrary ellipses seem to have two more degrees of freedom than circles, because in addition to a circle's radius and center there is the angle of rotation as well as the scaling of the ratio of the major axis to the minor axis (which circles do not have). So, I would have to check but I believe it is 5 distinct points that are needed to uniquely define an ellipse passing through those points, possibly subject to some restrictions on the points, such as no point is in the convex hull of any other 3. Anyway, let's say that in general you need a certain number of points to define a unique ellipse (probably 4 or 5 points) passing through them. Let's say we have more than that number of points. If we are dealing in IEEE floating point, say 64-bit, subject to round-off error in the mantissa, what is a robust way to determine whether the points MAY lie on a common ellipse, provided that round-off error in the points' manitissas can explain any discrepancy? Ideally I'd like to avoid extended precision arithmetic to get the answer, but if that is required, so be it, as long as its use and added running time can be minimized.

Answer 1

I would fit the best ellipse you can to your point set and then compute the error of points positions to it.

Finding/fitting ellipse:

1. If you have many points dispersed along whole circumference

   compute middle point point C (middle of bounding box) and find closest B and most far A points to it.

   • C - is the center
   • CA is major axis
   • CB is semi major axis

2. if you have just few points or not covering whole circumference

   ellipse equation with use of basis vectors a,b:

   x(t)=x0+ax*cos(t)+bx*sin(t);
   y(t)=y0+ay*cos(t)+by*sin(t);
   

   • a=(ax,ay) is major axis vector
   • b=(bx,by) is semi major axis vector
   • (x0,y0) is center point
   • t=<0,2.0*PI> is the angular parameter
   • (x(t),y(t)) is point on circumference

   

   We know that a,b are perpendicular to each other so in 2D:

   bx= ay*|b|/|a|;
   by=-ax*|b|/|a|;
   

   in that case we need at least 5 points. In case of 3D we need 6 points or somehow exploit cross product. Anyway this system of equations would most probably lead to transcendental equations so I would not even try to solve this algebraically but instead use approximation search see:

   • solve transcendental system
   • Approximation of n points to the curve with the best fit

   so just nest for approximation for x0,y0,ax,ay,|b|/|a| and compute bx,by then compute cumulative error e=(sum of distances of points to ellipse) and that is it.

Computing the distance of point and ellipse can be tricky

I would construct homogenous transform matrix representing ellipse coordinate system:

• a will be the x axis
• b will be the y axis

• (x0,y0) will be the origin

  then if you need to compute smallest distance between point (xi,yi) and ellipse transform (xi,yi) to ellipse coordinate system (multiply by inverse matrix) and let call the result (Xi,Yi) This can be done also by dot products instead of inverse matrix ...

  Then compute its mean angle on circle t=atan2(Yi*|a|/|b|,Xi)=M

  

  then compute point on that angle in real ellipse x(t),y(t) and the distance is then d=|(xi-x(t),yi-y(t))|

Now the error evaluation of yours

after approximation search found best fitted ellipse compute min distance of each point to it and remember/compute min,max,average distances. From that you now can estimate the error properties you need for example if (max<=round_error) ...

[Notes]

• 5 or 6 nested approximation loops can take a lot of time O(log^6(n)) so try to set the solution bounds as close as you can. Do not use too many recursion layers.
• Also you can combine both approaches first find A,B,C from many points and then use that as start solution and use approx search near that initial solution
• In 2D you can also construct (ax,ay),(bx,by) from |a|,|b|,angle

Answer 2

This is a nontrivial task, and numerous methods have been studied in the past. The topic is called "ellipse fitting" or "robust ellipse fitting". In the second case, you expect that a proportion of the points are not aligned on the ellipse (they are called outliers) and you don't want them to influence the fitting.

For a first approach, I would recommend you the RANSAC method:

• choose five points at random; this is indeed what you need to define an ellipse

• form a 6x6 matrix by plugging the coordinates of the five given points plus an unknown point, with rows (x², xy, y²,x, y, 1)

Developing the determinant of this matrix (6 cofactors), you will find the implicit equation of the ellipse,

a x² + b xy + c y² + d x + e y + f = 0.

• Then a goodness-of-fit is given by the sum (or maximum) of distances of the remaining points to the ellipse.

The Euclidean distance to an ellipse is an uneasy problem as it involves a quartic equation. A reasonable approximation is obtained by plugging the point coordinates in the LHS of the implicit equation and dividing the absolute value by the norm of the gradient.

• Repeat the process with several point sets and keep the smallest global residue, this is your robust fit.

• You can now look at individual residues.