Rotations, Orientations, and their Representations

Image courtesy pixabay

Introduction

This is a story about Quaternions. It will also turn out to be a story about Rust (the programming language), but that will come later. Quaternions are an extension of complex numbers that turn out also to be closely related to rotations. If you want to describe how an object like a satellite or a coffee mug is rotating, quaternions are an excellent way of doing so. They are also closely related to spin ½ particles in quantum mechanics, but I’m not going to talk much about that.

I was enamored from the moment I first encountered quaternions in undergrad. At first glance they appear to be a fairly trivial extension of complex numbers, but unlike complex numbers, multiplication is not commutative: if $a$ and $b$ are quaternions, then in general, $a \cdot b \neq b \cdot a$. This was my first exposure to a non-abelian group, a way of modeling phenomena where the order in which things happen matters. Far from being a mathematical curiousity, many real-life phenomena are non-abelian. If I spill some tea on the floor and mop it up, the final result is very different than if I mop the floor and then spill some tea.

It turns out rotations are also non-abelian. If I rotate a coffee mug one way, and then another, the final orientation will in general be different than if I had performed the rotations in the reverse order. The space of rotations is called $\textrm{SO}(3)$. If we represent the action of performing two rotations in sequence as the multiplication of two quaternions, then quaternions of unit magnitude correspond to a double cover of $\textrm{SO}(3)$. That means every rotation can be represented as a quaternion (and except for certain rotations, there are actually two quaternions corresponding to the same rotation).

Having started my career in the Attitude Control Systems department at Boeing, I spent a lot of time thinking about rotations and how we represent them. As I moved from project to project over my seven years there, I ended up writing the same basic library in Matlab, C++, Java, and even Perl!

I discovered that not only were quaternions interesting from a mathematical perspective; they also posed an interesting exercise in polymorphism (the object oriented programming concept). That’s because quaternions are not the only way of representing rotations, and rotations are not the only thing that quaternions can represent! Other ways of representing rotations include Euler angles, direction cosine matrices (also known as rotation matrices), and Rodrigues parameters.1 Moreover, anything that can be used to represent a rotation (like a quaternion), can also be used to represent an orientation. The code I wrote needed to support multiple representations (quaternion vs rotation matrix), and multiple applications (rotation vs orientation).

It gets even more confusing! When I started looking at others’ code, the equations were different than what I had learned in college, and different from what was on Wikipedia! That’s because there are different conventions about how to interpret a rotation. For example, should the rotation be interpreted in an inertial frame, or in the frame of the thing being rotated? Should rotations be interpreted in a left-handed or right-handed sense? Even the way the data is stored in the computer (column major order vs row major order) influenced the decisions other programmers before me had made. Since different teams were using different conventions, the software had to be flexible enough to support all the possibilities.

Every time I implemented the library in a new language, I ended up learning another nuance. I tried to program in such a way that abstracted these details away from the programmer (myself) so that the programmer (myself) could be confident in the code using those libraries.

In compiled, statically-typed languages with polymorphism (Java), I discovered that the right class structure could prevent certain types of mistakes. The compiler itself could catch certain types of errors. Coming from a more physics-y background, this seemed like magic to me!

Recently I started working through the Rust tutorial.2 Rust has a pretty interesting take on polymorphism. Rust has generics: we can write a single function that operates on different data types. For example, we can write an “add” function that operates on integers or floats. Instead of interfaces, Rust has traits, a way of defining a collection of behaviors certain types of entities should implement. To practice these ideas, I decided to implement a rotation and orientation library.

Quaternions and Rust

Before we get started, I should mention I’m developing on a Macbook, using the 2018 (1.31.1) version of Rust. Rust claims to be cross-platform, so hopefully the code in this post works regardless of your computer. Rust is under active development, so this code might not work in earlier or later versions of Rust. The ideas should be fairly stable, however. In the code examples, lines that begin with dollar signs indicate commands sent to the shell. Otherwise, code is Rust.

I created a new library project using cargo:

$ cargo new --bin orientations

That created a lib.rs file in a src/ folder, which is where our code will go. The code below shows our initial implementation.

struct Vector3d {
    data: [f64; 3]
}

struct Quaternion {
    real_part: f64,
    imaginary_part: Vector3d
}

trait Rotation {
}

trait Orientation {
}

I created a Quaternion struct consisting of real and imaginary parts. The real part is a single floating point number. The imaginary part is three floating point numbers, which I represented as a vector. (There are of course vector data types in Rust, but I wanted to create the functions I would need so that it would be more logically self-contained.) Then I created Rotation and Orientation traits which will contain all the functions rotation and orientation objects should support.

Then I created some convenience functions for creating vectors and quaternions, shown below.

impl Vector3d {
    fn new(data: [f64; 3]) -> Vector3d {
        Vector3d{ data }
    }

