Component of a quaternion rotation around an axis

Question

I'm having trouble finding any good information on this topic. Basically I want to find the component of a quaternion rotation, that is around a given axis (not necessarily X, Y or Z - any arbitrary unit vector). Sort of like projecting a quaternion onto a vector. So if I was to ask for the rotation around some axis parallel to the quaternion's axis, I'd get the same quaternion back out. If I was to ask for the rotation around an axis orthogonal to the quaternion's axis, I'd get out an identity quaternion. And in-between... well, that's what I'd like to know how to work out :)

Solution 1

I spent the other day trying to find the exact same thing for an animation editor; here is how I did it:

1. Take the axis you want to find the rotation around, and find an orthogonal vector to it.
2. Rotate this new vector using your quaternion.
3. Project this rotated vector onto the plane the normal of which is your axis

4. The acos of the dot product of this projected vector and the original orthogonal is your angle.

   public static float FindQuaternionTwist(Quaternion q, Vector3 axis)
   {
       axis.Normalize();
   
       // Get the plane the axis is a normal of
       Vector3 orthonormal1, orthonormal2;
       ExMath.FindOrthonormals(axis, out orthonormal1, out orthonormal2);
   
       Vector3 transformed = Vector3.Transform(orthonormal1, q);
   
       // Project transformed vector onto plane
       Vector3 flattened = transformed - (Vector3.Dot(transformed, axis) * axis);
       flattened.Normalize();
   
       // Get angle between original vector and projected transform to get angle around normal
       float a = (float)Math.Acos((double)Vector3.Dot(orthonormal1, flattened));
   
       return a;
   }
   

Here is the code to find the orthonormals however you can probably do much better if you only want the one for the above method:

private static Matrix OrthoX = Matrix.CreateRotationX(MathHelper.ToRadians(90));
private static Matrix OrthoY = Matrix.CreateRotationY(MathHelper.ToRadians(90));

public static void FindOrthonormals(Vector3 normal, out Vector3 orthonormal1, out Vector3 orthonormal2)
{
    Vector3 w = Vector3.Transform(normal, OrthoX);
    float dot = Vector3.Dot(normal, w);
    if (Math.Abs(dot) > 0.6)
    {
        w = Vector3.Transform(normal, OrthoY);
    }
    w.Normalize();

    orthonormal1 = Vector3.Cross(normal, w);
    orthonormal1.Normalize();
    orthonormal2 = Vector3.Cross(normal, orthonormal1);
    orthonormal2.Normalize();
}

Though the above works you may find it doesn't behave as you'd expect. For example, if your quaternion rotates a vector 90 deg. around X and 90 deg. around Y you'll find if you decompose the rotation around Z it will be 90 deg. as well. If you imagine a vector making these rotations then this makes perfect sense but depending on your application it may not be desired behaviour. For my application - constraining skeleton joints - I ended up with a hybrid system. Matrices/Quats used throughout but when it came to the method to constrain the joints I used euler angles internally, decomposing the rotation quat to rotations around X, Y, Z each time.

Good luck, Hope that helped.

Solution 2

There is an elegant solution for this problem, specially suited for quaternions. It is known as the "swing twist decomposition":

in pseudocode

/**
   Decompose the rotation on to 2 parts.
   1. Twist - rotation around the "direction" vector
   2. Swing - rotation around axis that is perpendicular to "direction" vector
   The rotation can be composed back by 
   rotation = swing * twist

   has singularity in case of swing_rotation close to 180 degrees rotation.
   if the input quaternion is of non-unit length, the outputs are non-unit as well
   otherwise, outputs are both unit
*/
inline void swing_twist_decomposition( const xxquaternion& rotation,
                                       const vector3&      direction,
                                       xxquaternion&       swing,
                                       xxquaternion&       twist)
{
    vector3 ra( rotation.x, rotation.y, rotation.z ); // rotation axis
    vector3 p = projection( ra, direction ); // return projection v1 on to v2  (parallel component)
    twist.set( p.x, p.y, p.z, rotation.w );
    twist.normalize();
    swing = rotation * twist.conjugated();
}

And the long answer and derivation of this code can be found here http://www.euclideanspace.com/maths/geometry/rotations/for/decomposition/

Solution 3

minorlogic's answer pointing out the swing-twist decomposition is the best answer so far, but it is missing an important step. The issue (using minorlogic's pseudocode symbols) is that if the dot product of ra and p is negative, then you need to negate all four components of twist, so that the resulting rotation axis points in the same direction as direction. Otherwise if you are trying to measure the angle of rotation (ignoring the axis), you'll confusingly get a mix of correct rotations and the reverse of the correct rotations, depending on whether or not the rotation axis direction happened to flip when calling projection(ra, direction). Note that projection(ra, direction) computes the dot product, so you should reuse it and not compute it twice.

Here's my own version of the swing-twist projection (using different variable names in some cases instead of minorlogic's variable names), with the dot product correction in place. Code is for the JOML JDK library, e.g. v.mul(a, new Vector3d()) computes a * v, and stores it in a new vector, which is then returned.

/**
 * Use the swing-twist decomposition to get the component of a rotation
 * around the given axis.
 *
 * N.B. assumes direction is normalized (to save work in calculating projection).
 * 
 * @param rotation  The rotation.
 * @param direction The axis.
 * @return The component of rotation about the axis.
 */
private static Quaterniond getRotationComponentAboutAxis(
            Quaterniond rotation, Vector3d direction) {
    Vector3d rotationAxis = new Vector3d(rotation.x, rotation.y, rotation.z);
    double dotProd = direction.dot(rotationAxis);
    // Shortcut calculation of `projection` requires `direction` to be normalized
    Vector3d projection = direction.mul(dotProd, new Vector3d());
    Quaterniond twist = new Quaterniond(
            projection.x, projection.y, projection.z, rotation.w).normalize();
    if (dotProd < 0.0) {
        // Ensure `twist` points towards `direction`
        twist.x = -twist.x;
        twist.y = -twist.y;
        twist.z = -twist.z;
        twist.w = -twist.w;
        // Rotation angle `twist.angle()` is now reliable
    }
    return twist;
}

Solution 4

Code for Unity3d

// We have some given data
Quaternion rotation = ...;
Vector3 directionAxis = ...;

// Transform quaternion to angle-axis form
rotation.ToAngleAxis(out float angle, out Vector3 rotationAxis);

// Projection magnitude is what we found - a component of a quaternion rotation around an axis to some direction axis
float proj = Vector3.Project(rotationAxis.normalized, directionAxis.normalized).magnitude;

View more solutions

Comments

[Ben Hymers]

I'm at work at the moment, I'll try this later and if it works I'll be somewhat chuffed :) (meant to put that on a separate line in the previous comment but apparently 'enter' means 'submit' in this particular text box)

[sebf]

Thank you (though check the first question I asked and its answer and see if you still think so! ;)). Just in case, if you need to decompose a quaternion to eulers (which makes it much easier to recreate), ed022 has a really good implementation here: forums.create.msdn.com/forums/p/4574/62520.aspx (17th post). Good luck with your app!

[angularsen]

If axis of observation equals axis of rotation then that would produce a zero-length vector, so this case would have to be checked.

[João Portela]

there is a small problem, the returned value is always positive. which means that if I have a rotation of -PI/4 (or 7*PI/4) I get PI/4 as a result. Am I missing something here?

[sebf]

@João, you are not missing anything but maybe misunderstand slightly how this works :) The core of this method is the acos of the dot product of the resulting vectors - to put it another way - this method does not operate 'directly' on the quaternion but rather 'observes the results' of the application of that quaternion. Therefore, the angle returned will always be the smallest between the two vectors, and will be limited to 0-360 deg. You can recover whether the angle of transformation was negative or positive however using the cross product. Good luck!

[Ben Hymers]

This looks like a much better solution with fewer edge cases, and much more efficient... I like it! When I get time I'll compare with the answer I accepted.

[minorlogic]

Please note paper "SWING-TWIST DECOMPOSITION IN CLIFFORD ALGEBRA" with general derivation

[Admin]

@minorlogic this looks great, but how do we deal with the singularity? I probably misunderstand something, but if this solution breaks some small percentage of the time, then surely it's not very useful?

[minorlogic]

singularity can be easy catched checking "twist " length is not zero. it is not present in code above (to simplify code).

[Admin]

How would I check for the singularity and where would it appear?

[minorlogic]

Singularity appears if rotation close to 180 deg, and rotation axis is orthogonal to "direction" vector. There is a infinity decompositions in this case. It can be checked if magnitude of unnormalized twist close to zero. Than you can select one of possible , valid twist.

[Jespertheend]

@sebf Could you elaborate? Which vectors do you need to apply cross product to in order to find out?

[Jespertheend]

ah, found it. You have to do the cross product between orthonormal1 and flattened

[Tim von Känel]

If i print the eulerAngles of my twist rotation it's (0, 0, x) which is expected. However if i print the eulerAngles of my swing rotation it rotates around the z-axis. Is this expected behaviour?

[minorlogic]

Swing twist is NOT EulerAngels decomposition.

[Luke Hutchison]

This is the best answer, however there is one step missing: if the dot product of ra and p is negative, then you should negate all four components of twist, so that the resulting rotation axis points in the same direction as direction. Otherwise if you are trying to measure the angle of rotation (ignoring the axis), you'll confusingly get a mix of rotations and their inverses, depending on whether or not the rotation axis direction flipped during projection(ra, direction). Note that projection(ra, direction) computes the dot product, so you should reuse it and not compute it twice.

[Luke Hutchison]

I added an answer with my correction for the direction of twist.

[Luke Hutchison]

See my answer, it may solve your problem about result repeatability.

[Aidenhjj]

here's a full implementation of the approach in @minorlogic's answer github.com/dtecta/motion-toolkit/blob/master/jointlimits/…