fixfix
2014 dxdy logo

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

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


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


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



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


29/09/06
4552
Я додумаю в вскр. (щас типа на работу надолго). Пока зафиксирую кое-что из памяти.
Алексей К. в сообщении #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$ --- наверняка эксцентриситет и фокальный параметр).

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


31/10/08
1244
Tur
Tur в сообщении #1367777 писал(а):
так что скажете по моему заходу? Нельзя полагать равным единице ни один из коэф-тов: а11 а22 а12?

:facepalm:
Выпишем уравнение эллипса в параметрической форме:
$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}$ легко найти угол вращения.
И далее вычислить центр эллипса.

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


29/09/06
4552
Pavia в сообщении #1367820 писал(а):
Возведём в квадрат и подставим в основное тригометрическое тожество.

Кнут, cepesh, Корн, Бронштейн, gris, Семендяев --- эти и многие другие люди приучили меня, что основное тригометрическое тожество (ОТТ) есть $$\cos^2 t +\sin^2 t=1.$$ С Вашим, Pavia, вариантом ОТТ, типа $$cos(t)^2 + sin(t)^2=1,$$я согласиться никак не могу.

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


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

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


18/05/10
75
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 = ...

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


31/10/08
1244
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

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


18/05/10
75
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

(Оффтоп)


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

(Оффтоп)


x y дуги

(Оффтоп)


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


31/10/08
1244
Код:
(X'*X)
- здесь переполнение для чисел double(вернее исчерпание разрядов мантиссы). Поэтому и срывается.

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


18/05/10
75
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}$
Может в чем то другом переполнение или оно не видно?

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

(Оффтоп)



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

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

(Оффтоп)


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


31/10/08
1244
Tur
Вот и у меня так же.
Изображение

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


18/05/10
75
Pavia в сообщении #1368433 писал(а):
Tur
Вот и у меня так же.
Изображение

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

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

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


31/10/08
1244
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}$ Вычисляется.
Изображение

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


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

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

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

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

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


07/03/16

3167
Fizykochemik в сообщении #1368546 писал(а):
Правильно ли я понимаю, что Вам нужен эллипс

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

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


18/05/10
75
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])));
Ни слова не понял :oops:
Цитата:
Изредка косяки. Но в нуль не уходит.
В старом варианте косяки один на тысячу и это неприемлимо. Ведь каждый день тесятки тысяч таблеток, сотни тысяч дуг и это только по одному каналу.

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, все верно

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

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



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

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


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

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