Work Product Summary

So after three months of working, it is now time to gather my thoughts and look back on what I’ve achieved. This brief post will also serve as a work product submission link.

First things first, all my merged changes live here as a new scipy.spatial.transform submodule. As of writing, the package include the following functionality:

  1. The Rotation class:
    1. Initialize a rotation from quaternion, direction cosine matrix, euler angles or rotation vector.
    2. Allows a user to convert freely between any of the above representations.
    3. Supports composition and inversion as operations on rotations.
    4. Supports application of rotations on points / vectors.
    5. Also contains a method of estimating the rotation between two frames of observation.
  2. The Slerp class:
    1. Spherical linear interpolation of rotations. Essentially the shortest path of constant angular velocity between keyframes rotations.

There is a single outstanding pull request for the RotationSpline class, which will implement cubic spline interpolation of rotations and rename Slerp to RotationSlerp to logically separate them from the generic spline and slerp methods.

Once this pull request has been finalized, I will also write a tutorial for the whole spatial.transform package as an introduction on how to work with rotations.

Do go through the Rotation class documentation to get a general idea of the usage of the class. (Documentation links are only finalised after every release of scipy, so take a look at the developer documentation for now).

Interpolations, matching and more

Now that a major chunk of the work has been done, I am making pull requests directly into the main scipy repository. This next phase of the project consitsts of rotation estimation and interpolations.

Wahba’s problem

Wahba’s problem is figuring out the rotation that maps one frame to another. Given a set of vector observations in two different frames, the new `Rotation.match_vectors` function returns the best estimate of the rotation that maps one to the other. This is often used to determine the orientation of spacecraft which only have stationary distant stars as reference.

The implementation will be completed within the next week and merged into the main repository.

Error bound based on input noise

While the estimation itself was pretty simple (the SVD method has been used for estimation) the major chunk of the work involved reading papers and figuring out how to give the user and estimate of the ‘correctness’ of the returned rotation. For this purpose, we also return a sensitivity matrix which is the [scaled] covariance matrix of the three component error vector expressed as a rotation of the initial frame.

I realise that’s a mouthful. In simple terms, the attitude error along each axes is treated as a random variable and the covariance matrix of these 3 random variables is computed.

Interpolations

Another useful aspect of having a rotation module is the ability to interpolate between orientation without having to specify all of them. This can easily be useful for 3D animation and simulations.

Two kinds of interpolations are supported:

Spherical Linear Interpolation

This is the ‘shortest-path’ equivalent along a 3D sphere. This leads to an orientation change with constant velocity between two successive keyframes.

Unfortunately, this often leads to a jerky motion at the keyframe endpoints when the angular velocity suddenly changes and the alternative is often preferred.

Spline Interpolation

Instead of using a linear function for interpolation, a spline (a polynomial of degree 3) is used. This has the additional benefits that we can choose the interpolation polynomial so as to ensure a continuous change in the angular velocity and acceleration at all times during the interpolation. This does not lead to the shortest path, but leads to an overall smooth orientation change.

This PR will also be merged in the coming week.

Next steps

As we head into these last few days of GSoC I plan to complete the two remaining pull requests and write a tutorial show casing the capabilities of this new module. Stay tuned!

Contribution Made

PR Merged

It’s been 7 weeks since GSoC has started and I’m proud to announce that the first phase of the work has been merged into scipy!

The module is production ready and completely agnostic of the representation used. The following are supported:

  • Quaternions
  • Direction Cosine Matrices
  • Rotation Vectors
  • Euler Angles

The following operations on rotations are also supported:

  • Composition
  • Inversion
  • Application on vectors

Since it is possible to work with multiple rotations in one object, we also support numpy-like indexing and slicing on rotations.

Documentation

Most of last week was spent in documenting and making minor changes to code. While I admit it was quite tedious, I honestly have a new-found appreciation for scipy (and numpy’s) extensive docs and the effort that goes into them. Without those, even this phase would have taken a much longer time.


It’s not over yet…

There’s still a lot of work to be done and exciting new features to be added:

  • Random sampling of rotations: While it seems simple in theory, the need to ensure that the rotations are uniformly distributed in 3D space makes it a little interesting.
  • Spherical Linear Interpolation of Quaternions: Ever used any video editing software where you specify a certain number of keyframes and the motion in interpolated between them automatically? Think that, but for quaternions. The Spherical Linear imposes a further constraint of constant angular velocity.
  • Qspline: Another interpolation algorithm that uses splines instead of straight lines for interpolation. Leads to smoother rotations.
  • Wahba’s Problem: Given two sets of vectors, how do you best estimate the rotation that describes them? This is a useful and well studied problem in spacecraft dynamics with a lot of interesting theory behind it.

Stay tuned for more updates!

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.


Enhancements

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.

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!