Насильственная аппроксимация дуги заданным эллипсом

[Tur]

Алексей К. в сообщении #1374850 писал(а):

У меня не грипп , но я тоже вынужден отложить эту развлекуху на некоторое время.

Для меня это не развлекуха, а работа и вы в ней мне здорово помогли, большое бремя сняли с плеч. Спасибо!
Опишу как это работает.
Дана дуга (синие 111 точек) и ее первоначальная аппроксимация эллипсом (на рисунке изображен черным и белым его центр и фокусы). Получившиеся параметры эллипса а = 187 и b = 58

Задача состоит в том чтобы эту дугу аппроксимировать эллипсом с параметрами а = 79 и b = 37.5
Вот код исходной аппроксимации:

(Оффтоп)

function [x y Mo MT] = move_zero(x, y)
ell = old_param_ell(x, y);
drawell_by_aa(ell.aa); 
mm(ell.F1,'w'); mm(ell.F2,'w');
MT = ell.MT; Mo = ell.Mo; mm(Mo,'w')
[x y] = movearr(x, y, MT);
end

function ell = old_param_ell(x, y)
% ell - aa MT fi a b Mo F1 F2
mx = mean(x); my = mean(y);
xn = x - mx; yn = y - my;
[a11 a12 a13 a22 a23] = aquadric(xn, yn);
[D I d] = invariant(a11, a12, a13, a22, a23);
type = type_quadric(D, I, d); % ell.type = 1;
if type == 1 || type == 2
  [a b fi s c] = a_b_fi(a11, a22, a12, d);
  Mo = first_Mo_d(a11, a12, a13, a22, a23);
  [Mo MT] = Mo_ell(Mo, mx, my, s, c);
  [F1 F2] = foci(Mo, s, c, a, b);
else
  if type == -1 return; end
  if type == 3 disp('parabola'); return; end
end
ell.type = type;
% ell.faa = [a11 a12 a13 a22 a23 1];
ell.aa = complete_equation(a, b, Mo, fi);
ell.MT = MT;
ell.ab = [a b];
ell.fi = fi;
ell.Mo = Mo;
ell.F1 = F1;
ell.F2 = F2;
end
function type = type_quadric(D, I, d)
if d == 0
  type = -1; return 
end
if 0 < D && d*I < 0
  type = 1; return % ellipse
end
if D < 0
  type = 2; return % hyperbola
end
if D == 0
  type = 3; return % parabola
end
type = -1;
end
function Mo = first_Mo_d(a11, a12, a13, a22, a23)
D = a11*a22 - a12*a12;
Mo.X = (a12*a23 - a13*a22)/D;
Mo.Y = (a13*a12 - a11*a23)/D;
end
function [a b fi s c] = a_b_fi(a11, a22, a12, d)
tg = 2*a12/(a11-a22); fi = 0.5*atan(tg);
s = sin(fi); c = cos(fi); s2 = s*s; c2 = c*c;
lambda1 = a11*c2 + 2*a12*s*c + a22*s2;
lambda2 = a11*s2 - 2*a12*s*c + a22*c2;
a2 = -d/(lambda1*lambda1*lambda2); 
b2 = -d/(lambda1*lambda2*lambda2);
if a2 < 0 a2 = -a2; end
if b2 < 0 b2 = -b2; end
if a2 < b2 
  tmp = b2; b2 = a2; a2 = tmp; 
  % s = sin(fi+pi/2); c = cos(fi+pi/2);
  fi = fi + pi/2;
  s = sin(fi); c = cos(fi);
end
a = sqrt(a2); b = sqrt(b2);
end
function [Mo MT] = Mo_ell(Mo, mx, my, s, c)
MT1 = [1 0 -mx; 0 1 -my; 0 0 1];
Mper = [1 0 -Mo.X; 0 1 -Mo.Y; 0 0 1];
Mpov = [c s 0; -s c 0; 0 0 1];
MT2 = Mpov*Mper; 
MT = MT2*MT1; invM = inv(MT);
P.X = 0; P.Y = 0; Mo = mp(invM, P);
end
function [F1 F2] = foci(Mo, s, c, a, b)
tg = s/c; p = sqrt(a*a-b*b)/sqrt(1+tg*tg);
F1.X = Mo.X - p; F1.Y = Mo.Y - tg*p;
F2.X = Mo.X + p; F2.Y = Mo.Y + tg*p;
end

function aa = complete_equation(a, b, Mo, fi)
s = sin(fi); c = cos(fi);
s2 = s*s; c2 = c*c;
xo = Mo.X; yo = Mo.Y;
xo2 = xo*xo; yo2 = yo*yo;
a2 = a*a; b2 = b*b;
a11 = c2/a2 + s2/b2;
a22 = s2/a2 + c2/b2;
a12 = s*c*(1/a2 - 1/b2);
a13 = -xo*a11 - yo*a12;
a23 = -xo*a12 - yo*a22;
a33 = a11*xo2 + a22*yo2 + 2*a12*xo*yo - 1;
aa = [a11 a12 a13 a22 a23 a33]';
end

