Euler angles from matrix

From CGAFaq

As discussed in the article on Euler angles, a triple of angles can represent a 3D rotation as a composition of coordinate axis-aligned rotations in 24 different ways.

xyz_s, xyx_s, xzy_s, xzx_s, yzx_s, yzy_s, yxz_s, yxy_s, zxy_s, zxz_s, zyx_s, zyz_s

xyz_r, xyx_r, xzy_r, xzx_r, yzx_r, yzy_r, yxz_r, yxy_r, zxy_r, zxz_r, zyx_r, zyz_r

This would seem to require either choosing one convention only, or writing separate code for each of the special cases. Neither choice is very appealing, and neither is necessary. By choosing matrix indices parametrically, one small piece of code can extract angles according to any system. Likewise, little code is required to generate a matrix.

Contents 1 Symmetries 2 First merger 3 Givens rotation 4 Archetype 5 Consolidation 6 Code 7 References

One caution: Extraction differs critically from generation in that, while a triple of angles generates a unique rotation matrix for a given system, extraction of angles suffers from three kinds of ambiguity no matter which system is used. For example, in system xzx_s some equivalent triples are (135°,60°,−90°) and (−45°,−60°,90°); or (90°,0°,0°) and (30°,0°,60°); or (90°,45°,−90°) and (90°,−315°,270°). Applications wishing to use angle extraction must compensate. Of special interest here, we cannot test the path angles → matrix → angles, but only the path matrix → angles → matrix.

Symmetries

All systems are equivalent to any one, say xyx_s, by reordering and negation. Three independent choices are:

1. choice of frame, such as xzx_s compared to xzx_r
2. choice of last axis, such as xzx_s compared to xzy_s
3. permutation of all axes, classified as:
   • even permutations, such as xzx_s compared to zyz_s
   • odd permutations, such as xzx_s compared to xyx_s

Angle order reversal relates the choice of static axes to rotating axes. For example, when rotations around static (“world”) axes x, z, and y are applied in this order to a column vector, the product of matrices is written R_y(R_z(R_x v)). Compare this to rotations around rotating (“body”) axes y, z, and x applied in this order; again the product of matrices is R_y(R_z(R_x v)). In effect, the only difference between these two systems is that xzy_s lists the angles in the order (θ_x,θ_z,θ_y), while yzx_r lists the angles in the order (θ_y,θ_z,θ_x); otherwise, the relationship between angles and matrix is identical.

Now consider the relationship between xzy_s and zyx_s. These are the same except for a reordering of the coordinate axes by the permutation (x,y,z) ↦ (z,x,y). This is an even permutation; that is, it requires an even number of axis pair swaps, such as z ↔ y then z ↔ x. Let P be the permutation matrix, here

For an even permutation of the axes, the columns of P are also an even permutation of the identity matrix columns; hence P has determinant +1. For an odd permutation of the axes, such as (x,y,z) ↦ (z,y,x), achieved by a single swap x ↔ z, the columns of P are an odd permutation of the identity matrix columns,

hence P has determinant −1. The sign of the determinant indicates that a right-handed coordinate basis becomes left-handed under odd permutations. Thus a positive rotation around the z axis, which turns x to y and turns y to −x, permutes to a negative rotation around the x axis, which turns z to y and turns y to −z.

Common code can extract angles for any permutation by permuting indices; however, because odd permutations reverse the sense of rotation, the extracted angles must then be negated.

First merger

Consider the matrix for xzy_s rotation with angles (θ_x, θ_z, θ_y). Letting c_x =cos θ_x, s_x =sin θ_x, and so on, the matrix is

Similarly, the xzx_s rotation with angles (θ_x, θ_z, θ_x′) is

“Undoing” the left-most rotation (either x′ or y) produces the same effect on both: a zero is produced in row z, column x.

At this point a four-quadrant arctangent, such as the C++ library function atan2, can extract θ_z = atan2(s_z,c_z) and θ_x = atan2(s_x,c_x). Numerical analysis provides a robust way to compute the “undoing” as a Givens rotation, shown below.

Givens rotation

Numerically stable calculation of a Givens plane rotation, mapping (a,b) to (r,0), depends on careful choice of the route to finding c and s, the cosine and sine of the angle θ, from a and b. The following method, choosing either t = cot θ or t = tan θ as route, is a variation of that recommended by Golub and Van Loan, made continuous as suggested by Edward Anderson.

if (b = 0) then ⟨c, s, r⟩ := ⟨sng(a), 0, |a|⟩
else if (a = 0) then ⟨c, s, r⟩ := ⟨0, sgn(b), |b|⟩
else if (|b| > |a|) then
  t := a/b
  u := sgn(b)*√(1+t2)
  s := 1/u
  c := s*t
  r := b*u
else
  t := b/a
  u := sgn(a)*√(1+t2)
  c := 1/u
  s := c*t
  r := a*u

Here sgn(x) is the signum function, which is +1, 0, or −1 as x is positive, zero, or negative. (In practice, the IEEE 754 copysign function would be used instead of a multiply.) The calculated c and s satisfy the equations r = c a+s b, 0 = c b−s a, and (of course) 1 = c²+s². This is a rotation by the negative of angle atan2(s,c), as required to undo it. Negating c and s gives a second, equally valid, solution which differs from the chosen one by 180°; the effect of this choice will propagate to the other two angles being extracted.

