Страница 1 из 1

Определение радиуса кривизны по координатам траектории

Добавлено: Пт янв 21, 2011 9:45 am
haker_fox
Здравствуйте, уважаемые коллеги!
Передо мной стоит следующая задача.
Имеется набор координат (xi, yi), описывающий некую траекторию. Координаты заданы только в виде пар чисел, аналитического выражения нет. Траектория может быть абсолютно любой. Ну, к примеру, это траектория гуляющей собаки :D
Необходимо по этим точкам вычислить радиус кривизны в каждой точки трактории. Вот тут http://mathworld.wolfram.com/RadiusofCurvature.html была найдена формула, которая позволяет вычислить радиус кривизны для функции вида y = f(x), а это наш случай. Основная проблема, что требуется вычисление второй производной, причем численно. А это погрешности, причем в некоторых случаях могут быть солидные. Приращение координат неравномерное.
В общем, возникает вопрос, можно ли вычислить радиус кривизны с приемлимой точностью?
Да, и еще, алгоритм будет реализован на языке Си++, следовательно использование математических пакетов исключается.
Заранее благодарен за любую помощь!
Извините, если что-то важное не сообщил.

Re: Определение радиуса кривизны по координатам траектории

Добавлено: Пт янв 21, 2011 12:32 pm
Kitonum
haker_fox писал(а):Здравствуйте, уважаемые коллеги!
Передо мной стоит следующая задача.
Имеется набор координат (xi, yi), описывающий некую траекторию. Координаты заданы только в виде пар чисел, аналитического выражения нет. Траектория может быть абсолютно любой. Ну, к примеру, это траектория гуляющей собаки :D
Необходимо по этим точкам вычислить радиус кривизны в каждой точки трактории. Вот тут http://mathworld.wolfram.com/RadiusofCurvature.html была найдена формула, которая позволяет вычислить радиус кривизны для функции вида y = f(x), а это наш случай. Основная проблема, что требуется вычисление второй производной, причем численно. А это погрешности, причем в некоторых случаях могут быть солидные. Приращение координат неравномерное.
В общем, возникает вопрос, можно ли вычислить радиус кривизны с приемлимой точностью?
Да, и еще, алгоритм будет реализован на языке Си++, следовательно использование математических пакетов исключается.
Заранее благодарен за любую помощь!
Извините, если что-то важное не сообщил.

Те формулы, на которые Вы ссылаетесь, Вам вряд ли подойдут! Дело в том, что Ваша траектория, даже приближённо, не может быть задана параметрическими уравнениями, т.к. для этого нужны тройки чисел (ti,xi,yi), где первая координата - время, а вторая и третья - координаты точки на плоскости, соответствующие моменту времени ti. Нельзя задать её и явным уравнением y = f(x), т.к. для этого все все xi должны быть различными, а в Вашей задаче такого может не быть. Например, собака может бегать по спирали и т.п.
Думаю, наиболее простой и ясный подход к решению Вашей задачи следующий. Как известно, радиусом кривизны кривой в данной точке называется радиус окружности, проходящей через 3 бесконечно близкие точки этой кривой. Разумеется, в строгом определении нужно использовать предельный переход, когда одна точка (в которой вычисляется радиус кривизны) фиксируется), а 2 другие (слева и справа от неё) к ней неограниченно приближаются. Этим методом, например, можете подсчитать радиус кривизны параболы в её вершине (без всяких производных!).
Применительно к Вашей задаче, речь не идёт ни о каком предельном переходе, т.к. у Вас просто дискретный набор точек на плоскости. Просто перебираете эти точки, начиная со 2 точки, а в качестве радиуса кривизны траектории в данной точке берёте радиус окружности, проходящей через эту точку, предыдущую и последующую. Если эти 3 точки лежат на одной прямой, то берёте R=∞. Все эти шаги легко программируются. Чем гуще расположены точки - тем точнее получится результат!

Re: Определение радиуса кривизны по координатам траектории

Добавлено: Пн янв 24, 2011 6:43 am
haker_fox
Kitonum писал(а):
haker_fox писал(а):Здравствуйте, уважаемые коллеги!
Передо мной стоит следующая задача.
Имеется набор координат (xi, yi), описывающий некую траекторию. Координаты заданы только в виде пар чисел, аналитического выражения нет. Траектория может быть абсолютно любой. Ну, к примеру, это траектория гуляющей собаки :D
Необходимо по этим точкам вычислить радиус кривизны в каждой точки трактории. Вот тут http://mathworld.wolfram.com/RadiusofCurvature.html была найдена формула, которая позволяет вычислить радиус кривизны для функции вида y = f(x), а это наш случай. Основная проблема, что требуется вычисление второй производной, причем численно. А это погрешности, причем в некоторых случаях могут быть солидные. Приращение координат неравномерное.
В общем, возникает вопрос, можно ли вычислить радиус кривизны с приемлимой точностью?
Да, и еще, алгоритм будет реализован на языке Си++, следовательно использование математических пакетов исключается.
Заранее благодарен за любую помощь!
Извините, если что-то важное не сообщил.

