四元数和三维旋转
在最初我们想在程序中旋转物体都会使用欧拉角方法。欧拉角包括偏航角(yaw)ϕ,俯仰角(pitch)θ,滚转角(roll)Φ。如下图所示:
我们每次获取旋转时,需要分别对这三个角进行旋转来获得最终的旋转矩阵。但是欧拉角方法会有一个严重的问题,万向节死锁。产生这个问题的主要原因是,我们在旋转时会设置一个层级,会将这三个角排序,比如yaw->pitch->roll,那么我们旋转yaw会影响pitch和roll轴,旋转pitch会影响poll轴,旋转roll则只影响它自己。因此在旋转roll轴和yaw到同一个平面时,就会产生死锁问题,我们会发现旋转少了一个维度。如下图所示:
我们发现想将箭头向下旋转90度却产生了奇怪的弯曲路径。
四元数
我们为了避免这种情况,可以采用四元数方法来处理万向节死锁问题。四元数,简单了理解就是四维空间的数,其表达式为:
其中a为实数。i,j,k,为虚数。且有:
为什么说四元数的旋转可以避免万向节死锁呢。这是因为四元数旋转方法是采用的轴角旋转。
轴角:
其表示设定中心轴u,将点旋转至点。这样不需要像欧拉角旋转那样分别旋转三次,也不会产生万向节死锁的问题了。
但是为什么四元数这样一个四维的数可以表示三维旋转呢?我们可以从1维开始理解。
1维
假设我们身处一个一维的世界,我们能不能尝试去理解2维的事务呢。
如上图所示,我们可以将2维的单位圆,以-1为投射点,将单位圆上所有的点投至i轴,如上图右边所示。我们可以将单位圆上所有的点映射在一维区间()上,而-1正好为未定义点,即又是正无穷,又是负无穷。值得注意的是这里的i点和-i点在投影前后是不变的(即是数部分为0)。
我们还可以发现,复数和复数相乘,如,可以看作对复数w做相应的空间变换(旋转,这里w为(1,0)),如下图:
我们可以看到单位圆的旋转在一维空间中的相应变化。单位圆逆时针旋转90度。一维中-1的点变换到i的点上。
如上图所示,我们可以吧转化为。并且我们发现若用代替i乘上w,正好为旋转45度,在一维上的表达为在轴i上移动至点。
2维
当我们处在2维世界,想要理解三维的世界时:
如上图所示,我们利用处理2维的思路去映射3维的单位球体。从-1点投射射线,将3维单位球上的点投至ij面。我们发现ij平面上的圆在投影前后也是一样的(实数部分为0)。
和之前提到的两个复数相乘可以看作:将后一个复数做相应的空间变换(旋转)。
如上图所示,我们同样用i乘上3维单位球上的点。我们发现单位球绕j轴旋转了90度。同时2维上的映射也产生了相应的变化,在1的点变换至i处。
我们的i可以转化为。可以发现在2维的图中我们圆在轴(1,0)方向进行变化。
如上图,我们将复数的轴设置为0.78i+0.63j,变换角度时,我们能发现三维中点会在该轴和实数轴形成的平面上旋转,而我们的重点二维图上,我们的点会沿着该轴方向移动变换。
3维
看完了1维和2维的处理,我们应该就能大概了解4维超球映射到3维的情况了。四元数在a=0时,通过我们上述的投影,和三维的单位球上的点是一致的,而0<a<1,投影到单位球内,a>1投影到单位球之外,如下图:
那么两个四元数相乘,看作其中一个对另一个四元数做的空间变换(旋转)是这么样的呢?我们可以把该旋转变换看成是两个独立的二维旋转,它们相互垂直,同时发生,但只能在四维的空间中存在,如下图所示:
该图中,作为作用的复数为i,因此也印证了:
,
我们再取作用复数q为,观察映射成的三维球体:
我们可以把q=-0.71+0.71i+0.00j+0.00k变换为。我们发现我们映射后的球体正好沿着轴(1.0i+0.0j+0.0k)移动且放大,该轴正好是我们轴角的轴u,并以u为中心轴旋转了135度。我们这下便可以想到是不只要抵消了球的移动和放大就可以单纯的实现四元数的旋转了。因此我们可以使用去右乘该变换,最后的结果为,但是这里我们发现球体的平移和缩放被抵消了,但是缺旋转了两倍的角度,这也是为什么我们最后选择旋转为。如下图:
更详细的四元数可视化可以去看3blue1brown的视频讲解。
欧拉角2四元数
如果要使在欧拉角的情况下避免万向节死锁,我们可以利用欧拉角获取四元数的值,然后利用四元数构造旋转矩阵,避免死锁。
我们首先定义两个四元数,。其中:
,u为前文的中心轴。
写成矩阵形式:
我们便获得了左乘矩阵L(q),同理可以获得右乘矩阵:
因此我们的四元数选择变换可以写成L(q)R()p。我们旋转矩阵为:
因为四元数的第一行第一列为四维空间处理实数的变换,所以三维旋转只需要矩阵的右下角3*3部分:
该矩阵在得知轴角的情况下可以直接计算。但如果提供的是欧拉角信息。我们则需要将欧拉角转换为四元数,再计算该矩阵。
在得知偏航角(yaw)ϕ,俯仰角(pitch)θ,滚转角(roll)Φ,旋转顺序为Z(ϕ),Y(θ),X(Φ)时,我们可以计算:
,
,
,
,
注意公式顺序要与旋转顺序一致,但是该公式适用于笛卡尔坐标系,在directX和OpenGL中的坐标系并不适用,他们X轴变为Z轴,Y轴变为X轴,Z轴变为Y轴,因此有公式:
,
,
,
,
再代入旋转矩阵即可。
我们代码:
glm::mat4 culQuaternion(float pitch,float yaw,float roll){
float a=cos(roll/2.0f)*cos(pitch/2.0f)*cos(yaw/2.0f)+sin(roll/2.0f)*sin(pitch/2.0f)*sin(yaw/2.0f);
float b=cos(roll/2.0f)*cos(pitch/2.0f)*sin(yaw/2.0f)+sin(roll/2.0f)*sin(pitch/2.0f)*cos(yaw/2.0f);
float c=cos(roll/2.0f)*sin(pitch/2.0f)*cos(yaw/2.0f)-sin(roll/2.0f)*cos(pitch/2.0f)*sin(yaw/2.0f);
float d=sin(roll/2.0f)*cos(pitch/2.0f)*cos(yaw/2.0f)-cos(roll/2.0f)*sin(pitch/2.0f)*sin(yaw/2.0f);
glm::mat4 rotateT=glm::mat4(glm::vec4(1.0-2*c*c-2*d*d,2*b*c+2*a*d,2*b*d-2*a*c,0),
glm::vec4(2*b*c-2*a*d,1.0-2*b*b-2*d*d,2*a*b+2*c*d,0),
glm::vec4(2*a*c+2*b*d,2*c*d-2*a*b,1.0-2*b*b-2*c*c,0),
glm::vec4(0,0,0,1.0));
return rotateT;
}
对应的旋转代码为:
glm::mat4 rotO=glm::mat4(1.0),rotQ=glm::mat4(1.0);
rotO=glm::rotate(rotO, OLpitch, glm::vec3(0.0,1.0,0.0));
rotO=glm::rotate(rotO, OLyaw, glm::vec3(1.0,0.0,0.0));
rotO=glm::rotate(rotO, OLroll, glm::vec3(0.0,0.0,1.0));
rotQ=culQuaternion(OLpitch,OLyaw,OLroll);
完成效果如下图:
上面的箱子为4元数计算的旋转矩阵的效果,下面的箱子为欧拉旋转计算的矩阵产生的效果。当前效果完全一致,但如果出现万向节死锁的情况,下面箱子欧拉旋转产生的效果就会出现问题了。
四元数插值
如果类似动画创建关键帧,给定两个角度,怎么能使中间的帧自动插值出正确的旋转矩阵呢?四元数球面插值作为一种非常好的解决方案。如下图所示:
左图显示了线性插值的结果,明显是不均匀的,而四元数球面插值则有非常好的效果。Shoemake(1985)最初提出了四元数的球面线性插值,其中给出了两个四元数q1和q2,而t〜[0,1]是要在它们之间进行插值的参数值:
Blow(2004)提出了一种理解该函数的直观方法。我们观察下图:
我们可以利用已知的向量v1,v2求出和v1垂直的,然后利用和v1以及求出。公式如下:
我们代码就非常简单:
glm::vec4 Slerp(float t,const glm::vec4& q1,const glm::vec4& q2){
float cosTheta=glm::dot(q1, q2);
if (cosTheta>0.9995) {
return glm::normalize((1-t)*q1+t*q2);
}else{
float theta=glm::acos(glm::clamp(cosTheta, -1.0f, 1.0f));
float thetap=theta*t;
glm::vec4 qperp=glm::normalize(q2-q1*cosTheta);
return q1*glm::cos(thetap)+qperp*glm::sin(thetap);
}
}
这里值得注意的是,避免非常小,引起错误,我们令大于0.9995时,使用线性插值,下图实现了在q1是在pitch的45度,q2在pitch的75度之间自动插值的效果: