Углы Эйлера и кватернионы

Автор: | 18.12.2018

Введение
Определение положения СК  углами Эйлера
Axis Angle представление вращения
Кватернионы
Основные операции над кватернионами
Преобразования между кватернионом и Axis Angle представлением
Конвертирование кватерниона в матрицу поворота
Конвертирование матрицы поворота в кватернион
Полезные ссылки

Введение

При определении положения объекта в пространстве обычно рассматривают 2 системы координат (СК) – та, относительно которой определяется положение (глобальная СК), и та, которая меняет положение относительно нее (локальная СК).

Положение локальной СК при вращении обычно определяется одним из 2-х способов:

  • углами Эйлера — три угла поворота относительно осей координат;
  • кватернионом — вектор и один угол поворота в плоскости, которая перпендикулярна вектору.

Кватернион — не единственный способ определения вращения СК  вектором и углом. Кроме него есть  Axis Angle представление вращения.

Преимущества углов Эйлера:

  • Минимально-возможное количество параметров (только 3 угла).
  • Просты в понимании.

Недостатки:

  • Множество альтернативных представлений, зависящих от последовательности применения поворотов (XYZ, ZYX, …).
  • Проблема “блокировки кардана”.

Преимущества кватернионов:

  • Для получения конечного поворота из 2х необходимо меньше вычислительных ресурсов.
  • Интегрирование скорости поворота производится быстрее.
  • Не подвержен проблеме “блокировки кардана”.

Недостатки:

  • Неоптимальность представления с точки зрения количества параметров (4 вместо 3-х).
  • Неочевидность (абстрактность) представления.

Определение положения СК  углами Эйлера

Рассмотрим пример определения положения СК углами Эйлера. По-простому углы Эйлера – углы вращения относительно осей координат.

  

На самом деле немного сложнее. В качестве оси вращения могут быть выбраны любая из осей как глобальной, так и локальной СК. При этом возникают следующие вопросы:

  1. Какие выбрать оси вращения?
  2. В какой последовательности вращать?
  3. В каком направлении вращать?

От выбора осей и последовательности вращения зависит конечный результат. На рисунках отображена следующая последовательность вращения относительно осей ЛСК :

  • оси  Z  (угол alpha);
  • оси X  (угол beta);
  • оси Z  (угол gamma).

То же преобразование можно получить, выполняя вращение  относительно осей ГСК:

  • оси z (угол (gamma+pi/2));
  • оси y (угол угол beta);
  • оси z (угол (-alpha)).

 Последовательность преобразований берется не откуда-то свыше, а определяется чисто на здравом смысле, используя пространственное воображение. Каждый конкретный случай может иметь свои особенности. Чтобы избежать недоразумений руководствуйтесь простым правилом — преобразования не должны искажать натуральную величину параметров (углов).

Умение правильно выбирать последовательность элементарных преобразований помогает в решении множества, задач, связанных с аффинными преобразованиями (см. Примеры аффинных преобразований).

Определение углов Эйлера через вращение относительно осей ГСК  позволяет достаточно просто получить зависимости для пересчета координат из ЛСК в ГСК через матрицы аффинных преобразований.

Также легко можно перейти к зависимостям для пересчета координат из ГСК в ЛСК.

Это выполняется  через перемножение вектора координат точек на результирующую обратную матрицу. Такая матрица получается перемножением обратных матриц в обратном порядке по сравнению с прямым преобразованием. При этом каждая из обратных матриц вращения может быть получена заменой знака угла на противоположный.

Детально с  теоретическими основами аффинных преобразований (включая и вращение) можно ознакомиться в статье Геометрические преобразования в графических приложениях

Примеры аффинных преобразований рассмотрены в статьях:

Axis Angle представление вращения

Выбрав подходящую ось (англ. rotation axis) и угол (англ. rotation angle) можно задать любую ориентацию объекта.

Обычно хранят ось вращения в виде единичного вектора и угол поворота вокруг этой оси в радианах или градусах.

q = [ x, y, z, w ] = [ v, w ]

В некоторых случаях удобно хранить угол вращения и ось в одном векторе. Направление вектора в этом случае совпадает с направлением оси вращения, а его длина равна углу поворота. В физике, таким образом хранят угловую скорость.  Направление вектора совпадает с направлением оси вращения, а длина вектора равна скорости (в радианах в секунду).

Кватернионы

Кватернион (как это и видно по названию) представляет собой набор из четырёх чисел, которые определяют вектор и угол вращения вокруг этого вектора. По сути такое определение ничем не отличается от Axis Angle представления вращения. Отличия лишь в способе представления. Как же хранят вращение в кватернионе?

