Вращение и кватернионы. Сборник рецептов. (3 стр)
Автор: Михаил Норель
Кватернионы единичной длины
Кватернионы единичной длины (далее единичные кватернионы) обладают рядом полезных свойств. При умножении двух единичных кватернионов получается единичный кватернион. Инверсный единичному кватерниону, является ему сопряженный.
Как мы отмечали раньше длина векторного и скалярного компонентов равны соответственно синусу и косинусу половины угла поворота. Часто удобно рассматривать кватернион как единичный, умноженный на некий скаляр, тогда многие операции можно вывести для единичных кватернионов, а затем дополнить длиной. Некоторые операции выполняются с единичными кватернионами быстрее, а некоторые медленнее, например, конвертирование в матрицу единичного кватерниона быстрее, чем произвольного. Но конвертирование из матрицы в единичный кватернион более медленная операция, чем конвертирование в неединичный. Именно для единичных кватернионов существуют эффективные способы интерполяции. Многие математические библиотеки используют только единичные кватернионы. Возможно, удобно реализовать отдельно классы кватернионов и единичных кватернионов.
Интерполяция
Эффективным применением кватернионов является интерполяция ориентации объекта. Интерполяция вращения часто важна при анимации по ключевым кадрам, экстраполяции и интерполяции объекта по рассчитанным состояниям (например, компенсация задержки данных по сети). Почему бы нам ни воспользоваться для этого углами Эйлера? Они же описывают вращение непосредственно в углах и, соответственно, интерполируя линейно углы, мы получим линейную интерполяцию ориентации? Неприятность заключена в том, что при интерполяции углов Эйлера объект проходит путь, который не является кратчайшим между ориентациями. Кратчайший путь между ориентациями проходит по дуге (вращение вокруг некой фиксированной оси), а интерполяция углов Эйлера будет лежать на дуге только в частных случаях, например, если все вращения происходят по оси совпадающей с базовой (башня танка вращается в одной плоскости). У кватернионов же есть замечательное свойство, сумма двух кватернионов q*s и q'*s' (s ,s' – скаляры) даст кватернион, который как раз и принадлежит плоскости вращения по кратчайшей дуге между q и q'. В этой же плоскости лежит и самая длинная дуга. Для интерполяции нам будет интересна только кратчайшая дуга.
Отметим еще один важный момент – кватернион описывает "ориентацию" объекта в пространстве и "вращение", но вращение ограниченно ±180 градусами. Если мы захотим ключевыми кадрами задать вращение минутной стрелки часов в диапазоне суток, то мы спокойно опишем это вращение углами Эйлера или Axis Angle представлением как 0 градусов в начале диапазона - 360*24 градусов в конце диапазона, а с кватернионом так не получится (как впрочем, и с матрицей ).
Как же происходит интерполяция кватернионов? Промежуточная ориентация объекта между q и q' может быть описана кватернионом, полученным из s*q и s'*q', где s и s' некоторые скаляры. Возникает соблазн вместо s и s' подставить (1-t) и t, и использовать линейную интерполяцию. Но дело в том, что в кватернионе хранятся не углы, а их проекции (синус и косинус половины угла вращения), а линейная интерполяция синуса (косинуса) угла не приведет к линейной интерполяции самого угла. Приведу конечную формулу для интерполяции кватернионов единичной длины:
iq = (q*sin((1-t)*omega) + q'*sin(t*omega))/sin(omega),где cos(omega) = inner_product(q,q')
Такая интерполяция называется сферической линейной интерполяцией SLERP (spherical linear interpolation).
В результате SLERP получим кватернион единичной длины. SLERP интерполяция не проста, и требует много тригонометрических операций. Но другие формы представления вращений неприемлемы для интерполяции даже в таком виде. Например, интерполяция компонентов матрицы может дать матрицу, которая вырождает объект в плоскость.
Есть ли способ интерполировать вращение быстрее? Насколько я знаю, это самый быстрый способ. Но если мы рассматриваем практическое применение интерполяции, рассмотрим несколько практических приемов. Интерполировать вращение можно приближенно, с небольшой ошибкой незаметной для глаза.
Что же произойдет, если мы будем интерполировать кватернионы линейно? Во-первых, интерполированный кватернион будет лежать на той же кратчайшей дуге. Во-вторых, длина кватерниона уже не будет равна длине исходных кватернионов. Наконец, интерполяция останется гладкой. Но угол поворота будет интерполироваться не линейно, а по некоторой непрерывной кривой. То есть, если параметр интерполяции менять плавно, то вращение будет идти неравномерно. Оценим эту ошибку для практического применения.
Максимальная ошибка при использовании линейной интерполяции составляет примерно 4 градуса (в случае, если ближайший поворот между ориентациями равен 180 градусам). Ошибка плавно проявляется и исчезает в зависимости от параметра интерполяции. При параметре интерполяции 0, 0.5, 1 – ошибка нулевая и плавно возрастает примерно при значениях параметра 0.23 и 0.76. 4 градуса это немного для глаза, поэкспериментируйте, и вы увидите, как тяжело заметить разницу в ориентации на 4 градуса. Но, скорее всего, вам не придется иметь дело с ошибкой в 4 градуса. Дело в том, что при уменьшении угла между кватернионами сферическая интерполяция стремится к линейной.
Если угол между ориентациями не превышает 90 градусов (что обычно достаточно для анимации по ключам "key frame animation"), то максимальная ошибка составляет не 2 (как можно ожидать), а 0.5 градуса. Такую ошибку практически нельзя увидеть. Из этого можно сделать вывод, что для задач интерактивной графики можно пользоваться линейной интерполяцией в большинстве случаев. Обратите внимание, что в итоге мы сможем (приближенно) интерполировать вращение без единой тригонометрической функции или извлечения корня. Это почти так же эффективно как интерполирование векторов.
Если вам необходима более точная интерполяция вращений, вы можете воспользоваться приближенными вычислениями с очень малой погрешностью. То же самое можно сказать и о плавной интерполяции сплайнами. Если есть набор ключевых ориентаций можно плавно менять ориентацию без резких изменений в угловой скорости. Существует много различных схем для плавной интерполяции, которые мы не будем рассматривать в этой статье, их легко найти в Интернете.
При реализации интерполяции рекомендую разнести интерполяцию на два этапа – инициализация, которая может быть сделана один раз для пары кватернионов и, непосредственно, интерполяция параметра, это может заметно ускорить процесс. Инициализация использует достаточно ресурсоемкий acos, который и можно вычислить только раз для пары кватернионов. Рассмотрим примеры реализации линейной и сферической линейной интерполяции.
/** class perform linear interpolation (approximate rotation interpolation), result quaternion nonunit, its length lay between. sqrt(2)/2 and 1.0 */ class q_lerper { xxquaternion q1, q2; public: inline void setup_from_unit(const xxquaternion& q1_, const xxquaternion& q2_){ q1 = q1_; q2 = q2_; float inner = inner_product( q1, q2 ); if( inner < 0 ) q2 = -q2; q2 -= q1; } inline void setup( const xxquaternion& q1_, const xxquaternion& q2_){ q1 = q1_; q1.normalize( ); q2 = q2_; q2.normalize( ); setup_from_unit( q1, q2); } inline void interpolate( float t, xxquaternion& result_q){ result_q = q1 + q2*t; } }; //------------------------------------------------------------------- /** Perform Spherical Linear Interpolation of the quaternions, return unit length quaternion */ class q_slerper { xxquaternion q1, q2; float omega; public: inline void setup_from_unit( const xxquaternion& q1_, const xxquaternion& q2_){ q1 = q1_ ; q2 = q2_ ; float cos_omega = inner_product( q1, q2); if( cos_omega < 0 ) { cos_omega = -cos_omega; q2 = -q2; } if( cos_omega > 0.9999f ) cos_omega = 0.9999f; omega = ( float)acos( cos_omega ); float inv_sin_omega = ( float)( 1.0/sin( omega )); q1.scale( inv_sin_omega); q2.scale( inv_sin_omega); } inline void setup( const xxquaternion& q1_, const xxquaternion& q2_){ q1 = q1_; q1.normalize( ); q2 = q2_; q2.normalize( ); setup_from_unit( q1, q2); } inline void interpolate( float t, xxquaternion& result_q ){ result_q = q1*( float)sin( ( 1.0 - t)*omega ) + q2*( float)sin( t*omega ); } };
Преобразование из двух направлений
Очень интересной операцией является создание кватерниона поворота из одного направления в другое по дуге, образованной этими направлениями. Причем это будет кратчайшая дуга между направлениями, называемая "Shortest arc". Это как бы кратчайший путь между точками на сфере, например, из одного города в другой по поверхности земного шара. Еще эффективнее можно получить кватернион, вращающий на две таких дуги. Это кватернион, образованный из векторного и скалярного произведений этих векторов, или кватернион, получившийся из произведения двух кватернионов, образованных из векторов. Поворот на "Shortest arc" можно получить из смеси кватерниона двойного поворота и идентичности (нулевого вращения). Просто сложив кватернион двойного поворота и кватерниона идентичности, получим кватернион поворота по кратчайшей дуге (разумеется, они должны быть одинаковой длины).
Так, например, можно ориентировать камеру в направлении объекта или плоский объект по направлению к камере (Sprite), легко реализуется контроллер вращения с помощью компьютерной мыши (Track Ball).
С помощью "Shortest arc" можно ориентировать ракету в направлении полета, причем ее повороты будут выглядеть естественно (разворот по кратчайшей дуге). Алгоритм очень прост, на каждом шаге берем предыдущий вектор направления. Строим "Shortest arc" от него к текущему направлению и поворачиваем объект на получившийся кватернион. Если мы применим повороты с помощью "Shortest arc" при движении по непрерывной кривой (например, по сплайну), мы реализуем очень полезный метод "parallel transport frame". Этот метод дает нам как бы ориентацию каната протянутого по этой кривой и минимизирует скручивание "twist" каната. Это особенно полезно для создания объектов по заданному трафарету и профилю кривой.
Для решения задачи инверсной кинематики, когда по заданному направлению надо найти необходимый поворот "Shortest arc" придется как нельзя кстати.
/** the shortest arc quaternion that will rotate one vector to another. create rotation from -> to, for any length vectors */ inline void shortest_arc(const vector3& from, const vector3& to ) { vector3 c( from*to ); set( c.x, c.y, c.z, from%to ); normalize( ); // if "from" or "to" is not unit, normalize it w += 1.0f; // reducing angle to halfangle if( w <= TINY ) // angle close to PI { if( ( from.z*from.z ) > ( from.x*from.x ) ) set( 0, from.z, - from.y, w ); //from*vector3(1,0,0) else set( from.y, - from.x, 0, w ); //from*vector3(0,0,1) } normalize( ); }
Композиция вращений
Умножение кватернионов требует меньше операций и может быть эффективнее умножения матриц 3x3. Заметим, что перемножение матриц содержащих только поворот можно оптимизировать. Архитектура процессоров обычно оптимизирована для перемножения матриц и перемножение матриц может быть реализовано на GPU. Для оценки производительности композиции вращений подсчитаем количество арифметических действий (получим приблизительную оценку). Сравните несколько вариантов комбинирования вращения.
– Перемножение двух матриц 3x3 требует 27 умножений и 18 сложений.
– Перемножение двух кватернионов: 16 умножений и 12 сложений. Существует способ умножений кватернионов требующий больше сложений, но меньше умножений. Возможно, на некоторых системах этот способ будет эффективнее.
– Конвертирование матрицы в кватернион требует 7 сложений и 2 условных перехода.
– Конвертирование кватерниона в матрицу требует 16 умножений 15 сложений и 1 деление.
– Конвертирование матрицы вращения в кватернион, умножение на другой кватернион, а затем конвертирование обратно в матрицу потребует 32 умножений 34 сложений 2 условных перехода и 1 деление.
Как мы видим, резкого превосходства нет и все приведенные комбинации вращения сравнимы по скорости выполнения. Целесообразность применения кватернионов зависит от конкретной задачи. Неочевидное преимущество использования кватерниона для представления вращения в том, что кватернион легко избавить от накапливаемых ошибок связанных с неточностью машинного представления чисел. Для того чтобы матрица вращения после многократного перемножения на другие матрицы не накапливала ошибок ее придется ортонормировать. Ортонормирование дорогая операция. Если этого не делать матрица вращения постепенно будет накапливать смещение и масштаб. Такие ошибки можно часто наблюдать в реализации физики и камеры. Типичный пример такой ошибки можно наблюдать при реализации камеры. С кватернионом такого не произойдет если следить за тем чтобы длина кватерниона держалась в пределах точности используемых чисел. Последняя из приведенных комбинаций вращения не будет накапливать ошибок в ни в матрице поворота, ни в длине кватерниона. Приведу пример типичной ошибки. Единичный кватернион конвертируется в матрицу с ошибкой вызванной округлением. Затем матрица конвертируется в единичный кватернион с расчетом что матрица содержит только поворот. Если мы пользуемся предположением что матрица содержит только поворот и кватернион единичный, тогда циклически выполняемые преобразования вызовут резкое приумножение ошибки. Если вы используете математику только единичных кватернионов, обязательно следите за тем чтобы его длина действительно была равна единице.
Физика
Кватернионы очень удобны для решения многих задач связанных с физикой твердого тела. Предлагаю рассмотреть интегрирование вращения твердого тела с использованием кватернионов простейшим методом Эйлера. Скорость вращения тела удобно хранить в обычном векторе, который совпадает по направлению с осью вращения и длина которого пропорциональна угловой скорости (радианы в секунду). Точное решение потребует найти угол на который тело повернулось за отведенный промежуток времени, это длина вектора угловой скорости помноженная на промежуток времени. Затем создать кватернион из угла и оси вращения и перемножить с текущим кватернионом ориентации. Это потребует как минимум один квадратный корень и вычисление синуса и косинуса. Часто применяют приближенное интегрирование. Для этого используют факт что при интегрировании применяют небольшие промежутки времени, и соответственно углы поворота невелики (интегрирование больших углов нестабильно с точки зрения динамики). Синус и косинус аппроксимируют соответственно уравнениями прямой и единицей (первые члены ряда Тейлора) или аппроксимируют несколько членов ряда Тейлора.
q( nv*sin(0.5*a), cos(0.5*a) ) заменяют на q( v*0.5*t, 1 ), где nv - единичный вектор, совпадающий с осью вращения, a - угол на который повернулся объект, v - исходная угловая скорость, t - время интегрирования. Это позволяет обойтись простыми арифметическими действиями. После перемножения двух кватернионов длина результирующего может измениться. Для того чтобы длина кватерниона держалась в пределах точности используемых чисел можно нормировать результирующий кватернион. Можно просто удерживать длину кватерниона в нужных пределах. Например, поделив его на сумму модулей компонентов. Или циклически получать кватернион из матрицы и обратно матрицу из кватерниона. В этом случае мы получим интегрирование вращения без использования тригонометрических функций или квадратных корней.
Исходный код к статье: 20040603.zip.
Полезные применения, не описанные в статье.
Хранение масштаба в кватернионе, интерполяция сплайнами, ограничители сочленений (joint constraints ), инверсная кинематика и многое другое.
Рекомендуемые ресурсы:
1. "Quaternions", Ken Shoemake. (ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/quatut.ps.Z).
2. "Euler Angle Conversion", Ken Shoemake, in "Graphics Gems IV".
3. "Arcball Rotation Control", Ken Shoemake in "Graphics Gems IV".
4. "Fiber Bundle Twist Reduction" Ken Shoemake in "Graphics Gems IV".
5. "Hacking Quaternions", The Inner Product (March 2002) Jonathan Blow, Game Developer Magazine.
Автор выражает благодарность Nom@ad-у за помощь в работе над статьей.
4 июня 2004 (Обновление: 7 дек 2011)