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

[Алексей К.]

Я додумаю в вскр. (щас типа на работу надолго). Пока зафиксирую кое-что из памяти.

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

...я когда-то уже думал над выбором условия нормировки коэффициентов коники

$$Ax^2+2Bxy+Cy^2+2Dx+2Ey+F=0$$

(и вряд ли это было в плане МНК; не помню). В качестве такового у меня получилось

$$B^2=(1-A)(1-C),\text{~~~~или~~~~}B^2-AC+A+C=1.$$

И с ним никаких особенностей не должно возникнуть.

Эту нормировку я выбирал из тех соображений, чтобы при повороте-переносе получить неявное уравнение кривой 2-го порядка в каноническом (одном из канонических) виде

$$ y^2 + (1-e^2)x^2 + 2 d e^2x - d^2 e^2 = 0 $$

($d$ --- абсцисса (вертикальной) директрисы, фокус --- в начале координат, $e$ и $p=ed$ --- наверняка эксцентриситет и фокальный параметр).

[Pavia]

Tur

Tur в сообщении #1367777 писал(а):

так что скажете по моему заходу? Нельзя полагать равным единице ни один из коэф-тов: а11 а22 а12?

Выпишем уравнение эллипса в параметрической форме:
$u=a\cdot cos(t)$
$v=b\cdot sin(t)$

Добавим вращение и смещение центра.
$x=u \cdot cos(\phi)- v \cdot sin(\phi)+x0$
$y=u \cdot sin(\phi)+ v \cdot cos(\phi)+y0$


Выведем уравнение эллипса в канонической форме.
Для чего представим в матричной форме и возьмем обратную матрицу. Очевидно что она будет такой же но угол с минусом и перемещение с обратным знаком.
$u= x \cdot cos(\phi)+ y \cdot sin(\phi)-x0$
$v=- x \cdot sin(\phi)+ y \cdot cos(\phi)-y0$

$a \cdot cos(t)= x \cdot cos(\phi)+ y \cdot sin(\phi)-x0$
$b \cdot sin(t)=- x \cdot sin(\phi)+ y \cdot cos(\phi)-y0$

Возведём в квадрат и подставим в основное тригометрическое тожество.

$cos(t)=( x \cdot cos(\phi)+ y \cdot sin(\phi)-x0)/a$
$sin(t)=(- x \cdot sin(\phi)+ y \cdot cos(\phi)-y0)/b$

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

Раскроем скобки и сгруппируем параметры

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

И сделаем замену
$a_{11} \cdot x^2 + a_{22} \cdot y^2 +2\cdot a_{12} \cdot xy+2\cdot a_{13} \cdot x+2 \cdot a_{23} \cdot y+a_{33}=0$

При фиксированных $a$ и $b$ из коэффициента $a_{12}$ легко найти угол вращения.
И далее вычислить центр эллипса.

[Алексей К.]

Pavia в сообщении #1367820 писал(а):

Возведём в квадрат и подставим в основное тригометрическое тожество.

Кнут, cepesh, Корн, Бронштейн, gris, Семендяев --- эти и многие другие люди приучили меня, что основное тригометрическое тожество (ОТТ) есть

$$\cos^2 t +\sin^2 t=1.$$

С Вашим, Pavia, вариантом ОТТ, типа

$$cos(t)^2 + sin(t)^2=1,$$

я согласиться никак не могу.

[Pavia]

Алексей К.
Согласен. Бывает ошибаюсь, устал очень, не уследил. Кстати там ещё один квадрат возле $a$ потерял, в следующей формуле его восстановил.

[Tur]

Pavia, спасибо! Я думал в том же направлении только без параметрической формы, а сразу подстановкой
$x = x\cdot cos(\phi) + y\cdot sin(\phi) - x_0$
$y = -x\cdot sin(\phi) + y\cdot cos(\phi) - y_0$
в уравнение
$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$
с дальнейшей группировкой (не знаете есть ли где то в литературе такая подстановка приведена?) и понял что длины полуосей вкраплены в коэф-ты а11 а22..., которые получаю в матлабе так