    fn zero() -> Vector3d {
        Vector3d::new( [0.0, 0.0, 0.0] )
    }
}

impl Quaternion {
    fn new(real_part: f64, imaginary_part: Vector3d) -> Quaternion {
        Quaternion {
            real_part,
            imaginary_part
        }
    }
}

Next up is to start defining the methods associated with a Rotation. We should be able to rotate things, obviously. It turns out there are two things we can rotate: vectors and three dimensional objects, and the equation for rotating those two things is different.

We should be able to get the identity rotation: the “rotation”, $e$, corresponding to no rotation at all! We should be able to invert any rotation: for any rotation $r$ there exists another rotation $r^{-1}$ that “undoes” the first rotation. This inverse is dual to the original rotation: the inverse of the inverse is the original rotation: $(r^{-1})^{-1} = r$. Rotating first by $r$ and then by $r^{-1}$ is equivalent to no rotation at all. This is also true if we rotate first by $r^{-1}$ and then by $r$. These facts make good test cases, so we are documenting them as we go.

We should be able to compose rotations. Given two rotations, $p$ and $q$, there exists a third rotation, $r$ that is equivalent to rotating first by $p$ and then rotating by $q$. There is a fourth rotation $r^\prime$ that is equivalent to rotating first by $q$ and then by $p$, and in general, $r \neq r^\prime$.

We should also be able to get various representations of that rotation: regardless of how it is actually represented, we should be able to get the quaternion or the rotation matrix or the axis and angle associated with a rotation.

Identity and Inverse Rotations

Let’s start with getting the identity rotation. A generic quaternion $q$ is written like $q_r + q_i i + q_j j + q_k k$, where $q_r$ is the real part of the quaternion, $[q_i, q_j, q_k]$ is the imaginary part, and $i$, $j$, $k$ are imaginary numbers satisfying $i^2 = j^2 = k^2 = ijk = -1$. The identity rotation has real part $1$ and imaginary part $[0, 0, 0]$.

trait Rotation {
    type R: Rotation;
    fn identity() -> Self::R;
}

impl Rotation for Quaternion {
    type R = Quaternion;

    fn identity() -> Quaternion {
        Quaternion::new(1.0, Vector3d::zero())
    }
}

Here we are specifying an associated type, and the identity function returns something of that type. For Quaternions, the associated type is a Quaternion, and in general the associated type is whatever is implementing the rotation trait.

Next we’ll implement the inverse. To do so we introduce the conjugate of a quaternion $p$: the quaternion $p^\dagger$ having the same real part as $p$ but negative imaginary part. Specifically, if $p = p_r + p_i i + p_j j + p_k k$, then $p^\dagger = p_r - p_i i - p_j j - p_k k$. We also introduce the norm or magnitude of a quaternion p, which is simply $| p | = \sqrt{p_r^2 + p_i^2 + p_j^2 + p_k^2}$. Then the inverse of a quaternion p is the conjugate of $p$ divided by the norm squared: $p^{-1} = p^\dagger / |p|^2$.

If we are using quaternions to represent rotations, then in principle $|p|$ will always be 1 and $p^{-1} = p^\dagger$; however, floating point errors may result in deviations from unit magnitude as we are composing rotations. We should be able to normalize quaternions, but it actually is not necessary (and in fact, it’s inefficient) to normalize as part of every relevant operation. If $q$ is a quaternion and $\alpha$ is some non-zero scalar, then $q$ and $\alpha q$ correspond to the same rotation. If we set $\alpha = 1 / |q|$ then $\alpha q$ is a quaternion of unit magnitude, showing that we can always normalize a quaternion. Moreover, if $q$ already has unit magnitude, then so does $-q$. These two quaternions correspond to the same rotation, showing the “double cover” property of quaternions mentioned above. (Typically we simply choose the quaternion having non-negative real part as the sole representative of a particular rotation.) Because of the equivalency of $q$ and $\alpha q$, we will only normalize on an as-needed basis.

First we’ll add some convenience functions that aren’t especially interesting.

impl Vector3d {
    fn new(data: [f64; 3]) -> Vector3d {
        Vector3d{ data }
    }

    fn dot(&self, other: &Vector3d) -> f64 {
        let mut dot_product: f64 = 0.0;
        for i in 0..3 {
            dot_product += self.data[i] * other.data[i];
        }
        dot_product
    }

    fn norm_squared(&self) -> f64 {
        self.dot(&self)
    }

    fn norm(&self) -> f64 {
        self.norm_squared().sqrt()
    }

    fn scalar_multiple(&self, alpha: f64) -> Vector3d {
        Vector3d::new(
            [
                alpha * self.data[0],
                alpha * self.data[1],
                alpha * self.data[2]
            ]
        )
    }

