I want to do this with as few copies as possible
This should not be a concern while your function is incorrect (your "sign error"). But also: should this -
vec_a = (1, 0)
vec_b = (0, 1)
really have the opposite angle sign from this?
vec_a = (1, 0)
vec_b = (0, 2)
A vector's argument is invariant under uniform scaling, so I think the answer must be "no", regardless of whether you're running a generic cosine similarity or some exotic low-angle algorithm.
Other bits -
Don't use two calls to asarray. If you're making a library-like function that does array coercion, then you're better off issuing a single call to broadcast_vectors.
I really don't understand why your code is prepending a new axis. It shouldn't need to do that.
Your code at least superficially resembles your linked Heron's method, but if it doesn't work, then no amount of optimisation makes sense. I suggest performing a comparative test between it and a better-known method such as the cosine similarity demonstration below.
import numpy as np
def angle_between(vec_a: np.typing.ArrayLike, vec_b: np.typing.ArrayLike, *, axis: int = -1) -> np.ndarray:
"""
Computes the angle in signed radians from vec_a to vec_b in a right-handed frame.
In all cases the magnitude of the angle follows simple cosine similarity.
In the two-dimensional case, the angle sign is based on a simple vector-argument comparison.
In the three-dimensional case, consider the plane shared by vec_a and vec_b. The sign of
the plane's normal is chosen so that it minimises the angular deflection from an
all-ones vector to disambiguate chirality. The angle sign is then the same as the 2D case
for the two vectors projected to the plane oriented by that normal.
Signed rotation demands a well-defined plane of rotation to disambiguate chirality.
Since two vectors cannot fix a hyperplane in four-dimensional or higher space, those
higher cases are all returned with a positive sign.
"""
vec_a, vec_b = np.broadcast_arrays(vec_a, vec_b)
# Simple cosine similarity
abs_angle = np.acos(
np.vecdot(vec_a, vec_b, axis=axis)
/np.linalg.norm(vec_a, axis=axis)
/np.linalg.norm(vec_b, axis=axis)
)
dim = vec_a.shape[axis]
if 2 <= dim <= 3:
if dim == 3:
# normal vector orienting the plane intersecting a and b
normal = np.linalg.cross(vec_a, vec_b, axis=axis)
# coerce the normal to conventional orientation
ones_basis = np.full(shape=normal.shape, fill_value=np.sqrt(1 / dim)) # use the normal closest to this
normal_to_ones = np.acos(
np.vecdot(normal, ones_basis, axis=axis)
/ np.linalg.norm(normal, axis=axis)
# the ones basis is already a unit vector, so no second norm in denominator
)
normal *= 1 - 2*(normal_to_ones > 0.5*np.pi) # -1 if > 180 # now the normal has conventional orientation
# project to the plane
u_basis = vec_a
v_basis = np.linalg.cross(normal, vec_a, axis=axis)
ax = np.vecdot(vec_a, u_basis, axis=axis) # u-dimension in first vector
bx = np.vecdot(vec_b, u_basis, axis=axis) # u-dimension in second vector
ay = np.vecdot(vec_a, v_basis, axis=axis) # v-dimension in first vector
by = np.vecdot(vec_b, v_basis, axis=axis) # v-dimension in second vector
else:
indexer = [slice(None)] * len(vec_a.shape) # to start, select all axes
xidx = list(indexer) # copy
xidx[axis] = 0 # select first (x) position in axis of interest
xidx = tuple(xidx) # multi-axis indexer must be tuple
ax = vec_a[xidx] # x-dimension in first vector
bx = vec_b[xidx] # x-dimension in second vector
yidx = list(indexer) # copy
yidx[axis] = 1 # select second (y) position in axis of interest
yidx = tuple(yidx) # multi-axis indexer must be tuple
ay = vec_a[yidx] # y-dimension in first vector
by = vec_b[yidx] # y-dimension in second vector
ta = np.atan2(ay, ax) # argument, first vector
tb = np.atan2(by, bx) # argument, second vector
sign = np.ones_like(ta) # sign vector, preserving shape and dtype
sign[tb < ta] = -1
else:
return abs_angle # implied sign of 1
return sign*abs_angle
def unit_test() -> None:
vec_a = (1, 0, 0)
vec_b = (0, 1, 0)
result = angle_between(vec_a, vec_b)
assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
result = angle_between(vec_b, vec_a)
assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
vec_a = (-0.5, 0, 0)
vec_b = ( 0, 3, 0)
result = angle_between(vec_a, vec_b)
assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
result = angle_between(vec_b, vec_a)
assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
vec_a = (1, 0)
vec_b = (0, 1)
result = angle_between(vec_a, vec_b)
assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
result = angle_between(vec_b, vec_a)
assert np.isclose(result, -0.5*np.pi)
vec_a = (1, 0)
vec_b = (0, 2)
result = angle_between(vec_a, vec_b)
assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
result = angle_between(vec_b, vec_a)
assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
vec_a = (
(1, 0),
(2, 0),
(3, 0),
)
vec_b = (
(0, 2),
(0, 1),
(0, 0.5),
)
result = angle_between(vec_a, vec_b)
assert np.allclose(result, (0.5*np.pi, 0.5*np.pi, 0.5*np.pi), atol=0, rtol=1e-14)
result = angle_between(vec_b, vec_a)
assert np.allclose(result, (-0.5*np.pi, -0.5*np.pi, -0.5*np.pi), atol=0, rtol=1e-14)
if __name__ == '__main__':
unit_test()
