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.

Rotation Vectors and Euler Angles

This concludes another week of working with scipy. This week was spent in implementing the conversion to and from rotation vectors. In addition, work was also started on initialization from Euler angles. Pretty standard stuff, but there’s always something to learn.

Rotation Vectors

A rotation vector is a 3D vector whose direction gives the axis of rotation and whose norm gives the angle of rotation. Calculating the rotation quaternion from these two quantities is straightforward and is given by the formula here.

Taylor Series

Even a simple formula needs to be implemented carefully. Were it not for my mentor, I would never have thought of using Taylor expansion for small numbers to avoid floating point inaccuracies.

Euler Angles

This brings us to a much more controversial part of the module. This part required some discussion with regards to the API, the kind of rotations we decided to allow, and how best to implement the general cases.

The first task was to determine the kind of rotations that we supported.

Extrinsic Rotations

In extrinsic rotations, the coordinate axes do not rotate with the body and remain fixed. In essence, these are rotations with respect to a global or ground frame of reference in physics.

Intrinsic Rotations

In intrinsic rotations, the coordinate axes rotate with the body. These are rotations with respect to a body-centered frame of reference.

In theory any possible combination of intrinsic and extrinsic rotations is possible.  In addition, there are few a definitions of Euler angles depending on the possible sequence of axes of rotation. It took some discussion to narrow down the range of possibilities and decide on the API that we would implement. After a few rounds of suggestions, initialization from euler is finally ready.


That’s it for this week I guess.

The Code

Hi all!

So it appears that with everything else going on, I forgot to share the link to the public repository. Not to worry, the project lives here currently. It will be merged into the main scipy repo at a later date.

Guess that’s it. Back to work.

Quaternions and Direction Cosine Matrices

This last week was spent implementing the initialization from and conversion of quaternions to direction cosine matrices. They are orthogonal matrices which can be used as transformation matrices for rotation.

Research

One key part of initializing from matrices is to decide what to do when the inputs are not orthogonal. Normally, these are not considered valid rotations but in many applications these matrices are valid with errors accumulated due to noise in measuring instruments or rounding off. I spent some time going through this resource to figure out how to best approximate an orthogonal matrix in that case.

Numpy speedups

Writing vectorized code for this particular algorithm was a fun challenge, mainly because there are 4 cases to deal with and vectorized operations don’t easily agree with if-else conditions.

I was able to generalize 3 out of the 4 cases and with a little help was able to figure out how to best implement the function, which lives here.

There was some discussion about the efficiency of vectorized code when there aren’t many objects to deal with. I ran benchmarks on both scalar and vectorized approaches. Those were quite informative and I realized that there are performance penalties vectorized code has to pay when dealing with very few objects.

Ultimately, however, Donald Knuth

Premature optimization is the root of all evil

and maintenance costs convinced us to stick to the vectorized approach.

Continuous Integration and Testing

This week I also set up Travis-CI and Circle-CI on my repo.  They were ridiculuosly easy to set up and went a long way towards making me capable of writing production ready code.  I’ve seldom had reason to worry about testing or maintability issues but setting this up made me realize I was in the big leagues now.

Cool Tricks

So this week I learned about the Einstein summation convention. It’s a powerful and succinct representation of complicated matrix operations.

For example, matrix multiplication is just:

...ij, ...jk-->...ik

I admit that’s slightly cryptic and not very useful now that numpy has a matrix multiplication function, but it’s very handy when you want to support versions that do not have an inbuilt __matmul__.


All in all, it was an amazing week and I look forward to working and learning even more.

Introductions

Hi all,

My name is Aditya Bharti and I am extremely excited to be working with scipy over this summer. I am an undergraduate researcher in IIIT Hyderabad working on Computer Vision and Machine Learning. My work so far can be found here.

My project is to implement 3D Rotations in scipy with a representation-independent interface, in addition to some interesting and useful algorithms for rotation fitting. This introductory blog post will explain some key aspects of my project.

What is a coordinate frame?

A coordinate frame is a system which uses one or more numbers, known as coordinates, to uniquely determine a position in space.

For the purposes of this project I will restrict myself to the familiar 3D Euclidean world. Each observer is free to choose any point as origin (0,0,0) and any 3 directions as the principal axes. Rotations are used to change between frames of references when the coordinate frames have the same origin.

Understanding Rotations

Rotations as Vector Transformations

The rotation transform applied to a point is easiest to visualize as the circular motion around an axis of rotation. This view of rotations as point transforms is less powerful than the next one under consideration.

Rotations as Frame Transformations

The second way of viewing rotations is as relating coordinates in the final coordinate frame to the coordinates in the initial coordinate frame.

Consider the example of a 90 degree anti-clockwise rotation. Let C1 be the coordinate frame before rotation and let C2 be the coordinate frame after rotation. Consider this rotation applied to the vector (1, 0, 0) in frame C1.

  1. Application on points
    The rotation can be thought of as simply rotating (1, 0, 0) [frame C1] about (0, 0, 0) to reach the new position (0, 1, 0) [frame C1].
  2. Relating frames
    If we want to answer the question “what are the coordinates of point (1, 0, 0) [frame C2] in frame C1?” this approach to rotations is more powerful. The answer is obtained by simply applying the rotation to (1, 0, 0) to get point (0, 1, 0) in frame C1. Intuitively, this can be seen by rotating the coordinate axes by 90 degrees anti-clockwise and observing that the new x-axis coincides with the old y-axis.

Representing rotations

Rotations can be represented in many different ways numerically

  • Vectors: The new positions of the 3 axes are stored
  • Discrete Cosine Matrices: Since a rotation is a linear transform, it can be represented as a matrix
  • Euler Angles: Once the sequence of axes is fixed, three angles are all that are required to represent any rotation
  • Axis-angle: Stores the direction of the axis of rotation and the angle by which to rotate
  • Quaternions: Since the rotation matrices are orthogonal, they can be represented using just 4 numbers. This mathematical structure can be seen as an extension of a complex number and has many interesting properties.

First Post!!

Hi all,
I’m extremely excited to be working with scipy over the summer.

I’ve tried to make full use of the bonding period so far by getting to know the community and solving setup related issues.

Stay tuned for more updates!