Нулевое приближение это центр найденного эллипса Mo.X = 195 Mo.Y = -105 и угол наклона эллипса fi = 0.7442
При аппроксимации получена также матрица поворота и переноса МТ. С ее помощью преобразуем исходную дугу так чтобы центр найденного эллипса совпал с началом координат, а главная ось эллипса с осью Х, см рис ниже.

Обозначим новую полученную дугу xn, yn. Теперь задача состоит в том чтобы аппроксимировать ее заданным эллипсом.

[Tur]

Эту аппроксимацию осуществим с помощью алгоритма предложенного Алексеем (немного преобразую его)
$L(x,y;u,v,\xi)=\frac{x}{a^2}u+\frac{y}{b^2}v+\left(\frac1{a^2}-\frac1{b^2}\right)xy\,\xi-0.5\cdot(\frac{x^2}{a^2}+ \frac{y^2}{b^2}-1) = 0


Параметры голубого эллипса аппроксимирующего синюю дугу xn, yn
a = 79 b = 37.5 Mo.X = 76 Mo.Y = -4.6 fi = 0.0135
Вот ф-ция, осуществляющая эту аппроксимацию

(Оффтоп)

function [Mo fi MT] = get_vell(a,b,x,y)
a2 = a*a; b2 = b*b; A = x/a2; B = y/b2;
C = (1/a2 - 1/b2)*x.*y;
D = -(x.*x/a2 + y.*y/b2 - 1)/2;
M = [A B C D]; X = pere_slau(M);
Mo.X = -X(1); Mo.Y = -X(2); fi = X(3);

Mper = [1 0 -Mo.X; 0 1 -Mo.Y; 0 0 1];
s = sin(fi); c = cos(fi);
Mpov = [c s 0; -s c 0; 0 0 1];
MT = Mpov*Mper;
end

function X = pere_slau(M)
% решение переопределенной
% системы уравнений
% A1*x + B1*y = C1    M = [A1 B1 -C1;
% A2*x + B2*y = C2         A2 B2 -C2;
% A3*x + B3*y = C3         A3 B3 -C3;
% A4*x + B4*y = C4         A4 B4 -C4];
[~, c] = size(M);
A = M(:,1:c-1);
B = M(:,c);
C = A'*A;
D = A'*B;
X = C\D;
end

Эта функция дает и соответствующую матрицу МТ. С ее помощью можно перенести новый полученный эллипс в координаты исходной дуги

Как видим аппроксимация еще далека от завершения.
Подобно тому как мы перенесли с помощью первичной аппроксимации исходную дугу в координаты где центр найденного эллипса попадал в начало координат - так и сейчас перенесем координаты дуги xn yn с помощью только что полученной МТ и получим (белые точки - обозначим их xn2 yn2):

На следующем шаге будем аппроксимировать уже эти точки

-- Пн фев 11, 2019 15:38:22 --

Ниже приведен зеленый эллипс с параметрами a = 79 b = 37.5 аппроскимирующий белые точки дуги

Новый зеленый эллипс в координатах исходной дуги

Как видим эта аппроксимация уже ближе к истинной.
Сделаем еще одну итерацию: с помощью новой "зеленой" МТ получим новую дугу xn3 yn3 (желтую)

новая аппроксимация желтой дуги (красная линия). Уже по рис видно на сколько ее центр близок с началу координат

Новая, красная итерация в исходных координатах

[Tur]

Нулевым приближением полагались два параметра найденных при первичной аппроксимации дуги произвольным эллипсом: угол наклона его главной оси и его центр Мо1. Я изменил это условие и нашел точку пересечения (Р) оси эллипса с его дугой или с ним самим и затем по оси эллипса отложит отрезок равный а = 79. Т.о. получил новый центр Мо.

В этом случае требуются всего одна-две итерации. Вот первая:

А вот вторая желтая, почти совпадающая с предыдущей


А вот случай когда дуга гипербола (так показывает первичная ее аппроксимация кривой второго порядка)

Находим ось гиперболы, ее пересечение либо с дугой либо с самой гиперболой в ее вершине (Р). От этой точки Р откладываем отрезок а = 79 и находим Мо - начальное приближение эллипса. Затем первая итерация:

и вторая:

Т.о. аппроксимировали гиперболу эллипсом с заданными осями.

[wrest]

Кажись, полезное чтение для ТС:

Grbić, R., Grahovac, D., & Scitovski, R. (2016). A method for solving the multiple ellipses detection problem. Pattern Recognition, 60, 824–834.
https://doi.org/10.1016/j.patcog.2016.06.031

(Abstract)

