Дешевый алгоритм для определения угла между векторами

Найти угол между двумя векторами несложно с помощью правила косинуса. Однако, поскольку я программирую для платформы с очень ограниченными ресурсами, я бы хотел избежать таких вычислений, как sqrt и arccos. Даже простые деления следует максимально ограничить.

К счастью, мне не нужен угол как таковой, а нужно только некоторое значение, пропорциональное указанному углу.

Поэтому я ищу дешевый в вычислительном отношении алгоритм для вычисления величины, связанной с углом между двумя векторами. Пока я не нашел чего-то, что подходило бы под все требования, и сам не смог что-то придумать.


person Jeroen    schedule 15.09.2009    source источник
comment
хм: важный вопрос: хранятся ли компоненты векторов в формате с фиксированной или плавающей точкой?   -  person Jason S    schedule 15.09.2009
comment
Ни один. Поскольку рассматриваемые координаты являются пиксельными координатами, они всегда являются целочисленными значениями. Нет необходимости в числах с плавающей / фиксированной точкой. Так что, я думаю, вы могли бы сказать, что они фиксированная точка с множителем 1 :)   -  person Jeroen    schedule 16.09.2009


Ответы (9)


Вы пробовали алгоритм CORDIC? Это общая структура для решения полярных прямоугольных задач с использованием только таблицы сложения / вычитания / битового сдвига +, по существу, выполняющего вращение на углы формы tan -1 (2 -n). Вы можете найти компромисс между точностью и временем выполнения, изменив количество итераций.

В вашем случае возьмите один вектор в качестве фиксированной ссылки и скопируйте другой во временный вектор, который вы поворачиваете с использованием угловых углов по направлению к первому вектору (примерно пополам), пока не достигнете желаемой угловой точности.

(изменить: используйте знак скалярного произведения, чтобы на каждом шаге определять, следует ли вращать вперед или назад. Хотя, если умножения достаточно дешевы, чтобы можно было использовать скалярное произведение, не беспокойтесь о CORDIC, возможно, используйте таблица пар sin / cos для матриц вращения углов / 2 n для решения задачи с делением пополам.)

(edit: Мне нравится предложение Эрика Бейнвилля в комментариях: поверните оба вектора к нулю и отслеживайте разницу углов.)

person Jason S    schedule 15.09.2009
comment
+1. CORDIC - это то, что вы хотите. Возможно, будет быстрее повернуть два вектора одновременно по направлению к оси X, отслеживая разницу углов. - person Eric Bainville; 15.09.2009
comment
... и CORDIC бесплатно предоставит вам норму каждого вектора. - person Eric Bainville; 15.09.2009
comment
Отлично. CORDIC, кажется, сделает то, что мне нужно. Я еще не реализовал это, но все равно отмечу ваш ответ как принятый. Спасибо! - person Jeroen; 15.09.2009

Если вам не нужен фактический евклидов угол, а что-то, что вы можете использовать в качестве основы для сравнения углов, тогда выбор может быть изменен на геометрию такси, потому что вы можете отказаться от тригонометрии и ее медлительности при ПОДДЕРЖАНИИ ТОЧНОСТИ (или, по крайней мере, с очень незначительной потерей точности, см. ниже).

В основных современных браузерных движках коэффициент ускорения составляет от 1,44 до 15,2, а точность почти такая же, как в atan2. Расчет угла ромба в среднем в 5,01 раз быстрее, чем atan2, а с использованием встроенного кода в Firefox 18 ускорение достигает коэффициента 15,2. Сравнение скорости: http://jsperf.com/diamond-angle-vs-atan2/2 < / а>.

Код очень простой:

function DiamondAngle(y, x)
{
    if (y >= 0)
        return (x >= 0 ? y/(x+y) : 1-x/(-x+y)); 
    else
        return (x < 0 ? 2-y/(-x-y) : 3+x/(x-y)); 
}

