Quaternion to Euler angles conversion: A direct, general and computationally efficient method

Current methods of the conversion between a rotation quaternion and Euler angles are either a complicated set of multiple sequence-specific implementations, or a complicated method relying on multiple matrix multiplications. In this paper a general formula is presented for extracting the Euler angles in any desired sequence from a unit quaternion. This is a direct method, in that no intermediate conversion step is required (no quaternion-to-rotation matrix conversion, for example) and it is general because it works with all 12 possible sequences of rotations. A closed formula was first developed for extracting angles in any of the 12 possible sequences, both “Proper Euler angles” and “Tait-Bryan angles”. The resulting algorithm was compared with a popular implementation of the matrix-to-Euler angle algorithm, which involves a quaternion-to-matrix conversion in the first computational step. Lastly, a single-page pseudo-code implementation of this algorithm is presented, illustrating its conciseness and straightforward implementation. With an execution speed 30 times faster than the classical method, our algorithm can be of great interest in every aspect.


Introduction
When dealing with 3D orientation problems, many different formalisms can be used to describe a given rotation [1], each of which has its own set of advantages and shortcomings. Arguably the most direct representation of a 3D rotation is a matrix R 2 SO(3), where SO(3) is the group of invertible 3 × 3 matrices such that det(R) = 1 and RR T ¼ R T R ¼ I, where I is the identity matrix. These rotation matrices represent direct linear transformations such that, with v 2 R 3 : Apart from being simple to use, a rotation matrix also has the advantage of being continuous, and a simple matrix multiplication can be used to compose rotations: R = R 2 R 1 is the rotation matrix corresponding to a rotation by R 1 followed by a rotation by R 2 . 3D rotation matrices have some numerical shortcomings, however. For example, as many as 9 numbers (and 6 constraints) are required to represent a 3 degree of freedom rotation, and it can be difficult and computationally costly to orthogonalize a rotation matrix numerically [2] (i.e., to check that the matrix has its determinant equal to 1 and its inverse equal to its transpose, which is necessary to compensate for the accumulated floating point errors). However, it is possible to parametrize this rotation matrix with a smaller set of numbers [3]. One of the most usual set of parameters are the Euler angles. The approach consists in decomposing the 3D rotation matrix into the product of three rotations: Where R θ e is a rotation by the angle θ around the axis e, and the consecutive axes are orthogonal (e 1 � e 2 = e 2 � e 3 = 0). The advantages of using Euler angles include the fact that only three numbers have to be stored, and due to their familiarity, they can be more easily understood, which explains why they are still being so widely used, even in cases where other forms of representation may be more appropriate. The use of Euler angles also has several disadvantages. For example, they are discontinuous and it is difficult to directly compose two 3D rotations expressed in Euler angles. Euler angles are also affected by the phenomenon commonly called "gymbal lock": when two axes become aligned, making the system underdetermined, special care has to be taken. In addition, since there are 12 possible axis sequences (24, when considering the difference between "intrinsic" and "extrinsic" rotations), the correct sequence has to be checked in the case of each application. An arguably preferable parametrization are quaternions. A quaternion is a hypercomplex number defined by one real part and three distinct imaginary parts (which can also be regarded as the vector part). When the norm of a quaternion is equal to 1, quaternions are a useful and efficient representation of 3D orientation: they can be composed as easily as rotation matrices, they are continuous, and they are easily constructed from the axis-angle representation. In addition, quaternions can be normalized trivially, which is much more efficient than having to cope with the corresponding matrix orthogonalization problem. For these reasons, most 3D graphical applications and rotation engines carry quaternions under the hood. Besides these advantages, Euler angles are still being preferred by many authors: Euler angles are the most familiar concept to most engineers and researchers. In addition, in the case of many problems in which there exists only one degree of freedom, angles can suffice.
To be able to perform fast calculations with quaternions and at the same time analyze rotations using angles, it might be necessary to have an efficient method of converting the one set of parameters to the other. Calculating the corresponding quaternion (or rotation matrix) for a given set of Euler angles is trivial. Extracting the Euler angles is much harder, however. One of the following two methods has generally been used up to now. The first method consists in adopting a different set of formulas for each possible angle sequence [4], which is difficult to implement and debug. The second method is that described in [5]. SciPy [6], for example, a widely used scientific library for the Python programming language, implements this method. It converts rotation matrices into Euler angles and involves many different matrix multiplications, including the inverse trigonometric functions required, which are naturally computationally costly. In addition, if rotations are stored in the form of quaternions (as is usually the case in many of the 3D rendering software tools dealing with rotations), an additional conversion step from quaternions to rotation matrices is necessary.
Since many robotic, graphic and other high-level applications involve the use of quaternions (even if they are hidden from the user), it can be necessary to have a concise, efficient method for the conversion between quaternions and Euler angles. The direct conversion formula from quaternions to Euler angles presented here requires fewer computational steps and less expensive computational resources. Moreover, this conversion formula is much simpler to implement and debug, making it a great option for any new applications needing to implement this kind of conversions.

