Определение центра окружности или эллипса по его точкам

Всем привет.

Есть следующая задача.

Имеется несколько точек на плоскости, как правило от 6 до 8. Я точно знаю, что они лежат на некотором эллипсе или окружности (с некоторой погрешностью). Мне нужны координаты центра окружности/эллипса. Эта задача уже 100500 раз решена и мне совсем не хочется лепить очередной велосипед. Гугление выдало массу обсуждений задачи и к сожалению ни одной строчки кода.

Тыкните плиз пальцем, где я могу скачать код, который решает эту задачу, например методом наименьших квадратов.
Заранее премного благодарен.
  • Вопрос задан
  • 17453 просмотра
Пригласить эксперта
Ответы на вопрос 4
@mayorovp
Так в чем же проблема?

Уравнение эллипса известно — Ax2+ Bxy + Cy2 + Dx + Ey = 1.
Если выписать его для каждой точки, рассматривая x и y как известные величины, а A,B,C,D,E — как неизвестные, то получится обычная СЛАУ с пятью неизвестными и от 6 до 8 уравнениями.

Если бы координаты точек были заданы точно, это этого было бы достаточно для точного решения. Поскольку точность ограничена — надо найти больше точек и использовать МНК.

Если вы не знаете, как применять МНК к СЛАУ, то проще всего его запомнить в матричной форме:
есть уравнение вида M t = r
его решение методом МНК — это (t = MT M)-1 MT r
осталось реализовать умножение и обращение матриц.

Осталось определить центр эллипса. Для этого запишем уравнение в виде
A(x-x0)2 + B(x-x0)(y-y0) + C(y-y0)2 = F
и заметим, что
D = -(2Ax0 + By0),
E = -(2Cy0 + Bx0),
F — 1 = (Ax02 + Bx0y0 + Cy02
)

Из первых двух уравнений получаем центр, про третье — забываем.

Сложнее обстоит дело с погрешностью — точных формул я не помню, да и что считать несколько непонятно. Возьмите 25 точек, сделайте из них от 10 до 1000 случайных выборок по 12 точек, и решите задачу для этих выборок, после чего можно найти дисперсию распределения центров.
Ответ написан
Комментировать
@megalol
Отстойные, но довольно быстрые, способы — через решение систем уравнений влоб. Типа, окружность по трем точкам, эллипс по восьми (?) точкам, перебираем комбинации все/случайным образом и усредняем результат. То, что предлагает mayorovp вторым пунктом. Почему отстойные — устойчивость плохая. Именно так мы делаем у себя — способы действительно отстойные и требуют костылей, но действительно быстрые. Разница между реальным центром и найденным таким способом намного меньше пикселя, но костылей действительно много. Например, отбирая треугольник при поиске окружности по 3-м точкам, нужно, чтобы все три точки были из разных четвертей окружности.
Я не рекомендую — кода много, а экономия оправдывает себя только на железе 20 летней давности, когда нужен околорилтайм. Для эллипса по-моему малореально вообще что-либо нормальное получить.

Универсальный способ — метод наименьших квадратов. То есть заставить компьютер подобрать такие параметры эллипса чтобы разница между эллипсом и нашими точками была минимальной. Лучше критерий сложно придумать — хотя иногда получается не то, что ожидаешь, но это скорее говорит о том, что исходные данные неправильные.

Алгоритм такой:
1. Выбираем подходящее уравнение эллипса.
2. Пишeм целевую функцию (сумма квадратов уравнения для каждой точки исходных данных). По необходимости добавляем к целевой функции дополнительные условия, типа невырождаемости эллипса в гиперболу. Типа, если условие нарушается, устремляем ее в бесконечность.
3. Выбираем минимизатор (градиентный спуск, нормальное уравнение т. д.)
4. Если результаты не устраивают по скорости или качеству, возвращаемся к п. 3 или к п. 1.

По поводу третьего пункта. В России почему-то считается, что МНК — это когда решение идет с помощью нормального уравнения (составляем матрицы и решаем (AtA)^-1*At*x=b без итераций в один присест). Это не так. Минимизировать можно каким угодно способом. Хоть генетическими алгоритмами. Нормальное уравнение хорошо тогда, когда член (AtA)^-1 получается маленькой квадратной матрицей. То есть когда параметров мало, а точек очень много. Систему уравнений для эллипса вообще к виду A*x-b=0 привести служно. Искать соответствующие статьи как это делается мне влом, поэтому я пойду влоб, в качестве минимизатора взяв стандартный матлабовский.

Используя формулы mayorovp.


%точки
xs = [1 2 3 4 5 6 7 8 9 10];
ys = [11 10 9 7 5 4 3 2 1 1];

%сначала приблизим окружностью - чтобы не вырождалось в гиперболу
%уравнение окружности
circ = @(xs, ys, x0, y0, Rsq) (xs - x0).^2 + (ys - y0).^2 - Rsq;
%целевая функция МНК
fitfunc1 = @(X) sum( ( circ(xs, ys, X(1), X(2), X(3)) ) .^2 );

%начальное приближение
x0 = mean(xs);
y0 = mean(ys);
R = 1;
result1 = fminsearch(fitfunc1, [x0 y0 R], optimset('display', 'iter'));


%уравнение эллипса
ellipse = @(xs, ys, x0, y0, A, B, C, F) A*(xs-x0).^2 + B*(xs-x0).*(ys-y0) + C*(ys-y0).^2 - F;
%целевая функция МНК
fitfunc2 = @(X) sum( ( ellipse(xs, ys, X(1), X(2), X(3), X(4), X(5), X(6))) .^2 );

%начальное приближение - найденная окружность
x0 = result1(1);
y0 = result1(2);
A = 1;
B = 0;
C = 1;
F = result1(3);
 
%сумма квадратов целевой функции для каждой точки
result = fminsearch(fitfunc2, [x0 y0 A B C F], optimset('display', 'iter'));
 
figure
hold on;
%изначальные точки
scatter(xs, ys,'b');
%эллипс
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');


Ответ написан
VasG
@VasG
Ого, приятно, что мою статейку в комментариях вспомнили.
Защищая себя, скажу, что методом МНК ну совсем просто всё делается, поэтому я счёл «Рассуждения есть, кода нет» достаточными.
По ссылке ниже можно скачать образец моего класса примерно годичной давности. Вы уж простите, даже в старой версии уже не совсем МНК :) Для оптимизации скорости я немного забрался вглубь проблемы, так что теперь всё работает быстро. Для работы нужна библиотека Jama. (исп. Jama.matrix, можно прям с ней собрать проект, чтобы остальную часть Jama не таскать).
Тынц.
Ответ написан
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы