Implement direct conversions to and from quaternions and Euler angles

Feature request

Basically, I want to implement this myself, and would like help knowing how it would be the best way to implement it.

The idea is to implement two conversion methods into the Quaternion class:

1. Returning the angles in any possible sequence (XYX, XYZ, etc).
2. The inverse operation, returning the Quaternion from the angles.

This way, we don't have to do the extra step converting it to a Matrix, which would greatly improve the performance. As an extra, every possible sequence would be implemented.

Feature description

I recently published an article about a direct formula for the conversion between a quaternion variable to Euler angles in any sequence, which can be read here (Open Access).

Compared to either having 12 separate formulas (or 24, taking into account both intrinsic and extrinsic rotations) or using the well known quaternion-to-matrix-to-euler method, this has 3 main advantages:

1. Numerically, it is up to 30 times faster than the previous quaternion-to-matrix-to-euler method (used both originally in SciPy and in transforms3d, for example).
2. It is a lot simpler to implement, debug and maintain than both methods.
3. It provides the simple formulas that can be used for theorerical work.

Because of points 1 and 2, it has already been merged into Scipy. Because point 3, it has also been merged into SymPy.

I truly believe it's the best possible way to do this conversion, so I'm trying to contribute it to as many open source projects I can!

Implementation considerations

I have some questions on how to implement this:

1. How should the variable defining the sequence be? In Scipy, for example, a 3 letter string is used. If the string is uppercase, the sequence is Intrinsic. If the variable is lowercase, it is Extrinsic.
2. Should I reimplement void setEuler(angle1, angle2, angle3, sequence) as a general function, and reimplement void setEuler(yaw, pitch, roll) and void setRPY(yaw, pitch, roll) as calls to the 4 parameter version?
3. Should I implement a void getEuler(angle1, angle2, angle3, solution_number, sequence) and then implement two 3 parameter functions as calls to this function?
4. What is the difference between both solutions in the getEuler functions as defined in the Matrix class?

Start by reviewing the Quaternion class conversion entry points, including setEuler and setRPY, then compare them with the existing getEuler functions in the Matrix class. Resolve the sequence representation, intrinsic versus extrinsic behavior, and the meaning of the two solution values before defining completion as direct conversions supporting every requested sequence.