In this paper, the multiple ellipses detection problem on the basis of a data points set coming from a number of ellipses in the plane not known in advance is considered. An ellipse is considered as a Mahalanobis circle with some positive definite matrix. A very efficient method for solving this problem is proposed. This method very successfully combines the well-known direct least squares method and the RANSAC algorithm with a realistic statistical model of multiple ellipses in the plane. The method is illustrated and tested on numerous synthetic and real-world applications. The method was also compared with other similar methods. In the case when a data points set comes from a number of ellipses with clear edges, the proposed method gives results similar to other known methods. However, when a data points set comes from a number of ellipses with noisy edges, the proposed method performs significantly better than the other methods. We should emphasize the advantage and utility of the proposed methods in a variety of applications such as: medical image analysis and ultrasound image segmentation.

Там они детектируют эллипсы без предварительно заданных параметров (полуосей). И вродь получается хорошо.

Ну и библиография там есть.
Например:
[16] Y. Qiao, S. H. Ong, Arc-based evaluation and detection of ellipses, Pattern Recognition 40 (2007) 1990–2003 ( DOI: 10.1016/j.patcog.2006.10.009 )
картинка оттуда:

[17] C. Akinlar, C. Topal, Edcircles: A real-time circle detector with a false detection control, Pattern Recognition 46 (2013) 725–740. ( https://doi.org/10.1016/j.patcog.2012.09.020 )

Насчет как почитать. Идете на sci hub (напр. sci-hub.se ) и туда вводите doi (всю строку, т.е. напр. https://doi.org/10.1016/j.patcog.2016.06.031 )

[Tur]

wrest, вы эту тему не смотрели. С детектированием эллипса уже много лет нет никаких проблем. Грубо говоря в этой теме речь шла об аппроксимации гиперболы эллипсом. Даны точки, дуга и мы даже не выясняем эти точки эллипса или гиперболы или еще чего-то. Это безразлично. Эти точки теперь можно аппроксимировать эллипсом причем длины его осей мы можем выбирать какие нам заблагорассудится.

[Алексей К.]

Алексей К. в сообщении #1374850 писал(а):

Ежели из Вашего (прежнего) набора из 127 точек взять первые 60, то уже получим гиперболу в качестве первого (5-точечного) приближения.

Результаты возни с этой гиперболой (относительно ценный вывод -- в конце сообщения).

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

Что значит здесь "эллипс -- не_эллипс"?
Я беру 5 "опорных" точек $P_{1,3,5}$ --- первую, среднюю (допустим, по номеру) и последнюю. И две промежуточные $(P_{2,4})$.
По 5 точкам однозначно строю конику и смотрю: эллипс -- не_эллипс?

А ещё пытаюсь построить эллипс с фиксированными полуосями по трём точкам, $P_{1,3,5}$. Таковых несколько (не уверен, что все вытащил).
Промежуточные точки $P_{2,4}$ служат для того, чтобы из этих вариантов выбрать более лучший.

На той картинке, однако, взяты все 127 точек.
А если я ограничусь 60-ю, выберу из них $P_{1,3,5}$, то эллипс с заданными полуосями, проходящий через эти точки, будет не существовать. Потому что нарушено условие

$$\frac{b^2}a<r<\frac{a^2}b,$$

где $r$ -- радиус окружности, проходящей через эти три точки (а штуки слева и справа -- радиусы кривизны в вершинах эллипса). У 60-точечной кривой $b^2/a=17.8,\; r=172.6,\; a^2/b=166.4$). Неча и надеяться на сходимость кого-то к кому-то.

А указанное условие существования либо очевидно, либо почти очевидно (т.е. слегка требует доказательства).

[Алексей К.]

Tur,
правильно ли я понял, что у Вас всё типа получается?

У меня сложилось впечатление, что если бы научиться искать эллипс по трём точкам, то это было бы и надёжное первое приближение, и, возможно, вполне приличный окончательный результат (без всяких там МНКов и итераций).
Я хоть и как-то научился численно его (их) находить, но плохо и ненадёжно.
В смысле, я пока даже не придумал, как это сделать хорошо. А интересно.

Изготавливая из Вашей кривой более короткие дуги (которые, естественно, дают приближения похуже или сильно хуже), я сталкивался с тем, что 5-точечное приближение не срабатывает. А 3-точечное --- надёжно. Но, разумеется, выводы надо делать, обладая Вашей богатой (насколько я понял) статистикой и возможностью проделать тыщу тестов.

Вот что у меня получилось с 70-точечной дугой, маленькой, совсем не захватывающей острую вершину, но тем не менее вполне прилично:



Фуксиевый эллипс --- 5-точечное приближение.
Фуксиевый пунктирный --- куда оно сошлось.
Зелёные эллипсы --- два 3-точечных приближения. Выбран тот, что пожирнее. Так как он поближе к точкам 2 и 4.
Зелёный пунктирный --- куда оно сошлось.

[Tur]

Алексей, спасибо что продолжаете попытки усовершенствовать алгоритм.

Я написал функцию возвращающую центр дуги и тестирую ее. Пока что удается ее завалить.
Суть этой функции в этом уравнении

$L(x,y;u,v,\xi)=\frac{x}{a^2}u+\frac{y}{b^2}v+\left(\frac1{a^2}-\frac1{b^2}\right)xy\,\xi-0.5\cdot(\frac{x^2}{a^2}+ \frac{y^2}{b^2}-1) = 0 обозначим это уравнение - (!!)

вот код, реализующий это уравнение:

a2 = a*a; b2 = b*b; A = x/a2; B = y/b2;
C = (1/a2 - 1/b2)*x.*y;
D = -(x.*x/a2 + y.*y/b2 - 1)/2;
M = [A B C D]; X = pere_slau(M);
Mo.X = -X(1); Mo.Y = -X(2); fi = X(3);

Т.е. аппроксимация ведется (пока) по всем точкам дуги и нахожу ее центр Мо и угол наклона fi.
Затем получаю матрицу поворота и переноса

Mper = [1 0 -Mo.X; 0 1 -Mo.Y; 0 0 1];
s = sin(fi); c = cos(fi);
Mpov = [c s 0; -s c 0; 0 0 1];
MT = Mpov*Mper;

Затем с ее (МТ) помощью переношу все точки дуги и повторяю процедуру нахождения Мо и fi. Это совсем краткое описание (могу это развернуть и подробно по шагам описать)

Главное вывод: любую (с некоторыми пока до конца не ясными ограничениями) дугу можно аппроксимировать эллипсом с произвольно заданными длинами осей. Уравнение (!!) дает такую возможность (должно дать). Полученный результат дает примерно то же самое как и аппроксимация набора точек на плоскости окружностью

(Оффтоп)

function [Mo R] = approx_circle(x, y)
N = length(x); k = 1/N; k2 = 2*k;
sumx = sum(x); sumy = sum(y);
x2 = x.*x; x3 = x.*x2; sumx2 = sum(x2); sumx3 = sum(x3);
y2 = y.*y; y3 = y.*y2; sumy2 = sum(y2); sumy3 = sum(y3);
sumxy = sum(x.*y);
sumx2y2 = sum(x2+y2);
sumxy2 = sum(x.*y2);
sumx2y = sum(x2.*y);

sumX = sumx2y2*sumx;
sumY = sumx2y2*sumy;

a1 = k2*sumx*sumx - 2*sumx2;
b1 = k2*sumx*sumy - 2*sumxy; 
b2 = k2*sumy*sumy - 2*sumy2;
c1 = sumx3 + sumxy2 - k*sumX;
c2 = sumy3 + sumx2y - k*sumY;

xo = (b1*c2 - c1*b2)/(a1*b2 - b1*b1);
yo = (b1*c1 - a1*c2)/(a1*b2 - b1*b1); Mo.X = xo; Mo.Y = yo;
R = sqrt(xo*xo + yo*yo + k*sumx2y2 - k2*(xo*sumx + yo*sumy));
end

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

На самом деле ур-е (!!) применяется не к исходным координатам дуги, а к уже перенесенным за счет нулевого приближения. Это большой недостаток этого метода, т.к. не всегда просто найди это нулевое приближение. Может быть возможно было бы в самое исходное ур-е эллипса подставить скажем не $x_1 = x - x_0$, а $x_1 = x - x_0 - \delta x$ тоже самое с $y_1$ и с $\varphi$, т.е. вместо $\varphi$ подставить $\varphi - $\delta\varphi$ и получим не это:

$(\frac{cos(\varphi)^2}{a^2}+\frac{sin(\varphi)^2}{b^2}) \cdot x^2+(\frac{sin(\varphi)^2}{a^2}+\frac{cos(\varphi)^2}{b^2}) \cdot y^2-2\cdot(\frac{cos(\varphi)sin(\varphi)}{a^2}-\frac{cos(\varphi)sin(\varphi)}{b^2})\cdot xy +2 \cdot (-\frac{x0 cos(\varphi)}{a^2}+\frac{y0 sin(\varphi)}{b^2}) \cdot x+2 \cdot (\frac{x0 sin(\varphi)}{a^2}+\frac{y0 cos(\varphi)}{b^2}) \cdot y+(\frac{x0^2}{a^2}+\frac{y0^2}{b^2}-1)=0$

а нечто гораздо сложнее. Но может быть его можно упростить за счет отбрасывания малых членов $\delta x$ $ \delta y$ $\delta \varphi$
Если бы это удалось сделать, то тогда бы нулевое приближение было бы не нужным! Но пришлось бы с самого начала находить $x_0$ $\delta x$ $y_0$ $\delta y$ $\varphi $ $\delta \varphi$ Впрочем может быть это мои иллюзии ибо как можно с самого начала отделить $x_0$ и $\delta x$ и т.д. ?

[Tur]

Алексей К. в сообщении #1376944 писал(а):

Tur,
правильно ли я понял, что у Вас всё типа получается?

По вашему алгоритму получается почти всегда и его можно было бы довести почти до абсолюта, но я оставил это пока по причине что нашел лучший способ и на половину его сделал.

Цитата:

У меня сложилось впечатление, что если бы научиться искать эллипс по трём точкам, то это было бы и надёжное первое приближение, и, возможно, вполне приличный окончательный результат (без всяких там МНКов и итераций).
Я хоть и как-то научился численно его (их) находить, но плохо и ненадёжно.
В смысле, я пока даже не придумал, как это сделать хорошо. А интересно.

Наверное это было бы выше всяких похвал, но для меня на первом месте надежность.

Алексей К. в сообщении #1376723 писал(а):

А если я ограничусь 60-ю, выберу из них $P_{1,3,5}$, то эллипс с заданными полуосями, проходящий через эти точки, будет не существовать. Потому что нарушено условие

$$\frac{b^2}a<r<\frac{a^2}b,$$

где $r$ -- радиус окружности, проходящей через эти три точки (а штуки слева и справа -- радиусы кривизны в вершинах эллипса). У 60-точечной кривой $b^2/a=17.8,\; r=172.6,\; a^2/b=166.4$). Неча и надеяться на сходимость кого-то к кому-то.

