quaternionic 0.1.1

Interpret numpy arrays as quaternionic arrays with numba acceleration

Quaternions by way of numpy arrays

This module subclasses numpy's array type, interpreting the array as an array of quaternions, and accelerating the algebra using numba. There is also basic initial support for symbolic manipulation of quaternions by creating quaternionic arrays with sympy symbols as elements, though this is a work in progress.

This package has evolved from the quaternion package, which adds a quaternion dtype directly to numpy. In many ways, that is a much better approach because dtypes are built in to numpy, making it more robust than this package. However, that approach has its own limitations, including that it is harder to maintain, and requires much of the code to be written in C, which also makes it harder to distribute. This package is written entirely in python code, but should actually have comparable performance because it is compiled by numba. Moreover, because the core code is written in pure python, it is reusable for purposes other than the core purpose of this package, which is to provide the numeric array type.

Installation

Because this package is pure python code, it can be installed with the simplest tools. In particular, you can just run

pip install quaternionic

For development work, the best current option is poetry. From the top-level directory, run poetry install or just poetry run <some command>.

Usage

Basic construction

Any numpy array a with a last axis of size 4 (and dtype=float) can be reinterpreted as a quaternionic array with quaternionic.array(a):

import numpy as np
import quaternionic

a = 1.0 - np.random.rand(17, 11, 4)  # Just some random numbers; last dimension is 4
q1 = quaternionic.array(a)  # Reinterpret an existing array
q2 = quaternionic.array([1.2, 2.3, 3.4, 4.5])  # Create a new array

In this example, q1 is an array of 187 (17*11) quaternions, just to demonstrate that any number of dimensions may be used, as long as the final dimension has size 4.

Here, the original array a will still exist just as it was, and will behave just as a normal numpy array — including changing its values (which will change the values in q1), slicing, math, etc. However, q1 will be another "view" into the same data. Operations on q1 will be quaternionic. For example, whereas 1/a returns the element-wise inverse of each float in the array, 1/q1 returns the quaternionic inverse of each quaternion. Similarly, if you multiply two quaternionic arrays, their product will be computed with the usual quaternion multiplication, rather than element-wise multiplication of floats as numpy usually performs.

Algebra

All the usual quaternion operations are available, including

• Addition q1 + q2
• Subtraction q1 - q2
• Multiplication q1 * q2
• Division q1 / q2
• Scalar multiplication q1 * s == s * q1
• Scalar division q1 / s != s / q1
• Reciprocal np.reciprocal(q1) == 1/q1
• Exponential np.exp(q1)
• Logarithm np.log(q1)
• Square-root np.sqrt(q1)
• Conjugate np.conjugate(q1) == np.conj(q1)

All numpy ufuncs that make sense for quaternions are supported. When the arrays have different shapes, the usual numpy broadcasting rules take effect.

Attributes

In addition to the basic numpy array features, we also have a number of extra properties that are particularly useful for quaternions, including

• Methods to extract and/or set components
  • w, x, y, z
  • i, j, k (equivalent to x, y, z)
  • scalar, vector (equivalent to w, [x, y, z])
  • real, imag (equivalent to scalar, vector)
• Methods related to norms
  • abs (square-root of sum of squares of components)
  • norm (sum of squares of components)
  • modulus, magnitude (equal to abs)
  • absolute_square, abs2, mag2, squared_norm (equal to norm)
  • normalized
  • inverse
• Methods related to array infrastructure
  • ndarray (the numpy array underlying the quaternionic array)
  • flattened (all dimensions but last are flattened into one)
  • iterator (iterate over all quaternions)

Note that this package makes a distinction between abs and norm — the latter being the square of the former. This choice agrees with the Boost library's implementation of quaternions, as well as this package's forerunner quaternion. This also agrees with the corresponding functions on the C++ standard library's complex numbers. Because this may be confusing, a number of aliases are also provided that may be less confusing. For example, some people find the pair abs and abs2 to be more sensible.

Rotations

