Quaternions#

Rotations are important in defining a scanning strategy. Here, we present a short tutorial on quaternions and explain the facilities provided by LBS.

There are several choices to describe a rotation in 3D space: Euler angles, rotation matrices, quaternions, etc. Each of these systems has its share of advantages and disadvantages. For instance, rotation matrices are handy when you have a vector and want to rotate it, as it is just a matter of doing a matrix-vector multiplication. Quaternions are more complicated in this regard, but they offer a mathematical operation called slerp (shorthand for spherical linear interpolation) that is not available with other representations, like rotation matrices. We assume that the reader knows what quaternions are and their mathematical properties; if you are not, be sure to read the book Visualizing quaternions, by Andrew J. Hanson (Elsevier, 2006, ISBN-0080474772) and the provocative essay by Marc ten Bosch, Let’s remove Quaternions from every 3D engine.

The LiteBIRD simulation framework models quaternions using the convention \((v_x, v_y, v_z, w)\); be aware that some textbooks use the order \((w, v_x, v_y, v_z)\). As the framework uses quaternions only to model rotations, they all must obey the relation \(v_x^2 + v_y^2 + v_z^2 + w^2 = 1\) (normalized quaternions), which is a property satisfied by rotation quaternions.

The class RotQuaternion can model time-varying quaternions. It is enough to provide a list of quaternions, a starting time, and a sampling frequency, which is assumed to be constant:

import litebird_sim as lbs

time_varying_quaternion = lbs.RotQuaternion(
    # Three rotation quaternions
    quats=np.array(
        [
            [0.5, 0.0, 0.0, 0.8660254],
            [0.0, -0.38268343, 0.0, 0.92387953],
            [0.0, 0.0, 0.30901699, 0.95105652],
        ]
    ),
    start_time=0.0,
    sampling_rate_hz=1.0,
)

This example assumes that time_varying_quaternion describes a rotation that evolves with time, starting from t = 0 and lasting 3 seconds, as the sampling frequency is 1 Hz.

Rotation quaternions can be multiplied together; however, they must refer to the same starting time and have the same sampling frequency.

Python functions for quaternions#

The LiteBIRD Simulation Framework provides three functions, quat_rotation_x(), quat_rotation_y(), quat_rotation_z() to compute simple rotation quaternions; they return plain the normalized quaternion representing a rotation by an angle \(\theta\) around one of the three axis x, y, and z. These quaternions are plain NumPy arrays and can be passed to the parameter quats of the constructor for RotQuaternion:

import litebird_sim as lbs
import numpy as np

def print_quaternion(q):
    print("{:.3f} {:.3f} {:.3f} {:.3f}".format(*q))

print("Rotation by π/3 around x:")
print_quaternion(lbs.quat_rotation_x(theta_rad=np.pi/3))
print("Rotation by π/3 around y:")
print_quaternion(lbs.quat_rotation_y(theta_rad=np.pi/3))
print("Rotation by π/3 around z:")
print_quaternion(lbs.quat_rotation_z(theta_rad=np.pi/3))

Rotation by π/3 around x:
0.500 0.000 0.000 0.866
Rotation by π/3 around y:
0.000 0.500 0.000 0.866
Rotation by π/3 around z:
0.000 0.000 0.500 0.866

There are two functions that implement in-place multiplication of quaternions: quat_right_multiply() performs the calculation \(r \leftarrow r \times q\), and quat_left_multiply() performs the calculation \(r \leftarrow q \times r\) (where \(\leftarrow\) indicates the assignment operation):

quat = np.array(lbs.quat_rotation_x(np.pi / 3))
lbs.quat_left_multiply(quat, *lbs.quat_rotation_z(np.pi / 2))
print("Rotation by π/3 around x and then by π/2 around z:")
print_quaternion(quat)

Rotation by π/3 around x and then by π/2 around z:
0.354 0.354 0.612 0.612

Note the syntax for quat_left_multiply(): you are supposed to pass the four components of the quaternion \(q\) as separate arguments, and thus we need to prepend the call to lbs.quat_rotation_z with * to expand the result (a 4-element tuple) into the four parameters required by quat_left_multiply(). The reason for this weird syntax is efficiency, as Numba (which is used extensively in the code) can easily optimize this kind of function call.

