Аппроксимация эллипсом , Фрагмент практической задачи

[Tur]

На плоскости даны N (~100-1000) точек, расположенные на дуге эллипса (известно что они расположены именно на эллипсе)

Требуется аппроксимировать эти точки эллипсом, т.е. найти координаты его фокусов и а, т.е. нужна функция [F1 F2 a] = approx_ellipse(x, y), где x и y - исходные массивы координат точек, F1, F2 - фокусы, и а - параметр эллипса. Всего пять чисел. Точки расположены равномерно (они на фото-матрице с постоянным расстоянием между пикселями)

Я сделал это через симплекс с нулевым приближением, но работает очень медленно(3 секунды для 100 точек).

Как решить эту задачу без нулевого приближения и чтобы работало на обычном среднем компьютере не более 0.1 сек?
Существуют ли какие-нибудь методы прямой аппроксимации эллипсом наподобие аппроксимации полиномами?

[shadeofgray]

Это нелинейная задача аппроксимации, может быть решена Левенбергом-Марквардтом, через минимизацию функции ошибки. Ошибка - это расстояние от экспериментальной точки до эллипса, причем расстояние берется не по оси Y (или X), как в обычных задачах на МНК, а в смысле кратчайшего пути от экспериментальной точки до ближайшей к ней точки предполагаемого эллипса.

Крайне желательно иметь хорошее стартовое приближение (попробуйте придумать какой-нибудь трюк для подбора параметров). Ещё мне кажется, что в случае, если эллипс имеет сильный эксцентриситет, то могут быть сложности со сходимостью.

А точки лежат точно на эллипсе (с поправкой на погрешность матрицы), или могут быть более крупные отклонения?

[MBo]

Этот вопрос встречался в группе
http://groups.google.ru/group/comp.graphics.algorithms/topics?hl=ru&lnk=gschg
однако у меня что-то не в порядке с поиском там.
Скорее всего, искать нужно по запросу "ellipse approximation", и отвечал Eberly
его страница www.geomеtrictools.com
Однако велика вероятность того, что идеи будут те же, что и у shadeofgray - общее уравнение квадрики, и задача минимизации

[Romtek]

По английски это называется Curve fitting.
1.5 Fitting an ellipse by geometric fit
Further Five Point Fit Ellipse Fitting

И здесь ещё набор материалов для прочтения:
Ellipse Detection, Ellipse Fitting, Elliptical Shapes

[Dethlord]

Глянь тут
http://www.geometrictools.com/LibFoundation/Approximation/Approximation.html

[Tur]

Спасибо за ссылки.

Цитата

А точки лежат точно на эллипсе (с поправкой на погрешность матрицы), или могут быть более крупные отклонения?

Пока ничего не могу сказать. Я пока на стадии симуляции. Я понимаю что шумы в данном случае могут сильно влиять на конечный результат. Об этом потом.
Думаю для начала попробовать простой алгоритм. Уравнение эллипса x^2 + axy + by^2 + cx + dy + e = 0. Разделяем дугу на пять равных по числу точек отрезков. Берем первые пять точек начала каждого отрезка, подставляем их в ур-е эллипса, получаем систему пяти уравнений, находим a,b,c,d,e, по ним находим координаты фокусов и главную полуось. Повторяем это для следующих пяти точек и т.д. до конца отрезков. Получаем набор точек фокусов и полуосей. Усредняем. Добавляем шумы и проверяем в какой мере они влияют на результат. Для каждых десяти чисел (5 точек на контуре) находим 5 чисел (F1, F2, a). Недостаток этого метода в том, что у малой выборки большая дисперсия средних. Но работать этот алгоритм будет быстро я думаю. Попробую сначала это, а там видно будет. Что скажите?
Вот нашел в сети http://www.duskyrobin.com/tpu/2006-07-00004.pdf Пока мало что там понял. Центр эллипса там известен заранее. По первому методу его можно найти и дисперсию его проверить. Моя проблема это мой низкий уровень математики.

[Tur]

Вот решение задачи согласно этой теории. Дальнейшее это комментарий к программе на матлабе, приведенной ниже. Даны семь точек на плоскости расположенных вблизи дуги эллипса x, y.

