This post is about two transform pairs of quaternion that reduce the working domain which are not commonly discussed. Specifically they both reduce the implied angle in quaternion space. This allows storing a quaternion value with three elements more attractive, but more importantly reduces the range of the implied angle. This angle reduction can be useful for interpolation, quantization and shaping distributions for example. The usefulness of these maps for quantization and interpolation depend on the target error bound and maximum implied angle the system allows. There are relatively well known techniques which use unit magnitude and that $Q$ and $-Q$ are equivalent for alternate storage and/or quantization. I will later examine the case of interpolation and will focus on the basic properties and present some simple quantization.

Once we get past the preliminary junk we’ll be almost entirely talking about a 2D plane and complex numbers. I will also drop most of the formalism. I’m also going to skip on actually defining the sets and will simply note how I name them and they are hopefully obvious by context.

I’m using pseudo-GLSL for code snippets in an attempt at clarity. I will make very few nods to optimizations other than carrying through some of the derivations.


Preliminary stuff


We generally work with quaternions $\left(\mathbb{H}\right)$ represented by four real values:

\[\begin{equation} \label{qxform} Q = a + \bvec{x,~y,~z} \end{equation}\]


Where I will call “$a$” the scalar part1 and $\bvec{x,y,z}$ the bivector part. For our purposes here we can consider a 3D bivector to be equivalent to a 3D vector.

We will also want to express a quaternion in terms of some implied unit bivector $\mathbf{u}$ $\left( \mathbb{S}^{2} \right) $ with implied magnitude $b$:

\[\begin{equation} \label{eq:qsb} Q = a + b~\mathbf{u} \end{equation}\]


and in polar form2 with an implied angle (measurement in quaternion space) of $\theta$:

\[\begin{equation} \label{eq:qpolar} Q = \cos\left(\theta\right)+\sin\left(\theta\right)~\mathbf{u} \end{equation}\]


Since all non-zero scalar multiples of a quaternion represent the same rotation (projective line3 $sQ$, for unit quaternion $Q$ and non-zero scalar $s$) so we can restrict ourselves to unit quaternions. Equation $\eqref{eq:qpolar}$ already makes this assumption. This gives the four dimension unit sphere $\left( \mathbb{S}^{3} \right) $ as a working domain.

