Quaternion interpolation in 3D

Question

I'm a chemist lost in the captivating world of mathematics thus if you could keep your answers simple it would be awesome!

Here is my problem: I have two mobiles (A,B) in 3D. Ideally, I would like to obtain a set of quaternions discretizing the transformation from A to B. If I am right one can define an axis and an angle of rotation around this axis that defines any rotation in space from A to B.

To give an example, lets say I have a rotation of 90 deg for going from A to B around the axis e, how to I get the quaternion corresponding to 1/10th,2/10th,3/10th, ..., 1 of this transformation ?

I'm very willing to rephrase and try to explain my problem better if it's not clear enough! Just ask.

Thanks in advance for you help,

Nicolas

Answer 1

If I get the description on the wikipedia site right. First build a quaternion $${\bf q} = \underset{\text{"real" part}}{\underbrace{\cos(\theta/2)}} + (e_x{\bf i} + e_y{\bf j} + e_z{\bf k})\sin(\theta/2)$$ where vector for the axis of rotation is $[e_x,e_y,e_z]$.

Now stuff the point you want to rotate into another quaternion, like this:

$${\bf p} = 0 + (p_x {\bf i} + p_y{\bf j} + p_z{\bf k})$$

That is put x,y and z coordinates of the point into ${\bf i,j,k}$ components and let the "real part" stay 0.

Then to perform the rotation, do ${\bf qpq}^{-1}$ then just extract the coordinates (what is in front of ${\bf i,j,k}$). Please note that it is important with the order of multiplication for quaternions. You can not switch place of the factors left-to-right like you can with complex, real and integer numbers.

Here you can choose $\theta$ to be any angle, for example $90/10=9$ deg if you want.

---

If you already have a $\bf q$ for one angle and want "shorter" rotations you can use formula for multiple-angles for sine and cos or just arcsin / arccos. Extract $\cos(\theta/2)$ as the "real part" from your $\bf q$, do arc cosine, divide the angle by the factor you want and then do sine and cosine and adjust.

---

EDIT: Also This will be a particular root (can you guess which?) to the quaternion polynomial:

$${\bf q}^{10} = {\bf q}_0 \Leftrightarrow {\bf q}^{10} - {\bf q}_0 = {\bf 0}$$

If the answer to this question is affirmative, we can implement a root-finding algorithm by stuffing matrix representations of the quaternions into the Characteristic polynomial using Kronecker products in the same way as is done for the complex numbers in the question.

Answer 2

Let's say you have two unit quaternions (also known as versors), $$\mathbf{q}_0 = w_0 + x_0 \mathbf{i} + y_0 \mathbf{j} + z_0 \mathbf{k}$$ and $$\mathbf{q}_1 = w_1 + x_1 \mathbf{i} + y_1 \mathbf{j} + z_1 \mathbf{k}$$ that represent rotations. You can blend between them arbitrarily using $0 \le t \le 1$: $$\mathbf{q}' = (1-t)\mathbf{q}_0 + t\mathbf{q}_1\\ \mathbf{q} = \frac{\mathbf{q}'}{\lVert\mathbf{q}'\rVert}$$ In other words, summing (scaled) quaternions and normalizing the result to unit length produces the desired intermediate rotation quaternions.

The Wikipedia article on quaternions and spatial rotation contains a lot of useful information, although most of it is a bit too mathematical for non-mathematicians to use/need. In particular, the quaternion-derived rotation matrix is very useful.

Note that an unit quaternion also represents rotation of $\theta$ around an axis that passes through origin and point $(x,y,z)$, $$\mathbf{q} = \cos\frac{\theta}{2} + \left(\frac{x}{\sqrt{x^2 + y^2 + z^2}}\sin\frac{\theta}{2}\right)\mathbf{i} + \left(\frac{y}{\sqrt{x^2 + y^2 + z^2}}\sin\frac{\theta}{2}\right)\mathbf{j} + \left(\frac{z}{\sqrt{x^2 + y^2 + z^2}}\sin\frac{\theta}{2}\right)\mathbf{k}$$ assuming $x^2 + y^2 + z^2 \ne 0$, i.e. $(x,y,z)$ is not a null vector. However, the $(1,0,0,0) = 1$ represents no rotation at all, and should be used instead of a null quaternion -- that is, for both zero $\theta$ and null $(x,y,z)$ vector.

The above is often shown in terms of unit vector as axis, $x^2+y^2+z^2=1$, in which case the division by square root is just division by one, and therefore vanishes. The above is valid for any point $(x,y,z)$ except origin.

Because of the axis-angle representation, rotation quaternions $-\mathbf{q}$ and $\mathbf{q}$ represent the same rotation. However, if the real part is negative, the rotation is over 180°. If you want the rotations to always be through the shorter arc (of the full circle), just negate all components of the quaternion if the real part is negative.

In pseudocode, blending between two unit quaternions can be written as

function blend(q0, q1, t):

    if (t <= 0):
       return q0
    end if
    if (t >= 1):
       return q1
    end if

    q.w = (1 - t)*q0.w + t*q1.w
    q.x = (1 - t)*q0.x + t*q1.x
    q.y = (1 - t)*q0.y + t*q1.y
    q.z = (1 - t)*q0.z + t*q1.z

    nn = q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z

    if (nn > 0):
        if q.w < 0:
            n = -sqrt(nn)
        else:
            n = sqrt(nn)

        q.w = q.w / n
        q.x = q.x / n
        q.y = q.y / n
        q.z = q.z / n

    else:
        q.w = 1
        q.x = 0
        q.y = 0
        q.z = 0
    end if

    return q  

If you intend to calculate many $t$ between the same pair of rotation quaternions, it is sufficient to ensure that both have a positive real part. However, even then you still need to normalize the result to unit length. It is always possible that the vector part is zero (nn == 0 in the above pseudocode) for some specific $t$, and you should ensure you produce a unit quaternion $1$ then (one as the real component, zero for others). Otherwise, you cannot use it as a versor, or derive a rotation matrix from it, and expect sane results.

Finally, if the transition using linear $t$ feels too abrupt, you can instead use $$t' = 3t^2 - 2t^3$$ which starts and ends with first derivative zero - this corresponds to camera movement that starts and stops with zero velocity, with highest change (velocity) in the middle. An even smoother transition is obtained using $$t' = 6t^5 - 15t^4 + 10t^3$$ which starts and stops with zero velocity and zero acceleration.

You can even rotate a unit quaternion by another unit quaternion (the same as applying the rotations consecutively), by calculating their Hamilton product. See here for details. Although the result should be an unit quaternion, you can always safely normalize the result, to ensure rounding errors do not compound. (When multiplying enough rotation matrices together, the rounding errors become visible at some point. While it is possible to normalize a rotation matrix too, there are many different ways to do it, and each favors a specific axis or direction. Unit quaternions/versors do not suffer from that; normalization does not favour any 3D direction.)