Numerically stable extraction of Axis-Angle from Unit Quaternion

You cannot do better than the formula $$ \omega = \frac{\begin{bmatrix} x & y & z\end{bmatrix}^T}{\|\begin{bmatrix}x & y & z\end{bmatrix}\|}$$ since all the information about $\omega$ is encoded in the relative magnitudes of $x,y,z$. Assuming these are stored in floating point with each component having a relative accuracy $\epsilon$, then the formula should be able to recover $\omega$ to the same relative accuracy. All you need is a robust way to compute the norm in the denominator. See for example LAPACK's DLAPY3 routine. If the norm is exceedingly small, then this corresponds to an exceedingly small rotation angle, and so the axis is ill-conditioned and essentially arbitrary corresponding to not-a-rotation.

Ultimately, it's a question of what you're trying to do with the axis-angle pair. If you need the axis to high absolute accuracy, then you can't represent it using a unit quaternion without losing accuracy in some extreme cases. This should only happen if the sine is near underflow (within a few orders of magnitude of DBL_MIN; even near it, gradual underflow will still help you to some extent).

The standard mathematical equation to convert a quaternion into axis-angle form is numerically unstable when the vector part of the quaternion is small. The correct way to handle this calculation is to use the Taylor expansion of the equation when the vector part of the quaternion is small.

This numerically-stable approach is presented clearly in Neil Dantam's course notes for Georgia Tech. See equations 50-57. http://www.neil.dantam.name/note/dantam-quaternion.pdf