Приведенный выше код дает угол между 0 и 4, а atan2 дает угол между -PI и PI, как показано в следующей таблице:

введите описание изображения здесь

Обратите внимание, что угол ромба всегда положительный и находится в диапазоне 0-4, тогда как atan2 дает также отрицательные радианы. Так что угол ромба более нормализован. И еще одно замечание: atan2 дает немного более точный результат, потому что длина диапазона составляет 2 * pi (т.е. 6,283185307179586), а в углах ромба - 4. На практике это не очень важно, например. Рад 2.3000000000000001 и 2.3000000000000002 оба находятся в углах ромба 1.4718731421442295, но если мы снизим точность, отбросив один ноль, рад 2.300000000000001 и 2.300000000000002 даст разные углы ромба. Эта «потеря точности» в углах ромба настолько мала, что оказывает значительное влияние только на больших расстояниях. Вы можете поиграть с преобразованиями в http://jsbin.com/bewodonase/1/edit?output (Старая версия: http://jsbin.com/idoyon/1):

введите описание изображения здесь

Приведенного выше кода достаточно для быстрого сравнения углов, но во многих случаях необходимо преобразовать ромбовидный угол в радианы и наоборот. Если вы, например. имеют некоторый допуск в виде углов в радианах, а затем у вас есть 100000-кратный цикл, в котором этот допуск сравнивается с другими углами, неразумно проводить сравнения с использованием atan2. Вместо этого перед зацикливанием вы меняете допуск в радианах на допуск такси (ромбовидные углы) и выполняете внутрицикловые сравнения, используя ромбовидный допуск, и таким образом вам не нужно использовать медленные тригонометрические функции в критических по скорости частях кода (= в петли).

Код, выполняющий это преобразование, таков:

function RadiansToDiamondAngle(rad)
{
  var P = {"x": Math.cos(rad), "y": Math.sin(rad) };
  return DiamondAngle(P.y, P.x);
}  

Как вы заметили, есть cos и sin. Как вы знаете, они медленные, но вам не нужно выполнять преобразование в цикле, но до цикла и ускорение огромно.

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

function DiamondAngleToRadians(dia)
{
  var P = DiamondAngleToPoint(dia);
  return Math.atan2(P.y,P.x);
}

function DiamondAngleToPoint(dia)
{
  return {"x": (dia < 2 ? 1-dia : dia-3), 
  "y": (dia < 3 ? ((dia > 1) ? 2-dia : dia) : dia-4)};
}

Здесь вы используете atan2, что медленно, но идея состоит в том, чтобы использовать это вне любых циклов. Вы не можете преобразовать угол ромба в радианы, просто умножив его на какой-то коэффициент, а вместо этого найдите точку в геометрии такси, в которой угол ромба между этой точкой и положительной осью X является рассматриваемым углом ромба, и преобразовав эту точку в радианы с помощью atan2.

Этого должно быть достаточно для быстрого сравнения углов.

Конечно, есть и другие методы ускорения atan2 (например, CORDIC и таблицы поиска), но AFAIK все они теряют точность и все же могут быть медленнее, чем atan2.

ИСТОРИЯ: Я протестировал несколько методов: скалярные произведения, внутренние произведения, закон косинуса, единичные круги, справочные таблицы и т. Д., Но ничего не было достаточно в случае, когда важны и скорость, и точность. Наконец, я нашел страницу в http://www.freesteel.co.uk/wpblog/2009/06/encoding-2d-angles-without-trigonometry/, который имел желаемые функции и принципы.

Сначала я предположил, что расстояния такси можно также использовать для точного и быстрого сравнения расстояний, потому что большее расстояние в евклидовом формате больше и в такси. Я понял, что, в отличие от евклидовых расстояний, угол между начальной и конечной точкой влияет на расстояние такси. Только длины вертикальных и горизонтальных векторов могут быть легко и быстро преобразованы между евклидовым вектором и такси, но во всех остальных случаях вы должны учитывать угол, и тогда процесс будет слишком медленным (?).

