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/