q = [ V*sin(alpha/2), cos(alpha/2) ]

В кватернионе параметры единичного вектора умножается на синус половины угла поворота. Четвертый компонент — косинус половины угла поворота.

// Угол вращения RotationAngle задан в радианах
x = RotationAxis.x * sin(RotationAngle / 2)
y = RotationAxis.y * sin(RotationAngle / 2)
z = RotationAxis.z * sin(RotationAngle / 2)
w = cos(RotationAngle / 2)

Таблица с примерами значений кватернионов:

Представление вращения кватернионом кажется необычным по сравнению с Axis Angle представлением. Почему параметры вектора умножаются на синус половины угла вращения, четвертый параметр — косинус половины угла вращения, а не просто угол?

Откуда получено такое необычное представление кватерниона детально можно ознакомиться в статье Доступно о кватернионах и их преимуществах  Хотя программисту не обязательно знать эти детали, точно также как и знать, каким образом получены матрицы преобразования пространства. Достаточно лишь знать основные операции с кватернионами, их смысл и правила применения.

Основные операции над кватернионами

Кватернион удобно рассматривать как 4d вектор, и некоторые операции с ним выполняются как над векторами.

Сложение, вычитание и умножение на скаляр.

 q + q' = [ x + x', y + y', z + z', w + w' ]
 q – q' = [ x - x', y - y', z - z', w - w' ]
 qs = [ xs, ys, zs, ws ]

Смысл операции сложения можно описать как «смесь» вращений, т.е. мы получим вращение, которое находится между q и  q’.  Умножение на скаляр на вращении не отражается. Кватернион, умноженный на скаляр, представляет то же самое вращение, кроме случая умножения на 0. При умножении на 0 мы получим «неопределенное» вращение.

Умножение двух кватернионов

qq' = [ vv' + wv' + w'v, ww' – v•v' ]

где vv’ — векторное произведение, v•v’ — скалярное произведение векторов.

Одна из самых полезных операций, она аналогична умножению двух матриц поворота. Итоговый кватернион представляет собой комбинацию вращений — объект повернули на q’, а затем на q (если смотреть из глобальной системы координат).

Кватернионы не коммутативны по умножению, как и матрица вращения: qq’ не равно q’q (это свойство вращений, а не их представления).

Норма и модуль

Следует различать (а путают их часто) эти две операции:

Норма (norm)

N(q) = xx + yy + zz + ww

Модуль (magnitude), или как иногда говорят «длина» кватерниона:

|q| = sqrt( N(q) ) = sqrt( xx + yy + zz + ww )

Сопряжение ( conjugate )

conjugate(q) = [-v, w]

Задает вращение обратное данному. Обратное вращение можно также получить, просто поменяв знак скалярного «w» компонента на обратный.

Инверсный (inverse) кватернион

Существует такой кватернион, при умножении на который произведение дает нулевое вращение и соответствующее тождественному кватерниону (identity quaternion), и определяется как:

 q–1 = conjugate(q)/N(q)

Тождественный кватернион

Записывается как q[0, 0, 0, 1]. Он описывает нулевой поворот (по аналогии с единичной матрицей), и не изменяет другой кватернион при умножении.

Скалярное произведение

q•q' = x*x' + y*y' + z*z' + w*w'

Скалярное произведение полезно тем, что дает косинус половины угла между двумя кватернионами умноженную на их длину. Соответственно скалярное произведение двух единичных кватернионов даст косинус половины угла между двумя ориентациями. Угол между кватернионами это угол поворота из q  в  q’ (по кратчайшей дуге).

Вращение 3d вектора

Вращение 3d вектора v кватернионом q определяется как

V' = qvq–1

причем вектор конвертируется в кватернион как

q = [x, y, z, 0]

и кватернион обратно в вектор как

v = [x, y, z]

Преобразования между кватернионом и Axis Angle представлением

/**
 create a unit quaternion rotating by axis angle representation
 */
 inline void unit_from_axis_angle(const vector3 &axis, const float &angle)
 {
 vector3 v(axis);
 v.norm();
 float half_angle = angle*0.5f;
 float sin_a = (float)sin(half_angle);
 set(v.x*sin_a, v.y*sin_a, v.z*sin_a, (float)cos(half_angle));
 };
 //-----------------------------------
 /**
 convert a quaternion to axis angle representation, 
 preserve the axis direction and angle from -PI to +PI
 */
 inline void to_axis_angle(const vector3 &axis, float &angle) 
 {
 float vl = (float)sqrt( x*x + y*y + z*z );
 if( vl > TINY )
 {
 float ivl = 1.0f/vl;
 axis.set( x*ivl, y*ivl, z*ivl );
 if( w < 0 )
 angle = 2.0f*(float)atan2(-vl, -w); //-PI,0 
 else
 angle = 2.0f*(float)atan2( vl, w); //0,PI 
 }else{
 axis = vector3(0,0,0);
 angle = 0;
 }
 };

