python - rotating coordinate system via a quaternion

Question

We have a gazillion spatial coordinates (x, y and z) representing atoms in 3d space, and I'm constructing a function that will translate these points to a new coordinate system. Shifting the coordinates to an arbitrary origin is simple, but I can't wrap my head around the next step: 3d point rotation calculations. In other words, I'm trying to translate the points from (x, y, z) to (x', y', z'), where x', y' and z' are in terms of i', j' and k', the new axis vectors I'm making with the help of the euclid python module.

I think all I need is a euclid quaternion to do this, i.e.

>>> q * Vector3(x, y, z)
Vector3(x', y', z')

but to make THAT i believe I need a rotation axis vector and an angle of rotation. But I have no idea how to calculate these from i', j' and k'. This seems like a simple procedure to code from scratch, but I suspect something like this requires linear algebra to figure out on my own. Many thanks for a nudge in the right direction.

Answer 1

Using quaternions to represent rotation isn't difficult from an algebraic point of view. Personally, I find it hard to reason visually about quaternions, but the formulas involved in using them for rotations are not overly complicated. I'll provide a basic set of reference functions here.¹ (See also this lovely answer by hosolmaz, in which he packages these together to create a handy Quaternion class.)

You can think of a quaternion (for our purposes) as a scalar plus a 3-d vector -- abstractly, w + xi + yj + zk, here represented by an ordinary tuple (w, x, y, z). The space of 3-d rotations is represented in full by a sub-space of the quaternions, the space of unit quaternions, so you want to make sure that your quaternions are normalized. You can do so in just the way you would normalize any 4-vector (i.e. magnitude should be close to 1; if it isn't, scale down the values by the magnitude):

def normalize(v, tolerance=0.00001):
    mag2 = sum(n * n for n in v)
    if abs(mag2 - 1.0) > tolerance:
        mag = sqrt(mag2)
        v = tuple(n / mag for n in v)
    return v

Please note that for simplicity, the following functions assume that quaternion values are already normalized. In practice, you'll need to renormalize them from time to time, but the best way to deal with that will depend on the problem domain. These functions provide just the very basics, for reference purposes only.

Every rotation is represented by a unit quaternion, and concatenations of rotations correspond to multiplications of unit quaternions. The formula² for this is as follows:

def q_mult(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
    x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
    y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
    z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
    return w, x, y, z

To rotate a vector by a quaternion, you need the quaternion's conjugate too:

def q_conjugate(q):
    w, x, y, z = q
    return (w, -x, -y, -z)

Quaternion-vector multiplication is then a matter of converting your vector into a quaternion (by setting w = 0 and leaving x, y, and z the same) and then multiplying q * v * q_conjugate(q):

def qv_mult(q1, v1):
    q2 = (0.0,) + v1
    return q_mult(q_mult(q1, q2), q_conjugate(q1))[1:]
    

Finally, you need to know how to convert from axis-angle rotations to quaternions and back. This is also surprisingly straightforward. It makes sense to "sanitize" input and output here by calling normalize.

def axisangle_to_q(v, theta):
    v = normalize(v)
    x, y, z = v
    theta /= 2
    w = cos(theta)
    x = x * sin(theta)
    y = y * sin(theta)
    z = z * sin(theta)
    return w, x, y, z

And back:

def q_to_axisangle(q):
    w, v = q[0], q[1:]
    theta = acos(w) * 2.0
    return normalize(v), theta

Here's a quick usage example. A sequence of 90-degree rotations about the x, y, and z axes will return a vector on the y axis to its original position. This code performs those rotations:

x_axis_unit = (1, 0, 0)
y_axis_unit = (0, 1, 0)
z_axis_unit = (0, 0, 1)
r1 = axisangle_to_q(x_axis_unit, numpy.pi / 2)
r2 = axisangle_to_q(y_axis_unit, numpy.pi / 2)
r3 = axisangle_to_q(z_axis_unit, numpy.pi / 2)

v = qv_mult(r1, y_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (0.0, 1.0, 2.220446049250313e-16)

Keep in mind that this sequence of rotations won't return all vectors to the same position; for example, for a vector on the x axis, it will correspond to a 90 degree rotation about the y axis. (Think of the right-hand-rule here; a positive rotation about the y axis pushes a vector on the x axis into the negative z region.)

v = qv_mult(r1, x_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (4.930380657631324e-32, 2.220446049250313e-16, -1.0)

As always, please let me know if you find any problems here.

---

^{1. These are adapted from an OpenGL tutorial archived here.}

^{2. The quaternion multiplication formula looks like a horrible rat's nest at first, but the derivation is easy, albeit tedious. Working with pencil and paper, you can represent the two quaternions like so: w + xi + yj + zk. Then note that ii = jj = kk = -1; that ij = k, jk = i, ki = j; and that ji = -k, kj = -i, ik = -j. Finally, multiply the two quaternions, distributing out the terms and rearranging them based on the results of each of the 16 multiplications. This helps to illustrate why you can use quaternions to represent rotation; the last six identities follow the right-hand rule, creating bijections between rotations from i to j and rotations around k, and so on.}

^{If you do this, you'll see that the identity ii = jj = kk = -1 explains the last three terms in the w formula; that jk = i explains the third term in the x formula; and that kj = -i explains the fourth term in the x formula. The y and z formulas work the same way. The remaining terms are all examples of ordinary scalar multiplication.}