2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


Посмотреть правила форума



Начать новую тему Ответить на тему На страницу Пред.  1, 2, 3, 4
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение11.02.2019, 12:45 


18/05/10
75
Алексей К. в сообщении #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. Теперь задача состоит в том чтобы аппроксимировать ее заданным эллипсом.

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение11.02.2019, 14:06 


18/05/10
75
Эту аппроксимацию осуществим с помощью алгоритма предложенного Алексеем (немного преобразую его)
$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 (желтую)
Изображение
новая аппроксимация желтой дуги (красная линия). Уже по рис видно на сколько ее центр близок с началу координат
Изображение
Новая, красная итерация в исходных координатах
Изображение

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение11.02.2019, 17:06 


18/05/10
75
Нулевым приближением полагались два параметра найденных при первичной аппроксимации дуги произвольным эллипсом: угол наклона его главной оси и его центр Мо1. Я изменил это условие и нашел точку пересечения (Р) оси эллипса с его дугой или с ним самим и затем по оси эллипса отложит отрезок равный а = 79. Т.о. получил новый центр Мо.
Изображение
В этом случае требуются всего одна-две итерации. Вот первая:
Изображение
А вот вторая желтая, почти совпадающая с предыдущей
Изображение

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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение12.02.2019, 12:55 


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

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 )

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение12.02.2019, 13:30 


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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение17.02.2019, 22:10 


29/09/06
4552
Алексей К. в сообщении #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$). Неча и надеяться на сходимость кого-то к кому-то.

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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение18.02.2019, 18:13 


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

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

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

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

Изображение

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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение19.02.2019, 13:00 


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

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

$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$ и т.д. ?

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение20.02.2019, 17:58 


18/05/10
75
Алексей К. в сообщении #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 малые смещения центра, то за несколько итераций можно придти к почти верному центру (у меня это уже заработало). В этом деле радует то что не надо перетаскивать все точки дуги к центру координат. Затем когда первичный центр будет найден можно добавить аналогично смещению центра смещение угла. Раскрою это позднее, видимо завтра или на след неделе представлю результаты.

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение27.02.2019, 14:54 


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

Повернем исходный эллипс $\\\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$

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение27.02.2019, 16:29 


18/05/10
75
Рассмотрим случай (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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение27.02.2019, 20:35 


18/05/10
75
Случей 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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение28.02.2019, 13:02 


18/05/10
75
К последнему уравнению случая 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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение28.02.2019, 14:25 


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

написано: $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)$

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

 Профиль  
                  
 
 Re: Насильственная аппроксимация дуги заданным эллипсом
Сообщение03.03.2019, 17:12 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 60 ]  На страницу Пред.  1, 2, 3, 4

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group