超级无敌杀人王教你怎么做两种坐标系十二种顺序的旋转矩阵与欧拉角之间转换!

852 阅读8分钟

我是超级无敌杀人王啊啊啊!!!

大家好,我是超级无敌杀人王DEF --Darth_Eternalfaith。

  • 欧拉角是在计算机图形学中的一种很常用的表示旋转的方式。虽然大部分情况下我们只需要用欧拉角转换成旋转矩阵,但是也有将矩阵转换成欧拉角的情况。这篇文章将带你用较少的代码把旋转矩阵转换成欧拉角。
  • 注意: 大多数情况下软件并不用提供完整的12种旋转顺序和欧拉角之间的转换,而是提供固定的顺序和固定的算法,那样运算效率更高,跑得更快!本文只是抛砖引玉,如果有更好的算法非常欢迎提出。

本文中的矩阵表示为:

    [ m11,   m12,   m13 ]
    [ m21,   m22,   m23 ]
    [ m31,   m32,   m33 ]

实际数据结构为线性结构

/** 矩阵
 * 矩阵的数据类型为1维线性表:
 * ```
 * [1, 2, 3]
 * [4, 5, 6]  >>  [1,2,3,4,5,6,7,8,9]
 * [7, 8, 9]
 * ```
 */
class Matrix extends CONFIG.VALUE_TYPE{...}

欧拉角

欧拉角表示在3d空间坐标系中绕坐标轴旋转的量,而根据旋转顺序的不同,最终呈现的旋转变换也会不同。

本文的旋转轴向参照物为世界坐标系,轴向不会变动。旋转方向为正视坐标轴正方向的顺时针旋转。
将欧拉角转换为旋转矩阵非常简单,就是对应轴向的旋转矩阵按顺序做矩阵乘法。我就不贴代码出来了,因为确实太简单了。

3D对应轴向的旋转也可以看作是对应平面的2D旋转。作为矩阵也是如此,如2D旋转矩阵(如下)

此处2D旋转矩阵坐标系为右手坐标系
mat_rotate = [+cos, +sin]
             [-sin, +cos]

而转换为3D矩阵则是这样(以下3d坐标系使用左手坐标系,右手坐标系的旋转矩阵是这些的转置,也就是只有sin的符号不同)

下面的东西使用了缩写: c: cos, s: sin


Z(xy)=  [ +cz,      -sz,        0 ]
        [ +sz,      +cz,        0 ]
        [ 0,        0,          1 ]

X(yz)=  [ 1,        0,          0 ]
        [ 0,        +cx,      -sx ]
        [ 0,        +sx,      +cx ]

Y(xz)=  [ +cy,      0,        +sy ]
        [ 0,        1,          0 ]
        [ -sy,      0,        +cy ]

欧拉角旋转顺序

选择两个轴向

选择两个不同的轴旋转,有六种配对(zx,zy,xy,xz,yz,yx),做乘法后:


XY= [ cy,           0,              sy      ]
    [ sx*sy,        cx,             -sx*cy  ]
    [ -cx*sy,       sx,             cx*cy   ]

XZ= [ cz,           -sz,            0       ]
    [ cx*sz,        cx*cz,          -sx     ]
    [ sx*sz,        sx*cz,          cx      ]

YX= [ cy,           sy*sx,          sy*cx   ]
    [ 0,            cx,             -sx     ]
    [ -sy,          cy*sx,          cy*cx   ]
    
YZ= [ cy*cz,        -cy*sz,         sy      ]
    [ sz,           cz,             0       ]
    [ -sy*cz,       sy*sz,          cy      ]
    
ZX= [ cz,           -cx*sz,         sz*sx   ]
    [ sz,           cx*cz,          -cz*sx  ]
    [ 0,            +sx,            cx      ]

ZY= [ cz*cy,        -sz,            cz*sy   ]
    [ cy*sz,        cz,             sz*cy   ]
    [ -sy,          0,              cy      ]

选择三个轴向

而选择三个轴向的话将有十二种顺序;分别为泰勒布莱顿欧拉角(xyz, xzy, yxz, yzx, zxy, zyx)、典型欧拉角(xyx, xzx, yxy, yzy, zxz, zyz)

由于全部推导实在太长了,所以我只拿xyz和xyx来举例。而它们都是通用的!如果想看各种旋转矩阵的结果,请看文献: www.geometrictools.com/Documentati…

XYZ =   [  cy*cz,               -cy*sz,                     sy       ]
        [  cz*sx*sy+cx*sz,       cx*cz-sx*sy*sz,           -cy*sx    ]
        [ -cx*cz*sy+sx*sz,       cz*sx+cx*sy*sz,            cx*cy    ]


X0YX1 = [  cy,                  sy*sx1,                     sy*cx1                 ]
        [  sy*sx0,              cx0*cx1 - cy*sx0*sx1,      -cy*cx1*sx0 - cx0*sx1   ]
        [ -sy*cx0,              cx1*sx0 + cy*cx0*sx1,       cy*cx0*cx1 - sx0*sx1   ]