The most common application of quaternions is to representing rotations by means of unit quaternions. Note that this package does not restrict quaternions to have unit norms, since it is usually better for numerical purposes not to do so. For example, whereas rotation of a vector v by a quaternion is usually implemented as R * v * np.conjugate(R), it is generally better to drop the assumption that the quaternion has unit magnitude and implement rotation as R * v * np.reciprocal(R). That is what this package does by default whenever rotations are involved.

Although this package does not restrict to unit quaternions, there are several converters to and from other representations of rotations, including

• to/from_rotation_matrix
• to_transformation_matrix (for non-unit quaternions)
• to/from_axis_angle representation
• to/from_euler_angles (though using Euler angles is almost always a bad idea)
• to/from_spherical_coordinates
• to/from_angular_velocity

Note that the last item relates to quaternion-valued functions of time. Converting to an angular velocity requires differentiation, while converting from angular velocity requires integration (as explored in this paper).

For these converters, the "to" functions are properties on the individual arrays, whereas the "from" functions are "classmethod"s that take the corresponding objects as inputs. For example, we could write

q1 = quaternionic.array(np.random.rand(100, 4)).normalized
m = q1.to_rotation_matrix

to obtain the matrix m from a quaternionic array q1. (Here, m is actually a series of 100 3x3 matrices corresponding to the 100 quaternions in q1.) On the other hand, to obtain a quaternionic array from some matrix m, we would write

q2 = quaternionic.array.from_rotation_matrix(m)

Also note that, because the unit quaternions form a "double cover" of the rotation group (meaning that quaternions q and -q represent the same rotation), these functions are not perfect inverses of each other. In this case, for example, q1 and q2 may have opposite signs. We can, however, prove that these quaternions represent the same rotations by measuring the "distance" between the quaternions as rotations:

np.max(quaternionic.distance.rotation.intrinsic(q1, q2))  # Typically around 1e-15

Distance functions

The quaternionic.distance contains four distance functions:

• rotor.intrinsic
• rotor.chordal
• rotation.intrinsic
• rotation.chordal

The "rotor" distances do not account for possible differences in signs, meaning that rotor distances can be large even when they represent identical rotations; the "rotation" functions just return the smaller of the distance between q1 and q2 or the distance between q1 and -q2. So, for example, either "rotation" distance between q and -q is always zero, whereas neither "rotor" distance between q and -q will ever be zero (unless q=0). The "intrinsic" functions measure the geodesic distance within the manifold of unit quaternions, and is somewhat slower but may be more meaningful; the "chordal" functions measure the Euclidean distance in the (linear) space of all quaternions, and is faster but its precise value is not necessarily as meaningful.

These functions satisfy some important conditions. For each of these functions d, and for any nonzero quaternions q1, q2, a, and b, we have

• symmetry: d(q1, q2) = d(q2, q1)
• invariance: d(a*q1, a*q2) = d(q1, q2) = d(q1*b, q2*b)
• identity: d(q1, q1) = 0
• positive-definiteness:
  • For rotor functions d(q1, q2) > 0 whenever q1 ≠ q2
  • For rotation functions d(q1, q2) > 0 whenever q1 ≠ q2 and q1 ≠ -q2

Note that the rotation functions also satisfy d(q1, -q1) = 0.

See Moakher (2002) for a nice general discussion.

Interpolation

Finally, there are also capabilities related to interpolation

• slerp (spherical linear interpolation)
• squad (spherical quadratic interpolation)

Related packages

Packages with some quaternion features include

• quaternion (renamed numpy-quaternion on pypi due to name conflict)
• clifford (very powerful; more general geometric algebras)
• pyquaternion (many features; pure python; no acceleration)
• quaternions (basic pure python package; no acceleration; specialized for rotations only)
• rowan (similar approach to this package, but no acceleration or overloading)
• Quaternion (minimal capabilities; unmaintained)
• scipy.spatial.transform.Rotation.as_quat (quaternion output for Rotation object)
• mathutils (a Blender package with python bindings)

Also note that there is some capability to do symbolic manipulations of quaternions in these packages:

• galgebra
• sympy.algebras.quaternion