4

Is there a straightforward algorithm for fitting data to an ellipse or other conic section? The data generally only approximately fits a portion of the ellipse. I am looking for something that doesn't involve a complicated iterative search, since this has to run at interactive speeds for data sets on the order of 100s of data points. I can downsample or cluster the data if it goes beyond that.

I found an article: "Least-Squares Fitting of Circles and Ellipses" (Gander et al 1994) and it does seem to address my needs, but it uses a lot of mathematical machinery that I either don't understand or have a library for. I'm sure I can grok it given time, but this isn't something I would like to allocate a week to doing.

J. M. ain't a mathematician
  • 75,051
Gabe Johnson
  • 143
  • Least-Squares fitting is exactly what you need. If the ellipse has axes parallel to the coordinate axes, there are only 4 parameters to fit. – lhf Mar 01 '11 at 20:28

2 Answers2

3

It baffles me a bit that people rush to iterative methods, even though in this case there are very straightforward direct methods. A conic section will implicitly be described by a degree $2$ algebraic equation of the form

$$Ax^2+Bxy+Cy^2+Dx+Ey+F=0$$

so if you are given a few points $(x_1,y_1),(x_2,y_2),...$ for each point you get an equation

$$Ax_1^2+Bx_1y_1+Cy_1^2+Dx_1+Ey_1+F=0$$

$$Ax_2^2+Bx_2y_2+Cy_2^2+Dx_2+Ey_2+F=0$$

they are $\textit{linear}$ in the coefficients $A,B,C,...$ of the conic section. So given enough non-degenerate points (technical term "in general poisiton") this homogenous system will have a unique solution. 'Enough' here means five since we will only know $A,B,C,...$ up to a scale factor. If you have more than five points, the system of equations will become overdetermined.

As a joke i once implemented this as a java applet - and it's still up. Maybe this will convince you that this method runs at 'interactive speed' (for five points at least...).

Of course the hard bit then becomes to work with the conic section in the form of the coefficients $A,B,C,...$ - so maybe you have to convert it to the usual parametric form. For this i refer you to for example the ellipse page of math world (equation 15-22). But wikipedia also has more on this.

  • Since in fact, it takes five points to determine a conic, only five of the six coefficients are needed (i.e. you can rescale such that one of those coefficients is 1). The problem is when the data is very noisy, and thus the usual assumptions of least squares (e.g. normally distributed errors) no longer apply. Of course, one can use the results of your suggestion as a starting point for iterative methods. – J. M. ain't a mathematician Jul 17 '11 at 16:48
  • Secondly: if even though you have many points, those points cover only a small arc of the conic, you may well end up with an ill-conditioned problem (then again, if you're using SVD to solve the overdetermined system, well...). – J. M. ain't a mathematician Jul 17 '11 at 16:51
  • This is an elegant solution (and a nice applet). However, it's discussed in the paper linked to in the question; the image on p. 65 shows that the solution may differ quite dramatically from the one that minimizes the sum of squared distances. – joriki Jul 17 '11 at 17:22
  • If five points is all you have, neither LP nor least squares is needed; in fact there is a determinantal expression for a conic through five points. On the note of optimizing in the Manhattan (and Chebyshev as well) norm, there are especially adapted LP routines in the literature, though we might be getting too far here. – J. M. ain't a mathematician Jul 17 '11 at 18:04
  • 1
    I wouldn't say Golub "invented" SVD proper (that's to Beltrami's credit); what Gene did (along with Velvel Kahan and Christian Reinsch) is to restructure the usual tridiagonal QR algorithm for computing eigenvalues of a symmetric tridiagonal matrix to compute the singular values of a bidiagonal matrix. He made a practical algorithm out of a theoretical construct, if I have to be concise... – Jerry Jul 18 '11 at 06:12
0

[Note: As lhf quite rightly pointed out in the comments, this answer is not very useful.]

I think iterative methods should be able to yield results interactively for hundreds of data points, but if you want to optimize for speed and simplicity, it's important which objective function you choose. Two obvious candidates in the circular case are

$$f(x_0,y_0,r) = \sum_i \left(\sqrt{(x_i-x_0)^2+(y_i-y_0)^2} - r\right)^2$$

and

$$g(x_0,y_0,r) = \sum_i \left((x_i-x_0)^2+(y_i-y_0)^2 - r^2\right)^2\;.$$

One might argue that the first alternative captures more accurately what you want to minimize, but the second one is computationally much simpler. There's a solution of the optimization conditions for it here. I believe that you can get quite a simple system of equations if you treat the elliptic case in the same vein:

$$g(x_0,y_0,r,\mu,\nu) = \sum_i \left(\mu(x_i-x_0)^2+\nu(y_i-y_0)^2 - r^2\right)^2\;.$$

I haven't worked out all the details, but I believe the resulting equations are linear in $\mu$ and $\nu$ and you can solve them through a simple iteration. If you want to pursue this approach and need further assistance, let me know.

joriki
  • 238,052
  • 1
    You are assuming that the ellipse is aligned with the coordinate axes. Otherwise, you need to use a full quadratic equation, which has 6 parameters. BTW, for your ellipse case, you don't need the $r$: you can take it as $1$. – lhf Mar 01 '11 at 21:37
  • @lhf: You're right on both counts -- in fact you have to take $r=1$, since otherwise the solution would trivially be $\mu=\nu=r=0$. – joriki Mar 01 '11 at 22:10
  • Your first two objective functions seem to be for fitting equilateral hyperbolas, not circles... – J. M. ain't a mathematician May 01 '11 at 01:28
  • @J.M.: Thanks. This answer was bad enough as it was without that typo... – joriki May 01 '11 at 02:45
  • FWIW: You did more or less what Gander proposed, but he added some refinements. I'll summarize Gander's algorithm maybe later for everybody's convenience. – J. M. ain't a mathematician May 01 '11 at 04:25