Prelim
I going to run through some basic quaternion quantization and only considers the worst case situation in terms of error vs. size:
 We have no context or prediction
 We have no known probability distributations
 We only consider quantization/reconstruction (no entropy coding, performace issues, etc)
 Fixed sized bitpackage
Basically we will assume that the quaternions being quantized are from a uniform random source with no restrictions. Or another way they are from an IID point sequence on the 3sphere. Given this toy problem statement I will consider the important metric be minimal peak magnitude error.
Minimal notation summary: I will write a unit quaternion as:
\[\begin{eqnarray*} Q &=& a + \left(x,~y,~z\right) \\ &=& a + b~\mathbf{u} \\ &=& \cos\left(\theta\right) + \sin\left(\theta\right)~\mathbf{u} \end{eqnarray*}\]
So scalars are lowercase, I put the scalar part first, lowercase bold face for 3D (bi)vectors and I use radian angle measures in $\mathbb{H}$ and not $\mathbb{R}^3$ except all error values which will be in degrees in $\mathbb{R}^3$. I will also assume that the scalar is positive, if not the sign of it will need to be carried through.
The following table is a recap of the Halfangle/Cayley transform post^{1}:
Name  Equation  $d\left(\theta\right)$  $ d'\left(\theta\right) $  max RT error 

identity (projected)  $\left<Q\right>_2$  $\sin\left(\theta\right) $  $\cos\left(\theta\right) $  0.055939 
halfangle (projected)  $\left<\sqrt{Q}\right>_2$  $\sin\left(\frac{\theta}{2}\right) $  $\frac{1}{2} \cos\left(\frac{\theta}{2}\right) $  0.000075 
Cayley transform  $\left(Q1\right)\left(Q+1\right)^{1}$  $\tan\left(\frac{\theta}{2}\right) $  $\frac{1}{2} \cos\left(\frac{\theta}{2}\right)^{2} $  0.000045 
harmonic mean  $\left(\sqrt{Q}1\right)\left(\sqrt{Q}+1\right)^{1}$  $\tan\left(\frac{\theta}{4}\right) $  $\frac{1}{4} \cos\left(\frac{\theta}{4}\right)^{2} $  0.000070 
exponential map  $\log Q$  $\theta$  1  0.000087 
Where bivector extraction (rank2) is denoted and expressed algebraically:
All the transforms map a unit quaternion to a point in a 3D ball in the direction of the axis of rotation. I’ll refer to this type of transform as ball maps. The table entries are:
 Distance from the origin of $\theta$ is $d\left(\theta\right)$, which is the cumlative density.
 The derivative $ d'\left(\theta\right) $ is the density.
 max RT error is the maximum round trip error^{2} (without quantization) measured in $\mathbb{R}^3$ and in degrees.