Quaternion algebra summary
In this section, the key properties of quaternions are summarized. It is assumed in this work that we are dealing with the classical Hamilton quaternions. Since the definitions concerning quaternion algebra are not perfectly consistent in the literature [7], some of the main notations and definitions used in this study are then presented. Quaternions form a non-commutative division algebra denoted by H, which extends the complex numbers. A quaternion q 2 H consists of four components: Where q r ; q x ; q y ; q z 2 R. All the properties of quaternions can be obtained using its fundamental property, as given by Hamilton: Using the above properties, the product of two quaternions q and p can be expressed by the Hamilton product: For the sake of simplicity, quaternions will be written here as 4 × 1 vectors (with the scalar q r as the first element): is the the imaginary/vector part of q. The Hamilton product between two quaternions in 4-vector form will be denoted by: Defining the conjugate q � ¼ q r À q " # and the absolute value as jqj ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi the inverse q −1 of q is given by: And for any quaternion q: If q is a unit quaternion, which means that |q| = 1 and q −1 = q � , it can be used to represent the rotation between two reference frames. Denoting v A and v B a vector v in frames A and B, respectively, and q ¼ q B A the unit quaternion corresponding to the rotation from A to B: The equivalent rotation matrix is given by: And the equivalent quaternion for a rotation of an angle θ around an axis e is given by:

Formula development
In the section, the formula for the conversion between a quaternion and any of the 6 proper Euler angle sequences is derived, and then an adaptation for the 6 remaining Tait-Bryan sequences is demonstrated.

Singularities.
There are two different singularities in these expressions. When θ 2 = 0, we have s 2 = 0 and θ − is undefined. When θ 2 = π, we have c 2 = 0 and θ + is undefined. In both cases, one degree of freedom is lost and we can argue that θ 1 (or alternatively, θ 3 ) loses its geometrical meaning. We can then either set θ 1 to zero, or keep it fixed in its latest value (for example, when updating an estimator, for the sake of continuity). Defining: Takingŷ 1 to be some constant (zero or otherwise), we can calculate: y 3 ¼ 2 atan2ðb; aÞ Àŷ 1 ; when y 2 ¼ 0

General formula for θ 1 and θ 3 in the absence of singularities.
If θ 2 6 ¼ 0 and θ 2 6 ¼ π/ 2, multiplying z + and z − yields: On similar lines, multiplying z + and the conjugate of z − yields: The angles can then be obtained using: Or, more simply, from Eq 25: Or: It is worth noting that Eq 32 requires fewer operations than Eq 30: only 2 calls to atan2, one addition and one subtraction, but a final wrapping step may be required in order to either keep the angles either in (−π, π] or [0, 2π).

General formulas for calculating θ 2 .
From Eq 24, we know that: And we can use any of the following equivalent formulas obtained directly from the definition: Where the factor n 2 = a 2 + b 2 + c 2 + d 2 = |q| 2 can be ignored if the quaternion is already normalized. Using the properties of inverse trigonometric functions, we can also find the following formula, which avoids the need for a square root:

Case 2: Tait-Bryan angles
We now define: Where −π/2 < ϕ 2 < π/2. Again assuming that e, e 0 and e 00 are orthogonal unit vectors and e 00 = ε e × e 0 , where ε = (e × e 0 ) � e 00 = ±1, we define: We note that: Which gives: Where: Where y 0 2 ¼ y 2 þ p=2 and y 0 3 ¼ εy 3 , and: PLOS ONE θ 1 0 // For simplicity, we are settingŷ 1 ¼ 0 Many operations are required to convert a quaternion into a rotation matrix. Using the homogeneous formula from Eq 11, for example, if special care is taken in order not to repeat any operations, we have to perform at least 4 2 = 16 floating point multiplications (all the possible products between two different components of the quaternion, plus all the squares of each component), 4 × 3 = 12 multiplications by 2 and 3 × 3 + 6 = 15 additions/subtractions. This conversion step alone is more than enough to make an algorithm based on [5] much slower than the proposed method. In addition, multiple matrix multiplications also have to be computed. By comparison, our algorithm replaces all the conversions and matrix multiplications by a simple permutation of the quaternion elements and in the case of Tait-Bryan angles, only 5 additional additions/subtractions and possibly a sign change are required.

Results
In this section, a performance comparison between our method and the Shuster method is presented. We adapted the SciPy library in order to compile the algorithm as described in Section 4. A real data set comprising the orientation of a spinning object with 3284 data points was used to compare the efficiency of the two algorithms. The full implementation and data set can be downloaded from [8]. First we noted that both methods give the same results: adding the absolute value of the differences between the two methods in a whole data set gives an error of the order of 10 −12 . The execution times required in our tests for each sequence (and their ratios) are presented in the Table 1. From this test, it can be clearly seen that the method presented here is about 30 times faster.

Conclusion
The Euler angles are still a useful intuitive 3D orientation parametrization. A fast method of conversion to/from any other set of parameters can therefore be of great interest for displaying or analyzing data, for instance. In this study, we therefore developed a general formula for this conversion which is concise, easy to implement and easy to debug. In addition, the fact that our method is about 30 times faster than the method proposed by [5], which required an intermediate conversion into rotation matrices, we believe that our proposed method can be of great interest. This faster execution time also makes this method suitable for use in embedded real time applications such as inertial measurement units (IMUs). We propose that this method could be adopted as the new standard method for converting quaternions into Euler angles, and we are now planning to contributing to several scientific libraries accordingly. Moreover, a possible further development is to generalize this formula for the Davenport angles [9], a generalization of the Euler angles in which any set of distinct non-orthogonal axes are used.