    fn negate(&self) -> Vector3d {
        Vector3d::new(
            [
                -self.data[0],
                -self.data[1],
                -self.data[2]
            ]
        )
    }

    fn zero() -> Vector3d {
        Vector3d::new( [0.0, 0.0, 0.0] )
    }
}

Now we are ready to implement the inverse method for the Rotation trait.

const DBL_EPSILON: f64 = 2.2204460492503131e-16;

trait Rotation {
    type R: Rotation;
    fn identity() -> Self::R;
    fn inverse(&self) -> Self::R;
}

impl Quaternion {
    // snip
    fn conjugate(&self) -> Quaternion {
        Quaternion::new(self.real_part, self.imaginary_part.negate())
    }

    fn norm_squared(&self) -> f64 {
        self.real_part * self.real_part + self.imaginary_part.norm_squared()
    }
}

impl Rotation for Quaternion {
    type R = Quaternion;
    // snip

    fn inverse(&self) -> Quaternion {
        // Check that norm is > 0
        let norm_squared = self.norm_squared();
        if norm_squared < DBL_EPSILON {
            panic!("Quaternion close to zero; cannot invert.")
        }

        let inv_norm_squared = 1.0 / norm_squared;
        let c = self.conjugate();
        let real_part = c.real_part * inv_norm_squared;
        let imaginary_part = c.imaginary_part.scalar_multiple(inv_norm_squared);
        Quaternion::new(real_part, imaginary_part)
    }
}

Note that to avoid re-writing the same code over and over again, I am concealing some code using // snip, as is done in the Rust tutorial. As mentioned above, as we compose rotations, floating point errors may introduce small deviations from unit magnitude. In the extreme case, the magnitude might become close to zero. We would get a division-by-zero error in that case, which we are explicitly checking for and panicking. It might make more sense to return an error for the calling function to handle, depending on the application. In most applications, this should never happen but it is good to be safe.

Getting Particular Representations

Next up, we want to be able to get particular representations of rotations. No matter how a rotation is represented, we should be able to represent this rotation as a quaternion, or a rotation matrix. For now, we will only implement the quaternion representation, since we haven’t even implemented rotation matrices yet.

impl Vector3d {
    // snip
    fn clone(&self) -> Vector3d {
        Vector3d::new( self.data )
    }
}

trait Rotation {
    // snip
    fn as_quaternion(&self) -> Quaternion;
}

impl Rotation for Quaternion {
    // snip
    fn as_quaternion(&self) -> Quaternion {
        Quaternion::new(self.real_part, self.imaginary_part.clone())
    }
}

This isn’t especially interesting, we’re pretty much just cloning the quaternion. Converting a rotation matrix to a quaternion or vice versa is less trivial. Mostly we need this function for what follows.

Composing Rotations

Next let’s implement composing rotations. As mentioned above, composing rotations is equivalent to multiplying quaternions; however, confusingly, $p \cdot q$ typically represents rotating first by $q$ and then by $p$. That is because we typically think of vectors as being column vectors. When we represent rotations as rotation matrices, we would write $R \cdot v$ to represent the action of rotating a vector $v$ by a rotation $R$. If we were to do a sequence of rotations, rotating first by $R_1$, then by $R_2$, we would write this as $R_2 \cdot (R_1 \cdot v) = (R_2 \cdot R_1) \cdot v$. In this convention, we “read right-to-left”: the rotation corresponding to a rotation first by $R_1$ and then by $R_2$ is equal to $R_2 \cdot R_1$. For this reason, we also interpret quaternion multiplication right-to-left: $p \cdot q$ is the rotation corresponding to rotating first by $q$ and then by $p$. If this isn’t confusing enough, not everyone does it this way, and sometimes $p \cdot q$ is interpreted the other way around.

For this reason, we prefer to separate the concept of composing rotations from the concept of multiplication. We can choose function names that explicitly describe the order of rotations, as shown in the function names below.

trait Rotation {
    type R: Rotation;
    fn identity() -> Self::R;
    fn before<T: Rotation>(&self, r: &T) -> T;
    fn after<T: Rotation>(&self, r: &T) -> T;
    fn multiply<T: Rotation>(&self, r: &T) -> Self::R;
}

(There is actually a highly-not-obvious problem with this code we will come back to shortly.) Given two rotations $p$ and $q$, p.before(q) represents the rotation first by $p$ and then by $q$, whereas p.after(q) represents the rotation first by $q$ and then by $p$. Of course, p.after(q) is equivalent to q.before(p), but importantly, if $p$ is represented by a quaternion, and $q$ is represented by a rotation matrix, then p.after(q) will be a rotation matrix, but q.before(p) will be a quaternion. We of course could have made it the other way around, but I liked the idea of having the operand have the same type as the result.

