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

Форум для обсуждения вопросов математики

Модератор: Admin

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

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

Сообщение haker_fox » Пт янв 21, 2011 9:45 am

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

Kitonum
Сообщения: 2010
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

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

Сообщение Kitonum » Пт янв 21, 2011 12:32 pm

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=∞. Все эти шаги легко программируются. Чем гуще расположены точки - тем точнее получится результат!

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

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

Сообщение haker_fox » Пн янв 24, 2011 6:43 am

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 - радиус окружности или радиус кривизны кривой в искомой точке.
Решая полученную систему, мы находим координаты центра окружности (они нам не нужны) и радиус (он нам очень нужен)))

Kitonum
Сообщения: 2010
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

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

Сообщение Kitonum » Пн янв 24, 2011 11:10 am

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

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

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

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

Сообщение haker_fox » Вт янв 25, 2011 8:09 am

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

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

muerte
Сообщения: 1
Зарегистрирован: Вт окт 28, 2008 11:49 am
Контактная информация:

Сообщение muerte » Вт янв 25, 2011 1:32 pm

Здравствуйте! Получилось так, что сейчас передо мной стоит примерно такая же задача. Поиски в интернете натолкнули меня на эту тему. Позвольте рассказать еще один способ решения (ну или расширение к способу, предложенному 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 за идею! Я уж было собиралась делать через аналитику, благо в моей задаче время есть. Однако производных все же хотелось избежать. Хотя вероятно что реализую оба варианта и посмотрю, какой точнее :)

hirnyk
Сообщения: 438
Зарегистрирован: Пт апр 08, 2005 1:41 pm

radius

Сообщение hirnyk » Вт янв 25, 2011 10:14 pm

> 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)

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

Re: radius

Сообщение haker_fox » Ср янв 26, 2011 10:09 am

Спасибо, коллеги!
Я использовал формулы описанной вокруг треугольника окружности, на которую меня направил уважаемый 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

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

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

Сообщение haker_fox » Сб фев 05, 2011 11:54 am

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

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Сообщение алексей_алексей » Сб фев 05, 2011 1:01 pm

Не очень понятно, траектория уже имеется или создаётся по ходу. Если по ходу, то это и есть траектория, и какие тут могут быть отклонения? А если она уже есть, то просто контролируется положение текущей точки “слева” и “справа” от ломаной, по подстановке в уравнение направленной прямой, проходящей через две ближайшие точки ломаной. Если траектория замкнутая, то однозначно, в противном случае выбираете критерий… Если всё же траектория строится по ходу, и Вы успеваете приблизить её уравнением для предыдущих точек, то точно также подставляете в это уравнение новые координаты точки и смотрите на знак и величину отклонения… ,наверное.

IVVA
Сообщения: 1036
Зарегистрирован: Вт апр 05, 2005 6:44 pm

Сообщение IVVA » Сб фев 05, 2011 1:05 pm

Можно рассмотреть три вектора и затем - векторные
произведения отклоняющихся веторов от направления
с веторм направления.И по знаку векторного произведения определить - в какую сторону отклонилась точка от установленного направления.

Kitonum
Сообщения: 2010
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Сообщение Kitonum » Вс фев 06, 2011 12:04 am

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(α)| превышен, соответствующую точку С отбрасываете и вместо неё пишите "+" или "-" в соответствии с приведённым критерием, и так далее.

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

Сообщение haker_fox » Вс фев 06, 2011 2:25 pm

алексей_алексей писал(а):Не очень понятно, траектория уже имеется или создаётся по ходу.

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

haker_fox
Сообщения: 7
Зарегистрирован: Пт янв 21, 2011 9:27 am
Контактная информация:

Сообщение haker_fox » Вс фев 06, 2011 2:27 pm

Уважаемые алексей_алексей, IVVA и Kitonum! Искренне благодарю Вас за ответы! Буду прорабатывать теорию.

133880
Сообщения: 1
Зарегистрирован: Вт окт 17, 2017 11:04 am

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

Сообщение 133880 » Вт окт 17, 2017 11:52 am

Столкнулся с такой проблемой. Поиск готовых решений ничего не дал. Может кому поможет:
% поиск радиуса кривизны
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 - искомый радиус
% Радиус получилось найти только переместив центр окружности в начало координат.