Ликбез по кватернионам, часть 6: поворот по кратчайшему пути

Часть 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-го порядка
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Кватернион поворота по кратчайшему пути

Рассмотрим часто встречающуюся на практике задачу, которая красиво решается с помощью кватернионов.
У нас есть единичные векторы и . Требуется найти поворот, который переведёт в по кратчайшему пути.


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


Здесь мы можем принять φ ≥ 0, поскольку вместо отрицательного угла можно развернуть ось на 180°. - это единичный вектор. Можно выразить их следующим образом:


Можно также найти скалярное произведение векторов:


Поскольку φ ≥ 0, то именно cosφ задаст угол φ однозначно, тогда как один и тот же sinφ может соответствовать двум разным углам. Например, если


это может соответствовать φ = π/4 либо φ = 3π/4.
Кватернион поворота по кратчайшему пути примет вид:


Мы можем использовать тригонометрические формулы половинных углов:



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


Скажем, мы используем числа с плавающей точкой одинарной точности (32 бита, IEEE 754), что не так уж редко как для компьютерной графики и вычислений через видеокарту, так и в космических приложениях, где «железо» заметно отстаёт по производительности от «земного».
Из 32 бит, 23 используются для мантиссы. Ближайшее к единице число, которое может быть записано – это

2^-1 × 1,11111111111111111111111₂ ≈ 0,999999940395355,

что соответствует углу 345 мкрад = 1'11’’ – более 1/60 градуса. Такой точности недостаточно для многих приложений.
Попробуем преобразовать формулу (7.5), подставив туда (7.2), (7.3) и (7.4):


Наконец, из формулы косинуса двойного угла


получаем


и приходим к формуле:


Заметим, что нам не нужно считать множитель перед скобкой напрямую по тригонометрическим формулам, достаточно сначала посчитать ненормированный кватернион:


и отнормировать его:


Вычисления по формулам (7.7) и (7.8) более просты и элегантны, чем непосредственное использование (7.2), (7.4), (7.5).
Кроме того, в них естественным образом учитывается вариант (т.е направления изначально совпадают). Если при вычислениях по (7.2) мы бы наткнулись на деление на ноль (векторное произведение равно нулю), то здесь мы получаем


Λ=1

Вычисления также являются численно более устойчивыми.
Пример 1.



для представления векторов используются числа одинарной точности.
При непосредственном использовании (7.2), (7.4), (7.5) получаем






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

Стоит обратить внимание, что потеря точности произойдёт, даже если при вычислении (7.2), (7.4), (7.5) мы будем использовать двойную или расширенную (80 бит) точность для представления промежуточных и конечных результатов, ведь потеря произошла ещё при попытке представить единичный вектор через числа с одинарной точностью!
Посмотрим теперь результат при использовании (7.7) и (7.8):





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

Пример 2.


для представления векторов используются числа одинарной точности.
Применение (7.2), (7.4), (7.5) по-прежнему приводит к полной потере точности – к нулевому повороту.
При использовании (7.7) и (7.8) получаем


В этот раз ошибка составила 6 × 10^-10 радиан, или около 10^-4 угловых секунд. Вероятно, представленный случай близок к наихудшему, когда при вычислении векторного произведения происходит вычитание близких величин. Тем не менее, результат очень достойный.

Особый случай.
Но формулы (7.7) и (7.8) всё-таки имеют одну «особую точку», которую придётся рассмотреть отдельно. Это случай (поворот ровно на 180 градусов), что даст ненормированный кватернион


и деление на ноль при попытке его нормализовать.
Проблема в том, что разворот на 180 градусов можно сделать по любой оси, перпендикулярной , и любой из них будет «кратчайшим». Из-за теоремы о причесывании ежа ни одна линейная (и вообще непрерывная) формула наподобие (7.7) не способна в такой ситуации выбрать «произвольный» вектор, поэтому мы должны рассмотреть этот случай отдельно. Если мы приходим к нулевому кватерниону, то выбираем какой-нибудь единичный вектор, перпендикулярный , назовём его (как его найти, см. https://nabbla1.livejournal.com/96041.html)
После чего строим итоговый кватернион:


Теперь наш метод становится универсальным. Запишем его целиком:


Если нас интересует ответ в виде матрицы поворота, мы можем вычислить её из этого кватерниона - см. часть 5. Как ни странно, существенно упростить вычисление этой матрицы, не прибегая к кватернионам - не удастся, они очень органично вписываются в данную задачу.

Рассмотрим, где эта задача может быть применена. В моей практике это - "коррекция по одной звезде". На спутнике стоит опорно-поворотное устройство с объективом. Его наводят на звезду по априорным координатам. Если бы ориентация спутника была точной, звезда возникла бы в "перекрестье", но в действительности она оказывается смещена от центра. Одна звезда не позволяет полностью определить ориентацию в пространстве, но, построив кватернион поворота из расчетного направления в фактическое, мы можем ввести поправку, которая уберёт 2 компоненты ошибки из 3. Нескомпенсированным окажется вращение вокруг направления на звезду.

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