Finally, the framework provides the function rotate_vector(), which applies the rotation described by a normalized quaternion to a vector. There are faster versions in rotate_x_vector(), rotate_y_vector(), and rotate_z_vector() that rotate the three basis vectors (1, 0, 0), (0, 1, 0), and (0, 0, 1). The functions all_rotate_vectors(), all_rotate_x_vectors(), all_rotate_y_vectors(), and all_rotate_z_vectors() can be applied to arrays of quaternions and vectors.

API reference#

litebird_sim.quaternions.all_rotate_vectors(result_matrix, quat_matrix, vec_matrix)#

Rotate a set of vectors using quaternions

Prototype:

all_rotate_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
    vec_matrix: numpy.array[N, 3],
)

Assuming that result_matrix and vec_matrix are two NumPy arrays with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_vector() to each row.

litebird_sim.quaternions.all_rotate_x_vectors(result_matrix, quat_matrix)#

Rotate the vector [1, 0, 0] using quaternions

Prototype:

all_rotate_x_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_x_vector() to each row.

litebird_sim.quaternions.all_rotate_y_vectors(result_matrix, quat_matrix)#

Rotate the vector [0, 1, 0] using quaternions

Prototype:

all_rotate_y_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_y_vector() to each row.

litebird_sim.quaternions.all_rotate_z_vectors(result_matrix, quat_matrix)#

Rotate the vector [0, 0, 1] using quaternions

Prototype:

all_rotate_z_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_z_vector() to each row.

litebird_sim.quaternions.multiply_quaternions_list_x_list(array_a: ndarray[tuple[Any, ...], dtype[_ScalarT]], array_b: ndarray[tuple[Any, ...], dtype[_ScalarT]], result: ndarray[tuple[Any, ...], dtype[_ScalarT]]) → None#

Multiply two sets of quaternions together

All the matrices must have the same shape (N, 4). The result of a × b is saved into result.

litebird_sim.quaternions.multiply_quaternions_list_x_one(array_a: ndarray[tuple[Any, ...], dtype[_ScalarT]], single_b: ndarray[tuple[Any, ...], dtype[_ScalarT]], result: ndarray[tuple[Any, ...], dtype[_ScalarT]]) → None#

Multiply a matrix of quaternions by one quaternion: array_a × single_b.

The result of a × b is saved into result.

litebird_sim.quaternions.multiply_quaternions_one_x_list(single_a: ndarray[tuple[Any, ...], dtype[_ScalarT]], array_b: ndarray[tuple[Any, ...], dtype[_ScalarT]], result: ndarray[tuple[Any, ...], dtype[_ScalarT]]) → None#

Multiply one quaternion by a matrix of quaternions: single_a × array_b.

The result of a × b is saved into result.

litebird_sim.quaternions.normalize_quaternions(quat_matrix: ndarray[tuple[Any, ...], dtype[_ScalarT]]) → None#

Normalize all the quaternions in a matrix with shape (N, 4)

All the quaternions are normalized in place.

litebird_sim.quaternions.quat_left_multiply(result, other_v1, other_v2, other_v3, other_w)#

Perform a multiplication between two quaternions

Prototype:

quat_left_multiply(
    result: numpy.array[3],
    other_v1: float,
    other_v2: float,
    other_v3: float,
    other_w: float,
)

This function implements the computation \(r = q \times r\); see also quat_right_multiply() for the computation \(r = r\times q\).

It’s easy to use NumPy quaternions for q as well:

import numpy as np
r = np.array([1.0, 2.0, 3.0, 4.0])
q = np.array([0.1, 0.2, 0.3, 0.4])
quat_right_multiply(r, *q)  # Unpack "q"
print("Result:", r)

litebird_sim.quaternions.quat_right_multiply(result, other_v1, other_v2, other_v3, other_w)#

Perform a multiplication between two quaternions

Prototype:

quat_right_multiply(
    result: numpy.array[3],
    other_v1: float,
    other_v2: float,
    other_v3: float,
    other_w: float,
)

