"""Quaternion arithmetic module. In this module quaternions are represented by
numpy arrays with 4 entries."""
import numpy as np
[docs]def point_to_quaternion(x):
r"""
Converts a point in :math:`\mathbf{R}^3` to an imaginary quaternion by adding
a zero as the first coordinate.
Parameters
----------
x : np.ndarray of shape `(n, 3)` or `(3,)`
the input array representing points in :math:`\mathbf{R}^3`.
Either a single point as 1D array of shape (3,) or an array of
points with shape (n, 3)
Returns
-------
np.ndarray
the array of resulting imaginary quaternions for each point of the input array
"""
if x.ndim == 1:
return np.insert(x, 0, 0)
else:
zeros = np.zeros((len(x), 1), dtype=x.dtype)
return np.concatenate([zeros, x], axis=1)
[docs]def multiply(x, y):
r""" Multiplication of quaternions.
Notice if we write a quaternion in a vectorize manner we have
.. math::
x &= [q_1, q_2, q_3, q_4] \\
y &= [p_1, p_2, p_3, p_4] \\
\Im(x) &:= [q_2, q_3, q_4] \\
\Im(y) &:= [p_2, p_3, p_4] \\
\Re(xy) &= q_{1}p_1 - \langle \Im(x), \Im(y) \rangle \\
\Im(xy) &= \Im(x) \times \Im(y) + q_{1}\Im(y) + p_{1}\Im(x)
If one of the inputs is a 2D array of several quaternions and the other
is a 2D array of a single quaternion, the output is the multiplication of
each element of the first array with the second array
Parameters
----------
x : np.ndarray
y : np.ndarray
Returns
-------
np.ndarray
"""
if x.ndim != 1:
re_x = x[:, 0]
re_y = y[:, 0]
im_x = x[:, 1:]
im_y = y[:, 1:]
re_xy = (re_x * re_y - (im_y * im_x).sum(axis=1))[:, None]
im_xy = np.cross(im_x, im_y) + re_x[:, None] * im_y + re_y[:, None] * im_x
xy = np.concatenate([re_xy, im_xy], axis=1)
else:
re_x = x[0]
re_y = y[0]
im_x = x[1:]
im_y = y[1:]
re_xy = re_x * re_y - np.dot(im_x, im_y)
im_xy = np.cross(im_x, im_y) + re_x * im_y + re_y * im_x
xy = np.concatenate([[re_xy], im_xy])
return xy
[docs]def inverse(X):
r"""
Calculates the inverse of a quaternion.
We have:
.. math::
x &= [q_1, q_2, q_3, q_4] \\
\bar{x} &= [q_1, -q_2, -q_3, -q_4] \\
||x||^{2} &= q_{1}^2 + q_{2}^2 + q_{3}^2 + q_{4}^2 \\
x^{-1} &= \frac{\bar{x}}{||x||^2}
Parameters
----------
X : np.ndarray
the input array of quaternions
Returns
-------
np.ndarray
an array of inversed values for each element of the input
"""
if X.ndim != 1:
if not X.any(axis=1).all():
raise ValueError("Zero quaternion does not have an inverse")
rotation_array = -1 * np.ones(X.shape)
rotation_array[:, 0] = np.ones(X.shape[0])
abs_x = 1 / np.sum(X * X, axis=1)
else:
if not X.any():
raise ValueError("Zero quaternion does not have an inverse")
abs_x = 1 / np.dot(X, X)
rotation_array = np.array([1, -1, -1, -1])
return (abs_x * (rotation_array * X).T).T
[docs]def cr(X1, X2, X3, X4):
r"""
The cross ratio of four quaternions. For each four
quaternions, the cross ratio is calculated by the formula:
.. math:: cr(x_1, x_2, x_3, x_4) = \frac{(x_1-x_2)(x_3-x_4)}{(x_2-x_3)(x_4-x_1)}.
Parameters
----------
X1 : np.ndarray
First quaternion or array of quaternions
X2 : np.ndarray
Second quaternion or array of quaternions
X3 : np.ndarray
Third quaternion or array of quaternions
X4 : np.ndarray
Fourth quaternion or array of quaternions
Returns
-------
np.ndarray
cross ratio of the input quaternions
If four 2D arrays of quaternions are given as input, the resulting array
is the 2D array of cross ratios of corresponding quaternions.
"""
numerator = multiply(X1 - X2, X3 - X4)
denominator = multiply(X2 - X3, X4 - X1)
cr = multiply(numerator, inverse(denominator))
return cr
[docs]def imag(X):
"""
Finds the imaginary part of the quaternion.
Parameters
----------
X : np.ndarray
the input array of quaternions
Returns
-------
np.ndarray
Imaginary parts of an array of quaternions.
"""
if X.ndim != 1:
return X[:, 1:]
else:
return X[1:]