It looks like the rotation order is $y x z$, but angles negated for some reason. In other words, the rotation matrix for Euler angles $(\chi,\gamma,\theta)$ used is most likely
$$\mathbf{R}(\chi,\gamma,\theta) = \left[\begin{matrix}
\cos \gamma \cos \theta - \sin \chi \sin \gamma \sin \theta & \cos \gamma \sin \theta + \sin \chi \sin \gamma \cos \theta & -\cos \chi \sin \gamma \\
-\cos \chi \sin \theta & \cos \chi \cos \theta & \sin \chi \\
\sin \gamma \cos \theta + \sin \chi \cos \gamma \sin \theta & \sin \gamma \sin \theta - \sin \chi \cos \gamma \cos \theta & \cos \chi \cos \gamma
\end{matrix}\right]$$
The rotation matrix corresponding to versor aka unit quaternion $\mathbf{Q} = (w, x, y, z)$ is
$$\mathbf{R}(w,x,y,z) = \left[\begin{matrix}
1 - 2 (y^2 + z^2) & 2 (x y - z w) & 2 (x z + y w) \\
2 (x y + z w) & 1 - 2 (x^2 + z^2) & 2 (y z - x w) \\
2 (x z - y w) & 2 (y z + x w) & 1 - 2 (x^2 + y^2)
\end{matrix}\right]$$
The stated problem is to find the correspondence between the two.
Looking at $\mathbf{R}_{23}$ we can tell that
$$\chi = \arcsin(2 y z - 2 x w)$$
Since $\mathbf{R}_{21}/\mathbf{R}_{22} = -\tan\theta$ if $\mathbf{R}_{22} \ne 0$,
$$\theta = -\left\lbrace\begin{align}
& -\arctan\left(\frac{x y + z w}{1/2 - x^2 - z^2}\right),\; & x^2 + z^2 \ne 1/2 \\
& -\frac{\pi}{2}, \; & x y + z w \ge 0 \\
& \frac{\pi}{2}, \; & x y + z w \lt 0
\end{align}\right.$$
Similarly, $\mathbf{R}_{13}/\mathbf{R}_{33} = -\tan\gamma$ if $\mathbf{R}_{33} \ne 0$,
$$\gamma = \left\lbrace\begin{align}
& -\arctan\left(\frac{x z + y w}{1/2 - x^2 - y^2}\right),\; & x^2 + y^2 \ne 1/2 \\
& -\frac{\pi}{2}, \; & x z + y w \ge 0 \\
& \frac{\pi}{2}, \; & x z + y w \lt 0 \end{align}\right.$$
As a result, I suggest trying the following pseudocode:
epsilon = 0.00001
halfpi = 0.5 * 3.14159265358979323846
function euler_angles(w, x, y, z):
temp = 2*(y*z - x*w)
if (temp >= 1 - epsilon):
xa = halfpi
ya = -atan2(y, w)
za = -atan2(z, w)
else if (-temp >= 1 - epsilon):
xa = -halfpi
ya = -atan2(y, w)
za = -atan2(z, w)
else:
xa = asin(temp)
ya = -atan2(x*z + y*w, 0.5 - x*x - y*y)
za = -atan2(x*y + z*w, 0.5 - x*x - z*z)
return (xa, ya, za)
which utilizes the atan2(dy,dx) function available in most programming languages.
(Above pseudocode edited on 2017-09-30, to add the checks for gimbal lock, where xa = ±0.5 π
. The epsilon
represents the machine epsilon in this context, and depends on the precision of the variables; typically, it should be a compile-time constant that is easy to modify when porting or running the code on different architectures. The value shown should work well for visualization.)
The corresponding conversion from Euler angles to a quaternion is easy to compute, if one realizes that the three rotations as quaternions are easy to express: $$\begin{align}
&\mathbf{R}_x(\chi) = \left( \cos\frac{\chi}{2} + \sin\frac{\chi}{2} \mathbf{i} \right )\\
&\mathbf{R}_y(\gamma) = \left( \cos\frac{\gamma}{2} + \sin\frac{\gamma}{2} \mathbf{j} \right )\\
&\mathbf{R}_z(\theta) = \left( \cos\frac{\theta}{2} + \sin\frac{\theta}{2} \mathbf{k} \right )\end{align}$$
Their Hamilton product $\mathbf{Q} = \mathbf{R}_y(\gamma)\mathbf{R}_x(\chi)\mathbf{R}_z(\theta)$ gives the desired quaternion here.
In pseudocode,
function quaternion(xa, ya, za):
cx = cos(-0.5*xa)
sx = sin(-0.5*xa)
cy = cos(-0.5*ya)
sy = sin(-0.5*ya)
cz = cos(-0.5*za)
sz = sin(-0.5*za)
w = cx * cy * cz + sx * sy * sz
i = cx * sy * sz + sx * cy * cz
j = cx * sy * cz - sx * cy * sz
k = cx * cy * sz - sx * sy * cz
return (w, i, j, k)
Remember that both of these functions depend on the rotation order, and that the order of rotations here is $y x z$ with negated angles.
Best Answer
I'll assume the convention used in this earlier answer and this earlier answer about the use of quaternions for rotation, that is, for an initial object vector $v$ and total rotation quaternion $q_w$ you get the rotated object vector by multiplying $q_w v q_w^{-1}.$
Part 1
Suppose you have a "parent" rotation with quaternion $q_p,$ and you want to combine this with a "child" rotation with quaternion $q_c$ that will rotate an object relative to the object's orientation after the "parent" rotation was performed.
Multiplication by quaternions in the conventional way transforms the world coordinates of whatever object you apply it to. In order to achieve the effect of first performing the parent rotation, then performing the child rotation relative to the rotated coordinate system that resulted from the "parent" rotation, you should apply the "child" rotation first in world coordinates.
Using the convention $q_w v q_w^{-1},$ whichever quaternion is on the right when we compute $q_w$ is the rotation that is performed first, so you want $q_w = q_p q_c.$ That way, when you compute $q_w v q_w^{-1},$ it is actually $$ q_w v q_w^{-1} = (q_p q_c) v (q_p q_c)^{-1} = (q_p q_c) v (q_c^{-1} q_p^{-1}) = q_p (q_c v q_c^{-1}) q_p^{-1}. $$
Part 2
If you have the total rotation quaternion $q_w$ and you also know the parent rotation quaternion $q_p,$ to recover the child rotation quaternion you can multiply on the left by $q_p^{-1},$ because then $$ q_p^{-1} q_w = q_p^{-1} (q_p q_c) = (q_p^{-1} q_p) q_c = q_c. $$
Part 3
If you have the total rotation quaternion $q_w$ and you also know the child rotation quaternion $q_c,$ to recover the parent rotation quaternion you can multiply on the right by $q_c^{-1},$ because then $$ q_w q_c^{-1} = (q_p q_c) q_c^{-1} = q_p (q_c^{-1} q_c) = q_p. $$
Checking the algorithm
After constructing a set of formulas like this, it's a good policy to verify that you put them together properly by taking some examples of "object" vectors and applying some examples of rotations to them.
A nice simple example is a parent rotation of $\frac\pi2$ ($90$ degrees) around the $z$ axis that takes points on the positive $x$ axis to points on the positive $y$ axis, and a child rotation of $\frac\pi2$ ($90$ degrees) around the $x$ axis that takes points on the positive $y$ axis to points on the positive $x$ axis. If you perform the parent rotation first, then perform the child rotation relative to the result of the parent rotation, the parent rotation will leave the $z$ axis where it was, but the child rotation will be about the $y$ axis in world coordinates (because that's the rotated position of the original $x$ axis), so it will take points on the positive $z$ axis to the positive $x$ axis. If you apply these rotations to a point on the positive $z$ axis and it ends up somewhere other than the positive $x$ axis, you will know you need to fix your implementation of the rotations somehow.
You can look at the rotated positions of some vectors after any one of the rotations in question (the parent rotation alone, the child rotation alone, or the combined parent-child rotation) to verify that each one is doing what you expect.