推导

6CC7FD5757C8DC338F5ED0081F3FDA28.jpg 由于我推导的时候过程全是草稿,所以就文章里干脆不写过程了。

结论 和 代码实现

可以看到旋转顺序为xyz的旋转矩阵的一次项为m13,旋转顺序为xyx的矩阵一次项为m11。
这个值的位置是如何确定的呢?其实非常简单!
它是第一个旋转矩阵(记作M1)和第三个矩阵(记作M3)决定的。
把M1的空省行记做v,M3的空省列记作u,则m[v][u]即为第二个旋转顺序的三角函数值。

判断值是cos还是sin也很简单,只需要判断u==v即可。而判断负号项位置的可以用未组合时矩阵的-sin位置,此时可以建立映射表表示左手坐标系和右手坐标系的符号位置和空行的uv

    /** @type {int[]} 对应轴向[z,x,y] (BPH)的旋转矩阵的空行/列的uv*/
    const _MAPPING__EulerAngles_NULL_UV_MATRIX=new Int8Array([2,0,1]);
    
    /** @type {int[]} 左手坐标系的 对应轴向[z,x,y] (BPH)的旋转矩阵的-sin的下标*/
    const _MAPPING__EulerAngles_INVERSE_INDEX_MATRIX__LEFT=new Int8Array([3,7,2]);
    
    /** @type {int[]} 右手坐标系的 对应轴向[z,x,y] (BPH)的旋转矩阵的-sin的下标*/
    const _MAPPING__EulerAngles_INVERSE_INDEX_MATRIX__RIGHT=new Int8Array([1,5,6]);

顺便把获取cos和sin的函数封装一下


    /** 获取atan2时使用的数据 
     * @param {List_Value} out         输出对象
     * @param {int}        op_static   静指针(下标)
     * @param {int}        op_move     动指针(下标)
     * @param {List_Value} mat         矩阵数据
     * @param {int[]}      map_inverse 坐标系-sin映射表
     * @param {boolean}    uv_o_vu     静指针-动指针 是u-v(true)还是v-u(false)
     * @param {int}        cos_uv      cos的uv
     * @param {boolean}    f_inverse   是否使用置负
     */
    function load_MK__EulerAngles_setup_Matrix(out,op_static,op_move,mat,map_inverse,uv_o_vu,cos_uv,f_inverse){
        var p=op_move,i,j;
        --p
        do{ // out[1]:cos, out[0]:sin
            i=uv_o_vu?Matrix.get_Index(3,op_static,p):Matrix.get_Index(3,p,op_static);
            j=p===cos_uv?0:1;
            out[j]=((~map_inverse.indexOf(i))&&(f_inverse))?-mat[i]:mat[i];
            p===0?p=2:--p;
        }while(p!==op_move);
    }

修正符号之后执行反三角函数asin或acos即可得到第二旋转顺序的旋转弧度。

而二次项是在m[v][u]所在的 行/列 上,记作 行m[v][^u] 和 列m[^v][u]。
第一个旋转顺序在列上,第三旋转顺序在行上,而它们都有一个共同的系数(第二旋转的cos值)(记作cosR2);此时我们只需要找出cos和sin再修正符号,就可以执行反三角函数atan2计算出旋转弧度。

对于泰勒布莱顿欧拉角的情况,我们可以继续使用u==v找到cos值,而负号项位置的正好就是未组合时矩阵的-sin位置;

但是对于典型欧拉角就不能如此了,对此我也没有找到可以用算式计算出cos值的uv的办法;另外符号也会由于两个-sin相遇变成正号;只好照着六种顺序又建立了映射表:

/** @type {int[]} 典型欧拉角的cos值uv值映射 */
const _MAPPING__PROPER_EULER_ANGLE_COS=new Int8Array([
    -1,
    1, // [0,1,0] 01
    0, // [0,2,0] 10
    1, // [1,0,1] 10
    -1,
    2, // [1,2,1] 01
    0, // [2,0,2] 01
    2  // [2,1,2] 10
]);
/** @type {int[]} 典型欧拉角是否使用置负的uv值映射 */
const _MAPPING__PROPER_EULER_ANGLE_INVERSE=new Int8Array([-1,0,1,1,-1,0,0,1]);

万向锁

众所周知,万向锁是由于第二旋转顺序正好为特殊值导致第一旋转顺序与第三旋转顺序重合导致的。只需要判断第二旋转顺序对应的矩阵上的值是为1即可,因为泰勒布莱顿欧拉角的万向锁条件是第二顺序旋转值为直角,而经典欧拉角的万向锁的对应条件是平角;正好对应cos和sin的值为1!

这时候我们只需要套用两个旋转轴的矩阵即可,第三顺序或第一顺序的旋转量直接视作0,或者计算完成后对半分

