Difference between two quaternions I am currently writing up a program that needs to know the difference in degrees for each axis between two quaternions.
Now, suppose I have a starting quaternion Qs and I need to compute at each step the difference between my current orientation represented by the quaternion Qc.
I do that: Q = Qc^-1 * Qs
Now I would like to get values from the quaternion Q that would represent for each axis the number of degrees that differ from Qc to Qs. Is that possible?
Edit: If that's not clear here are other explanations
I know I can get the angle value from a quaternion by doing 2 * cos-1(q.w). Now I would like to get an idea on how to project this angle onto each axis. The expected output is something like:

From Qs to Qc there is a rotation of 30 degrees on the x axis, 20 degrees on the y axis and 0 degrees on the z axis

 A: It sounds like you are just looking for a way to convert your quaternion $Q$ to Euler angles, and there are numerous webpages on the internet that cover this very well already. Just google "quaternions to Euler angles" and you get things like these:
http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToAngle/
https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
http://www.chrobotics.com/library/understanding-quaternions
A: So if I interpret your question right
you already have the axis in ${\bf i,j,k}$ components (if you disregard the factor $\sin(\theta/2)$), so just scrap the real part (pretend it's zero).
First let $\bf p_1,p_2$ be the quaternions after they have had a real part stripped away (and the rest renormalized so that square sum = 1).
We then have (from Wikipedia):
$${\bf p_1' = qp_1q^{-1} \Leftrightarrow p_1'q=qp_1 \hspace{1cm}} (1)$$
where $\bf q$ is one of the quaternions 
$${\bf q} = \underset{a}{\underbrace{\cos(\alpha/2)}} + \underset{b}{\underbrace{\sin(\alpha/2)}}  (b_1{\bf i}+ b_2 {\bf j}+ b_3{\bf k})$$
for the set of binary vector $${\bf b} = (b_1,b_2,b_3)^T \text{ such that } \cases {b_k \in \{0,1\}\\ \sum_{\forall i} b_i = 1}$$
or just simply $${\bf b} \in \{(1,0,0)^T,(0,1,0)^T,(0,0,1)^T\}$$
We now want to find the $\alpha$s for these so that the resulting rotation is as close as possible
$$\min_{\bf q} \left\|\bf qp_1q^{-1}-p_2 \right\| \Leftrightarrow \min_{\bf q} \left\|\bf qp_1-p_2q \right\|\hspace{1cm} (2)$$
is linear in the coefficients. Let us call them $a = \sin(\alpha), b = \cos(\alpha)$. However, we will at some point need to add some functionality to ensure $a^2+b^2 = 1$ (the trig 1, right?) and then it becomes non-linear.
A reasonably good approach would probably be to disregard the constraint, solve the linear (least squares) system then impose the constraint for the solution found by doing a renormalization update $$\cases{a = \frac{a}{a^2+b^2}\\b=\frac{b}{a^2+b^2}}$$ and then if needed iterate from the new position to refine the solution found.
A: Perhaps worth mentionning here, the angular difference between two rotations represented by unit norm quaternion is
$$d\left(\mathbf{q}_1, \mathbf{q}_2\right) = 2 \arccos\left(\mathbf{q}_1^T\mathbf{q}_2\right)$$
