zz :https://www.cnblogs.com/21207-iHome/p/6894128.html
- RPY角與Z-Y-X歐拉角
描述坐標系{B}相對於參考坐標系{A}的姿態有兩種方式。第一種是繞固定（參考）坐標軸旋轉：假設開始兩個坐標系重合，先將{B}繞{A}的X軸旋轉$\gamma$，然后繞{A}的Y軸旋轉$\beta$，最后繞{A}的Z軸旋轉$\alpha$，就能旋轉到當前姿態。可以稱其為X-Y-Z fixed angles或RPY角(Roll, Pitch, Yaw)。
Roll:橫滾
Pitch: 俯仰
Yaw: 偏航（航向）
由於是繞固定坐標系旋轉，則旋轉矩陣為（$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}$$
另一種姿態描述方式是繞自身坐標軸旋轉：假設開始兩個坐標系重合，先將{B}繞自身的Z軸旋轉$\alpha$，然后繞Y軸旋轉$\beta$，最后繞X軸旋轉$\gamma$，就能旋轉到當前姿態。稱其為Z-Y-X歐拉角，由於是繞自身坐標軸進行旋轉，則旋轉矩陣為：
$$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}$$
可以發現這兩種描述方式得到的旋轉矩陣是一樣的，即繞固定坐標軸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.
- Axis-Angle與四元數
繞坐標軸的多次旋轉可以等效為繞某一轉軸旋轉一定的角度。假設等效旋轉軸方向向量為$\vec{K}=[k_x,k_y,k_z]^T$，等效旋轉角為$\theta$，則四元數$q=(x,y,z,w)$，其中：
$$\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};
View Code
上面的代碼存在一個問題，即奇異性沒有考慮。下面看一種特殊的情況（參考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中有許多與剛體姿態變換相關的函數可以直接調用：
- 四元數乘法：XMQuaternionMultiply method --Computes the product of two quaternions.
- 旋轉矩陣轉四元數：XMQuaternionRotationMatrix method --Computes a rotation quaternion from a rotation matrix.
- 四元數轉旋轉矩陣：XMMatrixRotationQuaternion method -- Builds a rotation matrix from a quaternion.
- 歐拉角轉四元數：XMQuaternionRotationRollPitchYaw method --Computes a rotation quaternion based on the pitch, yaw, and roll (Euler angles).
- 四元數轉Axis-Angle：XMQuaternionToAxisAngle method --Computes an axis and angle of rotation about that axis for a given quaternion.
- 歐拉角轉旋轉矩陣：XMMatrixRotationRollPitchYaw method --Builds a rotation matrix based on a given pitch, yaw, and roll (Euler angles).
- Axis-Angle轉旋轉矩陣：XMMatrixRotationAxis method --Builds a matrix that rotates around an arbitrary axis.
- 構造繞X/Y/Z軸的旋轉矩陣：XMMatrixRotationX method --Builds a matrix that rotates around the x-axis.(Angles are measured clockwise when looking along the rotation axis toward the origin)
下面的代碼中坐標系繞X軸旋轉90°（注意這里不是按照右手定則的方向，而是沿着坐標軸向原點看過去以順時針方式旋轉，因此與傳統的右手定則剛好方向相反），來進行變換：
View Code
結果如下圖所示:
參考：
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