Итак, в качестве вывода я думаю, что в приложениях, критичных к скорости, где есть цикл или рекурсия нескольких сравнений углов и / или расстояний, углы быстрее сравниваются в пространстве такси и расстояния в евклидовом (в квадрате, без использования sqrt) пространстве.

person Timo Kähkönen    schedule 03.02.2013
comment
Это именно то, что мне было нужно - спасибо. Но я назвал параметры dX и dY, поскольку они на самом деле дельты осей, а не ординаты точек. Думаю, так понятнее. - person Michael Scheper; 11.02.2015

В те времена, когда было несколько килобайт оперативной памяти и машины с ограниченными математическими возможностями, я использовал таблицы поиска и линейную интерполяцию. Основная идея проста: создайте массив с таким разрешением, которое вам нужно (большее количество элементов уменьшает ошибку, создаваемую интерполяцией). Затем выполните интерполяцию между значениями поиска.

Вот пример в процессе обработки (исходная неработающая ссылка).

Вы также можете сделать это с помощью других триггерных функций. На процессоре 6502 это позволило вычислить полную трехмерную каркасную графику с увеличением скорости на порядок.

person Godeke    schedule 15.09.2009
comment
+1: интерполяция в таблицы - быстрая и маленькая. Не так точно, но часто достаточно точно. 256 различных значений косинуса часто бывает достаточно. - person S.Lott; 15.09.2009
comment
Спасибо за предложение, но это все равно избавит от arccos. для определения остальных элементов правила косинуса в этом случае (длины всех сторон треугольника) все равно потребуется несколько sqrt. - person Jeroen; 15.09.2009
comment
Меньшего может быть достаточно. Если вы хотите сделать небольшой дополнительный расчет, вам нужны только значения между нулем и половиной пи / 90 градусов, и, конечно же, вам нужно определить необходимую детализацию и точность. Трудно сказать что-либо определенное о деталях без учета конкретного случая. - person David Thornley; 15.09.2009
comment
SQRT легко заменить на цикл аппроксимации, который сходится быстро и с большой точностью. И в этом нет деления, просто умножение. - person S.Lott; 15.09.2009
comment
CORDIC - отличное решение, если у вас нет множителя, но если он у вас есть в вашей системе, я подозреваю, что интерполяция будет быстрее (понимая, что вы можете построить такую ​​таблицу для SQRT также и в своем домене). CORDIC является итеративным, в то время как поиск / интерполяция имеет фиксированную стоимость. Однако я специализируюсь на 3D на слабом оборудовании, так что, возможно, ваше оборудование сочтет CORDIC приемлемым. - person Godeke; 15.09.2009

Здесь, на SO, у меня все еще нет права комментировать (хотя у меня есть на math.se), так что на самом деле это ответ на сообщение Тимо об углах ромба.

Сама концепция углов ромба, основанная на норме L1, является наиболее интересной, и если бы это было просто сравнение того, какой вектор имеет большее / меньшее значение по сравнению с положительной оси X было бы достаточно. Однако OP упомянул угол между двумя общими векторами, и я предполагаю, что OP хочет сравнить его с некоторым допуском для определения состояния гладкости / угла или чего-то подобного, но, к сожалению, кажется, что только с формулами, предоставленными на jsperf.com или freesteel.co.uk (ссылки выше) кажется, что это невозможно сделать, используя ромбовидные углы.

Обратите внимание на следующий результат моей реализации формул в Asymptote:

Vectors : 50,20 and -40,40
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.21429
Convert that to degrees       : 105.255
Rotate same vectors by 30 deg.
Vectors : 33.3013,42.3205 and -54.641,14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.22904
Convert that to degrees       : 106.546
Rotate same vectors by 30 deg.
Vectors : 7.67949,53.3013 and -54.641,-14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.33726
Convert that to degrees       : 116.971

