Quaternions and Rotations 2: Euler angles etc.

Let me create a simple example. You may want to have read Quaternions and Rotations 1, and the introductory post as well.

I will specify an Euler angle sequence, “ZYX”, with angles:

The notation “ZYX” means that we first rotate about the z-axis, then about y, then about x; and I choose to let the angles be in the same order. This is typically used in aerospace applications. The other Euler angle sequence common in my experience is “ZXZ” for orbits. Altogether, I count 12 kinds: three choices for the first axis, two for the second, and two for the third. There is no point in choosing consecutive axes the same, because consecutive rotations about the same axis are equivalent to one rotation, through the sum of the angles, about that axis.

The obvious first step is to construct the combined rotation matrix for that Euler angle sequence, but I’m going to go in a different order. I’m going to test all 8 transformations among rotation matrices, angle/axis, quaternion, and Euler angle sequences.

What the heck? Pictures are good.

Let me first do Euler angle sequence to quaternion.

We know how to construct the quaternion representation for each of those rotations, and hence for the combined rotation. We construct 3 quaternions and multiply them together. Here’s the one for the x-axis rotation; since the axis is (1,0,0), the imaginary part is just (SX, 0, 0), where SX is the sine of the half-angle:

Then the composite rotation is given by the quaternion product (note the order, qz on the left):

This is some of the most complicated code I wrote for this stuff — not for the mathematics, but for the string and list manipulations! I’m sure they can be done far more elegantly than I did.

The algorithm is straight-forward, and looks rather free from peril to me. Worst case, I think, is that one or two of the angles are zero — but we still know the rotation axes in the Euler angle sequence. Still, one or two of the quaternions would collapse to (1, 0, 0, 0), because the imaginary part is multiplied by the sine of 0 (i.e. by 0), but that seems ok; it corresponds to an identity matrix. In addition, a rotation of 180° should not be an issue: the sine of half the angle is the sine of 90°, i.e. 1.

To get from the quaternion to the angle/axis representation, we take the ArcCos of the real part… and double it:

Since it is less than Pi/2, I’ll leave it alone. If it were greater than Pi/2, I would choose to replace it by an equivalent negative angle, by subtracting 2 Pi. (More on that in just a few lines further on.)

Having the angle, now make a unit vector out of the imaginary part, to get a unit axis of rotation:

For reference, here is my code. Don’t operate heavy machinery with it; my primary motivation in displaying it is to be crystal clear about what I did.

qToAa[q_Quaternion] := Module[{angle, ha, a1, unit, cha},
unit = Sign[q];
cha = unit[[1]];
ha = ArcCos[cha] // N;
a1 = 2 ha;
angle = If[a1 > Pi, a1 – 2 Pi, a1];
{angle, unitv[{unit[[2]], unit[[3]], unit[[4]]}]}
]

Sign[q] gets me a unit quaternion…
the real part, its first entry, is the cosine of the half angle…
get the half angle…
double it…
but if it’s more than 180°, write it as a negative angle…
finally, as we return, make a unit vector of the imaginary part of the quaternion…

You may wish to note that there are two choices in that code. First, the ArcCos returns an angle in the half-open interval [0, 180°). (Maybe we shouldn’t consider that a choice.) After I double that angle, I’m in the half-open interval [0, 360°). The real choice is that I personally want an angle in the half-open interval (-180°, 180°].

I’ll remark that half of my angle is now guaranteed to have a positive cosine.

To get from the angle/axis to the rotation matrix, we construct a skew-symmetric matrix N from the unit-vector axis of rotation…

Here’s the code.

toSkew[v_List] := {{0, v[[3]], -v[[2]]},
{-v[[3]], 0, v[[1]]},
{v[[2]], -v[[1]], 0}}

aaToRot[{angle_, axis_List}] := Module[{unit, n},
unit = unitv[axis];
n = toSkew[unit];
IdentityMatrix[3] + Sin[angle] n + (1 – Cos[angle]) n.n
]

Notice that I will construct a unit vector out of the given axis. Because my example starts with an Euler angle sequence, my subsequent transformations to angle/axis are automatically called with unit vectors. But I don’t ever want to have to compute a unit vector first in order to start with an angle/axis representation. As I mentioned before, I want to be able to say, “Do a 120° rotation (of the coordinate system) about the axis (1,1,1)” — and that’s not a unit vector.

Of course, that sequence of transformations is not how we would go from an Euler angle sequence to the rotation matrix. We can do that directly.

To construct the matrix directly from its Euler angle sequence, we just multiply three rotation matrices together, Rx Ry Rz with the corresponding angles:

