Ликбез по кватернионам, часть 5 5/8 - метрика ненормированных кватернионов
Оглавление "Ликбеза по кватернионам":
[Spoiler (click to open)]
Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
часть 5 5/8 - метрика ненормированных кватернионов
часть 5 11/16 - красивая псевдометрика произвольных кватернионов
Часть 5 3/4 - исследуем "пространство поворотов"
Часть 5 7/8 - почти изотропный ёжик
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 6 7/8 - "уполовинивание угла" на плоскости
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Часть 17 - лидарная задача
Задачка к части 17
Дэвенпорт VS Мортари-Маркли
Мортари-Маркли берут реванш!
Дэвенпорт VS Мортари-Маркли, раунд 3
В предыдущей главе мы вводили метрику ПОВОРОТОВ, и соответствующую ей ПСЕВДОМЕТРИКУ кватернионов. ПСЕВДО означает, что она может равняться нулю для двух разных кватернионов, нам это очень кстати, ведь один и тот же поворот всегда может выражаться двумя разными кватернионами, отличающимися только знаком.
Но мы ограничились рассмотрением кватернионов с единичной нормой (именно они выражают чистый поворот, без масштабирования), из-за чего получили несколько пессимистичную оценку относительно точности представления кватернионов 16-битными числами.
Сейчас мы рассмотрим более общий случай и ещё поразмышляем об оптимальной числовой системе для описания кватернионов.
Напомним: мы брали два кватерниона с единичной нормой, Λ и Μ, и вводили метрику следующим образом:
=2|\Lambda-M|\sqrt{1-\frac{|\Lambda-M|^2}{4}})
Корень необходим, чтобы метрика поворота корректно считалась, когда Λ≈-Μ. Если же мы рассматриваем Λ≈Μ, корнем можно пренебречь и в первом приближении записать:
=2|\Lambda-M|)
(p от слова proximity - близость)
Теперь уберём условие единичных норм кватернионов. В таком случае, метрику проще всего переписать следующим образом:
=2|\frac{\Lambda}{|\Lambda|}-\frac{M}{|M|}|)
Надо лишь потребовать, чтобы Λ≠0 и Μ≠0: нулевой кватернион не может выражать поворот, его невозможно отнормировать.
Это не единственный способ задать метрику. Сейчас нас не волнует, что при использовании формулы поворота из главы 3 вектор может не только повернуться, но и изменить длину. Можно ввести другую метрику, которая учитывает и паразитное изменение длины вектора, но пока мы считаем, что всегда сможем отнормировать либо кватернион, либо вектор, либо нас длина вектора в принципе не интересует.
Теперь применим такую метрику для оценки ошибки представления кватерниона в виде тех или иных чисел. Истинный кватернион с единичной нормой снова назовём Λ. Но из-за ограниченной точности компьютерной арифметики, в действительности мы сможем записать лишь кватернион Μ=Λ+ΔΛ, норма которого, скорее всего, будет отличаться от единицы.
Найдём, насколько поворот, выражаемый таким кватернионом, отличается от истинного:
=2|\Lambda-\frac{\Lambda+\Delta\Lambda}{|\Lambda+\Delta\Lambda|}|=2|\frac{\Lambda+\Delta\Lambda}{\sqrt{|\Lambda+\Delta\Lambda|^2}}-\Lambda|)
В предыдущей главе мы ввели скалярное произведение кватернионов по обычной формуле, сумма попарных произведений элементов, т.е
=\lambda_0\mu_0+\lambda_1\mu_1+\lambda_2\mu_2+\lambda_3\mu_3)
Мало того, тогда можно при возведении нормы кватерниона в квадрат раскрывать скобки по обычным правилам:
+|\Delta\Lambda|^2\approx1+2(\Lambda,\Delta\Lambda))
(мы учли, что кватернион Λ имеет единичную норму, а также отбросили слагаемое второго порядка малости)
Итак:
\approx2|\frac{\Lambda+\Delta\Lambda}{\sqrt{1+2(\Lambda,\Delta\Lambda)}}-\Lambda|)
В знаменателе, под корнем у нас обычное действительное число. Опять, считая ΔΛ достаточно малым, разложим выражение в ряд Тейлора, оставив первых два слагаемых:
\approx2|(1-(\Lambda,\Delta\Lambda))(\Lambda+\Delta\Lambda)-\Lambda|\approx2|\cancel{\Lambda}+\Delta\Lambda-(\Lambda,\Delta\Lambda)\Lambda-\cancel{\Lambda}|=2|\Delta\Lambda-(\Lambda,\Delta\Lambda)\Lambda|)
Полученное выражение имеет ясный физический смысл: мы взяли кватернион разности, ΔΛ, и вычли из него проекцию на кватернион Λ, поскольку данная составляющая ΔΛ только изменяет длину (норму) исходного кватерниона, но он продолжает выражать тот же самый поворот, а поэтому данное изменение не должно включаться в нашу метрику.
Заменим норму на корень из квадрата нормы, который можно раскрыть:
\Lambda|^2=|\Delta\Lambda|^2-2(\Lambda,\Delta\Lambda)^2+(\Lambda,\Delta\Lambda)^2=|\Delta\Lambda|^2-(\Lambda,\Delta\Lambda)^2)
Выпишем кватернионы покомпонентно:
Тогда получим покомпонентное выражение для метрики:
=2\sqrt{\Delta\lambda_0^2+\Delta\lambda_1^2+\Delta\lambda_2^2+\Delta\lambda_3^3-(\Delta\lambda_0\lambda_0+\Delta\lambda_1\lambda_1+\Delta\lambda_2\lambda_2+\Delta\lambda_3\lambda_3)^2})
Вряд ли эту формулу можно существенно упростить, при произвольных значениях Δλn. Скобку можно раскрыть, что даст 4 квадрата и 6 перекрёстных слагаемых ("каждый с каждым").
Но мы можем с хорошей точностью утверждать, что если ΔΛ - это ошибки представления кватерниона Λ, то его отдельные компоненты - это независимые случайные величины с нулевыми средними, т.е
=0,n{\ne}k)
Тогда мы можем найти математическое ожидание метрики dp:
)\approx2\sqrt{E(\Delta\lambda_0^2)+E(\Delta\lambda_1^2)+E(\Delta\lambda_2^2)+E(\Delta\lambda_3^2)-\lambda_0^2E(\Delta\lambda_0^2)-\lambda_1^2E(\Delta\lambda_1^2)-\lambda_2^2E(\Delta\lambda_2^2)-\lambda_3^2E(\Delta\lambda_3^2)})
Равенство лишь приближённое, потому что, вообще говоря, математическое ожидание нельзя "заносить под знак корня". Простейший пример: пусть у нас будет случайная величина ξ которая может принимать два значения, 0 или 1, с вероятностями по 50%. Её математическое ожидание составит 0,5, и корень из него: около 0,707. С другой стороны, если извлекать корень из самой случайной величины, у нас снова будет получаться 0 или 1, с вероятностями по 50%. И математическое ожидание от корня будет ровно 0,5. Но когда в подкоренном выражении участвует несколько случайных величин с непрерывными распределениями, больших проблем не будет.
Математическое ожидание от квадрата случайной величины с нулевым средним - это не что иное, как её дисперсия, или квадрат от среднеквадратичного отклонения:
=\sigma_n^2)
для n=0..3. Мы не хотим ограничиваться числами с фиксированной запятой, поэтому сейчас допускаем, что σn могут быть все разными.
Упростим выражение:
)=2\sqrt{\sigma_0^2(1-\lambda_0^2)+\sigma_1^2(1-\lambda_1^2)+\sigma_2^2(1-\lambda_2^2)+\sigma_3^2(1-\lambda_3^2)})
Это ответ.
Если кватернион Λ имеет все нулевые компоненты, кроме одной, то оставшаяся компонента должна равняться ±1, и тогда соответствующая скобка в выражении обратится в ноль, позволяя задавать единичное значение с огромными погрешностями - покуда нас интересует только поворот, это не имеет никакого значения! В частности, поэтому мы спокойно отнеслись к формату Q1.15, позволяющему записывать числа от -1 до 1-2-15≈0,999969 - мы понимали, что единица может получиться только в том случае, когда все остальные компоненты нулевые, а значит, направление кватерниона верное!
Рассмотрим ещё один частный случай: σ0=σ1=σ2=σ3=σ, т.е все компоненты записаны с одинаковой погрешностью. Тогда сигму вынесем за скобки и за знак корня:
)=2\sigma\sqrt{1-\lambda_0^2+1-\lambda_1^2+1-\lambda_2^2+1-\lambda_3^2}=2\sigma\sqrt{3})
(снова вспоминаем, что Λ имеет единичную норму, поэтому сумма квадратов всех компонентов обращается в 1)
Это даёт нам более корректную оценку точности представления кватернионов в числах с фиксированной запятой. В предыдущей главе мы умножали ошибку в одном компоненте на корень из четырёх, т.е на 2, поскольку не учли, что из 4 направлений, нам вредят только 3, а последнее меняет норму кватерниона, что не влияет на представление поворота. Кроме того, там мы по-простому записали ошибку как половинку младшего разряда. Если уж погружаться в случайные величины, то лучше представить ошибку как равномерно распределённую по [-lsb/2, lsb/2] случайную величину (lsb - это least significant bit, т.е цена младшего разряда). Такая ошибка имеет нулевое математическое ожидание и дисперсию, равную 1/12. При таких вводных, получим следующую метрику ошибки:
)=2\sigma\sqrt{3}=\frac{2\sqrt{3}lsb}{\sqrt{12}}=lsb)
Очень простой результат: среднеквадратичная ошибка представления кватерниона, в радианах, просто равна цене младшего разряда! В нашем случае Q1.15 (16-битное число со знаком, 1 бит перед запятой) цена младшего разряда составляет примерно 3·10-5, такова и будет ошибка в радианах, или 6 угловых секунд. С вероятностью 70% ошибка будет укладываться в диапазон [-6'', +6'']. А с вероятностью 95%, ошибка будет укладываться в диапазон [-12'', +12'']. Но в прошлой главе мы пришли к выводу, что диапазон [-12'', +12''] - это вообще максимальная ошибка, и для этого даже не пришлось "вычленять" часть ошибки, меняющую норму кватерниона. Можно ли дать более жёсткую оценку сверху - пока сказать не могу, для этого надо с математического ожидания вернуться к честной формуле, и отыскать наихудший случай.
Но вот эта формула:
намекает на ещё одну возможность, использовать нелинейную шкалу для компонентов кватерниона, более плотную вблизи нуля, и более редкую по мере приближения к единице!
Кажется правильным решением сделать шаг, пропорциональный
Если шаг вблизи нуля обозначить за единицу, то шаг вблизи 0,5 должен возрасти до 1,15 (на 15%), вблизи 0,9 - до 2,29. Надо сказать, это достаточно медленное увеличение. Это точно не числа с плавающей точкой, по крайней мере, не стандартные IEEE.
Аналитического выражения нужная нам кривая не имеет, она получается решением уравнения наподобие
Так что решал я его численно, тупо в экселе. Пока рассмотрел "игрушечную задачу", где каждый компонент кватерниона мы пытаемся вместить в 1 байт. Тогда, учитывая знак, у нас остаётся 128 отсчётов на интервал [0;1). Можно считать, мы заводим LUT (LookUp Table) из 128 значений, и каждое значение байта соответствует одному из чисел этой таблицы. График выглядит так, нас интересует центральная кривая.
Когда я только увидел выражение, мне показалось, что кривая должна являться арксинусом, но впечатление было обманчивым, видно, что это не арксинус.
Уменьшенный наклон нашей кривой вблизи нуля вроде как позволяет расставить там отсчёты более часто, тем самым повысив точность при той же самой разрядности.
Главный вопрос - НАСКОЛЬКО? Наклон вблизи нуля составляет примерно 0,8 от "линейного" - значит ли это, что мы повысили точность на 20%? Нет, не значит: у нас в кватернионе обязан присутствовать "крупный" компонент, и то, что мы "сэкономим" вблизи нуля, частично придётся искупить повышенной ошибкой в этом крупном компоненте.
Реальный выигрыш составляет 9%. Прямо скажем, негусто, столько возни ради этого!
Так что остаюсь при своём мнении: оптимальным представлением компонент кватерниона являются числа с фиксированной запятой, в диапазоне -1..≈+1.
[Spoiler (click to open)]
Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
часть 5 5/8 - метрика ненормированных кватернионов
часть 5 11/16 - красивая псевдометрика произвольных кватернионов
Часть 5 3/4 - исследуем "пространство поворотов"
Часть 5 7/8 - почти изотропный ёжик
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 6 7/8 - "уполовинивание угла" на плоскости
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Часть 17 - лидарная задача
Задачка к части 17
Дэвенпорт VS Мортари-Маркли
Мортари-Маркли берут реванш!
Дэвенпорт VS Мортари-Маркли, раунд 3
В предыдущей главе мы вводили метрику ПОВОРОТОВ, и соответствующую ей ПСЕВДОМЕТРИКУ кватернионов. ПСЕВДО означает, что она может равняться нулю для двух разных кватернионов, нам это очень кстати, ведь один и тот же поворот всегда может выражаться двумя разными кватернионами, отличающимися только знаком.
Но мы ограничились рассмотрением кватернионов с единичной нормой (именно они выражают чистый поворот, без масштабирования), из-за чего получили несколько пессимистичную оценку относительно точности представления кватернионов 16-битными числами.
Сейчас мы рассмотрим более общий случай и ещё поразмышляем об оптимальной числовой системе для описания кватернионов.
Напомним: мы брали два кватерниона с единичной нормой, Λ и Μ, и вводили метрику следующим образом:
Корень необходим, чтобы метрика поворота корректно считалась, когда Λ≈-Μ. Если же мы рассматриваем Λ≈Μ, корнем можно пренебречь и в первом приближении записать:
(p от слова proximity - близость)
Теперь уберём условие единичных норм кватернионов. В таком случае, метрику проще всего переписать следующим образом:
Надо лишь потребовать, чтобы Λ≠0 и Μ≠0: нулевой кватернион не может выражать поворот, его невозможно отнормировать.
Это не единственный способ задать метрику. Сейчас нас не волнует, что при использовании формулы поворота из главы 3 вектор может не только повернуться, но и изменить длину. Можно ввести другую метрику, которая учитывает и паразитное изменение длины вектора, но пока мы считаем, что всегда сможем отнормировать либо кватернион, либо вектор, либо нас длина вектора в принципе не интересует.
Теперь применим такую метрику для оценки ошибки представления кватерниона в виде тех или иных чисел. Истинный кватернион с единичной нормой снова назовём Λ. Но из-за ограниченной точности компьютерной арифметики, в действительности мы сможем записать лишь кватернион Μ=Λ+ΔΛ, норма которого, скорее всего, будет отличаться от единицы.
Найдём, насколько поворот, выражаемый таким кватернионом, отличается от истинного:
В предыдущей главе мы ввели скалярное произведение кватернионов по обычной формуле, сумма попарных произведений элементов, т.е
Мало того, тогда можно при возведении нормы кватерниона в квадрат раскрывать скобки по обычным правилам:
(мы учли, что кватернион Λ имеет единичную норму, а также отбросили слагаемое второго порядка малости)
Итак:
В знаменателе, под корнем у нас обычное действительное число. Опять, считая ΔΛ достаточно малым, разложим выражение в ряд Тейлора, оставив первых два слагаемых:
Полученное выражение имеет ясный физический смысл: мы взяли кватернион разности, ΔΛ, и вычли из него проекцию на кватернион Λ, поскольку данная составляющая ΔΛ только изменяет длину (норму) исходного кватерниона, но он продолжает выражать тот же самый поворот, а поэтому данное изменение не должно включаться в нашу метрику.
Заменим норму на корень из квадрата нормы, который можно раскрыть:
Выпишем кватернионы покомпонентно:
Тогда получим покомпонентное выражение для метрики:
Вряд ли эту формулу можно существенно упростить, при произвольных значениях Δλn. Скобку можно раскрыть, что даст 4 квадрата и 6 перекрёстных слагаемых ("каждый с каждым").
Но мы можем с хорошей точностью утверждать, что если ΔΛ - это ошибки представления кватерниона Λ, то его отдельные компоненты - это независимые случайные величины с нулевыми средними, т.е
Тогда мы можем найти математическое ожидание метрики dp:
Равенство лишь приближённое, потому что, вообще говоря, математическое ожидание нельзя "заносить под знак корня". Простейший пример: пусть у нас будет случайная величина ξ которая может принимать два значения, 0 или 1, с вероятностями по 50%. Её математическое ожидание составит 0,5, и корень из него: около 0,707. С другой стороны, если извлекать корень из самой случайной величины, у нас снова будет получаться 0 или 1, с вероятностями по 50%. И математическое ожидание от корня будет ровно 0,5. Но когда в подкоренном выражении участвует несколько случайных величин с непрерывными распределениями, больших проблем не будет.
Математическое ожидание от квадрата случайной величины с нулевым средним - это не что иное, как её дисперсия, или квадрат от среднеквадратичного отклонения:
для n=0..3. Мы не хотим ограничиваться числами с фиксированной запятой, поэтому сейчас допускаем, что σn могут быть все разными.
Упростим выражение:
Это ответ.
Если кватернион Λ имеет все нулевые компоненты, кроме одной, то оставшаяся компонента должна равняться ±1, и тогда соответствующая скобка в выражении обратится в ноль, позволяя задавать единичное значение с огромными погрешностями - покуда нас интересует только поворот, это не имеет никакого значения! В частности, поэтому мы спокойно отнеслись к формату Q1.15, позволяющему записывать числа от -1 до 1-2-15≈0,999969 - мы понимали, что единица может получиться только в том случае, когда все остальные компоненты нулевые, а значит, направление кватерниона верное!
Рассмотрим ещё один частный случай: σ0=σ1=σ2=σ3=σ, т.е все компоненты записаны с одинаковой погрешностью. Тогда сигму вынесем за скобки и за знак корня:
(снова вспоминаем, что Λ имеет единичную норму, поэтому сумма квадратов всех компонентов обращается в 1)
Это даёт нам более корректную оценку точности представления кватернионов в числах с фиксированной запятой. В предыдущей главе мы умножали ошибку в одном компоненте на корень из четырёх, т.е на 2, поскольку не учли, что из 4 направлений, нам вредят только 3, а последнее меняет норму кватерниона, что не влияет на представление поворота. Кроме того, там мы по-простому записали ошибку как половинку младшего разряда. Если уж погружаться в случайные величины, то лучше представить ошибку как равномерно распределённую по [-lsb/2, lsb/2] случайную величину (lsb - это least significant bit, т.е цена младшего разряда). Такая ошибка имеет нулевое математическое ожидание и дисперсию, равную 1/12. При таких вводных, получим следующую метрику ошибки:
Очень простой результат: среднеквадратичная ошибка представления кватерниона, в радианах, просто равна цене младшего разряда! В нашем случае Q1.15 (16-битное число со знаком, 1 бит перед запятой) цена младшего разряда составляет примерно 3·10-5, такова и будет ошибка в радианах, или 6 угловых секунд. С вероятностью 70% ошибка будет укладываться в диапазон [-6'', +6'']. А с вероятностью 95%, ошибка будет укладываться в диапазон [-12'', +12'']. Но в прошлой главе мы пришли к выводу, что диапазон [-12'', +12''] - это вообще максимальная ошибка, и для этого даже не пришлось "вычленять" часть ошибки, меняющую норму кватерниона. Можно ли дать более жёсткую оценку сверху - пока сказать не могу, для этого надо с математического ожидания вернуться к честной формуле, и отыскать наихудший случай.
Но вот эта формула:
намекает на ещё одну возможность, использовать нелинейную шкалу для компонентов кватерниона, более плотную вблизи нуля, и более редкую по мере приближения к единице!
Кажется правильным решением сделать шаг, пропорциональный
Если шаг вблизи нуля обозначить за единицу, то шаг вблизи 0,5 должен возрасти до 1,15 (на 15%), вблизи 0,9 - до 2,29. Надо сказать, это достаточно медленное увеличение. Это точно не числа с плавающей точкой, по крайней мере, не стандартные IEEE.
Аналитического выражения нужная нам кривая не имеет, она получается решением уравнения наподобие
Так что решал я его численно, тупо в экселе. Пока рассмотрел "игрушечную задачу", где каждый компонент кватерниона мы пытаемся вместить в 1 байт. Тогда, учитывая знак, у нас остаётся 128 отсчётов на интервал [0;1). Можно считать, мы заводим LUT (LookUp Table) из 128 значений, и каждое значение байта соответствует одному из чисел этой таблицы. График выглядит так, нас интересует центральная кривая.
Когда я только увидел выражение, мне показалось, что кривая должна являться арксинусом, но впечатление было обманчивым, видно, что это не арксинус.
Уменьшенный наклон нашей кривой вблизи нуля вроде как позволяет расставить там отсчёты более часто, тем самым повысив точность при той же самой разрядности.
Главный вопрос - НАСКОЛЬКО? Наклон вблизи нуля составляет примерно 0,8 от "линейного" - значит ли это, что мы повысили точность на 20%? Нет, не значит: у нас в кватернионе обязан присутствовать "крупный" компонент, и то, что мы "сэкономим" вблизи нуля, частично придётся искупить повышенной ошибкой в этом крупном компоненте.
Реальный выигрыш составляет 9%. Прямо скажем, негусто, столько возни ради этого!
Так что остаюсь при своём мнении: оптимальным представлением компонент кватерниона являются числа с фиксированной запятой, в диапазоне -1..≈+1.