nabbla (nabbla1) wrote,
nabbla
nabbla1

Category:

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

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


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

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

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

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

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

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

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

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

2-1 × 1,111111111111111111111112 ≈ 0,999999940395355,

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

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

получаем
image219.png

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

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

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

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

image224.png
Λ=1

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

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

image229.png

image230.png

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

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

image231.png

image232.png

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

Пример 2.
eq_7_8ex1.png

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

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

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

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

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

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

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

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

Recent Posts from This Journal

  • Лестница для самых жадных

    В эти выходные побывал на даче, после 3-недельной "самоизоляции". Забавно, как будто зима началась! Особенно грязные галоши остались на улице, в…

  • Возвращаемся к макету

    Очень давно макетом видеоизмерителя параметров сближения не занимался: сначала "громко думал" по поводу измерения его положения на аппарате, а потом…

  • Минутка живописи

    В процессе разгребания содержимого квартиры (после нескольких ремонтов) дошёл, наконец, и до картин. В кои-то веки их повесил. Куда их вешать -…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 4 comments