Let me call it m1 to distinquish it conceptually from the alternative calculation that led to m0; they are the same (to within computational noise, about 10^-16).

Anyway, my code matches that answer:

And, again, the algorithm is simple, but my string and list manipulations are probably unneccessarily complicated.

To get the angle/axis from the rotation matrix, I find the eigenvector corresponding to eigenvalue 1.

{$\lambda\$,PT} = Eigensystem[m1];

The eigenvalues are

$\lambda\$ = {1., 0.41839+0.908267 i, 0.41839-0.908267 i}

It’s the first one we want, and the corresponding eigenvector is… and it is already a unit vector:

The cosine of the angle can be found from the trace of the matrix:
trace = 1 + 2 cos $\theta\$,

As we have seen before, that difference is the numerical value of the symbolic matrix

$m0 - m0^T = \left(\begin{array}{ccc} 0 & 2 c \sin (\theta ) & -2 b \sin (\theta ) \\ -2 c \sin (\theta ) & 0 & 2 a \sin (\theta ) \\ 2 b \sin (\theta ) & -2 a \sin (\theta ) & 0\end{array}\right)$

All three off-diagonal terms are non-zero, so we can use any one of them to get the sine of the angle, sa. (My code will test all three off-diagonal terms if necessary, since two of them could be zero for a rotation about a coordinate axis.)

By the way, I’ve decided that using the eigenvector is safer than trying to extract a unit vector (a, b, c) from that matrix, since the off-diagonal terms could be very small, for small values of sin $\theta\$. And they could even be zero, legitimately, for a rotation of 180°.

Unfortunately, it’s not literally the same answer as we got before, but its negative:

a0 = {1.13912, {0.460404, 0.852003, 0.249237}}

We have the negative angle and and negative axis. The fortunes of war: these are two descriptions of exactly the same rotation.

To get from angle/axis to quaternion, we compute the cosine and the sine of half the angle, use the cosine for the real part, and the sine times the axis of rotation for the imaginary part:

And my code:

aaToQ[{angle_, v_List}] := Module[{unit},
unit = unitv[v];
Quaternion[Cos[angle/2],
Sin[angle/2] unit[[1]] ,
Sin[angle/2] unit[[2]],
Sin[angle/2] unit[[3]] ] // N
]

Note, like aaToRot, this code will always make a unit vector out of the input vector; I can call this code without creating a unit axis of rotation for it myself.

OK, we’ve come to the new part: getting from a quaternion to an Euler angle sequence. We have the quaternion representation of the rotation m1 (= m0):

q1 = Quaternion[0.842137,0.248279,0.459454,0.134404]

We can write the three quaternions for the specified Euler angle sequence “ZYX”:

Now just equate the imaginary parts, and solve. Numerically. And I replaced cosines by positive square roots involving sines. Looking for half-angles with positive cosines is consistent with my transformation from quaternion to angle/axis.

And I will give it starting values in the first quadrant. This warrants further investigation, but it will do for a start.

I get three equations in three unknowns. I call FindRoot with initial sines = 1/2.

That’s it, these are the sines of the half angles…

sha = {0.34202, 0.382683, 0.422618}

(I had to change the order of the answers to match the ZYX. More list manipulations in order to automate the process….)

So we get the half angles…

has = {0.349066, 0.392699, 0.436332}

and we double them:

2 has = {0.698132, 0.785398, 0.872665}

and those are, indeed, the angles we started with:

as = {0.698132,0.785398,0.872665}

Of course, I have code for that:

That system of equations to be solved was constructed by my code just as I constructed it here. And, as with a few other things, the messy part is the string and list manipulations to build the appropriate pieces.

Three equations in three unknowns is a lot simpler than the nine equations in three unknowns we would have if we tried to decompose a matrix instead of a quaternion. Among other things, it may be far easier to search for other solutions among three equations than among nine.

OK, those are seven of the eight transformations. As I said in the last post, I can find the Euler angle sequence decomposition of a matrix by finding the Euler angle sequence of its quaternion representation.

I have code for that, and it’s very simple:

rotToE[s_String, mat_List] := Module[{t1, q1},
t1 = rotToAa[mat];
q1 = aaToQ[t1];
qToE[s, q1]
]

It does exactly what I promised: it goes from matrix to angle/axis to quaternion to an Euler angle sequence. It traces out the three blue arrows in the picture way back at the beginning, but it looks like the red arrow, a direct decomposition of a rotation matrix into an Euler angle sequence.

Now I can play with Euler angle sequences.