Вот поэтому я и иду путем МНК и итераций.

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

[Tur]

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

Повернем исходный эллипс $\\\frac{x^2}{a^2}+\frac{y^2}{b^2} = 1$ на $\\\varphi$

$\frac{(x\cdot cos(\varphi)+y\cdot sin(\varphi))^2}{a^2}+\frac{(-x\cdot sin(\varphi)+y\cdot cos(\varphi))^2}{b^2} = 1$

после преобразований отсюда получим ( обозначим $e = \frac{1}{a^2}-\frac{1}{b^2}$ )

$(\frac{cos^2(\varphi)}{a^2}+\frac{sin^2(\varphi)}{b^2})\cdot x^2+(\frac{sin^2(\varphi)}{a^2}+\frac{cos^2(\varphi)}{b^2})\cdot y^2+2\cdot sin(\varphi)\cdot cos(\varphi)\cdot e\cdot x\cdot y - 1 = 0$

обозначим

$a_{11} = \frac{cos^2(\varphi)}{a^2}+\frac{sin^2(\varphi)}{b^2}$

$a_{22} = \frac{sin^2(\varphi)}{a^2}+\frac{cos^2(\varphi)}{b^2}$

$a_{12} = 2\cdot sin(\varphi)\cdot cos(\varphi)\cdot e$

