"Ликбез" - часть 5 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
Предлагаю приближённый метод построения единичных векторов, распределённых равномерно ПО ВСЕМ НАПРАВЛЕНИЯМ, то что 2born назвал изотропным ежом. Метод довольно простой концептуально и в вычислительном плане (один квадратный корень на каждый вектор, плюс немного сложений/умножений), хотя полной изотропности всё равно не получается - можно разглядеть, где полюса, а где "экватор".
"Плотность" векторов при таком методе должна отклоняться от средней не более чем на 20% в каждую сторону, что для многих применений должно быть достаточно.
Метод следующий:
1. Начинаем с точек (b,c), равномерно распределённых по квадрату
я выбирал число N - количество точек по каждой координате, и делил на (N-1) интервалов, что даёт N2 точек. Именно это N приведено в картинках выше,
2. Если
то отбрасываем данную точку, иначе - продолжаем,
3. находим корректирующий множитель:
(b^2+c^2))
4. помножаем точку на него:
5. находим значение a:
6. находим пару векторов (x,y,z) и (-x,y,z):
,)
Вот эти векторы, посчитанные для каждой пары (b,c), и будут распределены почти равномерно.
Есть более "грубый" способ, но ещё более простой.
Под катом объясним, откуда всё это взялось, и некоторые характеристики. По большому счёту, мы попытались превратить гиперсферу в шар, и в некоторой степени это удаётся, но с оговорками...
Идея примерно следующая: компьютеру гораздо проще равномерно выбрать точки на сплошной фигуре, например, на диске. Или, если нужно выбирать 3 точки - то на шаре. Можно просто идти с постоянным шагом, а "краевых эффектов" не так уж много, и с ростом количества точек они становятся всё менее и менее заметны. Куда как хуже попытаться "пройти по кривой" с равным шагом, а особенно "пройти по сфере" - тут у нас ОДНИ СПЛОШНЫЕ КРАЕВЫЕ ЭФФЕКТЫ!
Все возможные повороты можно представить в виде шара. В каждую точку шара проводим радиус-вектор. Его длина будет неким образом выражать величину поворота (можно напрямую в радианах, тогда радиус шара будет равен π, можно sin(φ/2), это будет векторная часть кватерниона и радиус шара 1, и пр.), а направление - ось поворота.
Тогда, кажется, остаётся лишь пройтись по всем точкам этого шара, затем повернуть вектор (1;0;0) с помощью каждого из этих поворотов - и результирующие векторы будто бы должны равномерно распределиться по сфере.
Между шаром и половинкой гиперсферы a2+x2+y2+z2=1 вроде как можно провести взаимно-однозначное соответствие. Каждой точке шара соответствует точка на гиперсфере. Но МЕТРИКА У НИХ РАЗНАЯ, что особенно заметно на поверхности этого шара, соответствующей всем поворотам на 180°.
Возьмём для примера кватернион Λ = i, это поворот на 180° вокруг оси X. Чтобы изобразить его на шаре, просто убираем скалярную часть, получается всё то же самое i. Возьмём "соседнюю" точку на этом шаре, 0,99i. С точки зрения метрики шара, расстояние между этими точками: 0,01. Но кватернион, соответствующий этой точке: 0,141+0,99i. Длина "кватерниона разности" составляет примерно 0,141, т.е в 14 раз выше! Ещё хуже будет, если отодвинуться на 0,001, тогда длина "кватерниона разности" получится в 42 раза больше расстояния между точками на шаре. Там производная стремится к бесконечности по мере приближения к поверхности шара.
Поэтому, если равномерно выбрать точки на этом шаре, вблизи его поверхности получатся огромнейшие прорехи!
Наверное, единственный способ избавиться от этой бесконечной производной на границе - это перейти к вектору Родрига:
Там, куда не глянь, с производной всё хорошо, но ВОТ НЕЗАДАЧА - повороты на 180° УБЕЖАЛИ В БЕСКОНЕЧНОСТЬ!
Перебрать всё пространство от минус до плюс бесконечности - задача тоже не самая приятная.
Но мы можем очень сильно облегчить себе жизнь, если ограничимся поворотами лишь до 90°. Тем самым, взяв вектор (1;0;0), мы сможем получить половину сферы. Но если мы распределим по ней точки равномерно, то можем спокойненько её "отразить" - и получить уже полную сферу!
Если мы возьмём шар в виде векторной части кватерниона:
то теперь мы дальше радиуса
не пойдём.
Возьмём точку, находящуюся на этом расстоянии от центра,
если мы отодвинемся на 0,01i к центру, то по метрике шара мы сместились на 0,01, логично. А когда мы достроим скалярную часть кватерниона и посмотрим длину кватерниона разности, мы получим 0,01404, т.е в корень из двух больше. Как видим, здесь различия уже не так существенны, не в 14 и не в 42 раза, а лишь в 1,4.
При этом, если смещаться не радиально (к центру), а тангенциально, т.е добавить 0,01j, то расстояние вообще совпадёт для шара и для гиперсферы. Имеем такую вот анизотропию.
Так что можно ожидать, что если мы равномерно расставим точки по шару с радиусом
, то разница в "плотности векторов" будет составлять около 40%, причём проблема возникнет как раз на поворотах близких к 90°, т.е "на экваторе". Можем даже сказать, что здесь будут "прорехи", поскольку мы-то думали, что выбираем кватернионы равномерно, но на самом деле они разнесены на 40% дальше.
И ещё можно пойти на существенное упрощение, убрав одну из векторных компонент кватерниона, например, по оси X. Дело в том, что если мы будем вращать конкретно вектор (1;0;0), то поворот вокруг оси X не оказывает на него никакого влияния. Как мы показывали в части 16, разложение кватерниона на повороты, можно поворот вокруг оси X "вычленить", и останется кватернион без компоненты при i. Таким образом, из всего шара нам понадобится лишь диск!
И вот уже вырисовывается более-менее рабочий метод:
1. Начинаем с точек (λy,λz), равномерно распределённых по квадрату
2. Если
то отбрасываем данную точку (она не входит в интересующий нас диск), иначе - продолжаем,
3. Достраиваем λy,λz до кватерниона Λ, задав ещё 2 компонента:
4. Поворачиваем вектор (1;0;0) с помощью данного кватерниона. Можно сделать это "честно", а можно воспользоваться 1-й строкой матрицы поворота, полученной из кватерниона, см. часть 5 - практическая реализация поворота, откуда получится:
,)
И чтобы получить полную сферу, возьмём не только вектор (x,y,z), но и (-x,y,z).
Покажем, что получается при квадратной сетке 11×11:
Слева изображены координаты (y,z) всех точек, затем (x,y), а справа - классическая изометрическая проекция, синим изображены оси координат. Видно, что вблизи (1;0;0) (центр первого рисунка) всё выглядит неплохо. То, что точки стоят реже всего посередине, а затем скучиваются к краям - это всё правильно, поскольку это нам "с такого ракурса" кажется, что там совсем небольшое пространство, куда этих точек накинуто выше всякой меры! На самом деле, мы видим проекцию той части поверхности под огромным углом, из-за чего она и ужимается.
А вот посередине (в координатах x,y) вы видим, где основная прореха. Точки явно идут сильно, сильно реже, чем должны. К краям (вверх и вниз) как будто бы всё выравнивается, но это опять впечатление обманчивое - из-за "проекции" они точки должны идти более кучно на краях, а они идут равномерно. Тоже неправильно.
Но есть возможность немного улучшить равномерность. Совсем "выправить" метрику мы не можем при всё желании - как ни бейся, шар и гиперсфера это разные вещи! Чтобы "вместить больше кватернионов", можно шар расширить, но при этом метрика "по касательной", которая до того была верна, станет чрезмерной! Но если сейчас мы имеем правильное значение "по касательной" и различие на 40% "по радиусу", то можно хотя бы прийти к Соломонову решению - пусть "по касательной" метрика будет на 20% выше, а "по радиусу" - на 20% ниже, чем надо!
Ровно этот "манёвр" мы и провели в самом начале. Этот корректирующий множитель,
равен единице в центре диска, и стремится к
на его краях. А диск мы также расширили в корень четвёртой степени из двух.
Теперь, надеюсь, становится понятно то колдунство в начале поста. К чему это приводит на примере N=11, также можно посмотреть в начале поста. Можно заметить, что точки немного "стягиваются" к экватору. Теперь прорехи как таковой почти не остаётся, но чувствуется анизотропия - в этом месте "вдоль экватора" точки идут реже чем "поперёк". Боюсь, с этой анизотропией ничего особенно не сделаешь. Но дело здесь не только в преобразовании диска в сферу (или шара в гиперсферу), но и в самом способе расставить точки равномерно по диску. Начав с квадратной решётки и попросту выкидывая все точки, выходящие за пределы диска, мы получаем "рваные края", а именно они и соответствют "экватору".
Ещё приведём немного картинок. "Вырожденный случай", N=3:
Получается ровно 6 направлений, по граням кубика. Логично...
N=5:
Надо сказать, неплохо, хотя при таком количестве точек ещё не очень понятна их "изотропность".
PS. Код этого алгоритма на Delphi :)
[Spoiler (click to open)]
thr := 1/sqrt(2)+epsilon;
initVal := -1/sqrt(sqrt(2));
step := -2*initVal / (N-1);
mul := sqrt(2)-sqrt(sqrt(2));
index := 0;
for j := 0 to N-1 do begin
z := initVal + step * j;
for i := 0 to N-1 do begin
y := initVal + step * i;
len := y*y + z*z;
if len > thr then continue;
len := 1 - mul * len;
y1 := y*len;
z1 := z*len;
len := y1*y1 + z1*z1;
a := sqrt(1-len);
Hedgehogs[index].X := 1-2*len;
Hedgehogs[index].Y := 2*a*z1;
Hedgehogs[index].Z := -2*a*y1;
inc(index);
Hedgehogs[index].X := -Hedgehogs[index-1].X;
Hedgehogs[index].Y := Hedgehogs[index-1].Y;
Hedgehogs[index].Z := Hedgehogs[index-1].Z;
inc(index);
end;
end;
Здесь epsilon - маленькая добавочка, я взял 10-7, иначе из-за ошибок плавающей арифметики у меня исчезают наиболее простые точки (0;1;0), (0;0;1) и пр., в общем, граница диска ему не принадлежит.
Hedgehogs - массив векторов, который мы последовательно заполняем.
[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
Предлагаю приближённый метод построения единичных векторов, распределённых равномерно ПО ВСЕМ НАПРАВЛЕНИЯМ, то что 2born назвал изотропным ежом. Метод довольно простой концептуально и в вычислительном плане (один квадратный корень на каждый вектор, плюс немного сложений/умножений), хотя полной изотропности всё равно не получается - можно разглядеть, где полюса, а где "экватор".
"Плотность" векторов при таком методе должна отклоняться от средней не более чем на 20% в каждую сторону, что для многих применений должно быть достаточно.
Метод следующий:
1. Начинаем с точек (b,c), равномерно распределённых по квадрату
я выбирал число N - количество точек по каждой координате, и делил на (N-1) интервалов, что даёт N2 точек. Именно это N приведено в картинках выше,
2. Если
то отбрасываем данную точку, иначе - продолжаем,
3. находим корректирующий множитель:
4. помножаем точку на него:
5. находим значение a:
6. находим пару векторов (x,y,z) и (-x,y,z):
Вот эти векторы, посчитанные для каждой пары (b,c), и будут распределены почти равномерно.
Есть более "грубый" способ, но ещё более простой.
Под катом объясним, откуда всё это взялось, и некоторые характеристики. По большому счёту, мы попытались превратить гиперсферу в шар, и в некоторой степени это удаётся, но с оговорками...
Идея примерно следующая: компьютеру гораздо проще равномерно выбрать точки на сплошной фигуре, например, на диске. Или, если нужно выбирать 3 точки - то на шаре. Можно просто идти с постоянным шагом, а "краевых эффектов" не так уж много, и с ростом количества точек они становятся всё менее и менее заметны. Куда как хуже попытаться "пройти по кривой" с равным шагом, а особенно "пройти по сфере" - тут у нас ОДНИ СПЛОШНЫЕ КРАЕВЫЕ ЭФФЕКТЫ!
Все возможные повороты можно представить в виде шара. В каждую точку шара проводим радиус-вектор. Его длина будет неким образом выражать величину поворота (можно напрямую в радианах, тогда радиус шара будет равен π, можно sin(φ/2), это будет векторная часть кватерниона и радиус шара 1, и пр.), а направление - ось поворота.
Тогда, кажется, остаётся лишь пройтись по всем точкам этого шара, затем повернуть вектор (1;0;0) с помощью каждого из этих поворотов - и результирующие векторы будто бы должны равномерно распределиться по сфере.
Между шаром и половинкой гиперсферы a2+x2+y2+z2=1 вроде как можно провести взаимно-однозначное соответствие. Каждой точке шара соответствует точка на гиперсфере. Но МЕТРИКА У НИХ РАЗНАЯ, что особенно заметно на поверхности этого шара, соответствующей всем поворотам на 180°.
Возьмём для примера кватернион Λ = i, это поворот на 180° вокруг оси X. Чтобы изобразить его на шаре, просто убираем скалярную часть, получается всё то же самое i. Возьмём "соседнюю" точку на этом шаре, 0,99i. С точки зрения метрики шара, расстояние между этими точками: 0,01. Но кватернион, соответствующий этой точке: 0,141+0,99i. Длина "кватерниона разности" составляет примерно 0,141, т.е в 14 раз выше! Ещё хуже будет, если отодвинуться на 0,001, тогда длина "кватерниона разности" получится в 42 раза больше расстояния между точками на шаре. Там производная стремится к бесконечности по мере приближения к поверхности шара.
Поэтому, если равномерно выбрать точки на этом шаре, вблизи его поверхности получатся огромнейшие прорехи!
Наверное, единственный способ избавиться от этой бесконечной производной на границе - это перейти к вектору Родрига:
Там, куда не глянь, с производной всё хорошо, но ВОТ НЕЗАДАЧА - повороты на 180° УБЕЖАЛИ В БЕСКОНЕЧНОСТЬ!
Перебрать всё пространство от минус до плюс бесконечности - задача тоже не самая приятная.
Но мы можем очень сильно облегчить себе жизнь, если ограничимся поворотами лишь до 90°. Тем самым, взяв вектор (1;0;0), мы сможем получить половину сферы. Но если мы распределим по ней точки равномерно, то можем спокойненько её "отразить" - и получить уже полную сферу!
Если мы возьмём шар в виде векторной части кватерниона:
то теперь мы дальше радиуса
не пойдём.
Возьмём точку, находящуюся на этом расстоянии от центра,
если мы отодвинемся на 0,01i к центру, то по метрике шара мы сместились на 0,01, логично. А когда мы достроим скалярную часть кватерниона и посмотрим длину кватерниона разности, мы получим 0,01404, т.е в корень из двух больше. Как видим, здесь различия уже не так существенны, не в 14 и не в 42 раза, а лишь в 1,4.
При этом, если смещаться не радиально (к центру), а тангенциально, т.е добавить 0,01j, то расстояние вообще совпадёт для шара и для гиперсферы. Имеем такую вот анизотропию.
Так что можно ожидать, что если мы равномерно расставим точки по шару с радиусом
, то разница в "плотности векторов" будет составлять около 40%, причём проблема возникнет как раз на поворотах близких к 90°, т.е "на экваторе". Можем даже сказать, что здесь будут "прорехи", поскольку мы-то думали, что выбираем кватернионы равномерно, но на самом деле они разнесены на 40% дальше.
И ещё можно пойти на существенное упрощение, убрав одну из векторных компонент кватерниона, например, по оси X. Дело в том, что если мы будем вращать конкретно вектор (1;0;0), то поворот вокруг оси X не оказывает на него никакого влияния. Как мы показывали в части 16, разложение кватерниона на повороты, можно поворот вокруг оси X "вычленить", и останется кватернион без компоненты при i. Таким образом, из всего шара нам понадобится лишь диск!
И вот уже вырисовывается более-менее рабочий метод:
1. Начинаем с точек (λy,λz), равномерно распределённых по квадрату
2. Если
то отбрасываем данную точку (она не входит в интересующий нас диск), иначе - продолжаем,
3. Достраиваем λy,λz до кватерниона Λ, задав ещё 2 компонента:
4. Поворачиваем вектор (1;0;0) с помощью данного кватерниона. Можно сделать это "честно", а можно воспользоваться 1-й строкой матрицы поворота, полученной из кватерниона, см. часть 5 - практическая реализация поворота, откуда получится:
И чтобы получить полную сферу, возьмём не только вектор (x,y,z), но и (-x,y,z).
Покажем, что получается при квадратной сетке 11×11:
Слева изображены координаты (y,z) всех точек, затем (x,y), а справа - классическая изометрическая проекция, синим изображены оси координат. Видно, что вблизи (1;0;0) (центр первого рисунка) всё выглядит неплохо. То, что точки стоят реже всего посередине, а затем скучиваются к краям - это всё правильно, поскольку это нам "с такого ракурса" кажется, что там совсем небольшое пространство, куда этих точек накинуто выше всякой меры! На самом деле, мы видим проекцию той части поверхности под огромным углом, из-за чего она и ужимается.
А вот посередине (в координатах x,y) вы видим, где основная прореха. Точки явно идут сильно, сильно реже, чем должны. К краям (вверх и вниз) как будто бы всё выравнивается, но это опять впечатление обманчивое - из-за "проекции" они точки должны идти более кучно на краях, а они идут равномерно. Тоже неправильно.
Но есть возможность немного улучшить равномерность. Совсем "выправить" метрику мы не можем при всё желании - как ни бейся, шар и гиперсфера это разные вещи! Чтобы "вместить больше кватернионов", можно шар расширить, но при этом метрика "по касательной", которая до того была верна, станет чрезмерной! Но если сейчас мы имеем правильное значение "по касательной" и различие на 40% "по радиусу", то можно хотя бы прийти к Соломонову решению - пусть "по касательной" метрика будет на 20% выше, а "по радиусу" - на 20% ниже, чем надо!
Ровно этот "манёвр" мы и провели в самом начале. Этот корректирующий множитель,
равен единице в центре диска, и стремится к
на его краях. А диск мы также расширили в корень четвёртой степени из двух.
Теперь, надеюсь, становится понятно то колдунство в начале поста. К чему это приводит на примере N=11, также можно посмотреть в начале поста. Можно заметить, что точки немного "стягиваются" к экватору. Теперь прорехи как таковой почти не остаётся, но чувствуется анизотропия - в этом месте "вдоль экватора" точки идут реже чем "поперёк". Боюсь, с этой анизотропией ничего особенно не сделаешь. Но дело здесь не только в преобразовании диска в сферу (или шара в гиперсферу), но и в самом способе расставить точки равномерно по диску. Начав с квадратной решётки и попросту выкидывая все точки, выходящие за пределы диска, мы получаем "рваные края", а именно они и соответствют "экватору".
Ещё приведём немного картинок. "Вырожденный случай", N=3:
Получается ровно 6 направлений, по граням кубика. Логично...
N=5:
Надо сказать, неплохо, хотя при таком количестве точек ещё не очень понятна их "изотропность".
PS. Код этого алгоритма на Delphi :)
[Spoiler (click to open)]
thr := 1/sqrt(2)+epsilon;
initVal := -1/sqrt(sqrt(2));
step := -2*initVal / (N-1);
mul := sqrt(2)-sqrt(sqrt(2));
index := 0;
for j := 0 to N-1 do begin
z := initVal + step * j;
for i := 0 to N-1 do begin
y := initVal + step * i;
len := y*y + z*z;
if len > thr then continue;
len := 1 - mul * len;
y1 := y*len;
z1 := z*len;
len := y1*y1 + z1*z1;
a := sqrt(1-len);
Hedgehogs[index].X := 1-2*len;
Hedgehogs[index].Y := 2*a*z1;
Hedgehogs[index].Z := -2*a*y1;
inc(index);
Hedgehogs[index].X := -Hedgehogs[index-1].X;
Hedgehogs[index].Y := Hedgehogs[index-1].Y;
Hedgehogs[index].Z := Hedgehogs[index-1].Z;
inc(index);
end;
end;
Здесь epsilon - маленькая добавочка, я взял 10-7, иначе из-за ошибок плавающей арифметики у меня исчезают наиболее простые точки (0;1;0), (0;0;1) и пр., в общем, граница диска ему не принадлежит.
Hedgehogs - массив векторов, который мы последовательно заполняем.