描述座標系{B}相對於參考座標系{A}的姿態有兩種方式。第一種是繞固定(參考)座標軸旋轉:假設開始兩個座標系重合,先將{B}繞{A}的X軸旋轉$\gamma$,而後繞{A}的Y軸旋轉$\beta$,最後繞{A}的Z軸旋轉$\alpha$,就能旋轉到當前姿態。能夠稱其爲X-Y-Z fixed angles或RPY角(Roll, Pitch, Yaw)。html
Roll:橫滾ios
Pitch: 俯仰windows
Yaw: 偏航(航向)ide
因爲是繞固定座標系旋轉,則旋轉矩陣爲($c\alpha$ is shorthand for $\cos\alpha$, $s\alpha$ is shorthand for $\sin\alpha$,and so on.)函數
$$R_{XYZ}(\gamma,\beta,\alpha)=R_Z(\alpha)R_Y(\beta)R_X(\gamma)=\begin{bmatrix}
c\alpha c\beta & c\alpha s\beta s\gamma-s\alpha c\gamma & c\alpha s\beta c\gamma+s\alpha s\gamma\\
s\alpha c\beta & s\alpha s\beta s\gamma+c\alpha c\gamma & s\alpha s\beta c\gamma-c\alpha s\gamma\\
-s\beta& c\beta s\gamma & c\beta c\gamma
\end{bmatrix}$$ui
另外一種姿態描述方式是繞自身座標軸旋轉:假設開始兩個座標系重合,先將{B}繞自身的Z軸旋轉$\alpha$,而後繞Y軸旋轉$\beta$,最後繞X軸旋轉$\gamma$,就能旋轉到當前姿態。稱其爲Z-Y-X歐拉角,因爲是繞自身座標軸進行旋轉,則旋轉矩陣爲:url
$$R_{Z'Y'X'}(\alpha,\beta,\gamma)=R_Z(\alpha)R_Y(\beta)R_X(\gamma)=\begin{bmatrix}
c\alpha c\beta & c\alpha s\beta s\gamma-s\alpha c\gamma & c\alpha s\beta c\gamma+s\alpha s\gamma\\
s\alpha c\beta & s\alpha s\beta s\gamma+c\alpha c\gamma & s\alpha s\beta c\gamma-c\alpha s\gamma\\
-s\beta& c\beta s\gamma & c\beta c\gamma
\end{bmatrix}$$idea
能夠發現這兩種描述方式獲得的旋轉矩陣是同樣的,即繞固定座標軸X-Y-Z旋轉$(\gamma,\beta,\alpha)$和繞自身座標軸Z-Y-X旋轉$(\alpha,\beta,\gamma)$的最終結果同樣,只是描述的方法有差異而已。In gerenal: three rotations taken about fixed axes yield the same final orientation as the same three rotations taken in opposite order about the axes of the moving frame.spa
繞座標軸的屢次旋轉能夠等效爲繞某一轉軸旋轉必定的角度。假設等效旋轉軸方向向量爲$\vec{K}=[k_x,k_y,k_z]^T$,等效旋轉角爲$\theta$,則四元數$q=(x,y,z,w)$,其中:3d
$$\begin{align*}
x &= k_x \cdot sin \frac{\theta}{2}\\
y &= k_y \cdot sin \frac{\theta}{2}\\
z &= k_z \cdot sin \frac{\theta}{2}\\
w &= cos \frac{\theta}{2}
\end{align*}$$
且有$x^2+y^2+z^2+w^2=1$
即四元數存儲了旋轉軸和旋轉角的信息,它能方便的描述剛體繞任意軸的旋轉。
四元數轉換爲旋轉矩陣:
$$R=\begin{bmatrix}
1-2y^2-2z^2 & 2(xy-zw) & 2(xz+yw)\\
2(xy+zw) & 1-2x^2-2z^2 & 2(yz-xw)\\
2(xz-yw)& 2(yz+xw) & 1-2x^2-2y^2
\end{bmatrix}$$
已知旋轉矩陣爲:
則對應的四元數爲:
定義兩個四元數:
四元數加法:
四元數乘法:
四元數的乘法的意義相似於矩陣的乘法,能夠表示旋轉的合成。當有屢次旋轉操做時,使用四元數能夠得到更高的計算效率。
<<Quaternions` (* This loads the package *) Quaternion[2, 1, 1, 3] ** Quaternion[2, 1, 1, 0] ** Quaternion[1, 1, 1, 1] (* Be sure to use ** rather than * when multiplying quaternions *)
計算結果爲:Quaternion[-12, 4, 14, 2]
$$q=\begin{bmatrix}\cos\frac{\gamma}{2}\\ 0\\ 0\\ \sin\frac{\gamma}{2}\end{bmatrix} \begin{bmatrix}\cos\frac{\beta}{2}\\ 0\\ \sin\frac{\beta}{2}\\ 0\end{bmatrix} \begin{bmatrix}\cos\frac{\alpha}{2}\\ \sin \frac{\alpha}{2}\\ 0\\ 0\end{bmatrix}=\begin{bmatrix}
\cos\frac{\alpha}{2}\cos\frac{\beta}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\sin\frac{\gamma}{2}\\
\sin\frac{\alpha}{2}\cos\frac{\beta}{2}\cos\frac{\gamma}{2}-\cos\frac{\alpha}{2}\sin\frac{\beta}{2}\sin\frac{\gamma}{2}\\ \cos\frac{\alpha}{2}\sin\frac{\beta}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\cos\frac{\beta}{2}\sin\frac{\gamma}{2}
\\ \cos\frac{\alpha}{2}\cos\frac{\beta}{2}\sin\frac{\gamma}{2}-\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\cos\frac{\gamma}{2}
\end{bmatrix}$$
$$\begin{bmatrix}\alpha\\ \beta\\ \gamma\end{bmatrix} = \begin{bmatrix}
\arctan\frac{2(q_0q_1+q_2q_3)}{1-2(q_1^2+q_2^2)}\\
\arcsin(2(q_0q_2-q_1q_3))\\
\arctan\frac{2(q_0q_3+q_1q_2)}{1-2(q_2^2+q_3^2)}
\end{bmatrix}$$
$$\begin{bmatrix}\alpha\\ \beta\\ \gamma\end{bmatrix} = \begin{bmatrix}
atan2(2(q_0q_1+q_2q_3),1-2(q_1^2+q_2^2))\\
\arcsin(2(q_0q_2-q_1q_3))\\
atan2(2(q_0 q_3+q_1 q_2),1-2(q_2^2+q_3^2))
\end{bmatrix}$$
θ = ATan(y / x)求出的θ取值範圍是[-PI/2, PI/2];
θ = ATan2(y, x)求出的θ取值範圍是[-PI, PI]。
當 (x, y) 在第一象限, 0 < θ < PI/2
當 (x, y) 在第二象限 PI/2 < θ≤PI
當 (x, y) 在第三象限, -PI < θ < -PI/2
當 (x, y) 在第四象限, -PI/2 < θ < 0
enum RotSeq{zyx, zyz, zxy, zxz, yxz, yxy, yzx, yzy, xyz, xyx, xzy,xzx};
// COMPILE: g++ -o quat2EulerTest quat2EulerTest.cpp #include <iostream> #include <cmath> using namespace std; /////////////////////////////// // Quaternion struct // Simple incomplete quaternion struct for demo purpose /////////////////////////////// struct Quaternion{ Quaternion():x(0), y(0), z(0), w(1){}; Quaternion(double x, double y, double z, double w):x(x), y(y), z(z), w(w){}; void normalize(){ double norm = std::sqrt(x*x + y*y + z*z + w*w); x /= norm; y /= norm; z /= norm; w /= norm; } double norm(){ return std::sqrt(x*x + y*y + z*z + w*w); } double x; double y; double z; double w; }; /////////////////////////////// // Quaternion to Euler /////////////////////////////// enum RotSeq{zyx, zyz, zxy, zxz, yxz, yxy, yzx, yzy, xyz, xyx, xzy,xzx}; void twoaxisrot(double r11, double r12, double r21, double r31, double r32, double res[]){ res[0] = atan2( r11, r12 ); res[1] = acos ( r21 ); res[2] = atan2( r31, r32 ); } void threeaxisrot(double r11, double r12, double r21, double r31, double r32, double res[]){ res[0] = atan2( r31, r32 ); res[1] = asin ( r21 ); res[2] = atan2( r11, r12 ); } void quaternion2Euler(const Quaternion& q, double res[], RotSeq rotSeq) { switch(rotSeq){ case zyx: threeaxisrot( 2*(q.x*q.y + q.w*q.z), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, -2*(q.x*q.z - q.w*q.y), 2*(q.y*q.z + q.w*q.x), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, res); break; case zyz: twoaxisrot( 2*(q.y*q.z - q.w*q.x), 2*(q.x*q.z + q.w*q.y), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, 2*(q.y*q.z + q.w*q.x), -2*(q.x*q.z - q.w*q.y), res); break; case zxy: threeaxisrot( -2*(q.x*q.y - q.w*q.z), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, 2*(q.y*q.z + q.w*q.x), -2*(q.x*q.z - q.w*q.y), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, res); break; case zxz: twoaxisrot( 2*(q.x*q.z + q.w*q.y), -2*(q.y*q.z - q.w*q.x), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, 2*(q.x*q.z - q.w*q.y), 2*(q.y*q.z + q.w*q.x), res); break; case yxz: threeaxisrot( 2*(q.x*q.z + q.w*q.y), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, -2*(q.y*q.z - q.w*q.x), 2*(q.x*q.y + q.w*q.z), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, res); break; case yxy: twoaxisrot( 2*(q.x*q.y - q.w*q.z), 2*(q.y*q.z + q.w*q.x), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, 2*(q.x*q.y + q.w*q.z), -2*(q.y*q.z - q.w*q.x), res); break; case yzx: threeaxisrot( -2*(q.x*q.z - q.w*q.y), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, 2*(q.x*q.y + q.w*q.z), -2*(q.y*q.z - q.w*q.x), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, res); break; case yzy: twoaxisrot( 2*(q.y*q.z + q.w*q.x), -2*(q.x*q.y - q.w*q.z), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, 2*(q.y*q.z - q.w*q.x), 2*(q.x*q.y + q.w*q.z), res); break; case xyz: threeaxisrot( -2*(q.y*q.z - q.w*q.x), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z, 2*(q.x*q.z + q.w*q.y), -2*(q.x*q.y - q.w*q.z), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, res); break; case xyx: twoaxisrot( 2*(q.x*q.y + q.w*q.z), -2*(q.x*q.z - q.w*q.y), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, 2*(q.x*q.y - q.w*q.z), 2*(q.x*q.z + q.w*q.y), res); break; case xzy: threeaxisrot( 2*(q.y*q.z + q.w*q.x), q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z, -2*(q.x*q.y - q.w*q.z), 2*(q.x*q.z + q.w*q.y), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, res); break; case xzx: twoaxisrot( 2*(q.x*q.z - q.w*q.y), 2*(q.x*q.y + q.w*q.z), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z, 2*(q.x*q.z + q.w*q.y), -2*(q.x*q.y - q.w*q.z), res); break; default: std::cout << "Unknown rotation sequence" << std::endl; break; } } /////////////////////////////// // Helper functions /////////////////////////////// Quaternion operator*(Quaternion& q1, Quaternion& q2){ Quaternion q; q.w = q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z; q.x = q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y; q.y = q1.w*q2.y - q1.x*q2.z + q1.y*q2.w + q1.z*q2.x; q.z = q1.w*q2.z + q1.x*q2.y - q1.y*q2.x + q1.z*q2.w; return q; } ostream& operator <<(std::ostream& stream, const Quaternion& q) { cout << q.w << " "<< showpos << q.x << "i " << q.y << "j " << q.z << "k"; cout << noshowpos; } double rad2deg(double rad){ return rad*180.0/M_PI; } /////////////////////////////// // Main /////////////////////////////// int main(){ Quaternion q; // x,y,z,w Quaternion qx45(sin(M_PI/8), 0,0, cos(M_PI/8) ); Quaternion qy45(0, sin(M_PI/8), 0, cos(M_PI/8)); Quaternion qz45(0, 0, sin(M_PI/8), cos(M_PI/8)); Quaternion qx90(sin(M_PI/4), 0,0, cos(M_PI/4) ); Quaternion qy90(0, sin(M_PI/4), 0, cos(M_PI/4)); Quaternion qz90(0, 0, sin(M_PI/4), cos(M_PI/4)); double res[3]; q = qz45*qx45; q.normalize(); quaternion2Euler(q, res, zyx); cout << "Rotation sequence: X->Y->Z" << endl; cout << "x45 -> z45" << endl; cout << "q: " << q << endl; cout << "x: " << rad2deg(res[0]) << " y: " << rad2deg(res[1]) << " z: " << rad2deg(res[2]) << endl << endl; q = qz90*qx90; q.normalize(); quaternion2Euler(q, res, zyx); cout << "Rotation sequence: X->Y->Z" << endl; cout << "x90 -> z90" << endl; cout << "q: " << q << endl; cout << "x: " << rad2deg(res[0]) << " y: " << rad2deg(res[1]) << " z: " << rad2deg(res[2]) << endl << endl; q = qx90*qz90; q.normalize(); quaternion2Euler(q, res, xyz); cout << "Rotation sequence: Z->Y->X" << endl; cout << "z90 -> x90" << endl; cout << "q: " << q << endl; cout << "x: " << rad2deg(res[0]) << " y: " << rad2deg(res[1]) << " z: " << rad2deg(res[2]) << endl; }
上面的代碼存在一個問題,即奇異性沒有考慮。下面看一種特殊的狀況(參考Maths - Conversion Quaternion to Euler):假設一架飛機繞Y軸旋轉了90°(俯仰角pitch=90),機頭垂直向上,此時如何計算航向角和橫滾角?
這時會發生自由度丟失的狀況,即Yaw和Roll會變爲一個自由度。此時再使用上面的公式根據四元數計算歐拉角會出現問題:
$\arcsin(2(q_0q_2-q_1q_3))$的定義域爲$[-1,1]$,所以$(q_0q_2-q_1q_3)\in[-0.5, 0.5]$,當$q_0q_2-q_1q_3=0.5$時(在程序中浮點數不能直接進行等於判斷,要使用合理的閾值),俯仰角$\beta$爲90°,將其帶入正向公式計算出四元數$(q_0,q_1,q_2,q_3)$,而後能夠發現逆向公式中atan2函數中的參數所有爲0,即出現了$\frac{0}{0}$的狀況!沒法計算。
$\beta=\pi/2$時,$\sin\frac{\beta}{2}=\cos\frac{\beta}{2}=0.707$,將其帶入公式中有
$$q=\begin{bmatrix}w\\ x\\ y\\ z\end{bmatrix}
\begin{bmatrix}
0.707(\cos\frac{\alpha}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\sin\frac{\gamma}{2})\\
0.707(\sin\frac{\alpha}{2}\cos\frac{\gamma}{2}-\cos\frac{\alpha}{2}\sin\frac{\gamma}{2})\\
0.707(\cos\frac{\alpha}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\sin\frac{\gamma}{2})\\
0.707(\cos\frac{\alpha}{2}\sin\frac{\gamma}{2}-\sin\frac{\alpha}{2}\cos\frac{\gamma}{2})
\end{bmatrix}=
\begin{bmatrix}
0.707\cos\frac{\alpha-\gamma}{2}\\
0.707\sin\frac{\alpha-\gamma}{2}\\
0.707\cos\frac{\alpha-\gamma}{2}\\
0.707\sin\frac{\alpha-\gamma}{2}
\end{bmatrix}$$
則$\frac{x}{w}=\frac{z}{y}=\tan\frac{\alpha-\gamma}{2}$,因而有
$$\alpha-\gamma = 2\cdot atan2(x,w)$$
一般令$\alpha=0$,這時$\gamma = -2\cdot atan2(x,w)$。能夠進行驗證:當四元數爲(w,x,y,z)=(0.653,-0.271,0.653,0.271)時,根據這些規則計算出來的ZYX歐拉角爲α=0°,β=90°,γ=45°
當俯仰角爲-90°,即機頭豎直向下時的狀況也與之相似,能夠推導出奇異姿態時的計算公式。比較完整的四元數轉歐拉角(Z-Y-X order)的代碼以下:
CameraSpacePoint QuaternionToEuler(Vector4 q) // Z-Y-X Euler angles { CameraSpacePoint euler = { 0 }; const double Epsilon = 0.0009765625f; const double Threshold = 0.5f - Epsilon; double TEST = q.w*q.y - q.x*q.z; if (TEST < -Threshold || TEST > Threshold) // 奇異姿態,俯仰角爲±90° { int sign = Sign(TEST); euler.Z = -2 * sign * (double)atan2(q.x, q.w); // yaw euler.Y = sign * (PI / 2.0); // pitch euler.X = 0; // roll } else { euler.X = atan2(2 * (q.y*q.z + q.w*q.x), q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z); euler.Y = asin(-2 * (q.x*q.z - q.w*q.y)); euler.Z = atan2(2 * (q.x*q.y + q.w*q.z), q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z); } return euler; }
在DirectXMath Library中有許多與剛體姿態變換相關的函數能夠直接調用:
下面的代碼中座標系繞X軸旋轉90°(注意這裏不是按照右手定則的方向,而是沿着座標軸向原點看過去以順時針方式旋轉,所以與傳統的右手定則恰好方向相反),來進行變換:
#include "stdafx.h" #include<iostream> #include <DirectXMath.h> using namespace DirectX; #define PI 3.1415926 int _tmain(int argc, _TCHAR* argv[]) { //-------------------Computes the product of two quaternions. XMVECTOR q1 = XMVectorSet(1, 1, 3, 2); XMVECTOR q2 = XMVectorSet(1, 1, 0, 2); XMVECTOR q3 = XMVectorSet(1, 1, 1, 1); XMVECTOR result = XMQuaternionMultiply(XMQuaternionMultiply(q3, q2), q1); // Returns the product of two quaternions as q1*q2*q3 std::cout << "Quaternion Multiply:" << std::endl; std::cout << XMVectorGetX(result) << "," << XMVectorGetY(result) << "," << XMVectorGetZ(result) << "," << XMVectorGetW(result) << std::endl << std::endl; //------------------Computes a rotation quaternion based on the pitch, yaw, and roll (Euler angles). float pitch = 90.0 * PI / 180.0; // Angle of rotation around the x-axis, in radians. float yaw = 0; // Angle of rotation around the y-axis, in radians. float roll = 0; // Angle of rotation around the z - axis, in radians. result = XMQuaternionRotationRollPitchYaw(pitch, yaw, roll); std::cout << "RPY/Euler angles to Quaternion:" << std::endl; std::cout << XMVectorGetX(result) << "," << XMVectorGetY(result) << "," << XMVectorGetZ(result) << "," << XMVectorGetW(result) << std::endl << std::endl; //-----------------Computes a rotation quaternion from a rotation matrix. float matrix[16] = { 1, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 1 }; XMMATRIX trans(matrix); // Initializes a new instance of the XMMATRIX structure from a sixteen element float array. result = XMQuaternionRotationMatrix(trans); // This function only uses the upper 3x3 portion of the XMMATRIX. std::cout << "Matrix to Quaternion:" << std::endl; std::cout << XMVectorGetX(result) << "," << XMVectorGetY(result) << "," << XMVectorGetZ(result) << "," << XMVectorGetW(result) << std::endl << std::endl; //-----------------Builds a rotation matrix from a quaternion. trans = XMMatrixRotationQuaternion(result); XMFLOAT3X3 fView; XMStoreFloat3x3(&fView, trans); // Stores an XMMATRIX in an XMFLOAT3X3 std::cout << "Quaternion to Matrix:" << std::endl; std::cout << fView._11 << "," << fView._12 << "," << fView._13 << std::endl << fView._21 << "," << fView._22 << "," << fView._23 << std::endl << fView._31 << "," << fView._32 << "," << fView._33 << std::endl << std::endl; //-----------------Computes an axis and angle of rotation about that axis for a given quaternion. float Angle = 0; XMVECTOR Axis; XMQuaternionToAxisAngle(&Axis, &Angle, result); Axis = XMVector3Normalize(Axis); // Returns the normalized version of a 3D vector std::cout << "Quaternion to Axis-Angle:" << std::endl; std::cout << "Axis: " << XMVectorGetX(Axis) << "," << XMVectorGetY(Axis) << "," << XMVectorGetZ(Axis) << std::endl; std::cout << "Angle: " << Angle*180.0 / PI << std::endl << std::endl; //-----------------Builds a matrix that rotates around an arbitrary axis. Angle = 90.0 * PI / 180.0; trans = XMMatrixRotationAxis(Axis, Angle); XMStoreFloat3x3(&fView, trans); // Stores an XMMATRIX in an XMFLOAT3X3 std::cout << "Axis-Angle to Matrix:" << std::endl; std::cout << fView._11 << "," << fView._12 << "," << fView._13 << std::endl << fView._21 << "," << fView._22 << "," << fView._23 << std::endl << fView._31 << "," << fView._32 << "," << fView._33 << std::endl << std::endl; //-----------------Builds a rotation matrix based on a given pitch, yaw, and roll(Euler angles). trans = XMMatrixRotationRollPitchYaw(pitch, yaw, roll); XMStoreFloat3x3(&fView, trans); // Stores an XMMATRIX in an XMFLOAT3X3 std::cout << "RPY/Euler angles to Matrix:" << std::endl; std::cout << fView._11 << "," << fView._12 << "," << fView._13 << std::endl << fView._21 << "," << fView._22 << "," << fView._23 << std::endl << fView._31 << "," << fView._32 << "," << fView._33 << std::endl << std::endl; //-----------------Builds a matrix that rotates around the x - axis. trans = XMMatrixRotationX(Angle); // Angles are measured clockwise when looking along the rotation axis toward the origin. XMStoreFloat3x3(&fView, trans); // Stores an XMMATRIX in an XMFLOAT3X3 std::cout << "Builds a matrix that rotates around the x-axis.:" << std::endl; std::cout << fView._11 << "," << fView._12 << "," << fView._13 << std::endl << fView._21 << "," << fView._22 << "," << fView._23 << std::endl << fView._31 << "," << fView._32 << "," << fView._33 << std::endl << std::endl; return 0; }
結果以下圖所示:
參考:
DirectXMath Library Quaternion Functions
Convert quaternion to euler rotations
Conversion between quaternions and Euler angles
Maths - Conversion Quaternion to Euler
Coordinate Transformations in Robotics—MATLAB
Introduction to Robotics - Mechanics and Control. Chapter 2 Spatial descriptions and transformations