This function implements the computation \(r = r \times q\), where r is the parameter result (a 3-element NumPy array) and q is the set of parameters other_v1, other_v2, other_v3, other_w. The reason why the elements of quaternion q are passed one by one is efficiency: in this way, the caller does not have to allocate a numpy.array for simple quaternions (like the ones returned by quat_rotation_x(), quat_rotation_y(), quat_rotation_z()).

It’s easy to use NumPy quaternions for q as well:

import numpy as np
r = np.array([1.0, 2.0, 3.0, 4.0])
q = np.array([0.1, 0.2, 0.3, 0.4])
quat_right_multiply(r, *q)  # Unpack "q"
print("Result:", r)

See also quat_left_multiply() for the computation \(r = q imes r\).

litebird_sim.quaternions.quat_rotation(theta_rad, x, y, z)#: This function rotates a quaternion by theta_rad about a specific axis. :param theta_rad: The angle to rotate the quaternion by in radians. :type theta_rad: float :param x: The x element of vector to rotate the quaternion about. :type x: float :param y: The y element of vector to rotate the quaternion about. :type y: float :param z: The z element of vector to rotate the quaternion about. :type z: float

litebird_sim.quaternions.quat_rotation_brdcast(theta_rad, vec_matrix)#: This function rotates a quaternion by theta_rad about a specific axis. :param theta_rad: The angle to rotate the quaternion by in radians. :type theta_rad: float :param vec_matrix: The vectors to rotate the quaternion about. :type vec_matrix: numpy.array[N, 3]

litebird_sim.quaternions.quat_rotation_x(theta_rad)#

Return a quaternion representing a rotation around the x axis

Prototype:

quat_rotation_x(theta_rad: float) -> tuple[float, float, float, float]

The angle theta_rad must be expressed in radians. The return value is the quaternion, using the order (v_x, v_y, v_z, w); it is returned as a 4-element tuple.

The fact that the result is a tuple instead of a NumPy array is because of speed: it helps in preventing unnecessary allocations in performance-critical code.

litebird_sim.quaternions.quat_rotation_x_brdcast(theta_rad_array)#

Return a quaternion representing a rotation around the x axis

Prototype:

quat_rotation_x(theta_rad: np.ndarray) -> np.ndarray

The angles theta_rad must be expressed in radians. The return value is the quaternion, using the order (v_x, v_y, v_z, w); it is returned as a 2D array with shape (n, 4), where n is the number of angles.

litebird_sim.quaternions.quat_rotation_y(theta_rad)#

Return a quaternion representing a rotation around the y axis

Prototype:

quat_rotation_y(theta_rad: float) -> tuple[float, float, float, float]

litebird_sim.quaternions.quat_rotation_y_brdcast(theta_rad_array)#

Return a quaternion representing a rotation around the y axis

Prototype:

quat_rotation_y(theta_rad: np.ndarray) -> np.ndarray

litebird_sim.quaternions.quat_rotation_z(theta_rad)#

Return a quaternion representing a rotation around the y axis

Prototype:

quat_rotation_z(theta_rad: float) -> tuple[float, float, float, float]

litebird_sim.quaternions.quat_rotation_z_brdcast(theta_rad_array)#

Return a quaternion representing a rotation around the z axis

Prototype:

quat_rotation_z(theta_rad: np.ndarray) -> np.ndarray

litebird_sim.quaternions.rotate_vector(result, vx, vy, vz, w, vect)#

Rotate a vector using a quaternion

Prototype:

rotate_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
    vect: numpy.array[3],
)

Applies a rotation, encoded through the quaternion vx, vy, vz, vw, to the vector vect (a 3-element NumPy array), storing the result in result (again a 3-element array).

Note: do not pass the same variable to vect and result!

The formula to rotate a vector v by a quaternion (q_v, w) is the following: \(v' = v + 2q_v ⨯ (q_v ⨯ v + w v)\), where q_v is the vector (vx, vy, vz).

litebird_sim.quaternions.rotate_x_vector(result, vx, vy, vz, w)#

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_x_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [1, 0, 0]), but it’s faster.

litebird_sim.quaternions.rotate_y_vector(result, vx, vy, vz, w)#

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_y_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [0, 1, 0]), but it’s faster.

litebird_sim.quaternions.rotate_z_vector(result, vx, vy, vz, w)#

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_z_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [0, 0, 1]), but it’s faster.