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

Автор: | 18.12.2018

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

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

Положение локальной СК при вращении может быть определено:

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

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

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

  

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

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

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

  • оси  Z  (угол alpha), одновременно и  z ГСК;
  • оси X  (угол beta);
  • оси Z  (угол gamma).

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

  • оси z (угол gamma);
  • оси y (угол угол beta);
  • оси z (угол alpha).

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

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

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

Примеры  реализации математических зависимостей для пересчета координат из локальной в глобальную СК и наоборот рассмотрены в статьях:

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

Углы Эйлера это комбинация (композиция) вращений вокруг базовых осей. Существует еще один, простой, способ задания вращения. Этот способ можно назвать «смесь» вращений вокруг базовых координатных осей, или просто вращение вокруг произвольной фиксированной оси. Три компоненты описывающие вращение образуют вектор, лежащий на оси, вокруг которой и поворачивается объект. Обычно хранят ось вращения в виде единичного вектора и угол поворота вокруг этой оси в радианах или градусах (Axis Angle). Выбрав подходящую ось и угол можно задать любую ориентацию объекта. В некоторых случаях удобно хранить угол вращения и ось в одном векторе. Направление вектора в этом случае совпадает с направлением оси вращения, а длина его равна углу поворота. В физике, таким образом, хранят угловую скорость. Вектор, совпадающий направлением с осью вращения и длиной представляющей скорость в радианах в секунду.

Кватернионы

В дополнение к Эйлеровым углам и матричным преобразованиям, существует ещё один, очень удобный способ представления вращений. Кватернион (как это и видно по названию) представляет собой набор из четырёх чисел, которые могут быть представлены как вектор (X,Y,Z) единичной длины  (англ. rotation axis) и угол вращения вокруг этого вектора (англ.rotation angle):

 

Кватернион моно рассматривать  как четыре действительных числа, например как 4d вектор или 3d вектор и скаляр.

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

Как же хранят вращение в кватернионе? Первые три компонента представляют вектор, лежащий на оси вращения, причем длина вектора зависит от угла поворота. Четвертый компонент зависит только от величины угла поворота. Зависимость довольно простая — если взять единичный вектор V за ось вращения и угол alpha за вращение вокруг этой оси, тогда кватернион представляющий это вращение можно записать как:

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)

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

Для понимания того, как хранит вращение кватернион, вспомним про двумерные вращения. Вращение в плоскости можно задать матрицей 2×2, в которой будут записаны косинусы и синусы угла поворота. Можно представить, что кватернион хранит комбинацию оси вращения и матрицы половины поворота вокруг этой оси.

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

Кватернион удобно рассматривать как 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]. Он описывает нулевой поворот (по аналогии с единичной матрицей), и не изменяет другой кватернион при умножении.

Скалярное произведение (inner product)

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]

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

Многие графические API в работают только с матрицами, некоторые используют и кватернионы. Мы рассмотрим в первую очередь  конвертирование кватерниона в матрицу поворота. Часто для задания вращений используют только кватернионы единичной длины, но это не обязательно и иногда даже не эффективно. Разница между конвертированием единичного и неединичного кватернионов составляет около 6-ти умножений и 3-х сложений, зато избавит во многих случаях от необходимости нормировать (приводить длину к 1) кватернион. Если кусок кода критичен по скорости и вы пользуетесь только кватернионами единичной длины тогда можно воспользоваться фактом что норма его равна 1.

    /**
      set the rotation to matrix
    */
    inline void to_matrix( matrix33& m  )const  {
        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] );
        }
    }

Преобразования между кватернионом и 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(vector3& axis, float& angle)const {
        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;
        }
    };

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

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

Вращение и кватернионы. Сборник рецептов

 

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

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

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