Конвертирование кватерниона в матрицу поворота

    /**
      set the rotation to matrix
    */
    inline void to_matrix(const matrix33 &m)  {
        float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
        float s  = 2.0f/norm();  // 4 mul 3 add 1 div
        x2 = x * s;    y2 = y * s;    z2 = z * s;
        xx = x * x2;   xy = x * y2;   xz = x * z2;
        yy = y * y2;   yz = y * z2;   zz = z * z2;
        wx = w * x2;   wy = w * y2;   wz = w * z2;

        m.m[0][0] = 1.0f - (yy + zz);
        m.m[1][0] = xy - wz;
        m.m[2][0] = xz + wy;

        m.m[0][1] = xy + wz;
        m.m[1][1] = 1.0f - (xx + zz);
        m.m[2][1] = yz - wx;

        m.m[0][2] = xz - wy;
        m.m[1][2] = yz + wx;
        m.m[2][2] = 1.0f - (xx + yy);
    }

Конвертирование матрицы поворота в кватернион

/**
 create a unit quaternion by rotation matrix
 martix must contain only rotation (not scale or shear)

For numerical stability we find first the greatest component of quaternion
 and than search others from this one
 */
 inline void unit_from_matrix( const matrix33 &mtx)
 {
 typedef float mtx_elm[3][3];
 const mtx_elm &m = mtx.m;
 float tr = m[0][0] + m[1][1] + m[2][2]; // trace of martix
 if (tr > 0.0f){ // if trace positive than "w" is biggest component
 set( m[1][2] - m[2][1], m[2][0] - m[0][2], m[0][1] - m[1][0], tr+1.0f );
 scale( 0.5f/(float)sqrt( w ) ); // "w" contain the "norm * 4"

}else // Some of vector components is bigger
 if( (m[0][0] > m[1][1] ) && ( m[0][0] > m[2][2]) ) {
 set( 1.0f + m[0][0] - m[1][1] - m[2][2], m[1][0] + m[0][1],
 m[2][0] + m[0][2], m[1][2] - m[2][1] );
 scale( 0.5f/(float)sqrt( x ) );

}else 
 if ( m[1][1] > m[2][2] ){
 set( m[1][0] + m[0][1], 1.0f + m[1][1] - m[0][0] - m[2][2],
 m[2][1] + m[1][2], m[2][0] - m[0][2] ); 
 scale( 0.5f/(float)sqrt( y ) );

}else{
 set( m[2][0] + m[0][2], m[2][1] + m[1][2],
 1.0f + m[2][2] - m[0][0] - m[1][1], m[0][1] - m[1][0] );
 scale( 0.5f/(float)sqrt( z ) );

}
 }

//----------------------------------------------------------------------------
 /**
 create a nonunit quaternion from rotation matrix 
 martix must contain only rotation (not scale or shear)
 the result quaternion length is numerical stable 
 */
 inline void from_matrix( const matrix33 &mtx)
 {
 typedef float mtx_elm[3][3];
 const mtx_elm &m = mtx.m;
 float tr = m[0][0] + m[1][1] + m[2][2]; // trace of martix
 if (tr > 0.0f){ // if trace positive than "w" is biggest component
 set( m[1][2] - m[2][1], m[2][0] - m[0][2], m[0][1] - m[1][0], tr + 1.0f );
 }else // Some of vector components is bigger
 if( (m[0][0] > m[1][1] ) && ( m[0][0] > m[2][2]) ) {
 set( 1.0f + m[0][0] - m[1][1] - m[2][2], m[1][0] + m[0][1],
 m[2][0] + m[0][2], m[1][2] - m[2][1] );
 }else 
 if ( m[1][1] > m[2][2] ){
 set( m[1][0] + m[0][1], 1.0f + m[1][1] - m[0][0] - m[2][2],
 m[2][1] + m[1][2], m[2][0] - m[0][2] ); 
 }else{
 set( m[2][0] + m[0][2], m[2][1] + m[1][2],
 1.0f + m[2][2] - m[0][0] - m[1][1], m[0][1] - m[1][0] );
 }
 }

Полезные ссылки:

 

 

Автор: Николай Свирневский

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *