Optimization, Maintainability and Major Milestones

The recent silence on the blog was a result of the debates of optimizations vs maintainability that were going on behind the scenes. Nevertheless, a major chunk of this project is now complete and is now ready to be merged into scipy.


Rotation Application

We can finally do what rotations were meant to do: rotate vectors. If an initial frame rotates to a final frame by rotation R, then the application of this rotation onto a set of vectors can be seen in two possible ways:

  1. As a projection of vector components of the original frame expressed in the basis of the final frame. In this respect, a rotation is simply a change of basis.
  2. As a physical rotation of a vector being glued to the original frame as it rotates. In this visualization, we stay in the initial reference frame and simply change coordinates.

There were a number of ways to approach this implementation, but we quickly settled on the simple approach of a matrix vector product.

The set of input output semantics to support, however, was infinitely more interesting and took some discussion. We now support applying N rotations on P points, where N and P follow numpy broadcasting rules: either one of them equals unity (in which case we broadcast the single object across the larger one), or they both equal each other (in which case we apply corresponding rotations on corresponding points). The output would be shaped (N  {or P}, 3) unless we had a single rotation being applied onto a single vectors.

Devising a simple implementation that would handle all these cases was challenging, but I ultimately found what I think is an elegant solution based on `numpy.einsum`.

Rotation Compositions

Compositions of a rotation were also added. There were quite a few implementations of this one. After some benchmarking and some discussion, we finally chose to compose quaternions using the Hamilton product. While not the fastest approach, it is easier to maintain in the long run and quickly catches up for large numbers of rotations.

Again, the number of rotations contained in two composable rotations R1 and R2 must follow the same rules as for the application case. Implementing this broadcasting was a much simpler task as numpy functions handled this implicitly for us.

Rotation Inversions and Indexing

We now support inverting rotations and extracting particular rotations from an object containing multiple rotations. This was one of the more straightforward and least controversial additions done since the last post.

Major Milestone

I’m happy to announce that the process of merging the Rotation class into the main scipy repository has been initiated and my work will soon be available to thousands of users across the globe.

This represents the completion of the first phase of this project. Now we need to focus on interpolations and estimations, all of which are interesting, mathematically and otherwise, in their own right. Stay tuned!

Euler angles, numerical stability and Gimbal Lock

This last week of work mainly dealt with the nitty-gritty issues of the Euler angle representation. Of all the representations so far this has to be the one most fraught with dangers and pitfalls, mainly because of one thing: Gimbal Lock

Gimbal Lock

Gimbal lock refers to the loss of a degree of freedom in a three-gimbal mechanism. This happens when to of the three axes are driven into a parallel configuration.

While implementing the `as_euler` function this led to a host of problems. Firstly, extening the algorithm mentioned in this paper to extrinsic rotations required careful thinking, and experimentation in itself.

Since we have the loss of a degree of freedom, it is not possible to uniquely determine all three angles. Additionally, the numerical instability means that making even tiny tweaks (such as not reversing the angle sequence returned by the algorithm for adapting to extrinsic rotations), which should work in theory, lead to failures by changing the sign of small quantities close to zero.

Angle Conventions

Since Euler angles are not unique, we had to decide a range of angles that we were allowed to return. Enforcing those restrictions led to another insight: the range of allowed Euler angles is determined by the choice of axes that we rotate about.

We chose to constrain the first and third angles in the range [-180, 180]. This covers the full circle and presented no issues.

Then we chose to constrain the second angle to [-90, 90] degrees. Again, this presented no issues in theory but we realised that this constraint cannot represent the full range of rotations when the first and third axes are parallel; that realization made us enforce a restriction within [0, 180] degrees in those cases.

All of the above insights are the result of a time consuming process of experimentation, trying to tease apart the algorithm, in order to extend it for our purposes.

This function will conclude the representation part of the class and will give way to cooler stuff like applications, compositions, and inversion in the coming weeks.