Те формулы, на которые Вы ссылаетесь, Вам вряд ли подойдут! Дело в том, что Ваша траектория, даже приближённо, не может быть задана параметрическими уравнениями, т.к. для этого нужны тройки чисел (ti,xi,yi), где первая координата - время, а вторая и третья - координаты точки на плоскости, соответствующие моменту времени ti. Нельзя задать её и явным уравнением y = f(x), т.к. для этого все все xi должны быть различными, а в Вашей задаче такого может не быть. Например, собака может бегать по спирали и т.п.
Думаю, наиболее простой и ясный подход к решению Вашей задачи следующий. Как известно, радиусом кривизны кривой в данной точке называется радиус окружности, проходящей через 3 бесконечно близкие точки этой кривой. Разумеется, в строгом определении нужно использовать предельный переход, когда одна точка (в которой вычисляется радиус кривизны) фиксируется), а 2 другие (слева и справа от неё) к ней неограниченно приближаются. Этим методом, например, можете подсчитать радиус кривизны параболы в её вершине (без всяких производных!).
Применительно к Вашей задаче, речь не идёт ни о каком предельном переходе, т.к. у Вас просто дискретный набор точек на плоскости. Просто перебираете эти точки, начиная со 2 точки, а в качестве радиуса кривизны траектории в данной точке берёте радиус окружности, проходящей через эту точку, предыдущую и последующую. Если эти 3 точки лежат на одной прямой, то берёте R=∞. Все эти шаги легко программируются. Чем гуще расположены точки - тем точнее получится результат!

Благодарю Вас за ответ!
Он очень нам помог! Возможность работы без производных сильно облегчает жизнь)
Можно небольшое уточнение?
Вы предлагаете составить систему трех уравнений окружности вида
(xi - a)^2 + (yi - b)^2 = r^2,
где xi, yi - координаты трех точек;
a, b - координаты центра окружности;
r - радиус окружности или радиус кривизны кривой в искомой точке.
Решая полученную систему, мы находим координаты центра окружности (они нам не нужны) и радиус (он нам очень нужен)))

Re: Определение радиуса кривизны по координатам траектории

Добавлено: Пн янв 24, 2011 11:10 am
Kitonum
haker_fox писал(а):Можно небольшое уточнение?
Вы предлагаете составить систему трех уравнений окружности вида
(xi - a)^2 + (yi - b)^2 = r^2,
где xi, yi - координаты трех точек;
a, b - координаты центра окружности;
r - радиус окружности или радиус кривизны кривой в искомой точке.
Решая полученную систему, мы находим координаты центра окружности (они нам не нужны) и радиус (он нам очень нужен)))

Да, можно и так, с помощью системы. А ещё проще прямым вычислением, используя известные формулы из аналитической геометрии! Используйте то, что центр описанной вокруг треугольника окружности можно получить как точку пересечения серединных перпендикуляров к любым двум сторонам треугольника. Так легко получить общую формулу для R для любой тройки точек, не лежащих на одной прямой, через их координаты, а затем уже использовать полученную формулу при составлении программы.

Re: Определение радиуса кривизны по координатам траектории

Добавлено: Вт янв 25, 2011 8:09 am
haker_fox
Kitonum писал(а):Да, можно и так, с помощью системы. А ещё проще прямым вычислением, используя известные формулы из аналитической геометрии! Используйте то, что центр описанной вокруг треугольника окружности можно получить как точку пересечения серединных перпендикуляров к любым двум сторонам треугольника. Так легко получить общую формулу для R для любой тройки точек, не лежащих на одной прямой, через их координаты, а затем уже использовать полученную формулу при составлении программы.

Уважаемый, Kitonum! Еще раз большое спасибо!
А то я начал численно решать систему нелинейных уравнений методом Ньютона и понял, что проблем еще больше стало: нужно выбрать начальные нулевые условия правильно, иначе можно и не решить систему; метод позволяет найти только каждый корень только одного знака(
Попробую использовать геометрию, хоть я в ней и не силен)
Еще раз спасибо!!!

