Ликбез по кватернионам, Дэвенпорт VS Мортари-Маркли!
Оглавление "Ликбеза по кватернионам":
[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
В прошлый раз мы решили " лидарную задачу", причём для нахождения поворота использовали линейный метод Мортари-Маркли. Но в оригинальной работе [Berthold K. P. Horn - Closed-form solution of absolute orientation using unit quaternions - 1987], по сути, предлагалось использовать метод Дэвенпорта.
Мы составили систему из 3 линейных уравнений с 3 неизвестными, решили её, получив вектор Родрига, и затем преобразовали его в кватернион. Но авторы предлагали составить матрицу 4х4, и найти собственный вектор, соответствующий МАКСИМАЛЬНОМУ собственному значению. Выглядит "посолиднее" - может, именно так и надо поступать? А мы пошли по упрощённому пути и получили ошибку в 0,59°. Вдруг дело не в зашумлённых входных данных, а именно в применённом методе, который не смог этот шум надёжно отфильтровать?
Давайте проверим...
Напомню, в нашей модельной задаче был набор из 8 "эталонных", отцентрированных векторов r'k (их сумма даст ноль):
r'1 (45; -104; -0,25),
r'2 (45; 104; -0,25),
r'3 (-15; -69; 42,75),
r'4 (-15; -84; -51,25),
r'5 (-15; -158; 8,75),
r'6 (-15; 158; 8,75),
r'7 (-15; 69; 42,75),
r'8 (-15; 84; -51,25)
и соответствующие им 8 "измеренных", отцентрированных векторов m'k:
m'1 (-104,125; -45,425; -0,9625),
m'2 (103,875; -44,225; 0,1375),
m'3 (-68,725; 13,875; 42,2375),
m'4 (-83,425; 15,675; -50,5625),
m'5 (-158,725; 15,175; 9,5375),
m'6 (158,875; 13,875; 7,9375),
m'7 (68,875; 15,375; 42,4375),
m'8 (83,375; 15,675; -50,7625)
Хочется найти такой поворот, который наилучшим образом переведёт r'k в m'k.
1. Посчитать матрицу B:
)
Для нашего примера получаем следующую матрицу:
)
Своеобразно... Один ОЧЕНЬ БОЛЬШОЙ коэффициент b21 и два просто БОЛЬШИХ: b12 и b33, остальное - мелочь пузатая.
2. Превращаем матрицу B в вектор Z и матрицу Дэвенпорта K:
)
)
В данном случае получается:
)
)
(решили числа округлить - когда коэффициент в сотню тысяч, можно уже без десятых долей обойтись)
3. Находим максимальное собственное значение матрицы Дэвенпорта и соответствующий собственный вектор M.
На самом деле, задача не настолько сурова, как кажется на первый взгляд. У матрицы Дэвенпорта важная особенность, она обладает нулевым следом! Это означает, что в характеристическом уравнении коэффициент при λ3 равен нулю. Дальнейшие изыскания, проведённые [здесь], показывают, что коэффициент при λ2 ВСЕГДА будет отрицательным, и это всё позволяет найти максимальное собственное значение не так тяжело, как "в общем случае".
Но поскольку мы пока не собираемся исполнять этот алгоритм "в железе", то не будем мучаться и снова напряжём WolframAlpha:
Как и должно быть в случае симметричной матрицы, все 4 собственных значения действительные, самое большое - первое. Можно увидеть, что WolframAlpha при построении собственных векторов просто делает последнее значение единичным, а все остальные уже выражаются однозначно.
4. Нормируем собственный вектор на единичную длину, и ещё домножим на -1, чтобы проще было сравнивать:
(индекс d указывает, что по методу Дэвенпорта)
В целом похож. Напомним, линейный метод Мортари-Маркли дал такой кватернион:
Расстояние между ними ( метрика) составляет 0,069°. Поворот, вычисленный по методу Мортари-Маркли, отличался от истинного на 0,59°. Наконец, поворот, вычисленный по методу Дэвенпорта, отличается на 0,52°, т.е всё-таки ошибка вышла на 11% меньше.
Результат немного неожиданный: я думал, что эти два метода вообще дадут практически идентичные кватернионы, но нет. И в этом конкретном примере метод Дэвенпорта дал ошибку на 11% меньше. Полагаю, что это anecdotal evidence: если поиграться подольше, то окажется: иногда один чуть точнее, иногда другой, тут уж как повезёт. По крайней мере, в статье Мортари-Маркли получается именно так.
Может, как-нибудь ещё с этим делом поковыряюсь, но пока моя точка зрения не изменилась: на практике вполне можно применять метод Мортари-Маркли, он простой до безобразия :) А в той статье использовали Дэвенпорта по простой причине, статья 1987 года (а отдана в редакцию и вовсе в 1986), а метод Мортари-Маркли появился в 2007 году!
PS Впрочем, ещё одну вещь надо проверить: не ухудшает ли работу метода Мортари-Маркли тот факт, что эталонный и измеренный векторы имеют разную длину? Что Дэвенпорту не мешают - никаких сомнений, это видно из структуры матрицы B (см. Дэвенпорт берёт след), где весовой коэффициент можно было бы внести в длину векторов, т.е "чем точнее мы знаем вектор - тем длиннее его сделаем". И пока речь шла о направлениях, т.е единичных векторах, всё хорошо. Но тут-то из-за ошибок измерения и длина может оказываться различной...
[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
В прошлый раз мы решили " лидарную задачу", причём для нахождения поворота использовали линейный метод Мортари-Маркли. Но в оригинальной работе [Berthold K. P. Horn - Closed-form solution of absolute orientation using unit quaternions - 1987], по сути, предлагалось использовать метод Дэвенпорта.
Мы составили систему из 3 линейных уравнений с 3 неизвестными, решили её, получив вектор Родрига, и затем преобразовали его в кватернион. Но авторы предлагали составить матрицу 4х4, и найти собственный вектор, соответствующий МАКСИМАЛЬНОМУ собственному значению. Выглядит "посолиднее" - может, именно так и надо поступать? А мы пошли по упрощённому пути и получили ошибку в 0,59°. Вдруг дело не в зашумлённых входных данных, а именно в применённом методе, который не смог этот шум надёжно отфильтровать?
Давайте проверим...
Напомню, в нашей модельной задаче был набор из 8 "эталонных", отцентрированных векторов r'k (их сумма даст ноль):
r'1 (45; -104; -0,25),
r'2 (45; 104; -0,25),
r'3 (-15; -69; 42,75),
r'4 (-15; -84; -51,25),
r'5 (-15; -158; 8,75),
r'6 (-15; 158; 8,75),
r'7 (-15; 69; 42,75),
r'8 (-15; 84; -51,25)
и соответствующие им 8 "измеренных", отцентрированных векторов m'k:
m'1 (-104,125; -45,425; -0,9625),
m'2 (103,875; -44,225; 0,1375),
m'3 (-68,725; 13,875; 42,2375),
m'4 (-83,425; 15,675; -50,5625),
m'5 (-158,725; 15,175; 9,5375),
m'6 (158,875; 13,875; 7,9375),
m'7 (68,875; 15,375; 42,4375),
m'8 (83,375; 15,675; -50,7625)
Хочется найти такой поворот, который наилучшим образом переведёт r'k в m'k.
1. Посчитать матрицу B:
Для нашего примера получаем следующую матрицу:
Своеобразно... Один ОЧЕНЬ БОЛЬШОЙ коэффициент b21 и два просто БОЛЬШИХ: b12 и b33, остальное - мелочь пузатая.
2. Превращаем матрицу B в вектор Z и матрицу Дэвенпорта K:
В данном случае получается:
(решили числа округлить - когда коэффициент в сотню тысяч, можно уже без десятых долей обойтись)
3. Находим максимальное собственное значение матрицы Дэвенпорта и соответствующий собственный вектор M.
На самом деле, задача не настолько сурова, как кажется на первый взгляд. У матрицы Дэвенпорта важная особенность, она обладает нулевым следом! Это означает, что в характеристическом уравнении коэффициент при λ3 равен нулю. Дальнейшие изыскания, проведённые [здесь], показывают, что коэффициент при λ2 ВСЕГДА будет отрицательным, и это всё позволяет найти максимальное собственное значение не так тяжело, как "в общем случае".
Но поскольку мы пока не собираемся исполнять этот алгоритм "в железе", то не будем мучаться и снова напряжём WolframAlpha:
Как и должно быть в случае симметричной матрицы, все 4 собственных значения действительные, самое большое - первое. Можно увидеть, что WolframAlpha при построении собственных векторов просто делает последнее значение единичным, а все остальные уже выражаются однозначно.
4. Нормируем собственный вектор на единичную длину, и ещё домножим на -1, чтобы проще было сравнивать:
(индекс d указывает, что по методу Дэвенпорта)
В целом похож. Напомним, линейный метод Мортари-Маркли дал такой кватернион:
Расстояние между ними ( метрика) составляет 0,069°. Поворот, вычисленный по методу Мортари-Маркли, отличался от истинного на 0,59°. Наконец, поворот, вычисленный по методу Дэвенпорта, отличается на 0,52°, т.е всё-таки ошибка вышла на 11% меньше.
Результат немного неожиданный: я думал, что эти два метода вообще дадут практически идентичные кватернионы, но нет. И в этом конкретном примере метод Дэвенпорта дал ошибку на 11% меньше. Полагаю, что это anecdotal evidence: если поиграться подольше, то окажется: иногда один чуть точнее, иногда другой, тут уж как повезёт. По крайней мере, в статье Мортари-Маркли получается именно так.
Может, как-нибудь ещё с этим делом поковыряюсь, но пока моя точка зрения не изменилась: на практике вполне можно применять метод Мортари-Маркли, он простой до безобразия :) А в той статье использовали Дэвенпорта по простой причине, статья 1987 года (а отдана в редакцию и вовсе в 1986), а метод Мортари-Маркли появился в 2007 году!
PS Впрочем, ещё одну вещь надо проверить: не ухудшает ли работу метода Мортари-Маркли тот факт, что эталонный и измеренный векторы имеют разную длину? Что Дэвенпорту не мешают - никаких сомнений, это видно из структуры матрицы B (см. Дэвенпорт берёт след), где весовой коэффициент можно было бы внести в длину векторов, т.е "чем точнее мы знаем вектор - тем длиннее его сделаем". И пока речь шла о направлениях, т.е единичных векторах, всё хорошо. Но тут-то из-за ошибок измерения и длина может оказываться различной...