перепишем исходное уравнение после поворота на $\varphi$ так: $ a_{11}\cdot x^2+a_{12}\cdot y^2+2\cdot a_{12}\cdot x\cdot y - 1 = 0$

теперь сдвинем на $x_0 \, y_0$

$a_{11}\cdot(x-x_0)^2+a_{22}\cdot(y-y_0)^2+2\cdot a_{12}\cdot(x-x_0)\cdot(y-y_0) - 1 = 0$

после преобразований получим:

$a_{11} x^2+a_{22} y^2+2 a_{12}\cdot x\cdot y-2 (a_{11}\cdot x_0+a_{12}\cdot y_0)\cdot x-2 (a_{12}\cdot x_0+a_{22}\cdot y_0)\cdot y+a_{33} = 0$

обозначим

$a_{13} = -2\cdot(a_{11}\cdot x_0+a_{12}\cdot y_0)$

$a_{23} = -2\cdot(a_{12}\cdot x_0+a_{22}\cdot y_0)$

$a_{33} = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$

[Tur]

Рассмотрим случай (1) при фиксированном $\varphi$
Коэффициенты $a_{11}\,\,\, a_{22}\,\,\, a_{12}$ есть функции только угла (не центра эллипса)

Преобразуем уравнение $a_{11} x^2+a_{22} y^2+2 a_{12} x y-2 (a_{11}\cdot x_0+a_{12}\cdot y_0)x-2 (a_{12}\cdot x_0+a_{22}\cdot y_0) y+a_{33} = 0$ к виду:

$a_{11}x^2+a_{22}y^2+2 a_{12}xy-2(a_{11}x+a_{12}y)x_0-2(a_{12}x+a_{22}y)y_0+a_{33} = 0$

где $\,\,\,a_{33} = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$

обозначим

$I = a_{11}x^2+a_{22}y^2+2 a_{12}xy$

$P = -2(a_{11}x+a_{12}y)$

$Q = -2(a_{12}x+a_{22}y)$

с учетом этого перепишем уравнение: $I + P\cdot x_0 + Q\cdot y_0 + a_{33} = 0$

Коэффициенты $I \,\,\, P \,\,\, Q $ есть функции только угла, т.е. не зависят от $x_0\,\,\,y_0$

Добавим небольшие изменения в значения координат центра $\Delta x \,\,\, \Delta y$
сначала для $a_{33}$

$a_{33} = a_{11}(x_0+\Delta x)^2+a_{22}(y_0+\Delta y)^2+2a_{12}(x_0+\Delta x)(y_0+\Delta y)-1$
после преобразований этого уравнения и отбрасывания малых членов получим:

$a_{33} = 2\cdot P_0 \cdot \Delta x + 2\cdot Q_0\cdot\Delta y + a_0$

где

$P_0 = 2\cdot(a_{11}\cdot x_0+a_{12}\cdot y_0)$

