ПрограммированиеСтатьиОбщее

Вращение и кватернионы. Сборник рецептов. (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;
        }
    };

Как видим, мы легко можем создавать вращения вокруг осей, в том числе и вокруг базовых. Перемножив три кватерниона, задающих вращение вокруг трех базовых осей, получим кватернион из углов Эйлера. Поскольку углы Эйлера могут интерпретироваться по-разному, в зависимости от очередности поворотов, приведу лишь один из вариантов конвертирования. Если вам понадобится другая очередность поворотов, рекомендую взять код из "Graphics Gems 4" к статье "Euler Angle Conversion, Ken Shoemake". В целом замечу, что конвертирование из Углов Эйлера и обратно неэффективная операция, которую стоит избегать. Конвертирование в Углы Эйлера для неединичных кватернионов так же эффективно, как и единичных.

Страницы: 1 2 3 Следующая »

#кватернионы, #математика

4 июня 2004 (Обновление: 7 дек 2011)

Комментарии [40]