Подставим в уравнение эллипса

эти значения хi, yi ( i = 1…7). Получим систему семи уравнений с пятью неизвестными (а11 – а23). Найдем сначала коэффициенты этой системы
data = [x2 y2 xy x y one];

Для удобства положим а33 = 1. Решим эту систему как указано например здесь. Решение системы Х = [0.02171 0.03839 -0.03186 -0.14497 -0.22382];

Найдем а11-а33 a11 = X(1); a22 = X(2); a12 = X(3)/2; a13 = X(4)/2; a23 = X(5)/2; a33 = 1;
Составим матрицу М = [a11 a12 a13; a12 a22 a23; a13 a23 a33];

и найдем инварианты




( справа от знака % в матлабе следует комментарий )
d = det(M); % d = -0.000152261
D = a11*a22 - a12*a12; % D = 0.0005797514
I = a11 + a22; % I = 0.0601

Находим центр эллипса

Mo.X = (a12*a23 - a13*a22)/D; % 7.8745
Mo.Y = (a13*a12 - a11*a23)/D; % 6.1822

Угол между положительным направлением оси Ox и каждым из двух главных направлений определяется формулой

tg = 2*a12/(a11-a22); % tg = 1.9099 здесь tg – имя переменной
Найдем угол альфа, выражая его тангенс из формулы двойного угла и решая соответствующее уравнение

k = (sqrt(tg*tg+1)-1)/tg; % k = 0.6052 – это тангенс угла наклона большей оси % эллипса к положительному направлению оси Х
Найдем сам угол наклона alpha = atan(k); % alpha = 0.5442 радиан

Чтобы найти оси эллипса найдем корни уравнения
λ2 − Iλ + D = 0
[issol lamda1 lamda2] = quadratic_equation(1, -I, D);

a = sqrt(-d/(lamda1*lamda2*lamda2)); % а = 4.6645 большая полуось
b = sqrt(-d/(lamda1*lamda1*lamda2)); % b = 2.3384 малая полуось
c = sqrt(abs(a*a-b*b)); % с = 4.0361 фокусное расстояние

Для того чтобы всегда правильно найти координаты фокусов учтем
if a < b
k = -1/k;
end

p = c/sqrt(1+k*k);
F1.X = Mo.X - p; F1.Y = Mo.Y - k*p; % F1.X = 4.4215 F1.Y = 4.0925
F2.X = Mo.X + p; F2.Y = Mo.Y + k*p; % F2.X = 11.3274 F2.Y = 8.2719

Складываем все найденные параметры в одну выходную структуру ell
ell.Mo = Mo; ell.F1 = F1; ell.F2 = F2;
ell.alpha = alpha; ell.alphagrad = 180/pi*alpha;
ell.a = a; ell.b = b; ell.c = c; ell.k = k; ell.M = M;

Вот эти функции

Цитата

function ell = detect_ellipse(x, y)
x2 = x.*x; y2 = y.*y; xy = x.*y; one = ones(length(x),1);
data = [x2 y2 xy x y one];
X = slau_mnk(data);
a11 = X(1); a22 = X(2); a12 = X(3)/2; a13 = X(4)/2; a23 = X(5)/2; a33 = 1;
M = [a11 a12 a13; a12 a22 a23; a13 a23 a33];
d = det(M);
D = a11*a22 - a12*a12;
I = a11 + a22;

Mo.X = (a12*a23 - a13*a22)/D;
Mo.Y = (a13*a12 - a11*a23)/D;
tg = 2*a12/(a11-a22);
k = (sqrt(tg*tg+1)-1)/tg;
alpha = atan(k);

[issol lamda1 lamda2] = quadratic_equation(1, -I, D);

a = sqrt(-d/(lamda1*lamda2*lamda2));
b = sqrt(-d/(lamda1*lamda1*lamda2));
c = sqrt(abs(a*a-b*b));

if a < b k = -1/k; end
p = c/sqrt(1+k*k);
F1.X = Mo.X - p; F1.Y = Mo.Y - k*p;
F2.X = Mo.X + p; F2.Y = Mo.Y + k*p;