$Q_0 = 2\cdot(a_{12}\cdot x_0+a_{22}\cdot y_0)$

$a_0 = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$

Добавим небольшие изменения в значения координат центра $\Delta x \,\,\, \Delta y$ в основное ур-е в этом случае:

$I + P(x_0+\Delta x)+Q(y_0+\Delta y) + a_{33} = 0$

после преобразований получим

$(P + P_0)\Delta x+(Q+Q_0)\Delta y+I+P_0+Q_0+a_0 = 0$

Это переопределенная система. Решаем ее и находим $\Delta x$ и $\Delta y$
Вот код

(Оффтоп)

function [dx dy] = dx_dy(fi, Mo, ab, x, y)
xo = Mo.X; xo2 = xo*xo;
yo = Mo.Y; yo2 = yo*yo;
x2 = x.*x; y2 = y.*y; xy = x.*y;
a = ab(1); a2 = a*a;
b = ab(2); b2 = b*b;
s = sin(fi); s2 = s*s;
c = cos(fi); c2 = c*c;
a11 = c2/a2 + s2/b2;
a22 = s2/a2 + c2/b2;
e = 1/a2 - 1/b2;
a12 = e*s*c;
I = a11*x2 + a22*y2 + 2*a12*xy;
P = -2*(a11*x + a12*y);
Po = 2*(a11*xo + a12*yo);
Q = -2*(a12*x + a22*y);
Qo = 2*(a12*xo + a22*yo);
ao = a11*xo2 + a22*yo2 + 2*a12*xo*yo - 1;
PP = P + Po; QQ = Q + Qo;
II = I + P*xo + Q*yo +ao;
M = [PP QQ II]; X = pere_slau(M);
dx = -X(1); dy = -X(2);
end

[Tur]

Случей 2. Будем искать изменение угла $\Delta \varphi$ при фиксированном центре.
исходное уравнение:
$a_{11} x^2+a_{22} y^2+2 a_{12}\cdot x\cdot y-2 (a_{11}\cdot x_0+a_{12}\cdot y_0)\cdot x-2 (a_{12}\cdot x_0+a_{22}\cdot y_0)\cdot y+a_{33} = 0$ ур-е (*)

здесь

$a_{11} = \frac{cos^2\varphi}{a^2}+\frac{sin^2\varphi}{b^2}$

$a_{22} = \frac{sin^2\varphi}{a^2}+\frac{cos^2\varphi}{b^2}$

$a_{12} = 2\cdot sin\varphi\cdot cos\varphi\cdot e$

$e = \frac{1}{a^2}-\frac{1}{b^2}$

Добавим небольшое изменение угла $\Delta \varphi$

$a_{11} = \frac{cos^2(\varphi + \Delta \varphi)}{a^2}+\frac{sin^2(\varphi + \Delta \varphi)}{b^2}$

$cos^2(\varphi + \Delta \varphi) = (cos\varphi\cdot cos\Delta \varphi-sin\varphi\cdot sin\Delta\varphi)^2 = (cos\varphi-\Delta\varphi\cdot sin\varphi)^2$

$cos^2(\varphi + \Delta \varphi) = cos^2\varphi-\Delta\varphi\cdot sin(2\varphi)$

$sin^2(\varphi+\Delta\varphi)=sin^2\varphi+\Delta\varphi\cdot sin(2\varphi)$

$a_{11} = \frac{cos^2\varphi-\Delta\varphi sin(2\varphi)}{a^2}+\frac{sin^2\varphi+\Delta\varphi sin(2\varphi)}{b^2}=\frac{cos^2\varphi}{a^2}+\frac{sin^2\varphi}{b^2}-\Delta\varphi\cdot e\cdot sin(2\varphi)$

$a_{11} = A - \Delta\varphi\cdot C$

$a_{22} = B + \Delta\varphi\cdot C$

$a_{12}= \frac{C}{2}+\Delta\varphi\cdot D$

$A =\frac{cos^2\varphi}{a^2}+\frac{sin^2\varphi}{b^2}$

$B =\frac{sin^2\varphi}{a^2}+\frac{cos^2\varphi}{b^2}$

$C = e\cdot sin(2\varphi)$

$D = e\cdot cos(2\varphi)$

далее преобразуем

$I = a_{11}x^2+a_{22}y^2+2 a_{12}xy$

$I=(A-\Delta\varphi C)x^2+(B+\Delta\varphi C)y^2+(C+2\Delta\varphi D)xy$

$I=\Delta\varphi (C(y^2-x^2)+2Dxy)+Ax^2+By^2+Cxy$

$I=\Delta\varphi \cdot I_1+I_2$

где $I_1=C(y^2-x^2)+2Dxy\,\,\,\,\,\,\,\, I_2=Ax^2+By^2+Cxy$

введем еще обозначение:

$J=-2(a_{11}x_0+a_{12}y_0)x-2(a_{12}x_0+a_{22}y_0)y$

Тогда уравнение (*) запишется так: $I+J+a_{33}=0$