The projective line intersects the unit-sphere at $Q$ and $-Q$, so we can further restrict ourself to quaternions with scalars part that are zero or positive $\left( a \geq 0 \right)$. In axis-angle think this is negating the angle and reversing the axis. The reduces the working domain to the half sphere $\left( \mathbb{S}^{3+} \right) $. $ \newcommand{\hcomplex}[1]{\mathbb{C}_{#1}} $ So far all standard fare.

If we take the set of all quaternion with fixed $\mathbf{u}$ definable in forms $\eqref{eq:qsb}~\eqref{eq:qpolar}$ they form a plane which spans the line-of-reals and $\mathbf{u}$. This plane forms a commutative subalgebra of $\mathbb{H}$ isomorphic to $\mathbb{C}$, we can denote this plane as $\hcomplex{u}$ or as the {1,Q}-plane. The set of all quaternions is the collection of all unique complex planes and they intersect at the line of reals.

The transform pairs which we will be considering are what I call (quaternion valued) complex functions which simply means if we have a function $f$ which given input in some complex plane then the result is in the same complex plane. In other words we are working in two dimensions using the standard set of complex numbers.

We are at the following in terms of $\eqref{eq:qsb}$:

  • effectively working with the complex number: $z = a + b~\mathbf{i} $
  • unit magnitude: $ a^2 + b^2 =1 $
  • positive scalar: $ a $ on $\left[0,1\right] $
  • $ b $ on $\left[0,1\right] $, since $b=\sqrt{x^2+y^2+z^2}$. Any sign of $b$ is contained by the implied $\mathbf{u}$
  • the implied angle $\theta$ is therefore on $\left[0,~\frac{\pi}{2}\right]$
  • the domain is the arc of unit circle in the first quadrant, axes included.

The presented forward transforms are also in the first quadrant, axes included.



Half/double angle


Quaternions naturally contain a half-angle/double-angle relation. To convert from an axis-angle representation we halve the angle of the rotation. If we have a rotation represented in axis-angle form with the axis represent by the unit vector $\hat{u}$ and an angle of rotation $\alpha$, then the conversion to quaternions can be expressed as:

\[\begin{equation*} \cos\left(\frac{\alpha}{2}\right)+\sin\left(\frac{\alpha}{2}\right)~\mathbf{u} \end{equation*}\]


where the unit vector $\hat{u}$ and unit bivector $u$ are the same $(x,y,z)$ values. I will from now on use $\alpha$ to indicate an angle measure in 3D and continue to use $\theta$ for in quaternion space with $\theta = \frac{\alpha}{2}$ and $\alpha = 2\theta$ in the standard (untransformed) representation.

In turn this half-angle representation is doubled when we apply the similarity transform4 (the rotation matrix conversion is the similarity transform applied to an orthonormal basis).

\[\begin{eqnarray*} A & = & QAQ^{-1} \\ & = & QAQ^{*} \end{eqnarray*}\]

We can introduce a second half-angle/double-angle transform pair. In complex number these translate to square root and squaring (or trig identities or geometry). The half angle transform is:

\[\begin{eqnarray*} \sqrt{z} & = & \sqrt{a + b~\mathbf{i}} \\ & = & \frac{1}{\sqrt{2+2a}} \left(1+a + b~\mathbf{i} \right) \end{eqnarray*}\]


And the double-angle transform:

\[\begin{eqnarray} z^2 & = & \left(a + b~\mathbf{i}\right)^2 \nonumber \\ & = & a^2 - b^2 + 2ab~\mathbf{i} \nonumber \\ & = & 1-2a^2 + 2ab~\mathbf{i} \label{zsq} \\ & = & 1-2b^2 + 2b\sqrt{1-b^2}~\mathbf{i} \label{zsqBv} \end{eqnarray}\]


Equation $\eqref{zsq}$ is for reconstructing from all four components and $\eqref{zsqBv}$ for when we only have the values of the bivector part.

After the half-angle transform we have:

  • $ a^2 + b^2 =1 $
  • $\theta \in \left[0,~\frac{\pi}{4}\right]$
  • $ a \in \left[\frac{1}{\sqrt{2}},1\right] $
  • $ b \in \left[0, \frac{1}{\sqrt{2}}\right] $
  • $ a \ge b $


For the complex plane as a whole we are halving the angle with respect to the line of reals and the magnitude is square-rooted. This maps the negative half-plane to the positive. Steven Wittens’ blog post “How to Fold a Julia Fractal” has some handy mathbox visualizations in “Like Hands on a Clock” section. From step 29->30 shows the square-root (with negative $x$ pointing down).

Remember in the plane all points along the projective line are equivalent. We can think of the projective line as being two half-lines on either size of the origin that have an angle of $\pi$ between them. After the half-angle transform these half lines have a relative angle of $\frac{\pi}{2}$ and since nothing has happened to the properties in the plane, so we now have two orthogonal projective lines all of which are equivalent representations. That means four unit quaternions are equivalent. (“Like Hands on a Clock”: step 30->31). When we square again these two images of the projective line become one again.

The half-angle transform does not effect angle measurements between bivectors so we cannot simply treat half-angle representations as the standard and expect everything to work. As an example we cannot simply compose rotations with a product. In this additional half-angle representation then the proper rotation is $ AABB $. If we attempted to compose with product and then square we would have $ ABAB $. I will say more about this in later posts.


// quaternion 'q' as vec4: q.w = scalar part 'a' , q.xyz = bivector part 'b'

// sqrt(a + bu)
vec4 q_usqrt(vec4 q) {
  float d = 1.0+q.w;           // 1+a
  float s = inversesqrt(d+d);  // 1/sqrt(2+2a)
  vec3  b = s*q.xyz;           // b/sqrt(2+2a) u
  return vec4(b, d*s);         // (1+a)/sqrt(2+2a) + b/sqrt(2+2a) u
}

// (a + bu)^2
vec4 q_upow2(vec4 q) {
  float s = q.w+q.w;           // 2a
  float w = s*q.w;             // 2a^2
  vec3  b = s*q.xyz;           // 2ab u
  return vec4(b, 1.0-w);       // (1-2a^2) + 2ab u
}

// (a + bu)^2  where we only have the 'bu' term
vec4 q_upow2b(vec3 v) {
  float d = dot(v,v);               // b^2
  float s = sqrt(1.0-d);            // sqrt(1-b^2)
  return vec4((s+s)*v, 1.0-(d+d));  // (1-2b^2) + 2b sqrt(1-b^2) u
}



Log/exp transform


Another pair of complex functions are $\log$ and $\exp$. For unit quaterions (complex numbers) $\log$ can be expressed in terms of the polar form \eqref{eq:qpolar} as:

\[\begin{eqnarray*} \log\left( \cos\left(\theta \right) + \sin\left(\theta \right)\mathbf{i} \right) & = & \theta~\mathbf{i} \\ \end{eqnarray*}\]


which maps a unit quaternion to a bivector. The reverse transform $\exp$ when given a bivector can be expressed as:

\[\begin{eqnarray*} \exp\left(\theta~\mathbf{i}\right) & = & \cos\left(\theta \right) + \sin\left(\theta \right)\mathbf{i} \end{eqnarray*}\]


which is simply another (and common) way to express the polar form. Notice that $\log$ is implicitly the conversion from an axis-angle representation and is sometime called Euler form5. Basically it is simply a connection between the algebra and trig operations. Since we are on a limited range for both transforms we could use a combination of trig identities and function approximation to come up with reasonably accurate and fast versions, however here we are only going to consider these as the “ideal target” transforms since $\log$ perfectly linearizes the angle where we have:

  • $ a = 0 $
  • $ b \in \left[0, \frac{\pi}{2}\right] $


For the space as a whole we are mapping our 4D unit half-sphere input to a 3D ball with a radius of $\frac{\pi}{2}$.



Cayley transform


The Cayley transform6 dates back to the 1840s7. It is a special case of a stereographic projection8 and it is a Möbius transformation9 (linear fractional transformation). It is also the connection between the Lie algebra and group and is related to tangent half-angle substitution10. Representing rotations as the Cayley transform of unit quaternions is called Cayley-Klein parameters5.

If we call the forward transform $f$ and the reverse $g$ then we have:

\[\begin{eqnarray*} f\left(z\right) & = & \frac{z-1}{z+1} \\ & = & \left(z-1\right)\left(z+1\right)^{-1} \\ & = & \frac{\left(z-1\right)\left(z^*+1\right)}{\left(z+1\right)\cdot\left(z+1\right)} \\ \\ g\left(z\right) & = & \frac{1+z}{1-z} = \frac{\left(1+z\right)\left(1-z^*\right)}{\left(1-z\right)\cdot\left(1-z\right)} \end{eqnarray*}\]


Expanding the forward transform with $\eqref{eq:qsb}$:

\[\begin{eqnarray*} f\left(a+b~\mathbf{i}\right) & = & \frac{\left(a-1+b~\mathbf{i}\right)\left(a+1-b~\mathbf{i}\right)}{\left(a+1+b~\mathbf{i}\right)\cdot\left(a+1+b~\mathbf{i}\right)} \\ & = & \frac{2b}{\left(a+1\right)^2+b^2}~\mathbf{i} \\ & = & \frac{2b}{2\left(1+a\right)}~\mathbf{i} \\ & = & \frac{b}{1+a}~\mathbf{i} \end{eqnarray*}\]


Like $\log$ the result is a bivector. In polar form we have:

\[f\left(\cos\left(\theta\right)+\sin\left(\theta\right)~\mathbf{i}\right) = \frac{\sin\left(\theta\right)}{1+\cos\left(\theta\right)}~\mathbf{i} = \tan\left(\frac{\theta}{2}\right)~\mathbf{i}\]


We only need the reverse transform for bivector input:

\[\begin{eqnarray*} g\left(b~\mathbf{i}\right) & = & \frac{\left(1+b~\mathbf{i}\right)\left(1+b~\mathbf{i}\right)}{\left(1-b~\mathbf{i}\right)\cdot\left(1-b~\mathbf{i}\right)} \\ & = & \frac{1-b^2+2b~\mathbf{i}}{1+b^2} \\ & = & \frac{2}{1+b^2} - 1 + \frac{2}{1+b^2}b~\mathbf{i} \end{eqnarray*}\]


NOTE: In general (since the quaternion product does not commute) division is an ill defined operation. It can mean either a product on the left or right by the inverse. However in this case we have a “complex function” so the product does commute and there is no ambiguity: $\left(Q-1\right)\left(Q+1\right)^{-1} = \left(Q+1\right)^{-1}\left(Q-1\right)$. Some text incorrectly state otherwise. I chose the division form since we should be thinking in complex numbers.

After the forward transform we have:

  • $ a = 0 $
  • $ b \in \left[0, 1\right] $

For the space as a whole we are mapping our 4D unit half-sphere input to the 3D unit ball. If we were to apply the forward transform above to the negative 4D half-sphere it maps to all of the 3D space outside of the unit ball. Inverting the original transforms is the opposite transform pair (negative to unit ball, positive to outside).


vec3 q_cayley(vec4 q) { return q.xyz*(1.0/(1.0+q.w)); }

vec4 q_cayley(vec3 v) {
  float s = 2.0/(1.0+dot(v,v)); // 2/(1+b^2)
  return vec4(s*v, s-1.0);
}



Recap and Cayley visualization


The following is a shadertoy to illistrate the Cayley transform and some basic concepts up to this point. It starts paused. Left-mouse down normalizes the direction to the pointer and set $Q$ to that point on the circle/sphere.


Recap

  • The plane of the figure is our complex plane $\hcomplex{u}$, with $u$ is in the positive $y$ direction and positive reals in the positive $x$ direction.
  • The unit circle is the intersection of the 4D unit sphere with this plane.
  • All rotations about $u$ are in this plane.
  • The two illustrated points on the circle are two unit quaternions $Q$ and $-Q$.
    • The point colored blue has a positive or zero scalar part. The one colored red has a negative scalar part.
    • They are initially $\pm 1$ so they are in all of the complex planes of the space.
    • The line between these points (except the origin) are all equvalient.
    • When the blue point is negative in $y$ we are logically working with $-u$ instead (our figure is upside down)



Cayley

The forward function $f$ is the stereographic projection from $(-1,0)$ onto line $x=0$ ($y$-axis). If we mouse-down and move in a counter-clockwise manner we have a blue line from $(-1,0)$ to our blue point. Where the line crosses the $y$ axis is a green point which is $f\left(Q\right)$. These two points meet at $\frac{\pi}{2}$. If we continue on the blue point changes to red and that green point exits the circle. If we repeat this process there is a second blue line through the red point which is $f\left(-Q\right)$. (If the two blue or magenta lines aren’t orthogonal it’s me being sloppy in the shader…move the mouse a bit). This second line does not appear at $\pm 1$, artistic choice (the red point is at infinity) or I’m being sloppy with the math in the shader, you choose. As we start to approach $\frac{\pi}{2}$ this second green point appears at the bottom of the screen and enters the circle when the red point changes to blue.

The magenta lines are the sterographic projection $(1,0)$ onto the $y$-axis, which are $\frac{1}{f\left(\pm Q\right)}$ respectively. As an aside we have: $f\left(Q\right) = \frac{1}{f\left( -Q\right) ^{*}}$



Comparison of mappings


These figures use plotly.js. You can click on the function in the legend to hide/show that plot.


The first figure is a simply comparison of how the magnitude of the bivector behaves with respect to the implied angle. The $x$-axis is the normalized angle, so in terms of the represented angle of rotation in 3D we would multiply by $\pi$. The $y$-axis is normalized magnitude of bivector (maximum value of $b$ divided out). If we were to put an upper limit on the implied angle then these plots would grow closer to that of $\log$.


The second figure is the same as the first with simply the linear part subtracted out. If we were to put an upper limit on the implied angle then the shape of plots stay pretty much the same.


We can think of the first figure as being plots of cumulative density where the $y$ value is “how many” and the corresponding $x$ is “where”. As an example if we are linearly quantizing ($N$ samples) and examine $x$ at 0.5 then the $y$ value of each plot is a multiplier $m$ and $mN$ is the number of samples on normalized angle range [0, 0.5]. Or examine the $y$ at 0.5 and then $x$ tells us how far along the angle range the first half of samples represents.

By taking the derivative we can examine localized density which gives a much better picture of how samples are distributed across the full range. Like the second figure the shape of the plots do not change much when adding a cap to maximum implied angle.

The density plots show how the representations behaves locally (in the neighboorhood) of a given angle. For quantization if we linearly quantize the magnitude of the bivector these show the distribution of samples.

Since $\log$ perfectly linearizes the angle we have a nice constant value of one and any samples are equidistributed across the range. The untransformed “standard” representation has a high density around zero and approaches zero at $\alpha = \pi$. The “half-angle” slight favors small angles and “Cayley” prefers large angles.

To be able to “eyeball” how density values change if you place a cap on maximum angle the following show how the untransformed range behaves. Notice the values at .5 vs the endpoints of “half-angle” above. And the “half/half” and “half/Cayley” plots are more to give a feel for evolution of density than to suggest actual usage.


The “take-away” from these should indicate that the “ideal” situation is not use a transform for quantization. The better solution in terms of minimizing quantization error is to have good context, prediction/dead-reckoning and store/transmit relative quaternions which are composed.



Bivector part


None of these transforms effect the implied unit bivector $u$. Accounting for symmetry we can consider the unit bivector $u = \bvec{x,y,z}$ where $x\geq y\geq z \geq 0$, then we have:

\[\begin{eqnarray*} x \in \left[\frac{1}{\sqrt{3}},~1\right] & \approx & \left[0.57735,~1\right] \\ y \in \left[0,\frac{1}{\sqrt{2}}\right] & \approx & \left[0,~.70711\right] \\ z \in \left[0,\frac{1}{\sqrt{3}}\right] & \approx & \left[0,~0.57735\right] \end{eqnarray*}\]


so the largest magnitude part of the bivector has the smallest range ($1-\frac{1}{\sqrt{3}} \approx 0.42265$)



References and Footnotes

  1. scalar locally just means a real value more generally its the working base field. 

  2. Some text prefer $e^{u\theta} $. 

  3. Projective line: Wikipedia 

  4. Quaternion rotation visualization: (local post

  5. “Euler form” and “Cayley-Klein parameters” are often grouped with “Euler angle” like representations, but the later are totally ad-hoc systems and the former are not.  2

  6. Cayley transform: Wikipedia 

  7. “Sur quelques propriétés des déterminants gauches.” Journal für die reine und angewandte Mathematik 32, 1846, Arthur Cayley. Reproduced in “The collected mathematical papers of Arthur Cayley, Volume 1”, 1889. 

  8. Stereographic projection: Wikipedia 

  9. Möbius transformation: Wikipedia 

  10. Tangent half-angle substitution: Wikipedia 



Comments





© 2023 Marc B. Reynolds
all original content is public domain under UNLICENSE