Under the hood, both before and after correspond to multiplications, so we will have another function that does just that. In a sense, before and after are syntactic sugar, but they make code easier to write correctly and reason about. So much so, that ideally multiply would be a private method, and external code wouldn’t even be able to call it directly. We’ll hold off on labeling functions as public and private so as not to clutter the code for now. Moreover, the multiply function has the same return type as the rotation outside the parentheses, not the operand. I had to do it this way to get the before and after behavior that we want, as we will see shortly.

impl Rotation for Quaternion {
    type R = Quaternion;
    // snip

    fn multiply<T: Rotation>(&self, r: &T) -> Quaternion {
        let rr = r.as_quaternion();
        let real_part = self.real_part * rr.real_part - self.imaginary_part.dot(&rr.imaginary_part);
        let q_i = self.imaginary_part.data[1] * rr.imaginary_part.data[2] - self.imaginary_part.data[2] * rr.imaginary_part.data[1];
        let q_j = self.imaginary_part.data[2] * rr.imaginary_part.data[0] - self.imaginary_part.data[0] * rr.imaginary_part.data[2];
        let q_k = self.imaginary_part.data[0] * rr.imaginary_part.data[1] - self.imaginary_part.data[1] * rr.imaginary_part.data[0];
        let imaginary_part = Vector3d::new( [q_i, q_j, q_k] );

        Quaternion::new(real_part, imaginary_part)
    }

    fn before<T: Rotation>(&self, r: &T) -> T {
        r.multiply(self)
    }

    fn after<T: Rotation>(&self, r: &T) -> T {
        r.inverse().multiply(&self.inverse()).inverse()
    }

}

The multiply function itself is the least interesting addition, at least from a Rust perspective. In the before function, we see that q.before(r) is equal to $r \cdot q$, reflecting the “right-to-left” interpretation discussed above. This is abstracted away from the user. The after function has a bizarre implementation: q.after(r) should simply be $q \cdot r$; however if we simply implemented this as self.multiply(r), the result would be a Quaternion, regardless of the type of r. If r is a rotation matrix, we want q.after(r) to be a rotation matrix, not a Quaternion, for the reasons discussed above. This issue did not come up in the before implementation since r.multiply(q) has the same type as r regardless of the type of q, as desired. We can use the fact that $q \cdot r = (r^{-1} \cdot q^{-1})^{-1}$ to get the desired behavior: r.inverse() has the same type as r, and therefore r.inverse.multiply(q.inverse()) has the same type as r. Finally, the inverse of this quantity has the same type as r. This is the first clue that domain knowledge of quaternions coupled with good software architecture instincts allows us to write code that abstracts the nuances of rotations away from the user!

If we try to compile this code, we will get an error.

$ cargo check
    Checking orientations v0.1.0 (/Users/rwilson4/Projects/orientations)
error[E0308]: mismatched types
   --> src/lib.rs:130:9
    |
129 |     fn before<T: Rotation>(&self, r: &T) -> T {
    |                                             - expected `T` because of return type
130 |         r.multiply(self)
    |         ^^^^^^^^^^^^^^^^ expected type parameter, found associated type
    |
    = note: expected type `T`
               found type `<T as Rotation>::R`

This function is supposed to return whatever r is: if r is a quaternion, it should return a quaternion for example. The multiply function returns the same type as r, however the compiler isn’t quite smart enough to realize that. We just need to tell it that not only is r of type T where T is a rotation, but T is a rotation having associated type R=T. All we need to change is the function signatures:

trait Rotation {
    type R: Rotation;
    // snip
    fn before<T: Rotation<R = T>>(&self, r: &T) -> T;
    fn after<T: Rotation<R = T>>(&self, r: &T) -> T;
    fn multiply<T: Rotation>(&self, r: &T) -> Self::R;
}

impl Rotation for Quaternion {
    // snip
    fn before<T: Rotation<R = T>>(&self, r: &T) -> T {
        r.multiply(self)
    }

    fn after<T: Rotation<R = T>>(&self, r: &T) -> T {
        r.inverse().multiply(&self.inverse()).inverse()
    }
}

Coda

This is an excellent breaking point, but there is still a lot to do. We haven’t implemented orientations, and so far quaternions are the only Rotation we have implemented. We also need to write test cases and benchmark performance (some of the code above is knowingly sub-optimal for simplicity as well as to give us something to optimize). Stay tuned for a follow up article!

Subscribe to Adventures in Why

* indicates required

  1. Wiley J. Larson and James R. Wertz (editors). Space Mission Analysis and Design, 3rd Edition. ↩︎

  2. https://doc.rust-lang.org/book/index.html ↩︎

Bob Wilson
Bob Wilson
Data Scientist

The views expressed on this blog are Bob’s alone and do not necessarily reflect the positions of current or previous employers.

Related