Добавлено: Вт янв 25, 2011 1:32 pm
muerte
Здравствуйте! Получилось так, что сейчас передо мной стоит примерно такая же задача. Поиски в интернете натолкнули меня на эту тему. Позвольте рассказать еще один способ решения (ну или расширение к способу, предложенному Kitonum)
Есть замечательная формула окружности, проходящей через 3 точки с координатами (x_i, y_i), паст из матлаб-скрипта:
m = [x^2 + y^2, x, y, 1;...
x1^2 + y1^2, x1, y1, 1;...
x2^2 + y2^2, x2, y2, 1;...
x3^2 + y3^2, x3, y3, 1];
где m = 0;
Раскрывая данный определитель и приводя его искуственно к стандартному виду:

(x-a)^2 + (y-b)^2 = R^2

Получаем R. Единственное, надо не забывать делать перед всем этим проверку на непринадлежность трех берущихся точек одной прямой, иначе там деление на ноль может получиться.
Спасибо Kitonum за идею! Я уж было собиралась делать через аналитику, благо в моей задаче время есть. Однако производных все же хотелось избежать. Хотя вероятно что реализую оба варианта и посмотрю, какой точнее :)

radius

Добавлено: Вт янв 25, 2011 10:14 pm
hirnyk
> with(geometry);
> point(A, x1, y1);

A
> point(B, x2, y2);

B
> point(C, x3, y3);

C
> assume(x1*y2-x1*y3+x2*y3-x2*y1+x3*y1-x3*y2 <> 0);
> circle(c, [A, B, C]);


c


> radius(c);

((x1-(1/2)*(x1^2*y2-x1^2*y3+y1^2*y2-y1^2*y3-y1*x2^2-y1*y2^2+y1*x3^2+y1*y3^2+y2^2*y3-y2*x3^2-y2*y3^2+x2^2*y3)/(x1*y2-x1*y3+x2*y3-x2*y1+x3*y1-x3*y2))^2 +(y1+(1/2)*(-x1*(x2^2+y2^2)+x1*(x3^2+y3^2)+x1^2*(-x3+x2)+y1^2*(-x3+x2)-x2*(x3^2+y3^2)+x3*(x2^2+y2^2))/(x1*y2-x1*y3+x2*y3-x2*y1+x3*y1-x3*y2))^2)^(1/2)

Re: radius

Добавлено: Ср янв 26, 2011 10:09 am
haker_fox
Спасибо, коллеги!
Я использовал формулы описанной вокруг треугольника окружности, на которую меня направил уважаемый Kitonum, ведь она же проходит через три заданные точки. Поскольку основное мое направление не математическое, пришлось потратить час на спокойное чтение википедии)
В итоге, я воспользовался формулами для нахождения длин векторов (сторон треугольника) по заданным координатам, формулой Герона для определения площади треугольника и формулой для радиуса описанной окружности. Получилось не совсем элегантное, но решение, не требующее итераций:
% Находим длины сторон треугольника
a = sqrt((x[i] - x[i+1])^2 + (y[i] - y[i+1])^)
% и так далее оставшиеся стороны
% По формуле Герона находим площать треугольника
p = (a + b + c) / 2
S = sqrt(p(p-a)(p-b)(p-c))
% Вычисляем радиус
R = a * b * c / 4 / S

Вот так) Наверно можно упростить, но для меня и это уже результат!

Добавлено: Сб фев 05, 2011 11:54 am
haker_fox
И снова здравствуйте, друзья!
В продолжение темы возникла следующая непреодолимая для меня задача. Есть траектория, по которой движется объект. Траектория - произвольная, объект - все та же старая добрая собака) Объект перемещается по траектории, но в силу различных причин (собака противоположной породы, порыв ветра и т.п.) он может отклоняться в одну или другую сторону с обязательным возвратом на траекторию, т.е. возможны кратковременные сходы и возвраты с/на траекторию. Необходимо определять знак этого отклонения. Предположим, что уход с траектории влево - это "плюс", а уход вправо - это "минус". Если такие отклонения анализирует человек, он легко определит знак отклонения, пройдя пальцем, взглядом или карандашом по траектории. Другое дело - компьютер. Какой бы алгоритм изобрести или найти готовый, позволяющий выполнить поставленную задачу. Еще раз напомню, что траектория задана в виде дискретных точек.