function [a11 a12 a13 a22 a23] = aquadric(x, y)
% входные x, y - координаты точек дуги
x2 = x.*x; y2 = y.*y; xy = x.*y;
one = ones(length(x),1); 
X = [x2 y2 xy x y]; Y = -one;
R = (X'*X)^(-1)*X'*Y;
a11 = R(1); a22 = R(2); a12 = R(3)/2; 
a13 = R(4)/2; a23 = R(5)/2; 
end

и далее c небольшими сокращениями код такой

M = [a11 a12 a13; a12 a22 a23; a13 a23 1];
d = det(M);
D = a11*a22 - a12*a12;
Mo.X = (a12*a23 - a13*a22)/D; 
Mo.Y = (a13*a12 - a11*a23)/D;
tg = 2*a12/(a11-a22); fi = 0.5*atan(tg);
s = sin(fi); c = cos(fi);
lambda1 = a11*c*c+2*a12*s*c+a22*s*s;
lambda2 = a11*s*s-2*a12*s*c+a22*c*c;
a2 = -d/(lambda1*lambda1*lambda2); 
b2 = -d/(lambda1*lambda2*lambda2);

т.е. полагаю а33 = 1 (не знаю как его найти в лоб, может быть вы знаете?)

Все эти коэф-ты а11 а22... получаются не нормативного эллипса. Поэтому у меня возникли сомнения а получу ли я тот же самый наклон из
уравнения
$a_{11} = \frac{cos^2(\phi)}{a^2}+\frac{sin^2(\phi)}{(b^2}$ или подобного ур-я а22 = ... или а12 = ...
Это в случае если каким то образом удастся найти а33, а если нет, то вообще неясно как находить наклон - так я думал и мысли повернулись в сторону НС. Дело в том, что у меня же таблетки не эллипсы, поэтому НС вероятно следующий шаг (я на него еще не решился, это же время...)

Если а33 найти не удастся, т.е. найденные коэф-ты а11 а22... на самом деле являются отношениями а11/а33 а22/а33 ... поэтому наклон можно будет найти из отношения а44 = а11/а22 т.е. при этом а33 сократиться, а44 будет известным числом и получится такая формула
$tg^2(\phi) = \frac {a44\cdot a^2-b^2}{a^2-a44\cdot b^2}$

Я понимаю, что пишу кратко и возможно неясно, это ввиду... не важно, если что то нужно раскрыть - с радостью. Просто боюсь не успеть до ночи.
Главная проблема думаю ясна:
1. а33 отсутствует
2. если а33 найти, то точно ли что из разных уравнений (замен на а11 а12...) получим теже самые углы?

Еще раз большое спасибо за большую проделанную работу!

У меня бывает, что аппроксимация дает не эллипс, а гиперболу - отсюда мои сомнения что наклон и центр можно однозначно найти из ур-я а12 = ...

[Pavia]

Tur в сообщении #1368081 писал(а):

а получу ли я тот же самый наклон из

Да в приделах погрешности.

Tur в сообщении #1368081 писал(а):

У меня бывает, что аппроксимация дает не эллипс, а гиперболу

Это у вас формула aquadric неправильная.

$a_{33}$ ищется так

function [a11 a12 a13 a22 a23 a33] = aquadric(x, y)
% входные x, y - координаты точек дуги
x2 = x.*x; y2 = y.*y; xy = x.*y;
one = ones(length(x),1); 
X = [x2 y2 xy x y 1]; Y = one; // По моему тут минус лишний
R = (X'*X)^(-1)*X'*Y;
a11 = R(1); a22 = R(2); a12 = R(3)/2; 
a13 = R(4)/2; a23 = R(5)/2; 
a33 = R(6) - 1;  
end

[Tur]

Pavia, спасибо за разъяснения
Почему то не получается т.о. верно находить параметры
слева по старому алгоритму, справа по исправленному (там и там эллипс произвольный)

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

На всякий случай приведу два кода drawell - старый, drawell2 - новый, x y дуги в самом конце. Эти два кода идентичны за исключением метода нахождения параметра а33, а именно:

старый вариант
Y = -ones(length(x),1); 
X = [x2 y2 xy x y];
R = (X'*X)^(-1)*X'*Y;
a11 = R(1); a22 = R(2); a12 = R(3)/2;
a13 = R(4)/2; a23 = R(5)/2;  a33 = 1;

новый вариант
Y = ones(length(x),1);
X = [x2 y2 xy x y Y];
R = (X'*X)^(-1)*X'*Y;
a11 = R(1); a22 = R(2); a12 = R(3)/2; 
a13 = R(4)/2; a23 = R(5)/2; 
a33 = R(6) - 1;

Кроме того старые а11... порядка от десятых до тысячных, а новые порядка $10^{-17}$
минус перед Y не влияет на порядок

старый вариант drawell

(Оффтоп)

function drawell(x, y)
[MT Mo fi a b F1 F2] = parameters_ell(x, y);
t = -0.1:0.1:2*pi; x1 = a*cos(t); y1 = b*sin(t);
invM = inv(MT); [x2 y2] = movearr(x1', y1', invM); 
plot(x2, y2, '-r'); plot(Mo.X, Mo.Y, '.w')
plot(F1.X, F1.Y, '.w'); plot(F2.X, F2.Y, '.w')
end
function [MT Mo fi a b F1 F2] = parameters_ell(x, y)
mx = mean(x); my = mean(y); xn = x - mx; yn = y - my;
MT1 = [1 0 -mx; 0 1 -my; 0 0 1];
[M Mo fi s c a2 b2] = param_ell(xn, yn);
if a2 < b2 
  tmp = b2; b2 = a2; a2 = tmp; 
  s = sin(fi+pi/2); c = cos(fi+pi/2);
end
a = sqrt(a2); b = sqrt(b2);
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);
tg = s/c; p = sqrt(a2-b2)/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 [M Mo fi s c a2 b2] = param_ell(x, y)
[a11 a12 a13 a22 a23] = aquadric(x, y);
M = [a11 a12 a13; a12 a22 a23; a13 a23 1];
d = det(M);
D = a11*a22 - a12*a12;
Mo.X = (a12*a23 - a13*a22)/D;
Mo.Y = (a13*a12 - a11*a23)/D;
tg = 2*a12/(a11-a22); fi = 0.5*atan(tg);
s = sin(fi); c = cos(fi);
lambda1 = a11*c*c+2*a12*s*c+a22*s*s;
lambda2 = a11*s*s-2*a12*s*c+a22*c*c;
a2 = -d/(lambda1*lambda1*lambda2); 
b2 = -d/(lambda1*lambda2*lambda2);
end
function [a11 a12 a13 a22 a23] = aquadric(x, y)
x2 = x.*x; y2 = y.*y; xy = x.*y;
Y = -ones(length(x),1); 
X = [x2 y2 xy x y];
R = (X'*X)^(-1)*X'*Y; % solution of system
a11 = R(1); a22 = R(2); a12 = R(3)/2;
a13 = R(4)/2; a23 = R(5)/2; 
end
function NewP = mp(M, P) 
% NewP = M*P matrix multiplied by the point
tmp = M*[P.X P.Y 1]'; NewP.X = tmp(1); NewP.Y = tmp(2);
end
function [xn yn] = movearr(x, y, M)
n = length(x); one = ones(n,1); xy = [x y one]; 
temp = M*xy'; xn = temp(1,:)'; yn = temp(2,:)';
end

новый вариант drawell2

(Оффтоп)

function drawell2(x, y)
[MT Mo fi a b F1 F2] = parameters_ell(x, y);
t = -0.1:0.1:2*pi; x1 = a*cos(t); y1 = b*sin(t);
invM = inv(MT); [x2 y2] = movearr(x1', y1', invM); 
plot(x2, y2, '-r'); plot(Mo.X, Mo.Y, '.w')
plot(F1.X, F1.Y, '.w'); plot(F2.X, F2.Y, '.w')
end
function [MT Mo fi a b F1 F2] = parameters_ell(x, y)
mx = mean(x); my = mean(y); xn = x - mx; yn = y - my;
MT1 = [1 0 -mx; 0 1 -my; 0 0 1];
[Mo fi s c a2 b2] = param_ell(xn, yn);
if a2 < b2 
  tmp = b2; b2 = a2; a2 = tmp; 
  s = sin(fi+pi/2); c = cos(fi+pi/2);
end
a = sqrt(a2); b = sqrt(b2);
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);
tg = s/c; p = sqrt(a2-b2)/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 [Mo fi s c a2 b2] = param_ell(x, y)
[a11 a12 a13 a22 a23 a33] = aquad(x, y);
M = [a11 a12 a13; a12 a22 a23; a13 a23 a33];
d = det(M);
D = a11*a22 - a12*a12;
Mo.X = (a12*a23 - a13*a22)/D;
Mo.Y = (a13*a12 - a11*a23)/D;
tg = 2*a12/(a11-a22); fi = 0.5*atan(tg);
s = sin(fi); c = cos(fi);
lambda1 = a11*c*c+2*a12*s*c+a22*s*s;
lambda2 = a11*s*s-2*a12*s*c+a22*c*c;
a2 = -d/(lambda1*lambda1*lambda2); 
b2 = -d/(lambda1*lambda2*lambda2);
end
function [a11 a12 a13 a22 a23 a33] = aquad(x, y)
x2 = x.*x; y2 = y.*y; xy = x.*y;
Y = ones(length(x),1);
X = [x2 y2 xy x y Y];
R = (X'*X)^(-1)*X'*Y;
a11 = R(1); a22 = R(2); a12 = R(3)/2; 
a13 = R(4)/2; a23 = R(5)/2; 
a33 = R(6) - 1;
end
function NewP = mp(M, P) 
% NewP = M*P matrix multiplied by the point
tmp = M*[P.X P.Y 1]'; NewP.X = tmp(1); NewP.Y = tmp(2);
end
function [xn yn] = movearr(x, y, M)
n = length(x); one = ones(n,1); xy = [x y one]; 
temp = M*xy'; xn = temp(1,:)'; yn = temp(2,:)';
end

x y дуги

(Оффтоп)

235 2
236 3
237 4
238 5
239 6
240 6
241 7
242 8
243 8
244 9
245 10
246 11
247 11
248 12
249 13
250 13
251 14
252 14
253 15
254 16
255 16
256 17
257 17
258 18
259 18
260 19
261 20
262 20
263 21
264 21
265 22
266 22
267 23
268 23
269 24
270 24
271 24
272 25
273 26
274 26
275 26
276 27
277 27
278 28
279 28
280 28
281 29
282 29
283 30
284 30
285 30
286 30
287 31
288 32
289 32
290 32
291 32
292 33
293 33
294 34
295 34
296 34
297 34
298 34
299 34
300 34
301 34
302 34
303 34
304 35
305 35
306 36
307 36
308 36
309 36
310 37
311 37
312 37
313 37
314 37
315 37
316 37
317 37
318 37
319 37
320 36
321 36
322 36
323 36
324 36
325 35
326 34
327 34
328 34
329 33
330 32
331 32
332 31
333 30
333 29
334 28
335 27
335 26
336 25
336 24
336 23
337 22
337 21
338 20
338 19
338 18
339 17
339 16
340 15
340 14
340 13
340 12
340 11
340 10
340 9
340 8
340 7
340 6
340 5
340 4
340 3
340 2

[Pavia]

(X'*X)

- здесь переполнение для чисел double(вернее исчерпание разрядов мантиссы). Поэтому и срывается.

[Tur]

Pavia в сообщении #1368355 писал(а):

(X'*X)

- здесь переполнение для чисел double(вернее исчерпание разрядов мантиссы). Поэтому и срывается.

Ок, спасибо, но я же вычитаю средние.
В drawell и drawell2 можно видеть что сначала я нахожу
xn = x - mean(x); yn = y - mean(y);
Именно эти xn yn входят в aquadric, для этих значений нахожу параметры эллипса и затем матрицу преобразования и с ее помощью точки x1 y1 переношу обратно в исходную область.
Может быть исходные точки x y следует сначала перевести в область от 0 до 1?
Как решать эту проблему переполнения? И почему ее нет в старом алгоритме?

Вот матрица М = X'*X для первой дуги

ее максимум 2.8116e+08, минимум (абсолютной величины) 1.1369e-13
Пределы double от $10^{-308 }$ до $10^{308}$
Может в чем то другом переполнение или оно не видно?

На всякий случай вот дуга побольше, но снова срыв

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

(Оффтоп)

215 2
216 3
216 4
217 5
218 6
219 7
219 8
220 9
221 10
222 11
222 12
223 13
224 14
225 15
225 16
226 17
227 18
228 19
229 20
230 21
231 22
232 23
233 24
234 25
235 26
235 27
236 27
237 28
238 29
239 30
240 31
241 32
242 33
243 34
244 35
245 35
246 36
247 37
248 38
249 38
250 39
251 40
252 41
253 41
254 42
255 43
256 43
257 44
258 45
259 45
260 46
261 47
262 47
263 48
264 48
265 49
266 50
267 50
268 51
269 51
270 52
271 52
272 53
273 53
274 54
275 54
276 55
277 55
278 56
279 56
280 57
281 57
282 57
283 58
284 58
285 59
286 59
287 60
288 60
289 60
290 61
291 61
292 61
293 62
294 62
295 62
296 62
297 63
298 63
299 64
300 64
301 64
302 64
303 64
304 64
305 64
306 65
307 66
308 66
309 66
310 66
311 66
312 66
313 66
314 66
315 66
316 66
317 66
318 66
319 66
320 66
321 65
322 65
323 64
324 64
325 64
326 63
327 62
328 62
329 61
330 60
331 60
332 59
333 58
334 57
334 56
335 55
336 54
336 53
337 52
338 51
338 50
338 49
339 48
339 47
340 46
340 45
340 44
340 43
340 42
340 41
340 40
340 39
340 38
340 37
340 36
340 35
340 34
340 33
340 32
339 31
339 30
338 29
338 28
337 27
337 26
336 25
336 24
335 23
334 22
334 21
333 20
333 19
332 18
331 17
330 16
330 15
329 14
328 13
328 12
327 11
326 10
326 9
325 8
324 7
323 6
322 5
322 4
321 3
320 2

-- Вс янв 13, 2019 20:59:33 --

Вот еще третья дуга 265 точек

ее матрица М = X'*X от 9 до -13 степени
просто чертовщина, неясно куда копать
третья дуга 265 точек

(Оффтоп)

98 2
97 3
96 3
95 4
94 4
93 4
92 4
91 5
90 5
89 6
88 6
87 6
86 7
85 8
84 8
83 9
82 9
81 10
80 11
79 12
78 13
77 14
76 15
75 16
75 17
74 18
74 19
73 20
72 21
72 22
72 23
72 24
71 25
71 26
71 27
71 28
71 29
71 30
71 31
71 32
71 33
71 34
71 35
72 36
72 37
72 38
72 39
72 40
73 41
73 42
74 43
74 44
75 45
76 46
77 47
77 48
78 49
79 50
80 50
81 51
82 52
83 52
84 53
85 54
86 54
87 55
88 55
89 56
90 56
91 57
92 57
93 58
94 58
95 58
96 59
97 59
98 60
99 60
100 60
101 61
102 61
103 62
104 62
105 62
106 63
107 63
108 63
109 64
110 64
111 64
112 64
113 65
114 65
115 66
116 66
117 66
118 66
119 66
120 67
121 67
122 67
123 67
124 67
125 68
126 68
127 68
128 68
129 68
130 68
131 68
132 68
133 69
134 69
135 69
136 69
137 70
138 70
139 70
140 70
141 70
142 70
143 70
144 70
145 70
146 70
147 70
148 70
149 70
150 70
151 70
152 70
153 70
154 70
155 70
156 70
157 70
158 70
159 70
160 70
161 70
162 70
163 70
164 70
165 70
166 70
167 70
168 70
169 70
170 70
171 69
172 69
173 69
174 69
175 69
176 69
177 69
178 68
179 68
180 68
181 68
182 68
183 68
184 68
185 67
186 67
187 66
188 66
189 66
190 66
191 66
192 66
193 65
194 65
195 65
196 65
197 64
198 64
199 64
200 64
201 63
202 63
203 62
204 62
205 62
206 61
207 61
208 60
209 60
210 60
211 59
212 58
213 58
214 57
215 56
216 56
217 55
218 54
219 53
220 53
220 52
221 51
222 50
222 49
222 48
223 47
223 46
224 45
224 44
224 43
224 42
225 41
225 40
225 39
225 38
226 37
226 36
226 35
226 34
226 33
226 32
226 31
225 30
225 29
225 28
224 27
224 26
224 25
223 24
223 23
222 22
222 21
221 20
220 19
220 18
219 17
218 16
217 15
216 14
215 14
214 13
213 13
212 12
211 12
210 11
209 10
208 10
207 10
206 9
205 8
204 8
203 8
202 7
201 7
200 6
199 6
198 6
197 5
196 5
195 4
194 4
193 4
192 4
191 4
190 3
189 2

[Pavia]

Tur
Вот и у меня так же.

[Tur]

Pavia в сообщении #1368433 писал(а):

Tur
Вот и у меня так же.

Спасибо что срочно ответили. У меня когда то такое уже было, но я так и бросил тогда, должен быть баланс между нуждой и временем. Я думаю что вы совершенно правы насчет вами предложенного решения той переопределенной системы, это же азы. Но это же не мистика, должно быть и объяснение и решение. Надеюсь что это только начало. Завтра попробую 1)перевести все в квадрат 0-1 и 2) найти центр через отношение а11/а22, как я уже писал, при а33 = 1

Сегодня старый новый год, будем считать эти срывы новогодним чудом

[Pavia]

Tur
Убрал фиксированные оси заменил на вычисления.

Delta:=Determinant(M3);
                
cent[0] := a[0]*cos(alfa)*cos(alfa)+2*a[2]*sin(alfa)*cos(alfa)+a[1]*sin(alfa)*sin(alfa);
cent[1] := a[0]*sin(alfa)*sin(alfa)-2*a[2]*sin(alfa)*cos(alfa)+a[1]*cos(alfa)*cos(alfa);

a_axis:= Sqrt(abs( -Delta/(cent[0]*cent[0]*cent[1])));
b_axis:= Sqrt(abs( -Delta/(cent[0]*cent[1]*cent[1])));

Изредка косяки. Но в нуль не уходит.

-- Пн янв 14, 2019 01:31:22 --

Tur в сообщении #1368376 писал(а):

Вот еще третья дуга 265 точек

После как исправил у меня норм. $a_{33}$ Вычисляется.

[Fizykochemik]

Мне кажется дискуссия слишком усложнилась и скорее дезориентирует.
Правильно ли я понимаю, что Вам нужен эллипс Mo, а выходят Mo1, ...Moi и т.д. для известных дуг i?
Дуга это точки расположенные рядом?

Вам известны фиксированные полуоси эллипса $(a,b)$, но координаты центра неизвестны?

Я бы сделал так (если я правильно понял).
Выражаете из ур-я эллипса $y$ через $x$, далее берете прогу для алгоритма нелинейных квадратов
Левенберга-Марквардта (ссылки или сами проги могу дать), и находите нужные неизвестные параметры.
Есть проблема связанная со сменой знака при корне в правой части, поэтому лучше перейти в обобщенные полярные координаты.
Более того, я бы попробовал аппроксимировать более сложной кривой, например 4го порядка.

Если неизвестна ориентация эллипса относительно системы координат. Это не проблема, так как в оптимизированных параметрах будут содержаться углы
поворота к канонической форме. Главное чтобы число параметров было $\leq$ числа точек.

[Emergency]

Fizykochemik в сообщении #1368546 писал(а):

Правильно ли я понимаю, что Вам нужен эллипс

Насколько я понял, нужен не эллипс, а его центр. То есть задача в принципе уже решена, хотя иногда врет с размерами.

[Tur]

Pavia в сообщении #1368469 писал(а):

Tur
Убрал фиксированные оси заменил на вычисления.

Delta:=Determinant(M3);
                
cent[0] := a[0]*cos(alfa)*cos(alfa)+2*a[2]*sin(alfa)*cos(alfa)+a[1]*sin(alfa)*sin(alfa);
cent[1] := a[0]*sin(alfa)*sin(alfa)-2*a[2]*sin(alfa)*cos(alfa)+a[1]*cos(alfa)*cos(alfa);

a_axis:= Sqrt(abs( -Delta/(cent[0]*cent[0]*cent[1])));
b_axis:= Sqrt(abs( -Delta/(cent[0]*cent[1]*cent[1])));

Ни слова не понял

Цитата:

Изредка косяки. Но в нуль не уходит.

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

Fizykochemik, спасибо за отклик. Вы понимаете правильно. Речь идет об анализе изображений, последовательности фреймов. Есть много контуров таблеток. Эти контуры аппроксимированы эллипсами. Найдены средние значения полуосей a и b. Их я называю нормативным эллипсом. Дуга это непрерывная последовательность пикселей изображения.
Для каждой дуги нужно найти центр тела к которому относится эта дуга. Поскольку нормативный эллипс известен, то если дугу аппроксимировать этим эллипсом, то его центр можно считать центром тела. Поэтому задача в том чтобы по координатам точек дуги и нормативным длинам полуосей a и b найти координаты центра нормативного эллипса для этой конкретной дуги.
Любой эллипс полностью задан если известны его полуоси a и b, его центр Мо и его наклон. Поэтому если в уравнение
$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$ подставить новые значения x и y смещенные на координаты центра xo yo и повернутые на угол $\varphi$
$x = x\cdot cos(\phi) + y\cdot sin(\phi) - x_0$
$y = -x\cdot sin(\phi) + y\cdot cos(\phi) - y_0$
то получим вместо уравнения
$a_{11} \cdot x^2 + a_{22} \cdot y^2 +2\cdot a_{12} \cdot xy+2\cdot a_{13} \cdot x+2 \cdot a_{23} \cdot y+a_{33}=0$ (*)
получим его же в преобразованном виде

Pavia в сообщении #1367820 писал(а):

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

т.е. здесь коэ-ты а11 а12.. выражены через полуоси a и b координаты центра xo yo и угол наклона $\varphi$
Из этого ур-я можно получить угол $\varphi$

Еще раз. Есть точки дуги x y. Подставляем их в ур-е (*), получаем переопределенную систему относительно а11 а22..., решаем ее и находим а11 а22... и из ур-я $a_{12} = \frac{cos(\phi)sin(\phi)}{a^2}-\frac{cos(\phi)sin(\phi)}{b^2}$ находим $\varphi$
Затем из системы уравнений :
$a_{13}=-\frac{x_0cos\varphi}{a^2}+\frac{y_0sin\varphi}{b^2}$
$a_{23}=-\frac{x_0sin\varphi}{a^2}+\frac{y_0cos\varphi}{b^2}$
находим центр $x_0$ $y_0$

Но оказалось что по чисто техническим причинам не удается всегда верно найти коэф-ты а11 а22... В старом варианте программы я находил отношения а11/а33 а22/а33..., поэтому как правило все работало почти всегда, но для эллипса вообще, а не для конкретного с осями а и b
Поэтому в принципе, т.е. с т. зр. математики эта задача уже решена, но технически пока не удается верно посчитать коэф-ты а11 а22..., т.к. как говорит Pavia происходит перегрузка double в строке (**) см ниже

function [a11 a12 a13 a22 a23 a33] = aquadric(x, y)
% входные x, y - координаты точек дуги
x2 = x.*x; y2 = y.*y; xy = x.*y;
one = ones(length(x),1); 
X = [x2 y2 xy x y 1]; Y = one;
R = (X'*X)^(-1)*X'*Y;       %  (**)
a11 = R(1); a22 = R(2); a12 = R(3)/2; 
a13 = R(4)/2; a23 = R(5)/2; 
a33 = R(6) - 1;  
end

Кстати, Pavia, я так и не понял почему следует вычитать единицу a33 = R(6) - 1 , а не просто уравнять а33 = R(6) ?


Emergency, все верно