ell.Mo = Mo; ell.F1 = F1; ell.F2 = F2;
ell.alpha = alpha; ell.alphagrad = 180/pi*alpha;
ell.a = a; ell.b = b; ell.c = c; ell.k = k; ell.M = M;
end

function x = slau_mnk(M)
[rows cols] = size(M); Z = [];
if cols-1 < rows
for col = 1:cols-1
N = [];
for i = 1:rows
N = [N; M(i,col)*M(i, : )];
end
Z = [Z; sum(N)];
end
end
b = -Z(:,end);
A = Z(1:end,1:end-1);
x = A\b; % решение приведенной системы средствами матлаба
end

function [issol x1 x2] = quadratic_equation(A, B, C)
D = B*B - 4*A*C; if D < 0 issol = 0; x1 = 0; x2 = 0; return; else issol = 1; end
x1 = (-B + sqrt(D))/(2*A); x2 = (-B - sqrt(D))/(2*A);
end

Этот код не всегда работает стабильно. Если angle меньше 180 градусов и отношение большой оси к малой более двух, то этот код может давать ложный эллипс (пример в конце поста).


Найдем приближенно и быстро расстояния от исходных точек до найденного эллипса. Сначала я думал соединить точку M (от которой ищется расстояние до эллипса) с фокусами эллипса и провести биссектрису этого угла F1MF2, найти точки пересечения этой биссектрисы с эллипсом и из них найти ближайшую точку к данной, и расстояние между ними считать расстоянием от данной точки до эллипса. Но затем взял не биссектрису, а нормаль. Уравнение нормали к эллипсу в точке (х1, y1)

Так как исходные точки близки к найденному эллипсу, то подставим вместо х1 y1 координаты М, получим нормаль, которая сечет найденный эллипс в двух точках. Найдем ближайшую из этих двух точек к М. Это и будем считать искомым расстоянием. Ошибки этих расстояний от истинных отклонений примерно на порядок-два меньше чем сами отклонения. Эти ошибки можно вычислить точно. Если точка вне эллипса, то ее расстояние до эллипса будем считать положительным, если точка внутри эллипса, то отрицательным. Вот расстояния от семи исходных точек до найденного эллипса.

Вот пример удачной аппроксимации. Это контур

Это аппроксимация каждой стороны контура своим эллипсом. Черные линии – это линии соединяющие фокусы. (Слева чуть виден другой график, не обращайте на него внимания)

Вот пример неудачной аппроксимации.

Почему это происходит и можно ли этого избежать мне до конца неясно. И если этого не избежать, то как приближенно, но быстро и устойчиво определить фокусы?
Странно, что в сети я не нашел подобного простого примера.

[Mikle]

Tur

Цитата

(~100-1000) точек, расположенные на дуге эллипса (известно что они расположены именно на эллипсе)

Раз известно, что ВСЕ они расположены именно на эллипсе, оставь из них четыре, остальные - избыточная информация.

Сообщение отредактировано: Mikle - 27.12.09, 08:53

[Pavia]

Mikle
Обычно точки если и лежат на чем-то, то имеют ошибку. Хотя бы потому что числа с плавающей точкой имеют ограниченную точность. А также обычное дело ошибка метода измерения.
Поэтому 4 точки брать нельзя.

Tur
Искать надо уметь. Проблема в том что есть еще платные вещи.
http://www.robots.ox.ac.uk/~awf/papers/ellipse-pami.pdf

[amk]

Mikle, через четыре точки можно провести более одного эллипса. Нужна пятая точка.

[05772]

Для чего так сложно? В уравнение для эллипса коеффициенты a,b,c,d,e входят линейно. Следовательно они находятся методом наименьших квадратов по стандартному алгоритму.

[amk]

05772, тоже хотел об этом написать, в прошлый раз забыл.
В каноническом уравнении
Ax² + Bxy + Cy² + Dx + Ey + F = 0
получается шесть неизвестных, но поскольку система получается однородная, необходимо пять точек. При большем МНК.

Правда в результате решения может получиться и не эллипс.

Все почему-то прицепились к поиску центра и угла поворота главных осей, вот и получилась нелинейная оптимизация.