Archetype

Take xzx_s as a convenient archetypal system. Using the undoing approach, code to extract angles θ_x, θ_z, and θ_x′ takes the following form.

⟨a, b⟩ := ⟨My,x, Mz,x⟩
⟨c, s, r⟩ := Givens(a, b)
sx := c*Mz,y - s*My,y
cx := c*Mz,z - s*My,z
θx := atan2(sx, cx)
θz := atan2(r, Mx,x)
θx′ := atan2(s, c)

Consolidation

General code replaces the fixed indices with parametric i, j, k, and h, and negates as necessary. For example, when the system is xzy_s, we should assign ⟨a,b⟩ := ⟨Mx,x,Mz,x⟩ and negate the last angle, θ_y := −atan2(s,c). The row containing Givens argument a is thus indexed by h, with (i, j, k) being a permutation of (x, y, z).

These can be derived from one axis, say the right-most in matrix multiplication order, and three binary symmetry choices. Let neg be 1 if (i, j, k) is an odd permutation of (x, y, z); let alt be 1 if the first and last system axis are different; and let rev be 1 if a rotating frame is used; otherwise each bit is 0.

Since the archetype is xzx_s, its tuple is ⟨i =x,neg=0,alt=0,rev=0⟩. Roll, pitch, and yaw of a camera with x axis (pitch) left, y axis (yaw) up, and z axis (roll) back are the system yxz_r, with tuple ⟨i =z,neg=1,alt=1,rev=1⟩. (A left-handed system with +z pointing forward would cause problems, though it is possible to work around them.)

A convenient way to encode the 4-tuple is as a single integer from 0 to 23, with neg,alt,rev as the low-order three bits and i above them. With this representation j, k, and h can be found by table lookup. For example, if EulNext is an array containing {y,z,x,y} (so that, say, EulNext[y]=z), and if ^ is a bitwise exclusive-or, then the following would suffice.

j := EulNext[i+neg]
k := (x+y+z)-i-j
h := EulNext[k+1^neg^alt]

Using this encoding, conventions translate to tuples as follows.

	i	n	a	r			Sys.		i	n	a	r			Sys.		i	n	a	r			Sys.		i	n	a	r			Sys.
⟨	1	0	0	0	⟩		xzxs	⟨	1	0	1	0	⟩		xzys	⟨	1	1	0	0	⟩		xyxs	⟨	1	1	1	0	⟩		xyzs
⟨	2	0	0	0	⟩		yxys	⟨	2	0	1	0	⟩		yxzs	⟨	2	1	0	0	⟩		yzys	⟨	2	1	1	0	⟩		yzxs
⟨	3	0	0	0	⟩		zyzs	⟨	3	0	1	0	⟩		zyxs	⟨	3	1	0	0	⟩		zxzs	⟨	3	1	1	0	⟩		zxys
⟨	1	0	0	1	⟩		xzxr	⟨	1	0	1	1	⟩		yzxr	⟨	1	1	0	1	⟩		xyxr	⟨	1	1	1	1	⟩		zyxr
⟨	2	0	0	1	⟩		yxyr	⟨	2	0	1	1	⟩		zxyr	⟨	2	1	0	1	⟩		yzyr	⟨	2	1	1	1	⟩		xzyr
⟨	3	0	0	1	⟩		zyzr	⟨	3	0	1	1	⟩		xyzr	⟨	3	1	0	1	⟩		zxzr	⟨	3	1	1	1	⟩		yxzr

Code

We thus arrive at the following code, which assumes column vectors, right-handed coordinates, right-handed rotations, and matrices stored according to C++ array conventions (row-major order). The input is a 3×3 rotation matrix, M, and a 4-tuple ⟨i,neg,alt,rev⟩ indicating the desired Euler angle system, and the output is the three angles in order, ⟨θ₁,θ₂,θ₃⟩.

⟨j, k, h⟩ := indices(i, neg, alt)
v∗ := M∗,i
⟨a, b⟩ := ⟨vh, vk⟩
⟨c, s, vh⟩ := Givens(a, b)
s1 := c*Mk,j - s*Mh,j
c1 := c*Mk,k - s*Mh,k
θ1 := atan2(s1, c1)
θ2 := atan2(vj, vi)
θ3 := atan2(s, c)
if (alt = 1) then θ3 := -θ3
if (neg = 1) then ⟨θ1, θ2, θ3⟩ := ⟨-θ1, -θ2, -θ3⟩
if (rev = 1) then ⟨θ1, θ2, θ3⟩ := ⟨θ3, θ2, θ1⟩

References

• Thomas Herter, Klaus Lott. Algorithms for decomposing 3-D orthogonal matrices into primitive rotations. Computers & Graphics. ISSN 0097-8493. Volume 17, Issue 5, September–October 1993, pp. 517–527.
• Ken Shoemake. Euler Angle Conversion. Graphics Gems IV. Paul Heckbert (ed.). Academic Press, 1994, ISBN: 0123361567. pp. 222–229. (Online code)