Skip to content

Commit

Permalink
Allow numpy to compute the angle of a quaternion (#38)
Browse files Browse the repository at this point in the history
If `q = np.exp(v̂ * θ/2)` for some unit vector `v̂` and an angle `θ` ∈[-2π,2π],
then `np.angle` returns `abs(θ)`.  This equals `2*abs(log(q))`, but is more
efficient.
  • Loading branch information
moble authored Mar 27, 2022
1 parent 7c4e853 commit 73e6ca5
Show file tree
Hide file tree
Showing 3 changed files with 52 additions and 4 deletions.
19 changes: 19 additions & 0 deletions quaternionic/algebra.py
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,25 @@ def log(q, qout):
qout[3] = f * q[3]


@attach_typelist_and_signature([(float64[:], float64[:])], '(n)->()')
def angle(q, qout):
"""Return angle (in radians) through which input quaternion rotates a vector
If `q = np.exp(v̂ * θ/2)` for some unit vector `v̂` and an angle `θ` ∈[-2π,2π],
then this function returns `abs(θ)`. This equals 2*abs(log(q)), but is more
efficient.
"""
b = np.sqrt(q[1]**2 + q[2]**2 + q[3]**2)

This comment has been minimized.

Copy link
@evbernardes

evbernardes Nov 21, 2022

Wouldn't it be better to calculate this as np.linalg.norm(q[1:]) instead of np.sqrt(q[1]**2 + q[2]**2 + q[3]**2)?

This comment has been minimized.

Copy link
@moble

moble Nov 21, 2022

Author Owner

The results should be identical, but there are two good reasons not to use that:

  1. Using np.linalg.norm(q[1:]) means that people who don't just know that function have to look it up to understand what's going on.
  2. The function np.linalg.norm is a very general and complicated piece of code, which makes it slow. In fact, it's about twice as slow as just taking the square-root of the sum of squares, even after numba compilation.

There are situations where we actually need the greater generality of np.linalg.norm, but this isn't one, because we know exactly what q is. So it's best to just stick with the most obvious approach.

This comment has been minimized.

Copy link
@evbernardes

evbernardes Nov 21, 2022

I see! Thanks :)

if b <= _quaternion_resolution * np.abs(q[0]):
if q[0] < 0.0:
qout[0] = 2*np.pi
else:
qout[0] = 0.0
else:
qout[0] = 2*np.abs(np.arctan2(b, q[0]))


@attach_typelist_and_signature([(float64[:], float64[:])], '(n)->(n)')
def sqrt(q, qout):
"""Return square-root of input quaternion √q.
Expand Down
16 changes: 12 additions & 4 deletions quaternionic/arrays.py
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,18 @@ def __array_finalize__(self, obj):
)

def __array_function__(self, func, types, args, kwargs):
output = super().__array_function__(func, types, args, kwargs)
if func in [np.ones_like]:
# Want the last dimension to equal [1,0,0,0] not [1,1,1,1]
output.vector = 0
from . import algebra_ufuncs as algebra

if func == np.angle:
output = np.zeros(args[0].shape[:-1], dtype=dtype)
algebra.angle(args[0].ndarray, output)
if kwargs.get("deg", False):
output *= 180/np.pi
else:
output = super().__array_function__(func, types, args, kwargs)
if func in [np.ones_like]:
# Want the last dimension to equal [1,0,0,0] not [1,1,1,1]
output.vector = 0
return output

def __array_ufunc__(self, ufunc, method, *args, **kwargs):
Expand Down
21 changes: 21 additions & 0 deletions tests/test_algebra.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,27 @@ def test_basis_multiplication():
assert i*j*k == -one


def test_array_function(array):
np.random.seed(1234)

# Test np.angle
assert np.angle(array(1, 0, 0, 0)) == 0.0
assert np.angle(array(-1, 0, 0, 0)) == 2*np.pi
assert np.angle(array(0, 1, 0, 0)) == np.pi
assert np.angle(array(0, -1, 0, 0)) == np.pi
assert np.angle(array(0, 0, 1, 0, 0)) == np.pi
assert np.angle(array(0, 0, -1, 0, 0)) == np.pi
assert np.angle(array(0, 0, 0, 1)) == np.pi
assert np.angle(array(0, 0, 0, -1)) == np.pi
angle_precision = 2e-15
for θ in np.random.uniform(-2*np.pi, 2*np.pi, size=100):
for _ in range(20):
= array.from_vector_part(np.random.normal(size=3)).normalized
R = np.exp( * θ/2)
assert abs(abs(θ) - np.angle(R)) < angle_precision
assert abs(abs(np.rad2deg(θ)) - np.angle(R, deg=True)) < np.rad2deg(angle_precision)


def test_array_ufunc(array):
np.random.seed(1234)

Expand Down

0 comments on commit 73e6ca5

Please sign in to comment.