Рассмотрим подробней J

$J=-2(a_{11}x_0x+a_{12}y_0x+a_{12}x_0y+a_{22}y_0y)$

$J=-2((A-\Delta\varphi\cdot C)x_0x+(\frac{C}{2}+\Delta\varphi\cdot D)y_0x+(\frac{C}{2}+\Delta\varphi\cdot D)x_0y+(B+\Delta\varphi\cdot C)y_0y)$

$J=-2(\Delta\varphi\cdot(C(y_0y-x_0x)+D(x_0y+y_0x))+Ax_0+By_0y+\frac{C}{2}(x_0y+y_0x))$

обозначим

$J_1 = \Delta\varphi\cdot(C(y_0y-x_0x)+D(x_0y+y_0x)$

$J_2 = Ax_0+By_0y+\frac{C}{2}(x_0y+y_0x)$

тогда $J=-2\cdot (\Delta\varphi\cdot J_1 + J_2)$

Подставляя все это в (*) получаем $\Delta\varphi\cdot I_1+I_2-2\cdot\Delta\varphi\cdot J_1-2\cdot J_2 + a_{33}=0$
или

$\Delta\varphi\cdot(I_1-2\cdot J_1)+I_2-2\cdot J_2+a_{33}=0$

Решаем эту переопределенную систему и находим $\Delta\varphi$

Вот код

(Оффтоп)

function dfi = d_fi(fi, Mo, ab, x, y)
xo = Mo.X; xo2 = xo*xo;
yo = Mo.Y; yo2 = yo*yo;
x2 = x.*x; y2 = y.*y; xy = x.*y;
a = ab(1); a2 = a*a;
b = ab(2); b2 = b*b;
s = sin(fi); s2 = s*s;
c = cos(fi); c2 = c*c;
A = c2/a2 + s2/b2;
B = s2/a2 + c2/b2;
e = 1/a2 - 1/b2;
C = 2*e*s*c;
D = e*(c2 - s2);
I1 = C*(y2 - x2) + 2*D*xy;
I2 = A*x2 + B*y2 + C*xy;
J1 = C*(yo*y - xo*x) + D*(xo*y + yo*x);
J2 = A*xo*x + B*yo*y + C/2*(xo*y + yo*x);
K = C*(yo2 - xo2) + 2*D*xo*yo;
ao = A*xo2 + B*yo2 + C*xo*yo - 1;
T = I1 - 2*J1 + K;
S = I2 - 2*J2 + ao;
M = [T S]; X = pere_slau(M); dfi = -X(1);
end

[Tur]

К последнему уравнению случая 2 забыл добавить раскрытие $a_{33}$

$a_{33} = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$
или
$a_{33} = (A-\Delta\varphi\cdot C)\cdot x_0^2+(B+\Delta\varphi\cdot C)\cdot y_0^2+2\cdot (\frac{C}{2}+\Delta\varphi\cdot D)x_0y_0-1$

$a_{33}=\Delta\varphi\cdot (C\cdot(y_0^2-x_0^2)+2\cdot D\cdot x_0\cdot y_0)+A\cdot x_0^2+B\cdot y_0^2+C\cdot x_0\cdot y_0-1$

обозначим $K =(C\cdot(y_0^2-x_0^2)+2\cdot D\cdot x_0\cdot y_0)\,\,\,\,\,\,a_0=A\cdot x_0^2+B\cdot y_0^2+C\cdot x_0\cdot y_0-1$

$a_{33}=\Delta\varphi\cdot K+a_0$ $\,\,\,\,\,\,\, .$ В коде случая 2 видно использование $K$ и $a_0$
=========================================================================================

Перехожу к случаю 3 когда ничто не фиксировано кроме начальных величин осей эллипса. Буду пользоваться прежними обозначениями.

$a_{11}\cdot x^2+a_{22}\cdot y^2+2 a_{12}\cdot x\cdot y-2\cdot(a_{11}\cdot x_0+a_{12}\cdot y_0)\cdot x-2\cdot(a_{12}\cdot x_0+a_{22}\cdot y_0)\cdot y+a_{33} = 0$ ур-е (*)

$I = a_{11}x^2+a_{22}y^2+2 a_{12}xy$

$J=-2(a_{11}x_0+a_{12}y_0)x-2(a_{12}x_0+a_{22}y_0)y$

$a_{33} = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$

уравнение (*) запишется так: $I+J+a_{33}=0$

С учетом введенных ранее обозначений запишем каждый элемент уравнения (*) как функцию $\,\,\, \Delta\varphi \,\,\, \Delta x\,\,\, \Delta y$

$I=\Delta\varphi \cdot I_1+I_2\,\,\,\,\,\,\,\,\, -$ функция только угла

$I_1=C(y^2-x^2)+2Dxy\,\,\,\,\,\,\,\, I_2=Ax^2+By^2+Cxy$

$J=-2\cdot(a_{11}\cdot x_0+a_{12}\cdot y_0)\cdot x-2\cdot(a_{12}\cdot x_0+a_{22}\cdot y_0)\cdot y$

Подставим в это ур-е значения

$a_{11} = A - \Delta\varphi\cdot C$

$a_{22} = B + \Delta\varphi\cdot C$

$a_{12}= \frac{C}{2}+\Delta\varphi\cdot D$

и вместо $x_0$ подставим $x_0+\Delat x$ а вместо $y_0$ подставим $y_0+\Delta y$

после преобразований получаем

$J=-2\cdot(\Delta\varphi\cdot J_1+\Delta x\cdot J_2+\Delta y\cdot J_3+J_4)$

$J_1=C\cdot(y_0\cdot y-x_0\cdot x)+D\cdot(y_0\cdot x+x_0\cdot y)(y_0+\Delta y)-1$

$J_2=A\cdot x+\frac{C}{2}\cdot y\,\,\,\,\,\,\,\,\,\,\, J_3=\frac{C}{2}\cdot x+B\cdot y\,\,\,\,\,\,\,\,\,\, J_4=Ax_0x+By_0y+\frac{C}{2}(y_0x+x_0y)$

Теперь в $a_{33} = a_{11}\cdot x_0^2+a_{22}\cdot y_0^2+2\cdot a_{12}\cdot x_0\cdot y_0-1$ подставляем значения коэффициентов а11 а22 а12 и смещаем центр

$a_{33}=(A-\Delta\varphi C)(x_0+\Delta x)^2+(B+\Delta\varphi C)(y_0+\Delta y)^2+2(\frac{C}{2}+\Delta\varphi D)(x_0+\Delta x)(y_0+\Delta y)-1$

после преобразований получаем

$a_{33}=\Delta\varphi \cdot E_1+\Delta x\cdot E_2+\Delta y\cdot E_3+E_4$

$E_1=C(y_0^2-x_0^2)+2Dx_0y_0\,\,\,\,\,\,\,\,\,\,\,E_2=2Ax_0+Cy_0$

$E_3=2By_0+Cx_0\,\,\,\,\,\,\,\,\,\,\,E_4=Ax_0^2+By_0^2+Cx_0y_0-1$

подставляем найденные значения в $I+J+a_{33}=0$ и после преобразований получаем

$\Delta\varphi\cdot P+\Delta x\cdot Q+\Delta y\cdot S+T=0\,\,\,\,\,\,\,\,\,\,\,\,\,(**)\,\,\,\, -$ Это и есть финальное, целевое уравнение.

$P=I_1-2\cdot J_1+E_1\,\,\,\,\,\,\,\,\,\,Q=-2\cdot J_2+E_2\,\,\,\,\,\,\,\,\,\,S=-2\cdot J_3+E_3\,\,\,\,\,\,\,\,\,\,T+I_2-2\cdot J_4+E_4$

Получили переопределенную систему с тремя неизвестными $\Delta\varphi \,\,\, \Delta x\,\,\, \Delta y\, .$ Решаем ее и находим их.
Вот код

(Оффтоп)

function [dfi dx dy] = ddd(fi, Mo, ab, x, y)
xo = Mo.X; xo2 = xo*xo;
yo = Mo.Y; yo2 = yo*yo;
x2 = x.*x; y2 = y.*y; xy = x.*y;
a = ab(1); a2 = a*a;
b = ab(2); b2 = b*b;
s = sin(fi); s2 = s*s;
c = cos(fi); c2 = c*c;
A = c2/a2 + s2/b2;
B = s2/a2 + c2/b2;
e = 1/a2 - 1/b2;
C = 2*e*s*c;
D = e*(c2 - s2);

I1 = C*(y2 - x2) +2*D*xy;
I2 = A*x2 + B*y2 + C*xy;
J1 = C*(yo*y - xo*x) + D*(yo*x + xo*y);
J2 = A*x + C/2*y; J3 = C/2*x + B*y;
J4 = A*xo*x + B*yo*y + C/2*(yo*x + xo*y);
E1 = C*(yo2 - xo2) + 2*D*xo*yo;
E2 = 2*A*xo + C*yo; E3 = 2*B*yo + C*xo;
E4 = A*xo2 + B*yo2 + C*xo*yo - 1;
P = I1 - 2*J1 + E1;
Q = -2*J2 + E2;
S = -2*J3 + E3;
T = I2 - 2*J4 + E4;
M = [P Q S T]; X = pere_slau(M);
dfi = -X(1); dx = -X(2); dy = -X(3);
end

[Tur]

Без опечаток никак не мог. Вот нашел одну в последнем посту.

написано: $J_1=C\cdot(y_0\cdot y-x_0\cdot x)+D\cdot(y_0\cdot x+x_0\cdot y)(y_0+\Delta y)-1$
следует читать: $J_1=C\cdot(y_0\cdot y-x_0\cdot x)+D\cdot(y_0\cdot x+x_0\cdot y)$

Вообще во всей теме полно опечаток и недоразумений. Но код работает. Поэтому следует сверяться с кодом.

[fan_of_algoritms]

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