Toy/proofofconcept C code: here
Aside: in memory compression 4:3
Other than identity (drop and reconstruct the scalar) any of the ball maps perform pretty well at round tripping between 3 and 4 component representations. Performance of runtime transform(s) and if the reverse can be folded with the next operation are probably more interesting considerations.
Smallest three
One existing method of quantizing quaternions is called smallest three. This does exactly what it says: determine the component with the maximum magnitude, note which one it is and store the other three components (negated if the largest is negative). This sideline information requires two bits of storage. For whatever it’s worth we can think of this geometically as performing one of four isoclinic rotations by $1, \mathbf{x}, \mathbf{y}$ and $\mathbf{z}$.
The shape generated by this mapping in our case is that of a ball of radius $\frac{\sqrt{2}}{2}$ when then the scalar is the largest component (implied angle is on $\left[0,\frac{\pi}{2}\right]$). Beyond that point the scalar becomes one of the encoded values and the largest bivector component is dropped. This in turn fills the ball and then “inflates” slightly beyond. I haven’t bothered to formulate an equation for the surface but roughly speaking if we put an axe aligned bounding box around the ball then slightly inflate within that contraint is roughly the surface. The approximate volume once we’ve scaled up to unit extent is 5.35284 which translates to using ~66.9% of the quantization space vs ~52.4% for the unit ball assuming uniform grid as we’ll start with.
TaitBryan ZYX and Euler ZXZ
In a prequel post^{3} I ran through the basic math and some computation issues wrt to error for both of these system.
An upside to both of these is they use all of the quantization space without any extra effort. The downside is they are computationally expensive for both encode and decode. Notiably (as presented) converting to ZYX requires two two parameter atans and one one parameter and produces:
\[\theta \in \left[\frac{\pi}{2},\frac{\pi}{2}\right] \\ \phi,\psi \in \left[\pi,\pi\right]\]Converting to ZXZ is slightly nicer only requiring one two parameter atan and two single parameter. However the formulation that I gave in that post does not produce minimal range angles, specifically returning:
\[\beta \in \left[0, \pi \right] \\ \alpha,\gamma \in \left[\frac{3\pi}{2},\frac{3\pi}{2}\right]\]So $\alpha$ and $\gamma$ need to be reduced to $\left[\pi,\pi\right]$. Reducing the values to range can be folded with the argument reduction of the two parameter atan, which I’ll ignore for the moment.
If either of these methods are determined to be of interest we can examine its error performance and potentially replace the arctangents and/or forward trig operations with lower accurancy approximations.
Skipahead
The next few sections are simply a recap of the transforms in the halfangle/Cayley post combined with scaling to the unit ball. So you might want to simply skip to “Example configuration and results” section.
Basic halfangle $\left<\sqrt{Q}\right>_2$
Recall if we take the squareroot of a unit quaternion ($Q$) we half its implied angle and we can reverse this by squaring. By choosing $Q$ with a zero or positive scalar this insures that the scalar ($w$) is the largest component so it can be dropped. The result is a radius $\sin\left(\frac{\theta_m}{2}\right)$ ball which in our case is $\frac{1}{\sqrt{2}}$. We can rescaling the equation so we are mapping to the unit ball.
The forward and inverse transform pairs then become:
\[Q = w + \left(x,~y,~z\right) \rightarrow \frac{1}{\sqrt{1+w}} \left(x,~y,~z\right) = p\] \[p = \left(x,~y,~z\right) \rightarrow \left(1p \cdot p \right) + \sqrt{2p \cdot p} \left(x,~y,~z\right) = Q\]
GLSL code snippets for the forward and inverse transforms (assuming $w$ is already positive)^{4}:
Basic Cayley $\left(Q1\right)\left(Q+1\right)^{1}$
The Cayley transform maps a unit quaternion with a zero/positive scalar to a ball with radius $\tan\left(\frac{\theta_m}{2}\right)$ which is the unit ball in the case we are considering. The forward and inverse transforms are:
GLSL code snippets for the forward and inverse transforms (assuming $w$ is already positive)^{4}:
Basic harmonic mean $\left(\sqrt{Q}1\right)\left(\sqrt{Q}+1\right)^{1}$
The harmonic mean is the compostion of halfangle followed by Cayley. This produces a radius $\tan\left(\frac{\theta_m}{4}\right)$ ball which in our case is: $\sqrt{2}1$.
Rescaling to the unit ball gives the forward and inverse transform pairs:
\[Q = w + \left(x,~y,~z\right) \rightarrow \frac{1+\sqrt{2}}{1+w+\sqrt{2+2w}} \left(x,~y,~z\right) = p\] \[p = \left(x,~y,~z\right) \rightarrow \frac{16a+a^2}{\left(1+a\right)^2} + \frac{4\left(\sqrt{2}1\right)\left(1a\right)}{\left(1+a\right)^2}\left(x,~y,~z\right) = Q\]where:
\[\begin{eqnarray*} a &=& \left(32\sqrt{2}\right)p \cdot p \end{eqnarray*}\]
GLSL code snippets for the forward and inverse transforms (assuming $w$ is already positive)^{4}:
Basic exponential map $\log Q$
Taking the logarithm produces a $\theta_m$ radius ball, so in our case we have $\frac{\pi}{2}$.
Rescaling to the unit ball gives the forward and inverse transform pairs:
\[Q = w + \left(x,~y,~z\right) \rightarrow \frac{\frac{2}{\pi}~ \text{atan}\left(\frac{\sqrt{1w^2}}{w}\right)}{\sqrt{1w^2}} \left(x,~y,~z\right) = p\] \[p = \left(x,~y,~z\right) \rightarrow \cos \left(\frac{\pi}{2} \sqrt{p \cdot p}\right) + \frac{\sin\left( \frac{\pi}{2} \sqrt{p \cdot p} \right)}{ \frac{\pi}{2} \sqrt{p \cdot p}} \left(x,~y,~z\right) = Q\]
GLSL code snippets for the forward and inverse transforms (assuming $w$ is already positive)^{4}:
Example configuration and results
Now let’s look at some basic quantization with:
 allowing an arbitrary number of bits in the fixed package (not requiring each component to have equal number of bits)
 linear quantized with round to negative infinity for quant and center reconstruct for dequant. This choice means we cannot exactly represent zero and therefore identity. I’m not considering this a problem since this post isn’t useful for relative/update information.
 on the standard integer lattice (Cartesian/uniform grid)
