Determining the center of a circle or ellipse from its points

Question

Algorithms

Hello.

There is the following task.

There are several points on the plane, usually from 6 to 8. I know for sure that they lie on some ellipse or circle (with some error). I need the coordinates of the center of a circle/ellipse. This problem has already been solved 100,500 times, and I don’t feel like sculpting another bicycle at all. Googling produced a lot of discussions of the problem and, unfortunately, not a single line of code.

Poke pliz with your finger where I can download the code that solves this problem, for example, by the least squares method.
Thank you very much in advance.

Answer 1

So what's the problem?
The ellipse equation is known - Ax ² + Bxy + Cy ² + Dx + Ey = 1 .
If you write it out for each point, considering x and y as known values, and A, B, C, D, E as unknowns, then you get the usual SLAE with five unknowns and 6 to 8 equations.
If the coordinates of the points were given exactly, this would be enough for an exact solution. Since the accuracy is limited, you need to find more points and use the least squares method.
If you don’t know how to apply LSM to SLAE, then the easiest way to remember it is in matrix form:
there is an equation of the form M t = r
, its solution by the LSM method is (t = M ^T M) ^-1 M^T r
it remains to implement the multiplication and inversion of matrices.
It remains to determine the center of the ellipse. To do this, we write the equation in the form
A(xx ₀ ) ² + B(xx ₀ )(yy ₀ ) + C(yy ₀ ) ² = F
and note that
D = -(2Ax ₀ + By ₀ ),
E = -( 2Cy ₀ + Bx ₀ ),
F - 1 = (Ax ₀ ² + Bx ₀ y ₀ + Cy ₀ ² )
From the first two equations we get the center, we forget about the third.
The situation with the error is more complicated - I don’t remember the exact formulas, and what to count is somewhat incomprehensible. Take 25 points, make from 10 to 1000 random samples of 12 points, and solve the problem for these samples, after which you can find the variance of the distribution of centers.

Answer 2

http://habrahabr.ru/post/135332/

Answer 3

Отстойные, но довольно быстрые, способы — через решение систем уравнений влоб. Типа, окружность по трем точкам, эллипс по восьми (?) точкам, перебираем комбинации все/случайным образом и усредняем результат. То, что предлагает mayorovp вторым пунктом. Почему отстойные — устойчивость плохая. Именно так мы делаем у себя — способы действительно отстойные и требуют костылей, но действительно быстрые. Разница между реальным центром и найденным таким способом намного меньше пикселя, но костылей действительно много. Например, отбирая треугольник при поиске окружности по 3-м точкам, нужно, чтобы все три точки были из разных четвертей окружности.
Я не рекомендую — кода много, а экономия оправдывает себя только на железе 20 летней давности, когда нужен околорилтайм. Для эллипса по-моему малореально вообще что-либо нормальное получить.
Универсальный способ — метод наименьших квадратов. То есть заставить компьютер подобрать такие параметры эллипса чтобы разница между эллипсом и нашими точками была минимальной. Лучше критерий сложно придумать — хотя иногда получается не то, что ожидаешь, но это скорее говорит о том, что исходные данные неправильные.
Алгоритм такой:
1. Выбираем подходящее уравнение эллипса.
2. Пишeм целевую функцию (сумма квадратов уравнения для каждой точки исходных данных). По необходимости добавляем к целевой функции дополнительные условия, типа невырождаемости эллипса в гиперболу. Типа, если условие нарушается, устремляем ее в бесконечность.
3. Выбираем минимизатор (градиентный спуск, нормальное уравнение т. д.)
4. Если результаты не устраивают по скорости или качеству, возвращаемся к п. 3 или к п. 1.
Regarding the third point. In Russia, for some reason, it is believed that the least squares is when the decision is made using a normal equation (we compose matrices and solve (AtA)^-1*At*x=b without iterations in one sitting). This is not true. You can minimize it in any way. Even genetic algorithms. The normal equation is good when the term (AtA)^-1 is obtained by a small square matrix. That is, when there are few parameters, and there are a lot of points. The system of equations for an ellipse in general can be brought to the form A*xb=0. Search for relevant articles on how this is done to me, so I'll go ahead, taking the standard Matlab as a minimizer.
Using formulas mayorovp.

% points
xs = [1 2 3 4 5 6 7 8 9 10];
ys = [11 10 9 7 5 4 3 2 1 1];

% first approximate by a circle - so as not to degenerate into a hyperbole
%circle equation
circ = @(xs, ys, x0, y0, Rsq) (xs - x0).^2 + (ys - y0).^2 - Rsq;
%objective function LSM
fitfunc1 = @(X) sum( ( circ(xs, ys, X(1), X(2), X(3)) ) .^2 );

% initial approximation
x0 = mean(xs);
y0 = mean(ys);
R=1;
result1 = fminsearch(fitfunc1, [x0 y0 R], optimset('display', 'iter'));


% ellipse equation
ellipse = @(xs, ys, x0, y0, A, B, C, F) A*(xs-x0).^2 + B*(xs-x0).*(ys-y0) + C*(ys -y0).^2 - F;
%objective function LSM
fitfunc2 = @(X) sum( ( ellipse(xs, ys, X(1), X(2), X(3), X(4), X(5), X(6))) .^2 ) ;

%initial approximation - found circle
x0 = result1(1);
y0 = result1(2);
A=1;
B=0;
C=1;
F = result1(3);
 
%sum of squared objective function for each point
result = fminsearch(fitfunc2, [x0 y0 ABCF], optimset('display', 'iter'));
 
figures
hold on;
%original points
scatter(xs, ys,'b');
%ellipse
xs1 = -2:0.005:20;
ys1 = -2:0.005:20;
[xss yss] = meshgrid(xs1, ys1);
ell = abs(ellipse(xss(:), yss(:), result(1), result(2), result(3), result(4), result(5), result(6)) ) < 0.0001;
scatter(xss(ell), yss(ell),'r');

Answer 4

Wow, nice that my article was remembered in the comments.
Defending myself, I’ll say that by the method of least squares, everything is done quite simply, so I considered “There are arguments, there is no code” sufficient.
You can download a sample of my class from about a year ago at the link below. Forgive me, even in the old version it's not quite MNC :) To optimize the speed, I climbed a little deeper into the problem, so now everything works quickly . To work, you need the Jama library. (Spanish Jama.matrix, you can build a project directly with it so that you don’t carry the rest of Jama).
Tynts .