Here is the deal. I have multiple points (X,Y) that form an 'ellipse like' shape.
I would like to evaluate/fit the 'best' ellipse possible and get its properties (a,b,F1,F2), or just the center of the ellipse.
Any ideas/leads would be appreciated.
Gilad.
There's a Matlab function fit_ellipse that can do the job. There's also this paper on methods for orthogonal distance fitting of ellipses. A web search for orthogonal ellipse fit will probably turn up a lot of other resources as well.
The ellipse fitting method proposed by:
Z. L. Szpak, W. Chojnacki, and A. van den Hengel. Guaranteed ellipse fitting with a confidence region and an uncertainty measure for centre, axes, and orientation. J. Math. Imaging Vision, 2015.
may be of interest to you. They provide estimates of both algebraic and geometric ellipse parameters, together with covariance matrices that express the uncertainty of the parameter estimates. They also provide a means of computing a planar 95% confidence region associated with the estimate that allows one to visualise the uncertainty in the ellipse fit.
A pre-print version of the paper is available on the authors websites (http://cs.adelaide.edu.au/~wojtek/publicationsWC.html). A MATLAB implementation of the method is also available for download: https://sites.google.com/site/szpakz/source-code/guaranteed-ellipse-fitting-with-a-confidence-region-and-an-uncertainty-measure-for-centre-axes-and-orientation
I will explain how I would approach the problem. I would suggest a hill climbing approach. First compute the gravity center of the points as a start point and choose two values for a and b in some way(probably arbitrary positive values will do). You need to have a fit function and I would suggest it to return the number of points (close enough to)lying on a given ellipse:
int fit(x, y, a, b)
int res := 0
for point in points
if point_almost_on_ellipse(x, y, a, b, point)
res = res + 1
end_if
end_for
return res
Now start with some step. I would choose a big enough value to be sure the best center of the elipse will never be more then step away from the first point. Choosing such a big value is not necessary, but the slowest part of the algorithm is the time it takes to get close to the best center so bigger value is better, I think.
So now we have some initial point(x, y), some initial values of a and b and an initial step. The algorithm iteratively chooses the best of the neighbours of the current point if there is any neighbour better then it, or decrease step twice otherwise. Here by 'best' I mean using the fit function. And also a position is defined by four values (x, y, a, b) and it's neighbours are 8: (x+-step, y, a, b),(x, y+-step, a, b), (x, y, a+-step, b), (x, y, a, b+-step)(if results are not good enough you can add more neighbours by also going by diagonal - for instance (x+-step, y+-step, a, b) and so on). Here is how you do that
neighbours = [[-1, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0], [0, 1, 0, 0],
[0, 0, -1, 0], [0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 1]]
iterate (cx, cy, ca, cb, step)
current_fit = fit(cx, cy, ca, cb)
best_neighbour = []
best_fit = current_fit
for neighbour in neighbours
tx = cx + neighbour[0]*step
ty = cx + neighbour[1]*step
ta = ca + neighbour[2]*step
tb = cb + neighbour[3]*step
tfit = fit(tx, ty, ta, tb)
if (tfit > best_fit)
best_fit = tfit
best_neighbour = [tx,ty,ta,tb]
endif
end_for
if best_neighbour.size == 4
cx := best_neighbour[0]
cy := best_neighbour[1]
ca := best_neighbour[2]
cb := best_neighbour[3]
else
step = step * 0.5
end_if
And you continue iterating until the value of step is smaller then a given threshold(for instance 1e-6). I have written everything in pseudo code as I am not sure which language do you want to use.
It is not guaranteed that the answer found this way will be optimal but I am pretty sure it will be good enough approximation.
Here is an article about hill climbing.
I think that Wild Magic library contains a function for ellipse fitting. There is article with method decription
The problem is to define "best". What is best in your case? The ellipse with the smallest area which contains n% of pointS?
If you define "best" in terms of probability, you can simply use the covariance matrix of your points, and compute the error ellipse.
An error ellipse for this "multivariate Gaussian distribution" would then contain the points corresponding to whatever confidence interval you decide.
Many computing packages can compute the covariance, with its corresponding eigenvalues and eigenvectors. The angle of the ellipse is the angle between the x axis and the eigenvector corresponding to the largest eigenvalue. The semi-axes are the reciprocal of the eigenvalues.
If your routine returns everything normalized (which it should), then you can decide by what factor to multiply everything to obtain an alpha-confidence interval.