Let’s examine how the maximum error globally behaves as we add bits to each transform, starting from 9 up to 32 bits. If we have a $n$ bit package and $n \bmod 3 = 1$ then the last component gets an extra bit and if $n \bmod 3 = 2$ then the last two get an extra bit each. As a log plot we get:
The overall is that this follows the reasonable expectation that the error should approximately linearly decrease in a logplot since adding a bit doubles the number of representable values. We can also see that ZYX and ZXZ are very close with ZXZ always the slightly better. Since quaternion to ZXZ can be formulated with lower runtime cost than ZYX we can stop considering ZYX. If we examine the methods with the lowest error at each bit size we will notice a pattern emerges. A table of the top three:
$n \bmod 3$  best  second  third 

0  expmap  harmonicmean  ZXZ 
1  expmap  ZXZ  harmonicmean 
2  ZXZ  expmap  harmonicmean 
Now let’s take the previous graph and view the relative effectiveness of adding a bit $\left(\frac{\text{error}\left(n+1\right)}{\text{error}\left(n\right)}\right)$:
This is a bit of a mess because halfangle isn’t well behaved here. Click on it in the legend to turn off that trace. Now the graph is pretty much two sets of traces: ZXZ vs. the ball maps and the $n \bmod 3$ behavior is clear.
$n \bmod 3$  maps  ZXZ 

0  0.7071  0.9220 
1  0.8660  0.7670 
2  0.8165  0.7071 
For the ball maps all three components are equally important so the largest improvement is when the points form a cube (mod is zero). ZXZ has one angle with half the range so it’s best is when mod is 2. Adding three bits halves the error for both types.
Although the toy problem statement only cares about global peak error let’s go ahead and examine error as a function of rotation angle. If we reexamine the density function from Halfangle/Cayley post we expect that higher density should result in lower errors:
Now a sequence of error vs rotation angle plots for each mod n behavior:
These plots still don’t show the whole picture. The part we are ignoring (for ball maps) is the impact of the direction of the bivector (axis of rotation). Some things to note:
 error plots of the ball maps are inversely related to density as expected and they behave the same regardless of bitsize
 ZYZ in it’s various position in the top three
 I’m showing identity map projected to ball which is useless for full range input. The point of including it is to compare with smallestthree. Specifically in the 32bit case there’s only mimimal impact to error on the range where they are the same transform. So sideline information can be effective espectially when they are replacing what would otherwise be less effective given the bit allocation.
Ideally I would now show some 3D density (heatmap) plots on the unit ball with a cutout for the various ball maps. I’m too lazy to do that ATM, so instead let’s look at best and worst case directions. The issue is our quantization space is currently a cube with a volume of eight, we’re mapping to a unit ball (volume $\frac{4}{3}\pi \approx 4.189$) so we are wasting $1\frac{\pi}{6}$ or approximately 47.6% of the space. If we have $n$ sample in one extent, then that’s how many we’ll pass for an axial aligned rotation and only $\frac{1}{\sqrt{3}}n$ if the axis is along a diagonal (and they’re farther apart). A 12bit encoding of the exponential map clearly shows:
For ball maps with 30 bit encoding:
Note the behaviour of the halfangle transform (click away the others). We can leave aligned quantization points alone and improve along the diagonal by adding radial stretching to the encode/decode. The global, along x and diagonal now look like:
The error along the diagonal has been overkill reduced. It seems likely some tweaks could produce a better global result by moving some of the point away from the diagonal toward higher peak error directions, but I’m not going to make the effort without looking at some density plots. Let’s look at halfangle plus radial stretch vs. the previous top performing methods at 30, 31 and 32bit packages:
So halfangle + radial stretch takes over the top spot for 30bit, is very close second at 31bits and loses it completely for 32bit packages. Recall that halfangle isn’t well behaved at $n \bmod 3$ effectiveness, so these different package size choice will vary.
Smallestthree encoded with 32bits shows that someone motivated could probably give it a nice bump in a similar manner.
Breaking point
That’s enough for one post. I do have some more bits and pieces on this topic before moving to context and prediction based.
References and Footnotes

“Quaternion half/double angle and Cayley transforms” (local post) ↩

Measured using the toy C implemenations. ↩

“Converting to Euler & TaitBryan” (local post) ↩

Otherwise replace $w$ by $\abs{w}$ and add $\text{sgn}\left(w\right)$ to the final scale value. In this case the inverse transform returns $Q$ when $w$ is negative. ↩ ↩^{2} ↩^{3} ↩^{4}