April 19th, 2017 This is a massively chopped down version of a post I haven’t been completing, so maybe breaking into parts will work better. The goals are as follows: Derive slerp, attempt to make it less black-box like and show it is a 2D problem Show that it is equivalent to the commonly seen function Form a reference version that only requires one forward trig operation Present analysis tools which require only reals Give basic analysis for lerp The following shadertoy is a viz for three interpolations formulated directly in 2D: slerp (blue) nlerp (magenta) half-angle+nlerp (red). Clicking down sets the initial point and the final point is the current mouse position. On button release the initial position resets to zero and the final remains. As a 2D slice of quaternions the maximum angle between the endpoints is $\frac{\pi}{2}$. From complex real powers to the algebraic form of slerp Since quaternions are an extension of complex numbers let’s first consider the complex case. Given a unit complex number $z$ both in explict and polar form: \[\begin{eqnarray} z & = & a + b\mathbf{i} \label{zr} \\ & = & \cos\left(\theta\right)+\sin\left(\theta\right)~\mathbf{i} \label{zp} \end{eqnarray}\] The real (principle) power can be expressed in terms of the polar form as $\left( t\in \mathbb{R} \right) $: \[\begin{equation} \label{power} z^t = \cos\left(t\theta\right)+\sin\left(t\theta\right)~\mathbf{i} \end{equation}\] This linearly parameterizes the angle between one and $z$. If $z$ respresents a rotation then this produces a continuous interpolation between the identity and the rotation represented by $z$ following the minimal magnitude path. To generalize to an arbitrary starting point $A$ and ending point $B$ we can simply transform the problem: \[\begin{equation} \label{zslerp} s \left( t \right) = \left(BA^*\right)^tA \end{equation}\] Which is the algebraic form of slerp in two-dimensions. SLERP the dreaded small angle circular arc Replacing $z$ with a unit quaternion: \[\begin{eqnarray} Q & = & a + \left(x,~y,~z\right) \label{qxform} \\ & = & a + b~\mathbf{u} \label{qiform} \\ & = & \cos\left(\theta\right)+\sin\left(\theta\right)~\mathbf{u} \label{qpform} \end{eqnarray}\] changes hardly anything. I’m assuming familarity with the contents of half/double angle and Cayley: preliminaries, the bullet points of which are: 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. Since the quaternion product doesn’t commute we can express slerp as: \[\begin{eqnarray} s\left( t \right) & = & \left(BA^*\right)^tA \label{slerp} \\ & = & A\left(A^*B\right)^t \label{slerp1} \\ & = & \left(AB^*\right)^{1-t}B \label{slerp2} \\ & = & B\left(B^*A\right)^{1-t} \label{slerp3} \end{eqnarray}\] Where $\eqref{slerp2}$ and $\eqref{slerp3}$ are simply reversing the points and parameterization (holds in complex as well). Using $\eqref{slerp}$ we have $Q=BA^*$ where $A \cdot B \geq 0$ and the action of slerp is $Q^t$, which maps a coordinate in the {1,$Q$}-plane to a coordinate in the {1,$Q$}-plane. Basic analysis tools The action of slerp is defined by the principle real power $\eqref{power}$ and is bound to a complex plane. So we can reason about it as walking a path1 from $\left(1,0\right) $ to $\left(a,b\right) = \left(\cos \theta,\sin \theta\right)$. If the path in the plane is described by: \[\begin{equation} \label{eq:fpoint} p\left(t\right) = \left(p_x\!\left(t\right),~p_y\!\left(t\right)\right) \end{equation}\] Then the angle at $t$ is: \[\begin{equation} \label{eq:fangle} \phi\left(t\right) =\text{atan2}\left(p_y\!\left(t\right), ~p_x\!\left(t\right)\right) \end{equation}\] From $\phi\left(t\right)$ we can compute other quantites such as angular velocity $\left(\omega=\frac{d\phi}{dt}\right)$ and acceleration $\left(\alpha=\frac{d\omega}{dt}\right)$. Taking the values for slerp $\eqref{power}$ and pluging them in we get the expected: \[\begin{eqnarray*} p\left(t\right) & = & \left( \cos\left(t\theta\right), ~\sin\left(t\theta\right)\right) \\ \phi\left(t\right) & = & \text{atan2}\left( \sin\left(t\theta\right) , ~\cos\left(t\theta\right)\right) = t \theta \\ \omega\left(t\right) & = & \theta \\ \alpha\left(t\right) & = & 0 \end{eqnarray*}\] This gives a simple framework for computing errors (excluding those from floating-point computation). Given some approximation of slerp we can compute the function $T$ to correct to constant angular velocity: \[\begin{eqnarray} \text{atan}\left(\frac{p_y\left(T(t)\right)}{p_x\left(T(t)\right)}\right) = t \theta \nonumber \\ \frac{p_y\left(T(t)\right)}{p_x\left(T(t)\right)} = \tan\left(t \theta\right) \label{reparam} \\ \end{eqnarray}\] by solving for $T$. Algebraic to geometric form the version you almost always see Rather than grind the math let’s directly place our 2D coordinate $\eqref{power}$ back into the original space: \[\begin{equation} \label{slerpxy} s(t) = \cos\left(t\theta\right)X + \sin\left(t\theta\right)Y \end{equation}\] Note that this holds for any number of dimensions. Using $\eqref{slerp}$ or $\eqref{slerp1}$ then the ${X}$ direction is $A$ and ${Y}$ is the direction orthogonal to $A$ that contains $B$ (normalized). Since the inputs are normalized we have: \[Y = \left\Vert B - (A \cdot B)A\right\Vert\] which brings us to: \[\begin{equation} \label{srecon} s(t) = \cos\left(t\theta\right)A + \sin\left(t\theta\right)\frac{B - \left(A\cdot B\right)A}{\sqrt{1 -\left(A\cdot B\right)^{2}}} \end{equation}\] If we were to swap $A$ and $B$ and replace $t$ by $1-t$ to instead use $\eqref{slerp2}$ or $\eqref{slerp3}$. Making some trig replacements: \[\begin{eqnarray*} \cos\left(\theta\right) & = & A\cdot B \\ \sin\left(\theta\right) & = & \sqrt{1 -\left(A\cdot B\right)^{2}} \end{eqnarray*}\] then we have: \[\begin{eqnarray} s(t) & = & \cos(t\theta) A + \sin(t\theta) \frac{B - A \cos(\theta)}{\sin(\theta)} \nonumber \\ & = & \left(\cos(t\theta) - \frac{\cos(\theta)\sin(t\theta)}{sin(\theta)}\right)A + \frac{\sin(t\theta)}{\sin(\theta)}B \nonumber \\ & = & \frac{\cos(t\theta)\sin(\theta) - \cos(\theta)\sin(t\theta)}{\sin(\theta)}A + \frac{\sin(t\theta)}{\sin(\theta)}B \nonumber \\ & = & \frac{\sin\left(\left(1-t\right)\theta\right)}{\sin\left(\theta\right)}A+\frac{\sin\left(t\theta\right)}{\sin\left(\theta\right)}B \label{gdavis} \end{eqnarray}\] Where $\eqref{gdavis}$ can be directly formulated from geometric reasoning and has somehow has be taken to be the true defination of slerp. Ken Shoemake’s 1985 paper2 gives $\eqref{slerp1}$ and $\eqref{gdavis}$. This equation has be covered very well, as an example J.M.P van Waveren3 2005. Reference formulation of slerp one forward trig op is enough Starting from $\eqref{srecon}$, replace $A\cdot B$ by $d$ and collecting terms gives: \[\begin{eqnarray*} s(t) & = & \cos\left(t\theta\right)A + \sin\left(t\theta\right)\frac{B - dA}{\sqrt{1 -d^2}} \\ & = & \left( \cos\left(t\theta\right)-d\frac{\sin\left(t\theta\right)}{\sqrt{1 -d^2}} \right) A + \frac{\sin\left(t\theta\right)}{\sqrt{1 -d^2}}B \end{eqnarray*}\] If we examine the terms in the limit: \[\begin{eqnarray*} \lim_{\theta \to 0} \ \ \cos\left(t\theta\right) & = 1 \\ \lim_{\theta \to 0} \ \ \frac{\sin\left(t\theta\right)}{\sin\left(\theta\right)} & = t \end{eqnarray*}\] which yields linear interpolation as $\theta$ approaches zero ($d$ approaching one): \[\begin{eqnarray*} \left(1-t \right) A + tB \end{eqnarray*}\] Or more simply the coord and arc approach being the same for small angles. The two forward trig ops can be converted into one and given the range of both are positive4. Showing both options as a pair of scale values where $s_t=\sin\left(t\theta\right)$ and $c_t=\cos\left(t\theta\right)$: \[\begin{eqnarray*} s(t) & = & \left\{\sqrt{1-s_t^2}-d \frac{s_t}{\sqrt{1-d^2}},~\frac{s_t}{\sqrt{1-d^2}}\right\} \\ & = & \left\{c_t-d \sqrt{\frac{1-c_t^2}{1-d^2}},~\sqrt{\frac{1-c_t^2}{1-d^2}}\right\} \end{eqnarray*}\] This leaves the computation of theta. I’ll paritially dismiss (directly) using $\text{acos}$ and $\text{asin}$ as being nearly hopeless functions5 and given the range of inputs we don’t need to yet distinguish between one and two parameter $\text{atan}$. There are a fair number of ways we can express $\theta$ via $\text{atan}$ and here are four of them: \[\begin{eqnarray} \theta & = & \text{atan}\left( \frac{\sqrt{1-d^2}}{d} \right) \label{acos} \\ & = & \frac{\pi}{2} - \text{ atan}\left( \frac{d}{\sqrt{1-d^2}} \right) \label{asin} \\ & = & 2 \text{ atan}\left( \frac{1-d}{\sqrt{1-d^2}} \right) \label{hacos} \\ & = & 2 \text{ atan}\left( \frac{\sqrt{1-d^2}}{1+d} \right) \label{hacos2} \end{eqnarray}\] where $\eqref{acos}$ and $\eqref{asin}$ are simply $\text{acos}$ and $\text{asin}$ in $\text{atan}$ form and $\eqref{hacos}$ $\eqref{hacos2}$ are applying half-angle to $\eqref{acos}$ in two forms. Now for some practical considerations. By reference we can go two ways: close enough to detect problems in actually used methods or minimal error with respect to representable so we can accurately measure errors. I’m going with the former. Also we have to decide the meaning of a negative dot product. It is an error so assert or spew a message. Can I haz minimal anyway plz? Go the long way around. I like memes as much as the next person so will go with option 2. Using C standard functions we could make a concise6 implemention: void slerp_ref_0(quat_t* R, quat_t* A, quat_t* B, float t) { float d = quat_dot(A,B); // cosine of angle between A & B = scalar part of BA^* float sd = sgn(d); // need sign of 'd' for choosing minimal path float s0,s1; d = fabs(d); // deals with degenerate case if (d < SLERP_CUT) { // slerp weights float s = sqrtf((1.f+d)*(1.f-d)); // sine of relative angle float a = atan2f(s,d); // acos->atan2. some number of divides here :( float c = cosf(t*a); // ta on [0, pi/2], range reduction :( s1 = sqrtf((1.f+c)*(1.f-c))/s; // atan2 might be computing 1/s :( s0 = c-d*s1; } else { // lerp s0 = 1.0f - t; s1 = t; } // weighted sum: R = s0*A + sd*s1*B quat_wsum(R, A, B, s0, sd*s1); } This leaves a fair number of thing we can know to the compiler to deduce and hopefully do something about. Let’s convert the “slerp weights” to single parameter $\text{atan}$ with input on $\left[ 0,1 \right]$ using $\eqref{hacos}$: // slerp weights t += t; // account for the 2 term float s2 = (1.f+d)*(1.f-d); float rs = rsqrt_nr(s2); // ~x^-.5 to at least 12-bit, followed by 1 NR step float y = 1-d; float a = atanf(y*rs); // still range reduction float s = sinf(t*a); // still range reduction float c = sqrtf((1.f+s)*(1.f-s)); // or cosf(ta). can't go through ~1/sqrt if max angle allowed s1 = s*rs; // sd*s1*B can be computed while c in progress s0 = c-d*s1; Showing the use of approximation plus one step of newton-raphson for $x^{-\frac{1}{2}}$ since it has minimal impact on accuracy. Good idea or not is a different story. Likewise for removing either of the forward trig ops. Assumming no hardware sin/cos then we could replace with approximations on $\left[0,\frac{\pi}{2}\right]$. Another option is move the angle doubling to after the approximation(s): // set-up for atan float s2 = (1.f+d)*(1.f-d); float rs = rsqrt_nr(s2); float y = 1-d; float a = atan_x(y*rs); // ~atan on [0,1] // forward trig approx (or root for one) float ta = t*a; // above we doubled t float c = cos_x(ta); // ~cos on [0, pi/4] float s = sin_x(ta); // ~sin on [0, pi/4] // double angle float tc = c+c; float cx = tc*c-1.f; // double angle (c^2-s^2 = 2c^2-1 = 1-2s^2) float sx = tc*s; // double angle s1 = sx*rs; s0 = cx-d*s1; Although the range of arctan is narrow enough to directly approximate it requires a fair number of terms to achieve small error. From the sum-of-tangents identity: \[\text{atan}\left(a\right) + \text{atan}\left(b\right) = \text{atan}\left(\frac{a+b}{1-ab}\right)\] To halve our angle range set $b=-1$ since $\text{atan}\left(-1\right)=-\frac{\pi}{4}$. Using this substitution allows narrowing the range to $\left[0, \sqrt{2}-1\right]$ Totally ignoring all the range reduction comments and toss in some lowest order static inline float rsqrtf_a(float x) { return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x))); } static inline float rsqrt_nr(float x) { float x0 = rsqrtf_a(x); return x0*(1.5f - x*0.5f*x0*x0); } static inline float p45(float x) { float x2=x*x; return x * (-2793121 * 0x1p-23f + x2 * (+13098501 * 0x1p-26f + x2 * (-4020455 * 0x1p-25f + x2 * ( 14980357 * 0x1p-28f + x2 * (-13750363 * 0x1p-30f))))); } static inline float atan_7(float x) { // 1.66290684247846942061144487572104929234984377645523e-6 float x2=x*x; return x * ( 0.99997723102569580078125f + x2 * (-0.332623004913330078125f + x2 * ( 0.19354121387004852294921875f + x2 * (-0.116428203880786895751953125f + x2 * ( 5.2648954093456268310546875e-2f + x2 * (-1.1719689704477787017822265625e-2f)))))); } // error = 6.77119023748673498630523681640625e-5 // zero = [|0.678780257701873779296875, 1.22571027278900146484375, 1.531032562255859375|] static inline float sin_2_4(float x) { float x2 = x*x; float s; s = x2*( +7.5143859721720218658447265625e-3f); s = x2*(s-0.1656731069087982177734375f); s = x *(s+0.999696791172027587890625f); return s; } #define SLERP_CUT 0x1.ffffdep-1 // max angle error ~= 0.097 degrees in 3D void slerp_sloppy_ex(quat_t* R, quat_t* A, quat_t* B, float t) { float d = quat_dot(A,B); float sd = sgn(d); float s0,s1; d = fabsf(d); if (d < SLERP_CUT) { t += t; float s2 = (1.f+d)*(1.f-d); float rs = rsqrt_nr(s2); float y = 1-d; float a = p45(y*rs); float ta = t*a; float s = sin_2_4(ta); float c = sqrtf((1.f+s)*(1.f-s)); s1 = s*rs; s0 = c-d*s1; } else { s0 = 1.0f - t; s1 = t; } quat_wsum(R, A, B, s0, sd*s1); } The reason for walking through some possible reductions of “reference form” are lost since I’ve ripped out the almost all the sections, but felt it was better to leave in place. Also losing forward trig operations is more a CPU side issues than GPU. LERP and nlerp Going back to our reference formula written with trig operations: \[\begin{eqnarray*} \left( \cos\left(t\theta\right)-\cos\left(\theta\right)\frac{\sin\left(t\theta\right)}{\sin\left(\theta\right)} \right) A + \frac{\sin\left(t\theta\right)}{\sin\left(\theta\right)}B \end{eqnarray*}\] Approximating with either lerp or nlerp produce the same angular error values. Errors introduced by using a non-unit quaternion depend on how the consumer of the output is written (the simplifications made using unit magnitude as a given). Plugging in the approximation we get: \[\begin{eqnarray*} p_l\left(t\right) & = & \left( 1+t\left(\cos\left(\theta\right)-1\right), ~t\sin\left(\theta\right)\right) \\ \phi_l\left(t\right) & = & \text{atan}\left(~t\sin\left(\theta\right), ~1+t\left(\cos\left(\theta\right)-1\right)\right) \\ \omega_l\left(t\right) & = & \frac{\sin\left(\theta\right)}{1-2t\left(\cos\left(\theta\right) \left(t-1\right) - \left(t-1\right) \right)} \\ \alpha_l\left(t\right) & = & \frac{2\sin\left(\theta\right)\left(\cos\left(\theta\right)-1\right)\left(2t-1\right)} {\left(1+2t\left(\cos\left(\theta\right) \left(t-1\right) - \left(t-1\right) \right)\right)^2} \end{eqnarray*}\] To get a quick feel let’s examine the angle and angular velocity errors for some sample angles: $ \frac{\pi}{2}, \frac{\pi}{4}, \frac{\pi}{8} $ (Rotations of: 180⁰, 90⁰, 45⁰). First the absolute error of the angle at ‘t’ (normalized to range X axis): \[\frac{1}{\pi} \left( t\theta - \phi_l \left(t \right) \right)\] As expected the result is exact at the end points and at $t=\frac{1}{2}$. The latter since a ray that bisects a coord bisects the arc. We can see the error decreases fast and by progressively turning off traces (click the legend) eye-ball it around a factor of 10 each time we halve the angle. Also the ‘t’ of the peak error shifts to smaller values and appears to be slowing toward some final point. Now the angular velocity absolute error. I’ve flipped the ordering so positive is too fast and negative is too slow: \[\omega_l \left(t \right) - \theta\] We can form an expression for the $t$ value of the maximum angle error as a function of $\theta$ by solving for $\omega_l\left(t\right)=\theta$: \[t_{\text{max}}\left(\theta\right) = \frac{1}{2} \left(1 - \sqrt{\frac{1 + \cos \theta - \frac{2}{\theta} \sin \theta}{ \cos \theta -1 }} \right)\] which is a very boring function. Here are some sample values: \[\begin{align*} t_{\text{max}}\left(\frac{\pi}{2} \right) & = \frac{1}{2}\left( 1-\sqrt{\frac{4}{\pi}-1} \right) & = \ \sim 0.238638 \\ t_{\text{max}}\left(\frac{\pi}{4} \right) & = & = \ \sim 0.217459 \\ t_{\text{max}}\left(\frac{\pi}{16} \right) & = & = \ \sim 0.211697 \\ t_{\text{max}}\left(0\right) & = \frac{1}{6}\left( 3-\sqrt{3} \right) & = \ \sim 0.211325 \end{align*}\] So the $t_{max}$ is rapidly approaching the fixed point $t_{max}(0)$. The max error for angular velocity is always at $t=\frac{1}{2}$. Choosing to examine the maximal angle range we have worst case error values of: \[\begin{align*} \phi_{max} & = \frac{\pi}{2} t_{max}\left( \frac{\pi}{2} \right) - \phi_l \left( t_{max}\left(\frac{\pi}{2} \right) \right) & = \ \sim 0.0711146 \\ \omega_{max} & = \omega_l \left(\frac{1}{2}\right) - \frac{\pi}{2} & = \ \sim 0.429204 \end{align*}\] Normalizing the abs error by these and varying theta (also normalized) we can examine max error vs angle: And the same plot as a log/log plot: Since these are approximately linear in log-log implies the error terms are approximately monomials. References and Footnotes Actually any scalar multiple of the end-points is mathematically correct. ↩ “Animating Rotation with Quaternion Curves”, Ken Shoemake, 1985. (PDF) ↩ “Slerping Clock Cycles”, van Waveren, 2005. ↩ no sign fix-up needed. ↩ arcos can be approximated by: $\sqrt{1-x}P\left(x\right)$ ↩ The abs and sgn computations aren’t needed with atan2. ↩ Comments quaternions (15) interpolation (5)