Ликбез по кватернионам, часть 6 7/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
Прямое продолжение части 6 3/4, "кватернион из синуса и косинуса угла". Там мы настаивали на произвольной точке (x,y), из которой надо вычленить лишь угол (наплевав на длину) и построить кватернион поворота на такой угол, желательно "малой кровью".
Но давайте задачу упростим: потребуем, чтобы
то есть нам, по сути, заданы косинус и синус угла, в виде пары (x,y).
Задача: получить новую пару (x', y'), которая будет выражать половинный угол (т.к в кватернионе используются половинные углы) и будет иметь длину максимально близкую к 1, чтобы потом нормализовать её минимальным количеством итераций. В идеале, если исходно мы имели
)
,
из этого мы должны получить
)
но так и быть, позволим этим двум значениям быть помноженным на небольшую константу, её мы убирать умеем...
Тогда алгоритм, описанный в предыдущей главе, сводится к следующему:
%20\rightarrow\left\{\begin{array}{rl}(1+x;y),&x\ge0\\(y;1-x),&x%3C0\end{array}\right.)
Такой "вектор" будут выражать В ТОЧНОСТИ половинный угол, но длина его будет варьироваться от
до 2. В реализации для QuatCore я делил эти выражения на 2, получая длину от ≈0,707 до 1, после чего проводил нормировку методом Ньютона. Так, после одной итерации длина будет отличаться от единицы не более чем на 11,6%.
Но чтобы длина с самого начала была наиболее близка к единичке, неплохо добавить "волшебную константу":
%20\rightarrow\left\{\begin{array}{rl}0{,}585(1+x;y),&x\ge0\\0{,}585(y;1-x),&x%3C0\end{array}\right.)
Теперь с самого начала длина будет отличаться от истинной не более чем на 17,3%, а после одной итерации метода Ньютона - не более, чем на 4,58%.
Но вчера я нашёл новый метод, также полностью ЛИНЕЙНЫЙ, который также даёт В ТОЧНОСТИ ПОЛОВИННЫЕ УГЛЫ, а длина отличается от единицы не более чем на 4%. А если к нему добавить одну итерацию метода Ньютона, отклонение уменьшается до 0,24% (по наихудшему случаю). Вот он:
%20\rightarrow\left\{\begin{array}{ll}0{,}5198(1+x;y),&-x%3Cy\le{x}\\0{,}3676(1+x+y;1-x+y),&-y%3Cx%3Cy\\0{,}5198(y;1-x),&x%3Cy\le{-x}\\0{,}3676(1+x-y;x+y-1),&y\le{x}%3C-y\end{array}\right.)
А найти его удалось после того, как старому методу была дана геометрическая интерпретация!
Первая формула, которая прибавляет единичку к иксу, объясняется элементарно:
У нас есть единичная окружность с центром в точке O. Нам задан вектор OB, координаты которого (x,y) = (cos α; sin α).
Обратим внимание на дугу AB. На неё опирается центральный угол AOB, равный α, а также вписанный угол AO1B, который, по соответствующей теореме, равен половине центрального, т.е α/2. Как же нам получить вектор O1B - да очень просто! Координату y надо оставить без внимания, а к координате x прибавить единичку, поскольку OO1=1, и этот отрезок строго горизонтален.
Такое построение будет работать для любого угла, за исключением 180° (π радиан): в этом случае точка B совпадёт с точкой O1, и мы построим вектор нулевой длины, который никакой угол выражать не может!
Поэтому необходима была вторая формула, которая сработает корректно вблизи особой точки. Её обосновать несколько труднее, но за рамки школьной геометрии выходить не придётся, просто чуть муторнее:
Это всё та же окружность радиусом 1, и снова мы имеем дело с вектором OB с координатами (x,y) = (cos α; sin α). Но теперь мы заходим в область отрицательных x, именно там "вотчина" второй формулы.
От вектора OB мы откладываем перпендикуляр по часовой стрелке, OC, также единичной длины. Вектор OC будет иметь координаты (sin α; -cos α).
Точку O1 мы в этот раз поместили внизу, и отложим от неё отрезок O1C. Мы утверждаем, что теперь он будет выражать нужный нам угол α/2 c горизонталью.
В противоположной части окружности отметим точку D, так что теперь мы имеем дугу CD, на которую опирается центральный угол COD и вписанный угол CO1D.
Угол COD равен π-α. Его можно найти, например, так: угол AOC равен α-π/2 (из полного угла α мы исключили 90 градусов, т.е угол COB), но они в сумме с углом COD должны составить 90°, или π/2. Поэтому
COD = π/2 - (α - π/2) = π - α
Вписанный угол CO1D составляет ровно половину от этого значения, или
π/2 - α/2
И наконец, угол с горизонталью составит
π/2 - (π/2 - α/2) = α/2.
Наконец, определяем координаты вектора O1C, для чего к вектору O1O (0; 1) прибавляем вектор OC (sin α; -cos α) и тут же приходим к нашей второй формуле, (y; 1-x).
Что и требовалось доказать...
В прошлой главе мы доказывали то же самое через тригонометрические выкладки. Здесь и вовсе можно было синусы и косинусы не привлекать, оставить просто x,y...
А можно ли какое-то построение придумать, в котором не будет особых точек? Судя по всему, нет! Требование получить половинный угол сразу же приводит нас к вписанным углам, т.е располагает точку аккурат на окружности! А значит, если мы будем произвольно крутиться по этой окружности, неизбежно мы в эту точку попадём, приходя к нулевому вектору. Злой ёжик не отпускает.
Эти рисунки наглядно показывают, как длина вектора, который мы вычисляем по этим формулам, варьируется от
до 2.
И ещё по этим рисункам можно догадаться, что выбор положения точки O1 на окружности, вообще-то, произволен! То есть эти две формулы не единственные.
Пока мы делали выкладки через тригонометрию, мы вроде как рассмотрели два самых очевидных варианта: домножали и делили наш кватернион поворота на синус и на косинус половинного угла.
На рисунках мы видим: в первом случае ничего крутить не пришлось, во втором мы получили поворот на 90°.
Так что возникло предчувствие: если мы будем умножать и делить наш кватернион на cos(α/2+φ) с произвольным углом φ, мы получим целое семейство таких формул.
Давайте попробуем!
Мы пытаемся получить кватернион
(учитывая, что у него используется только две компоненты из 4, здесь бы и комплексное число подошло, но для них нехарактерно использовать половинные углы. Так что уж пускай кватернион, это у нас ликбез по кватернионам или где!?)
Давайте домножим и поделим его на cos(α/2-π/4):
}\left(cos\frac{\alpha}{2}cos(\frac{\alpha}{2}-\frac{\pi}{4})+i\,sin\frac{\varphi}{2}cos(\frac{\alpha}{2}-\frac{\pi}{4})\right))
Раскроем этот косинус:
=\frac{1}{\sqrt{2}}(cos\frac{\alpha}{2}+sin\frac{\alpha}{2}))
И подставим это в наш кватернион. Множитель в начале просто обозначим как k, он не так уж интересен, отнормировать мы всегда успеем:
\right))
Каждое слагаемое выразим через двойной угол:
Двойки в знаменателе пока что тоже внесём в константу k, и получим:
\right))
И теперь, вспоминая, что x=cos α, y=sin α, и получаем вторую формулу для нашего улучшенного метода:
\right))
Надо только не забыть, что у нас дополнительно возникало деление на корень из двух, т.е выражение в скобках имеет бОльшую длину. Максимальное значение достигается для входного вектора (0;1), из него мы получаем (2;2), т.е вектор длиной
Аналогичные выкладки можно проделать, умножив и поделив кватернион на cos(α/2+π/4), и получим:
\right))
Теперь, имея 4 разные формулы, мы можем поделить окружность на 4 части, по 90 градусов каждая. Прелесть в том, что поскольку точка O1 вынесена далеко, длина вектора не успевает сильно измениться при движении на ±45°. Скажем, на первом рисунке у нас максимальная длина вектора O1B достигается при нулевом угле и составляет 2. А если повернуться в любую сторону на 45°, длина будет равняться
^2+\frac{1}{2}}\approx{1{,}8478)
что составляет ≈0,924 от максимальной длины, уменьшение менее чем на 8%. Так что если теперь ещё подобрать "волшебные константы", чтобы длина болталась в разные стороны вокруг среднего, мы и получим обещанные ±4%!
Это мы и сделали.
Лично для меня было откровением, что начав с полной окружности (с вообще любого угла), ограничившись рассмотрением 4 вариантов, чисто линейно, можно не только получить в точности половинный угол, но и практически "угадать" с длиной! Впрочем, не уверен, что данный метод найдёт применение в моих приборах: пока что у меня абсолютно нет проблем с временем вычислений, а вот память кода хочется экономить, чтобы вместиться в блоки внутренней памяти ПЛИС... Поэтому пока что для меня оптимально даже с "волшебными константами" не заморачиваться - ну проведу я пяток итераций метода Ньютона, вообще не почувствую.
Одно плохо: сейчас мы поделили плоскость на квадранты "по диагонали". Удобнее было бы проверять ещё более простые условия: знаки x, y. Что ж, вывод соответствующих формул оставляем читателю в качестве несложного упражнения :)
[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
Прямое продолжение части 6 3/4, "кватернион из синуса и косинуса угла". Там мы настаивали на произвольной точке (x,y), из которой надо вычленить лишь угол (наплевав на длину) и построить кватернион поворота на такой угол, желательно "малой кровью".
Но давайте задачу упростим: потребуем, чтобы
то есть нам, по сути, заданы косинус и синус угла, в виде пары (x,y).
Задача: получить новую пару (x', y'), которая будет выражать половинный угол (т.к в кватернионе используются половинные углы) и будет иметь длину максимально близкую к 1, чтобы потом нормализовать её минимальным количеством итераций. В идеале, если исходно мы имели
,
из этого мы должны получить
но так и быть, позволим этим двум значениям быть помноженным на небольшую константу, её мы убирать умеем...
Тогда алгоритм, описанный в предыдущей главе, сводится к следующему:
Такой "вектор" будут выражать В ТОЧНОСТИ половинный угол, но длина его будет варьироваться от
до 2. В реализации для QuatCore я делил эти выражения на 2, получая длину от ≈0,707 до 1, после чего проводил нормировку методом Ньютона. Так, после одной итерации длина будет отличаться от единицы не более чем на 11,6%.
Но чтобы длина с самого начала была наиболее близка к единичке, неплохо добавить "волшебную константу":
Теперь с самого начала длина будет отличаться от истинной не более чем на 17,3%, а после одной итерации метода Ньютона - не более, чем на 4,58%.
Но вчера я нашёл новый метод, также полностью ЛИНЕЙНЫЙ, который также даёт В ТОЧНОСТИ ПОЛОВИННЫЕ УГЛЫ, а длина отличается от единицы не более чем на 4%. А если к нему добавить одну итерацию метода Ньютона, отклонение уменьшается до 0,24% (по наихудшему случаю). Вот он:
А найти его удалось после того, как старому методу была дана геометрическая интерпретация!
Первая формула, которая прибавляет единичку к иксу, объясняется элементарно:
У нас есть единичная окружность с центром в точке O. Нам задан вектор OB, координаты которого (x,y) = (cos α; sin α).
Обратим внимание на дугу AB. На неё опирается центральный угол AOB, равный α, а также вписанный угол AO1B, который, по соответствующей теореме, равен половине центрального, т.е α/2. Как же нам получить вектор O1B - да очень просто! Координату y надо оставить без внимания, а к координате x прибавить единичку, поскольку OO1=1, и этот отрезок строго горизонтален.
Такое построение будет работать для любого угла, за исключением 180° (π радиан): в этом случае точка B совпадёт с точкой O1, и мы построим вектор нулевой длины, который никакой угол выражать не может!
Поэтому необходима была вторая формула, которая сработает корректно вблизи особой точки. Её обосновать несколько труднее, но за рамки школьной геометрии выходить не придётся, просто чуть муторнее:
Это всё та же окружность радиусом 1, и снова мы имеем дело с вектором OB с координатами (x,y) = (cos α; sin α). Но теперь мы заходим в область отрицательных x, именно там "вотчина" второй формулы.
От вектора OB мы откладываем перпендикуляр по часовой стрелке, OC, также единичной длины. Вектор OC будет иметь координаты (sin α; -cos α).
Точку O1 мы в этот раз поместили внизу, и отложим от неё отрезок O1C. Мы утверждаем, что теперь он будет выражать нужный нам угол α/2 c горизонталью.
В противоположной части окружности отметим точку D, так что теперь мы имеем дугу CD, на которую опирается центральный угол COD и вписанный угол CO1D.
Угол COD равен π-α. Его можно найти, например, так: угол AOC равен α-π/2 (из полного угла α мы исключили 90 градусов, т.е угол COB), но они в сумме с углом COD должны составить 90°, или π/2. Поэтому
COD = π/2 - (α - π/2) = π - α
Вписанный угол CO1D составляет ровно половину от этого значения, или
π/2 - α/2
И наконец, угол с горизонталью составит
π/2 - (π/2 - α/2) = α/2.
Наконец, определяем координаты вектора O1C, для чего к вектору O1O (0; 1) прибавляем вектор OC (sin α; -cos α) и тут же приходим к нашей второй формуле, (y; 1-x).
Что и требовалось доказать...
В прошлой главе мы доказывали то же самое через тригонометрические выкладки. Здесь и вовсе можно было синусы и косинусы не привлекать, оставить просто x,y...
А можно ли какое-то построение придумать, в котором не будет особых точек? Судя по всему, нет! Требование получить половинный угол сразу же приводит нас к вписанным углам, т.е располагает точку аккурат на окружности! А значит, если мы будем произвольно крутиться по этой окружности, неизбежно мы в эту точку попадём, приходя к нулевому вектору. Злой ёжик не отпускает.
Эти рисунки наглядно показывают, как длина вектора, который мы вычисляем по этим формулам, варьируется от
до 2.
И ещё по этим рисункам можно догадаться, что выбор положения точки O1 на окружности, вообще-то, произволен! То есть эти две формулы не единственные.
Пока мы делали выкладки через тригонометрию, мы вроде как рассмотрели два самых очевидных варианта: домножали и делили наш кватернион поворота на синус и на косинус половинного угла.
На рисунках мы видим: в первом случае ничего крутить не пришлось, во втором мы получили поворот на 90°.
Так что возникло предчувствие: если мы будем умножать и делить наш кватернион на cos(α/2+φ) с произвольным углом φ, мы получим целое семейство таких формул.
Давайте попробуем!
Мы пытаемся получить кватернион
(учитывая, что у него используется только две компоненты из 4, здесь бы и комплексное число подошло, но для них нехарактерно использовать половинные углы. Так что уж пускай кватернион, это у нас ликбез по кватернионам или где!?)
Давайте домножим и поделим его на cos(α/2-π/4):
Раскроем этот косинус:
И подставим это в наш кватернион. Множитель в начале просто обозначим как k, он не так уж интересен, отнормировать мы всегда успеем:
Каждое слагаемое выразим через двойной угол:
Двойки в знаменателе пока что тоже внесём в константу k, и получим:
И теперь, вспоминая, что x=cos α, y=sin α, и получаем вторую формулу для нашего улучшенного метода:
Надо только не забыть, что у нас дополнительно возникало деление на корень из двух, т.е выражение в скобках имеет бОльшую длину. Максимальное значение достигается для входного вектора (0;1), из него мы получаем (2;2), т.е вектор длиной
Аналогичные выкладки можно проделать, умножив и поделив кватернион на cos(α/2+π/4), и получим:
Теперь, имея 4 разные формулы, мы можем поделить окружность на 4 части, по 90 градусов каждая. Прелесть в том, что поскольку точка O1 вынесена далеко, длина вектора не успевает сильно измениться при движении на ±45°. Скажем, на первом рисунке у нас максимальная длина вектора O1B достигается при нулевом угле и составляет 2. А если повернуться в любую сторону на 45°, длина будет равняться
что составляет ≈0,924 от максимальной длины, уменьшение менее чем на 8%. Так что если теперь ещё подобрать "волшебные константы", чтобы длина болталась в разные стороны вокруг среднего, мы и получим обещанные ±4%!
Это мы и сделали.
Лично для меня было откровением, что начав с полной окружности (с вообще любого угла), ограничившись рассмотрением 4 вариантов, чисто линейно, можно не только получить в точности половинный угол, но и практически "угадать" с длиной! Впрочем, не уверен, что данный метод найдёт применение в моих приборах: пока что у меня абсолютно нет проблем с временем вычислений, а вот память кода хочется экономить, чтобы вместиться в блоки внутренней памяти ПЛИС... Поэтому пока что для меня оптимально даже с "волшебными константами" не заморачиваться - ну проведу я пяток итераций метода Ньютона, вообще не почувствую.
Одно плохо: сейчас мы поделили плоскость на квадранты "по диагонали". Удобнее было бы проверять ещё более простые условия: знаки x, y. Что ж, вывод соответствующих формул оставляем читателю в качестве несложного упражнения :)