完整的使用矩阵生成欧拉角的代码


    /** 使用矩阵生成欧拉角
     * @param  {Matrix_3}     mat         仅做过旋转变换的矩阵
     * @param  {int[]}       [_axis]    创建旋转矩阵时的乘法顺序 [z,x,y] 默认为 [0,1,2] (BPH)(zxy)
     * @param  {List_Value}  [_out]      接收数据的对象
     * @return {EulerAngles} 修改并返回 out, 或返回一个新的欧拉角
     */
    static create_EulerAngles__Matrix(mat,_axis,_out){
        var axis=_axis||[0,1,2];
        var u,v,index,cos_uv,flag_inverse;
        var acs=axis[0]===axis[2]?acos:asin;
        var sin_axis1;
        var sin_cos_axis0 = new Vector(2),
            sin_cos_axis2 = new Vector(2);
        var rtn=_out||new EulerAngles(3);
        var map__inverse=CONFIG.AXIS?  _MAPPING__EulerAngles_INVERSE_INDEX_MATRIX__LEFT:_MAPPING__EulerAngles_INVERSE_INDEX_MATRIX__RIGHT;
        
        v=_MAPPING__EulerAngles_NULL_UV_MATRIX[axis[0]];
        u=_MAPPING__EulerAngles_NULL_UV_MATRIX[axis[2]];
        index=Matrix.get_Index(3,u,v);
        sin_axis1= ~map__inverse.indexOf(index)?-mat[index]:mat[index];

        if(approximately(abs(sin_axis1),1)){ // Euler Angles Lock
            v=_MAPPING__EulerAngles_NULL_UV_MATRIX[axis[0]];
            u=_MAPPING__EulerAngles_NULL_UV_MATRIX[axis[1]];
            load_MK__EulerAngles_setup_Matrix(sin_cos_axis0,u,v,mat,map__inverse,true,u,true);
            rtn[0]=atan2(sin_cos_axis0[1],sin_cos_axis0[0]);
            rtn[1]=axis[0]===axis[2]?DEG_180:DEG_90;
            rtn[2]=0;
        }else{ // default
            if(axis[0]===axis[2]){
                index=axis[0]*3+axis[1];
                cos_uv=_MAPPING__PROPER_EULER_ANGLE_COS[index];
                flag_inverse=_MAPPING__PROPER_EULER_ANGLE_INVERSE[index];
                load_MK__EulerAngles_setup_Matrix(sin_cos_axis0,u,v,mat,map__inverse,true ,cos_uv,!flag_inverse);
                load_MK__EulerAngles_setup_Matrix(sin_cos_axis2,v,u,mat,map__inverse,false,cos_uv,flag_inverse);
            }else{
                load_MK__EulerAngles_setup_Matrix(sin_cos_axis0,u,v,mat,map__inverse,true ,u,true);
                load_MK__EulerAngles_setup_Matrix(sin_cos_axis2,v,u,mat,map__inverse,false,v,true);
            }
            rtn[0]=atan2(sin_cos_axis0[1],sin_cos_axis0[0]);
            rtn[1]=acs(sin_axis1);
            rtn[2]=atan2(sin_cos_axis2[1],sin_cos_axis2[0]);
        }
        return rtn;
    }
    

测试用例:

function test_rotate(){
    var f=false;
    var m;
    var axis=[2,1,2];
    var list_deg=[123,180,321];
    var vec=new Vector([1,1,1]);
    var temp_v1,
        temp_v2,
        temp_v3,
        temp_v4;

    console.log("顺序",axis);
    console.log("角度",list_deg);
    console.log("向量",vec);
    
    do{
        console.log("\n\n\n");
        if(NML.CONFIG.AXIS===0)console.log("当前为左手坐标系");
        if(NML.CONFIG.AXIS===1)console.log("当前为右手坐标系");

        console.log("欧拉角 生成 旋转矩阵");
        
        m=Matrix_3.create_Rotate__EulerAngles(list_deg.map(value=>value*DEG),axis);
        log_m3(m);
        console.log("v*m",temp_v1=Matrix.multiplication(vec,m,1,3,3));
        console.log("m*v",temp_v2=Matrix.multiplication(m,vec,3,3,1));
        
        var ea=EulerAngles.create_EulerAngles__Matrix(m,axis);
        console.log("旋转矩阵 生成 欧拉角");
        console.log(ea);
        
        console.log("验证矩阵 to 欧拉角结果");
        m=Matrix_3.create_Rotate__EulerAngles(ea,axis);
        console.log("v*m",temp_v3=Matrix.multiplication(vec,m,1,3,3));
        console.log("m*v",temp_v4=Matrix.multiplication(m,vec,3,3,1));

        console.log("判断相等 v*m",Matrix.check_Equal(temp_v1,temp_v3));
        console.log("判断相等 m*v",Matrix.check_Equal(temp_v2,temp_v4));
        
        f=!f;
        NML.CONFIG.AXIS=1;
    }while(f);
}

测试结果:

zxy万向锁情况

012万向锁.png

zxy正常情况

012正常.png

yxy万向锁情况

212万向锁.png

yxy正常情况

212正常.png

最后: 我是超级无敌杀人王啊啊啊!!!