nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Ликбез по кватернионам, часть 5: практическая реализация поворота


Поздравим с днём рождения Виктора Чапаева victor_chapaev, который сподвиг меня написать данный "опус".

Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты

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

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

Чтобы опять не зависнуть на пару месяцев, буду писать по маленьким кусочкам. Сегодня будет уникальный материал, который вы вряд ли найдёте где-то ещё: в библиотеках раскиданы ссылки на форумы и блоги, где товарищи выводили эффективные формулы, но все эти ссылки уже мертвы. Чтобы не работать с формулами, оставленными мудрыми предками (или не воспринимать их как магические заклинания, которые надо просто выучить), автор вывел их заново и убедился в их корректности.

Реализация поворота вектора с помощью кватерниона
У нас есть кватернион Λ с единичной нормой:

и некоторый вектор , который мы записываем как кватернион с нулевой скалярной (действительной) частью:

Рассмотрим, как повернуть вектор наиболее эффективным способом.

«По определению»
Как мы уже знаем, поворот вектора можно произвести по формуле (3.2):
(6.1)


Если не написать специализированной функции, учитывающей, что у вектора «нет скалярной части», это потребует 32 умножения и 24 сложения, тогда как обычное умножение матрицы 3х3 на вектор требует 9 умножений и 6 сложений. Далее в тексте, для удобочитаемости будем писать в виде: 32m + 24a, что означает «32 умножения и 24 сложения».
Можно слегка упростить выражения, учитывая нулевой скалярный компонент векторов и . Сначала вычислим промежуточный кватернион:







(6.2)


Выпишем выражения для компонентов кватерниона в явном виде:




(6.3)


Убрав умножение на нулевую скалярную часть, мы снижаем количество операций до 12m+8a. Далее, вычислим итоговый вектор, учитывая, что его скалярная часть заведомо равна нулю:






(6.4)


Выпишем выражения для компонентов вектора в явном виде:


(6.5)


Здесь мы обошлись 12m+9a. Итого, на поворот вектора требуется 24 умножения и 17 сложений – чуть лучше, чем до упрощения, но всё равно значительно более затратно, чем при умножении вектора на матрицу.

Магическая формула
Следующие формулы применены в нескольких библиотеках для работы с кватернионами, но ссылки на форум, где программист приводил вывод этих формул – давным-давно «мертвы». Для начала, сами формулы:

,

(6.6)


Как видно, используется 2 векторных умножения, 2 умножения вектора на скаляр (один из них – и вовсе число 2) и два сложения векторов. Векторное умножение требует 6m+3a, умножение на скаляр – 3m, сложение векторов – 3a.
В общей сложности, поворот вектора по (6.6) требует 18 умножений и 12 сложений, что хуже обычного матричного умножения ровно вдвое, во всяком случае, по числу операций. С другой стороны, 4 компоненты кватерниона замечательно умещаются в один регистр xmm, и при хорошей реализации на SSE полное время выполнения (включая передачу аргументов функции и возврат результата) может оказаться даже ниже.
Попробуем вывести (6.6) геометрически. Такое доказательство будет наиболее наглядным.



Мы хотим повернуть вектор вокруг оси на угол φ. Угол φ «закодирован» внутри кватерниона:


Мы можем разложить вектор на компоненту, параллельную и перпендикулярную (см. рисунок):



При повороте, параллельная компонента останется на том же месте, тогда как перпендикулярная повернётся, но останется в плоскости p, перпендикулярной . После поворота, этот вектор можно будет выразить следующим образом:



Вектор перпендикулярен как , так и , поэтому он должен быть сонаправлен их векторному произведению:

(мы можем вписать в векторное произведение полный вектор , поскольку параллельная компонента даст ноль)

Его длина составит


Мы хотим, чтобы , т.е вектор описывал бы круговое движение, сохраняющее его длину. Отсюда найдём:


Также нам необходимо выразить через исходные данные. Наиболее простой способ – через ещё одно векторное произведение:


Снова найдём длину:


Отсюда,


Наконец, выразим итоговый вектор:







Вводим вектор :


и получаем формулу (6.6).

Алгебраическое доказательство.
Ту же формулу можно вывести непосредственно из (6.1). Сначала выразим промежуточный кватернион:


Далее, найдём итоговый кватернион, который должен оказаться вектором:



Мы «выкинули» (говоря русским языком: похерили), поскольку векторное произведение ортогонально каждому из своих множителей, т.е имеет нулевое скалярное произведение с ними. Можно также сказать, что это слагаемое – не что иное, как смешанное произведение трёх векторов и равно нулю, если они линейно зависимы.
После упрощения имеем:
(6.7)


Заметим, что скалярная часть равна нулю, то есть перед нами действительно вектор.
Запишем формулу для двойного векторного произведения «бац минус цаб»:
(6.8)


но не будем убирать из (6.7) двойное векторное произведение, вместо этого выразим:


Подставляем его в (6.7):




что снова совпадает с формулой (6.6).

Вариации на тему
Может показаться, что два векторных произведения – слишком много для поворота вектора, надо обойтись одним. Да, это возможно. Подставим (6.8) в (6.7):




Первый множитель можно упростить, учитывая, что кватернион имеет единичную норму, а также вынести двойку:
(6.9)


Посчитаем количество операций в формуле (6.9).
Вычисление первого множителя требует 3m+1a. Затем 3m, чтобы помножить вектор на скаляр. Далее, 6m+3a для векторного произведения, 3m – умножение вектора на скаляр, 3m+2a – скалярное произведение, 3m – умножение вектора на скаляр, 3a – сложение векторов, 3a – их удвоение (либо 3m, не так уж важно в современных компьютерах), 3a – сложение векторов.
Итого: 21 умножение и 15 сложений – несколько хуже, чем по формуле (6.6).
Но как мы знаем, кроме количества арифметических операций, в современных процессорах очень важно, насколько эти операции «независимы» - требует ли каждая следующая результатов предыдущей, или можно начать выполнять их параллельно. В этом плане формула (6.9) и даже самые простые выкладки (6.3), (6.5) могут оказаться предпочтительнее.

Вычисление матрицы поворота
В компьютерной графике (и не только) часто встречается ситуация, когда необходимо повернуть множество векторов с помощью одного и того же кватерниона. В таком случае имеет смысл преобразовать кватернион в матрицу и производить поворот уже с её помощью.
Вывести матрицу поворота, соответствующую кватерниону Λ – несложно. Заметим, что все выведенные нами формулы линейны относительно входного вектора . Поэтому достаточно повернуть базисные векторы, после чего выразить поворот произвольного вектора как их взвешенную сумму.
Найдём, как повернётся вектор (1;0;0)T, он же i в кватернионной записи. Как ни странно, наиболее просто это делается «по определению». Здесь, когда используется всего один кватернион, для упрощения работы уберём индексы, записав его так:


Тогда:








Первый множитель можно упростить, учитывая единичную норму:



Повторим то же самое для вектора (0;1;0)T, он же j:









И, наконец, для вектора (0;0;1)T, он же k:









Наконец, из векторов можно составить матрицу:



На её подсчет требуется не более 13 умножений и 12 сложений. Такое количество операций достигается, если мы первоначально умножим наш кватернион на (потребует 4 умножения), благодаря чему все двойки из выражения уйдут. Дальше мы вычислим произведения x2, y2, z2, ax, ay, az, xy, xz, yz (ещё 9 умножений), и останется лишь сложить всё вместе.
Возможно, есть и более быстрые методы, но, как правило, вычисление матрицы занимает ничтожно мало времени по сравнению с последующим поворотом сотен векторов, поэтому вовсе не обязательно оптимизировать эту операцию.


Пока что кватернионы проигрывают в простоте и эффективности матрицам поворота (хотя и не значительно). Их полюбили в космической технике явно за другое… Продолжение следует.
Tags: кватернионы-это просто (том 1), математика, работа
Subscribe

  • Ремонт лыжных мостиков

    Вернулся с сегодняшнего субботника. Очень продуктивно: отремонтировали все ТРИ мостика! Правда, для этого надо было разделиться, благо народу…

  • Гетто-байк

    В субботу во время Великой Октябрьской резни бензопилой умудрился петуха сломать в велосипеде. По счастью, уже на следующий день удалось купить…

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 7 comments