Добавлено: Сб фев 05, 2011 1:01 pm
алексей_алексей
Не очень понятно, траектория уже имеется или создаётся по ходу. Если по ходу, то это и есть траектория, и какие тут могут быть отклонения? А если она уже есть, то просто контролируется положение текущей точки “слева” и “справа” от ломаной, по подстановке в уравнение направленной прямой, проходящей через две ближайшие точки ломаной. Если траектория замкнутая, то однозначно, в противном случае выбираете критерий… Если всё же траектория строится по ходу, и Вы успеваете приблизить её уравнением для предыдущих точек, то точно также подставляете в это уравнение новые координаты точки и смотрите на знак и величину отклонения… ,наверное.

Добавлено: Сб фев 05, 2011 1:05 pm
IVVA
Можно рассмотреть три вектора и затем - векторные
произведения отклоняющихся веторов от направления
с веторм направления.И по знаку векторного произведения определить - в какую сторону отклонилась точка от установленного направления.

Добавлено: Вс фев 06, 2011 12:04 am
Kitonum
haker_fox писал(а):И снова здравствуйте, друзья!
В продолжение темы возникла следующая непреодолимая для меня задача. Есть траектория, по которой движется объект. Траектория - произвольная, объект - все та же старая добрая собака) Объект перемещается по траектории, но в силу различных причин (собака противоположной породы, порыв ветра и т.п.) он может отклоняться в одну или другую сторону с обязательным возвратом на траекторию, т.е. возможны кратковременные сходы и возвраты с/на траекторию. Необходимо определять знак этого отклонения. Предположим, что уход с траектории влево - это "плюс", а уход вправо - это "минус". Если такие отклонения анализирует человек, он легко определит знак отклонения, пройдя пальцем, взглядом или карандашом по траектории. Другое дело - компьютер. Какой бы алгоритм изобрести или найти готовый, позволяющий выполнить поставленную задачу. Еще раз напомню, что траектория задана в виде дискретных точек.

1) Воспользуемся следующим известным фактом. Если на плоскости заданы 2 вектора (для наших целей считаем, что конец первого совпадает с началом второго) своими координатами AB={x1,y1} и BC={x2,y2}, то синус угла между ними вычисляется по формуле sin(α)=det(matrix([[x1,y1],[x2,y2]]))/(|AB|⋅|BC|), где det означает определитель соответствующей матрицы 2 порядка. Знак определителя будет "+", если кротчайший поворот 1 вектора до совмещением с направлением 2 вектора происходит против часовой стрелки (точка С относительно направления АВ сдвинута влево), и знак "-" в противном случае. Это всё следует из теории векторного произведения.

2) Теперь введём числовой критерий того, что точка С сошла с траектории. Будем считать что это произошло, если sin(α), вычисленный выше, по модулю больше наперёд заданного положительного числа, например 0.7. Это число Вы должны ввести сами, исходя из частоты расположения точек на траектории. Если траектория предполагается гладкой, то при неограниченном сближении точек sin(α) должен стремиться к 0.

3) Теперь можете писать свою программу. Просматриваете последовательно список точек и как только |sin(α)| превышен, соответствующую точку С отбрасываете и вместо неё пишите "+" или "-" в соответствии с приведённым критерием, и так далее.

Добавлено: Вс фев 06, 2011 2:25 pm
haker_fox
алексей_алексей писал(а):Не очень понятно, траектория уже имеется или создаётся по ходу.

Траектория имеется, она задана заранее в виде дискретных точек) Когда траектория подготовлена, по ней запускается объект. Он может во время движения отклоняться влево или вправо.

Добавлено: Вс фев 06, 2011 2:27 pm
haker_fox
Уважаемые алексей_алексей, IVVA и Kitonum! Искренне благодарю Вас за ответы! Буду прорабатывать теорию.

Re: Определение радиуса кривизны по координатам траектории

Добавлено: Вт окт 17, 2017 11:52 am
133880
Столкнулся с такой проблемой. Поиск готовых решений ничего не дал. Может кому поможет:
% поиск радиуса кривизны
P=max(size(X));
m0=ones(P,1);
H=[m0 2*X' 2*Y']; % X и Y - массивы координат траектории
Z1=X.*X+Y.*Y; Z1=Z1';
K=inv(H'*H)*H';
Xoc=K*Z1;
X=X-Xoc(2); % Xoc(2) - координата центра окружности
Y=Y-Xoc(3); % Xoc(3) - координата центра окружности
H=[m0 2*X' 2*Y'];
Z1=X.*X+Y.*Y; Z1=Z1';
K=inv(H'*H)*H';
Xoc=K*Z1;
R=sqrt(Xoc(1)); %R - искомый радиус
% Радиус получилось найти только переместив центр окружности в начало координат.