Итак, дело в том, что вы не можете сделать алмаз (альфа) -алмаз (бета) и сравнить его с некоторым допуском, в отличие от вывода atan2. Если все, что вы хотите сделать, это алмаз (альфа)> бриллиант (бета), то я полагаю, что с бриллиантом все в порядке.

person jamadagni    schedule 03.06.2013
comment
Да ты прав. Углы ромба можно только сравнивать, но нельзя складывать, вычитать, умножать или делить. Единственным исключением являются градусы 0, 45, 90, 135, 180, 225, 270, 315, 360. Вы можете безопасно переносить эти углы между такси и евклидовыми пространствами, используя разделение / умножение. Перенос любого другого угла вызывает ошибку, и эта ошибка варьируется от 0 до 4,074568596222775 градусов. Я попытался выяснить логику этой разницы и можно ли ее учесть для повышения точности передачи из космоса в космос. - person Timo Kähkönen; 30.08.2015

Перекрестное произведение пропорционально углу между двумя векторами, и когда векторы нормализованы, а угол мал, перекрестное произведение очень близко к фактическому углу в радианах из-за приближения малого угла.

конкретно:

I1Q2-I2Q1 пропорционален углу между I1Q1 и I2Q2.

person Dan Boschen    schedule 18.12.2010

Решение было бы тривиальным, если бы векторы были определены / сохранены с использованием полярных координат вместо декартовых координат (или «а также» с использованием декартовых координат).

person ChrisW    schedule 15.09.2009
comment
К сожалению, это не так. Единственный способ - преобразовать координаты в полярные в моем коде, но это приводит к той же проблеме, которую я пытаюсь решить в первую очередь. - person Jeroen; 15.09.2009

скалярное произведение двух векторов (x1, y1) и (x2, y2) равно

x1 * x2 + y1 * y2 

и эквивалентен произведению длин двух векторов на косинус угла между ними.

Итак, если вы сначала нормализуете два вектора (разделите координаты на длину)

Where length of V1 L1 = sqrt(x1^2 + y1^2),  
  and length of V2 L2 = sqrt(x2^2 + y2^2),

Тогда нормализованные векторы равны

(x1/L1, y1/L1),  and (x2/L2, y2/L2),  

И точечный продукт нормализованных векторов (который совпадает с косинусом угла между векторами) будет

 (x1*x2 + y1*y2)
 -----------------
     (L1*L2)

конечно, это может быть так же сложно в вычислительном отношении, как вычисление косинуса

person Charles Bretana    schedule 15.09.2009
comment
OP знает, что для скалярного произведения требуется arccos, чтобы получить угол - person Jason S; 15.09.2009
comment
@Jason, да, но точечный продукт сам по себе является величиной, которая связана с углом, который у вас нет; t обязательно нужно преобразовать его в радианы или градусы, чтобы использовать его ... - person Charles Bretana; 16.09.2009

если вам нужно вычислить квадратный корень, рассмотрите возможность использования хака invsqrt.

acos((x1*x2 + y1*y2) * invsqrt((x1*x1+y1*y1)*(x2*x2+y2*y2)));
person neoneye    schedule 15.09.2009
comment
Если нет проблем с переполнением, я бы объединил invsqrt в один вызов invsqrt ((x1 * x1 + y1 * y1) * (x2 * x2 + y2 * y2)); лучше один раз запустить сложный алгоритм, чем дважды. - person Jason S; 15.09.2009

В вашем случае может работать точечный продукт. Он не пропорционален углу, а «связан».

person Malte Clasen    schedule 15.09.2009
comment
OP знает, что для скалярного произведения требуется arccos, чтобы получить угол - person Jason S; 15.09.2009
comment
Или квадратные корни для нормализации длины. И то, и другое ОП пытается избежать. - person David Thornley; 15.09.2009