"""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 of shape (4,) or (n, 4)
the array of resulting imaginary quaternions for each point of the input array
Examples
--------
>>> import numpy as np
>>> from ddg.math.quaternion import point_to_quaternion
>>> point_to_quaternion(np.arange(3))
array([0, 0, 1, 2])
>>> x = np.arange(12).reshape((4, 3))
>>> x
array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11]])
>>> point_to_quaternion(x)
array([[ 0, 0, 1, 2],
[ 0, 3, 4, 5],
[ 0, 6, 7, 8],
[ 0, 9, 10, 11]])
"""
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 of shape (4,) or (n, 4)
y : np.ndarray of shape (4,) or (n, 4)
Returns
-------
np.ndarray of shape (4,) or (n, 4)
Examples
--------
>>> import numpy as np
>>> from ddg.math.quaternion import multiply
>>> x = np.arange(4)
>>> y = np.array([5, 6, 7, 8])
>>> multiply(x, y)
array([-44, 0, 20, 10])
>>> x = np.arange(8).reshape((2, 4))
>>> y = y.reshape((1, 4))
>>> multiply(x, y)
array([[ -44, 0, 20, 10],
[-108, 48, 60, 66]])
>>> multiply(x, x)
array([[-14, 0, 0, 0],
[-94, 40, 48, 56]])
"""
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 of shape (4,) or (n, 4)
the input array of quaternions
Returns
-------
np.ndarray of shape (4,) or (n, 4)
an array of inversed values for each element of the input
Raises
------
ValueError
If any of the quaternions is the zero quaternion.
If any of the quaternions doesn't have an inverse.
Examples
--------
>>> import numpy as np
>>> from ddg.math.quaternion import inverse
>>> x = np.arange(8).reshape((2, 4))
>>> inverse(x)
array([[ 0. , -0.07142857, -0.14285714, -0.21428571],
[ 0.03174603, -0.03968254, -0.04761905, -0.05555556]])
>>> inverse(np.array([10, 5, 0, 1]))
array([ 0.07936508, -0.03968254, 0. , -0.00793651])
"""
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 of shape (4,) or (n, 4)
First quaternion or array of quaternions
X2 : np.ndarray of shape (4,) or (n, 4)
Second quaternion or array of quaternions
X3 : np.ndarray of shape (4,) or (n, 4)
Third quaternion or array of quaternions
X4 : np.ndarray of shape (4,) or (n, 4)
Fourth quaternion or array of quaternions
Returns
-------
np.ndarray of shape (4,) or (n, 4)
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.quaternion import cr
>>> X = np.arange(16).reshape((4, 4))
>>> X
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> cr(X[0], X[1], X[2], X[3])
array([-0.33333333, 0. , 0. , 0. ])
"""
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 of shape (4,) or (n, 4)
the input array of quaternions
Returns
-------
np.ndarray of shape (3,) or (n, 3)
Imaginary parts of an array of quaternions.
Examples
--------
>>> import numpy as np
>>> from ddg.math.quaternion import imag
>>> imag(np.array([10, 15, 20, 30]))
array([15, 20, 30])
>>> x = np.arange(8).reshape((2, 4))
>>> x
array([[0, 1, 2, 3],
[4, 5, 6, 7]])
>>> imag(x)
array([[1, 2, 3],
[5, 6, 7]])
"""
if X.ndim != 1:
return X[:, 1:]
else:
return X[1:]