https://fgiesen.wordpress.com/2012/08/24/quaternion-differentiation/ Skip to content Follow: RSS Twitter The ryg blog When I grow up I'll be an inventor. * Home * About * Coding * Compression * Computer Architecture * Demoscene * Graphics Pipeline * Maths * Multimedia * Networking * Papers * Stories * Thoughts * Uncategorized Quaternion differentiation August 24, 2012 I wanted to point someone to a short explanation of this today and noticed, with some surprise, that I couldn't find something concise within 5 minutes of googling. So it seems worth writing up. I'm assuming you know what quaternions are and what they're used for. First, though, it seems important to point one thing out: Actually, there's nothing special at all about either integration or differentiation of quaternion-valued functions. If you have a quaternion-valued function of one variable q(t), then \displaystyle \dot{q}(t) = q'(t) = \frac{\mathrm{d}q}{\mathrm{d}t} = \lim_{h \rightarrow 0} \frac{q(t+h) - q(t)}{h} same as for any real- or complex-valued function. So what, then, is this post about? Simple: unit quaternions are commonly used to represent rotations (or orientation of rigid bodies), and rigid-body dynamics require integration of orientation over time. Almost all sources I could find just magically pull a solution out of thin air, and those that give a derivation tend to make it way more complicated than necessary. So let's just do this from first principles. I'm assuming you know that multiplying two unit quaternions quaternions q[1]q[0] gives a unit quaternion representing the composition of the two rotations. Now say we want to describe the orientation q(t) of a rigid body rotating at constant angular velocity. Then we can write q(0) = q_0 q(1) = q_\omega q_0 where q_\omega describes the rotation the body undergoes in one time step. Since we have constant angular velocity, we will have q(2) = q_ \omega q_\omega q_0 = q_\omega^2 q_0, and more generally q(k) = q_\ omega^k q_0 for all nonnegative integer k by induction. So for even more general t we'd expect something like q(t) = q_\omega^t q_0. Now, q[o] is a unit quaternion, which means it can be written in polar form q_\omega = \cos(\theta/2) + \sin(\theta/2) (\mathbf{n}_x i + \mathbf {n}_y j + \mathbf{n}_z k) where th is some angle and n is a unit vector denoting the axis of rotation. That part is usually mentioned in every quaternion tutorial. Embedding real 3-vectors as the corresponding pure imaginary quaternion, i.e. writing just \mathbf{n} for the quaternion \mathbf{n}_x i + \mathbf{n}_y j + \mathbf{n}_z k, is usually also mentioned somewhere. What usually isn't mentioned is the crucial piece of information that the polar form of a quaternion, in fact, just the quaternion version of Euler's formula: Any unit complex number z can be written as the complex exponential of a pure imaginary number z=\exp(it)=\cos(t) + i \sin(t), and similarly any unit-length quaternion (and q[o] in particular) can be written as the exponential of a pure imaginary quaternion q_\omega = \exp(\frac{\theta}{2}\mathbf{n}) which gives us a natural definition for q_\omega^t = \exp(t \frac{\theta}{2} \mathbf{n}). Now, what if we want to write a differential equation for the behavior of q(t) over time? Just compute the derivative of q(t) as you would for any other function of t. Using the chain and product rules we get: \dot{q}(t) = \frac{\mathrm{d}}{\mathrm{d}t} (q_\omega^t q_0) = \frac {\theta}{2} \exp(t \frac{\theta}{2} \mathbf{n}) \mathbf{n} q_0 Now t and th are real numbers, so the exponential is of a real number times n; exp is defined as a power series, and for arbitrary c real we have \displaystyle \exp(c \mathbf{n}) = \sum_{k=0}^{\infty} \frac{(c \ mathbf{n})^k}{k!} = \sum_{k=0}^{\infty} \frac{c^k \mathbf{n}^k}{k!} n commutes with powers of n, and real numbers commute with everything; therefore, n commutes with exp(cn) for arbitrary c, and thus \frac{\theta}{2} \exp(t \frac{\theta}{2} \mathbf{n}) \mathbf{n} q_0 = \frac{\theta}{2} \mathbf{n} \exp(t \frac{\theta}{2} \mathbf{n}) q_0 = \frac{\theta}{2} \mathbf{n} q(t) The vector thn is in fact just the angular velocity o, which yields the oft-cited but seldom-derived equation: \dot{q} = \frac{\mathrm{d}q}{\mathrm{d}t} = \frac{1}{2} \omega q This is usually quoted completely without context. In particular, it's not usually mentioned that q(t) describes the orientation of a body with constant angular velocity, and similar for the crucial link to the exponential function. Share this: * Facebook * X * Like Loading... Related From - Maths 20 Comments 1. [bc8b] Goran Andersson permalink Hi, It was with great satisfaction I found your deduction of the quarternion time derivate. I've been searching for this myself, but until now in vain. However, I don't understand when you say: "The vector thn is in fact just the angular velocity o,.." The angular velocity o have the units radians/s, but the quantity thn have the unit radians as n is dimensionless. Or have I missed something... Regards, Goran Reply + [3a11] fgiesen permalink In my derivation, t is just a real number; I'm not bothering with units throughout the post, which is what muddles things a bit. If you want to patch it up and use dimensional quantities, you need to do it throughout: The initial definition q(t) = q_omega^t * q_0 doesn't work as-is; you need a fudge factor k [1/s] so you can use a power function here (we need to nail down after what amount of time we've turned by "one omega"). So you get q(t) = q_omega^(k*t) * q_0 and for the differential equation, that gives us q'(t) = k*(theta/2)*n * exp(k*t*(theta/2)*n) * q_0 = k*(theta /2)*n * q(t) = 1/2 * omega * q(t) where omega = k*theta*n. k is still [1/s], theta is [rad] and n is a dimensionless vector. Hope that helps. Sorry for the confusion. Reply o [b022] Hongwei permalink Hi, your explanation about this question regarding unit is already very clear. Thanks a lot for that. But I'm still confused about the relationship. From my understanding, omega should be how much theta per second. Then it should be: omega = k*theta Why omega = k*theta*n? I think it doesn't make sense that n is representing the direction of angular rate, because it should be 'orthogonal' between n and angular rate, shouldn't it? Otherwise, what does n represent? Hope to get your answer soon! Regards, Hongwei. o [3a11] fgiesen permalink n is the normalized axis. theta is [radians], k is [1/s]. k*theta is radians per second but doesn't specify the axis of rotation. 2. [305c] Daniel Smith permalink Hello! When checking the derivation of partial derivatives of a function of a quaternion for another mathematician at work, I found a fairly obscure flaw, this is not mentioned in your article above. Since the complex part is a unit vector, the obvious differentiation does not work, because it effectivily wiggles the individual parts without considering that wiggling one part would change all the other parts to keep the sum 1. Mathematically speaking, the individuial terms are not independent. To get the right result you must differentiate the imaginary parts of quaternion divided by the length of the unit vector i.e. R + Im/sqrt(x^2 + y^2 + z^2), this breaks the dependency between the terms. Reply + [3a11] fgiesen permalink I'm not sure what you mean. The imaginary part of the quaternions involved is *not* a unit vector. A unit quaternion is one such that conj(q)*q = R^2 + x^2 + y^ 2 + z^2 (in your notation) = 1. I don't see why you would normalize just the imaginary part of a quaternion, or what it's intended to accomplish. I'm simply talking about the quaternion-valued function q(t) = q_omega^t * q_0 here, which has the derivative (by time!) given above. Furthermore, for unit quaternions (which represent rotations), we have |q(t)| = |q_omega^t| * |q_0| = |q_omega|^t * |q_0| = 1^t * 1 = 1; so if q_omega and q_0 are both unit quaternions, the same will hold for any quaternion returned by q(t), without ever enforcing that constraint explicitly. "Wiggling the individual parts" doesn't make sense in this context; both q_omega and q_0 are constants, and while I'm looking at a quaternion-valued function, it's a function of one real variable, t, which is the one I'm differentiating by. I don't get to "wiggle" my quaternions off the curve described by q(t), which as I've described above is entirely made up of unit quaternions. Reply 3. [pict] Matyas Czeman permalink Hello! You are my last hope. As u mentioned q1=qw*q0 and so on. So my question would be, what is qw? So to make it clear: I got an initail quaternion, and the angular velocities at 50Hz or so (doesnt matter) How can I update my orientation quaternion based on the w.x, w.y, and w.z (the 3 angulary velocity)? Reply + [3a11] fgiesen permalink q_w in polar form is given in the article, this is assuming you have a normalized axis (n) and an angle in radians per time step (theta). The angular velocity vector w (well, omega really...) is the axis vector times the angular speed in radians per second. If w=0, q_w is just the identity quaternion. Otherwise, use theta=len(w)*secs_per_timestep, n=normalize(w). Reply o [fa65] darus permalink I am a bit confused. If talking about linear velocity, then we can write: dp/dt = v where p - position vector, v - linear velocity vector. Then we can write: p(t + Dt) = p(t) + v*Dt One would expect the same principle work for quaternions. So, if we have dq/dt = 0.5*o*q then we can write: q(t + Dt) = q(t) + dq/dt * Dt = q(t) + 0.5*o*q(t)*Dt Can this scheme be used to integrate position in quaternion over time? o [3a11] fgiesen permalink First-order, yes, and that is indeed the most common use for this. However, position is linear; rotation fundamentally isn't. If the forces acting on a particle sum to 0, "p(t + Dt) = p(t) + v*Dt" is exact. The same is not true for rotation. Even at constant angular velocity, there are higher-order terms in the Taylor series q(t + Dt) = q(t) + dq/dt * Dt + d^2q/dt^2 * Dt^2 + ... and they are (in general) nonzero. The exact expression for constant angular velocity is already given in the article as part of the derivation: q(t + Dt) = q_omega^(Dt) * q(t) For comments on how to modify this when not working with unitless quantities, see my earlier comment. 4. [c1d1] Nick permalink This is exactly the missing piece I needed. I've been trying to get from angular velocity to quaternions and keep getting lost in a sea of trigonometry. Thanks! Reply 5. [c521] Revathi Ravula permalink Hi, The derivation was very helpful. This is when we have a function for a quaternion but what about differenciation when we have numbers instead of functions. Basically numerical differenciation of quaternions. I was wondering if you can suggest some ways to do it. Thanks Reply + [3a11] fgiesen permalink You do it exactly the same way you would differentiate real or complex numbers: (q(t + h) - q(t)) / h. Reply 6. [c521] revathi permalink thank you very much. it was very helpful :) Reply 7. [3663] Michael permalink Hello! By searching one cannot often find any practical formulae about quaternion space derivatives I'd like to point such formulae. The best and simplest way to compute a quaternionic derivative is as follows. We represent a quaternion argument p = x+y[?]i+z[?]j+u[?]k (i,j,k are basic quaternion units; "[?]" is the quaternion multiplication) and a quaternion-differentiable (holomorphic) function of that argument ps(p) = ps_1(x,y,z,u)+ps_2(x,y,z,u)[?]i+ps_3 (x,y,z,u)[?]j+ps_4(x,y,z,u)[?]k in the Cayley-Dickson doubling form: p = a+b[?]j, where a'= x+y[?]i ; b = z+u[?]i and ps(p)=ps(a,b)=Ph_1(a,b)+Ph_ (a,b)[?]j, where Ph_1(a,b)=ps_1(a,b)+ps_2(a,b)[?]i and Ph_2(a,b)=ps_3(a,b) +ps_4(a,b)[?]i. Each expression for ps(p) is initially to be obtained from a complex function of the same kind by means of the direct replacement of a complex variable with a quaternion variable in the expression for the complex function. For example, ps(p)=p^ (-1). Just as a complex- holomorphic function satisfies Cauchy-Riemann's equations in complex analysis, a quaternion- holomorphic function satisfies the following quaternionic generalization of Cauchy-Riemann's equations: (1) [?]Ph_1/[?]a = [?]Ph_2*/[?]b*, (2) [?]Ph_2/[?]a = - [?]Ph_1*/[?]b* (3) [?]Ph_1/[?]a = [?]Ph_2/[?]b, (4) [?]Ph_2/[?]a* = - [?]Ph_1/[?]b* after doing a = a* = x, where the complex conjugation is denoted by * and the partial differentiation with respect to some variable s by [?]/[?]s. For example, by [?]Ph_2*/[?]b* is denoted the partial derivative of the complex conjugate of a function Ph_2 with respect to the complex conjugate of a variable b. Firstly, we compute the partial derivatives of the functions Ph_1, Ph_2, Ph_1*, Ph_2* (with respect to the variables a, b, a*, b*); secondly, we put a = a* =x in the computed expressions of partial derivatives; and thirdly, we check whether equations (1) - (4) hold. One of the formulae to compute the first quaternionic derivative of the quaternion-holomorphic function is the following: ps(p)[1] = ([?]Ph_1/[?]a + [?]Ph_1/[?]a*) + ([?]Ph_2/[?]a + [?]Ph_2/[?]a*)[?]j . Higher derivatives of quaternion-holomorphic functions can be computed analogically and they are holomorphic like the first derivative. For details and examples I refer to http://vixra.org/abs/ 1609.0006 Reply 8. [6b6a] Lucas Kunc permalink I really like this explanation, although there's one thing I don't fully understand. When calculating the derivative, how do you get (theta/2) * n to the left of the exponential? When I apply the chain rule it ends up to the right [if f = exp and g(x) = x * n * theta/2 then (f.g)'(t) = exp(t * n * theta/2) * n * theta/2], and because of non-commutativity of quaternions I don't see how I can get to the desired result. Am I doing something wrong? Reply + [3a11] fgiesen permalink No, you're right, that derivation incorrectly assumes commutativity as written and doesn't actually work. :) Reply o [3a11] fgiesen permalink I'm about to fix this up in the post, but the short version is this: n commutes with powers n^k of n, and then you can use the definition of exp as a power series in n to show that n also commutes with exp(cn) with c= arbitrary real. 9. [e239] Kristians permalink Hi, Why you are not computing derivatives of n and th as well? They do change in general rotation. Reply + [3a11] fgiesen permalink Because the problem statement assumes they're constant, and I'm trying to cover the basic problems without introducing extraneous details! This is still differential calculus, once you know the derivatives of the elementary functions you're using, the rest is a matter of using the rules. If you need varying n or th, you can use the chain and product rules yourself. Reply Leave a comment Cancel reply [ ] [ ] [ ] [ ] [ ] [ ] [ ] D[ ] << Z-transform done quick Linear interpolation past, present and future >> * Recent Posts + Entropy decoding in Oodle Data: x86-64 6-stream Huffman decoders + Computational complexity of texture encoding + A very brief BitKnit retrospective + Notes on FFTs: for implementers + Notes on FFTs: for users + What's that magic computation in stb__RefineBlock? + On AlphaTensor's new matrix multiplication algorithms + Morton codes addendum + Entropy decoding in Oodle Data: x86-64 3-stream Huffman decoders + Entropy decoding in Oodle Data: Huffman decoding on the Jaguar * Categories + Coding + Compression + Computer Architecture + Demoscene + Graphics Pipeline + Maths + Multimedia + Networking + Papers + Stories + Thoughts + Uncategorized * Archives + October 2023 + July 2023 + May 2023 + March 2023 + November 2022 + October 2022 + September 2022 + April 2022 + October 2021 + August 2021 + July 2021 + July 2019 + April 2019 + February 2019 + December 2018 + September 2018 + March 2018 + February 2018 + January 2018 + December 2017 + November 2017 + September 2017 + August 2017 + April 2017 + October 2016 + August 2016 + April 2016 + March 2016 + February 2016 + January 2016 + December 2015 + October 2015 + September 2015 + July 2015 + May 2015 + February 2015 + December 2014 + October 2014 + August 2014 + July 2014 + June 2014 + May 2014 + March 2014 + February 2014 + December 2013 + November 2013 + October 2013 + September 2013 + August 2013 + July 2013 + June 2013 + May 2013 + March 2013 + February 2013 + January 2013 + August 2012 + July 2012 + June 2012 + April 2012 + March 2012 + February 2012 + November 2011 + October 2011 + September 2011 + August 2011 + July 2011 + May 2011 + February 2011 + January 2011 + December 2010 + November 2010 + October 2010 + September 2010 + August 2010 + March 2010 + December 2009 + October 2009 Blog at WordPress.com. * Comment * Subscribe Subscribed + [wpcom-] The ryg blog Join 476 other subscribers [ ] Sign me up + Already have a WordPress.com account? Log in now. * + [wpcom-] The ryg blog + Customize + Subscribe Subscribed + Sign up + Log in + Copy shortlink + Report this content + View post in Reader + Manage subscriptions + Collapse this bar %d [b]