Angles between two n-dimensional vectors in Python

Question

I need to determine the angle(s) between two n-dimensional vectors in Python. For example, the input can be two lists like the following: [1,2,3,4] and [6,7,8,9].

Answer 1

There are numerous solutions to this question, yet none employ a simple expression using numpy capable of handling n-dimensional arrays across any axis, as given by:

import numpy as np
def angle(vec_0, vec_1, axis):
    return np.arctan2(np.cross(vec_0, vec_1, axis = axis), np.sum(vec_0*vec_1, axis = axis))

---

Example

We select points on a grid and draw two vectors at each point: one in black and the other in grey. Then, we calculate the angle between the black and grey vectors, and plot the points colored according to this angle.

#Define grid of N*N
N = 4
rng = np.random.default_rng(10)

pos =  np.array(np.meshgrid(np.linspace(0,1,N),np.linspace(0,1,N), indexing = 'ij')).T
vec_black = rng.uniform(-1, 1, (N*N, 2))
vec_grey = rng.uniform(-1, 1, (N*N, 2))
theta_black_grey = angle(vec_black, vec_grey, axis=-1)

import matplotlib.pyplot as plt
width = 0.008
scale = 8
plt.figure()
plt.quiver(*pos.T, *vec_black.T, color = 'k', width = width, scale = scale)
plt.quiver(*pos.T, *vec_grey.T, color = 'grey', width = width, scale = scale)
plt.scatter(*pos.T, c=theta_black_grey,  cmap ='hsv', vmin=-np.pi, vmax = np.pi)
plt.xlim([-0.3,1.3])
plt.ylim([-0.3,1.3])
plt.colorbar(label='Angle (black arrow, grey arrow)')

Output:

Answer 2

Given two vectors vec(u) and vec(v), then it can be shown that the following formul is the most numerically stable solution:

2 * atan ( norm(norm(u) * v - norm(v) * u)/norm(norm(u) * v + norm(v) * u) )

The benefits of this formula over any of:

acos(dot_product(u,v)/(norm(u)*norm(v)))
asin(norm(cross_product(u,v))/(norm(u)*norm(v)))

are:

1. atan is numerically stable for small angles. The acos function is not as cos(th) ~ 1 - 0.5*th^2 for small angles
2. atan is numerically stable for angles around Pi/2. As the formula computes the half-angle, it corresponds to the angle Pi for the formula using the cross-product. Under these conditions, the cross-product computation is unstable and hence then formula using asin is unstable.

One could pre-normalize the vectors u and v and use the form:

2 * atan ( norm(v - u)/norm(v + u) )

however, a division actually looses numeric accuracy. In this case this loss of precision would appear in the normalizations and the final division.

Hence, implementation wise, we should use the classic atan2 function which avoids another division. Hence, we achieve the best numeric stable solution.

So, using numpy, the implementation is straightforward:

_nu = numpy.linalg.norm(u)
_nv = numpy.linalg.norm(v)
2.0 * numpy.arctan2( numpy.linalg.norm( u * _nv - v * _nu),
                     numpy.linalg.norm( u * _nv + v * _nu))

Answer 3

The OP was asking for n-dimensions n>2, but many people will end up here for a two dimensional problem, so I want to clarify the best solution for that special case.

The proposed solutions all use the arccosine function (other than MK83 and faken) and so, are all very inaccurate and prone to error for angles near 0 or 180 degrees. That is because very large changes in angle cause almost no change in the cosine of the angle at those values.

A second problem with arccos (for the two dimensional case) is that it can not distinguish between a positive angle and a negative angle. So the angle between (0,1) and (1,0) will be the same as that between (1,0) and (0,1) although the first angle should be 90 and the second -90 degrees.

faken has an excellent answer to to OPs multidimensional problem that avoids the use of arccos, and so is accurate over the whole range.

MK83 solves the 2 dimensional problem using atan2which is a function that is provided for this exact problem. The answers range from -180 degrees to 180 degrees. I propose a solution here only for two dimensions, which is simpler and faster than MK83

def angle(a, b, c=None):
    """ This function computes angle between vector A and vector B when C is None
        and the angle between AC and CB, when C is a vector as well.
        All vectors must be two dimensional.
    """
    if c is None:
        angles = np.arctan2([a[1], b[1]], [a[0], b[0]]])
    else:
        angles = np.arctan2([a[1]-c[1], b[1]-c[1]], [a[0]-c[0], b[0]-c[0]])
    return np.degrees(angles[1] - angles[0])

This is about three times faster than MK83's solution.

Answer 4

import math

ax, ay = input('Enter x and y of vector a: ').split()

ax, ay = float(ax), float(ay)

bx, by = input('Enter x and y of vector b: ').split()

bx, by = float(bx), float(by)

ab = ax * bx + ay * by

a = math.sqrt(ax * ax + ay * ay)

b = math.sqrt(bx * bx + by * by)

cos = ab / (a*b)

rad = math.acos(cos)

deg = math.degrees(rad)

print (f'θ = {deg}')

Answer 5

Using numpy (highly recommended), you would do:

from numpy import (array, dot, arccos, clip)
from numpy.linalg import norm

u = array([1.,2,3,4])
v = ...
c = dot(u,v)/norm(u)/norm(v) # -> cosine of the angle
angle = arccos(clip(c, -1, 1)) # if you really want the angle

Answer 6

The traditional approach to obtaining an angle between two vectors (i.e. arccos(dot(u, v) / (norm(u) * norm(v))), as presented in the other answers) suffers from numerical instability in several corner cases. The following code works for n-dimensions and in all corner cases (it doesn't check for zero length vectors, but that's easy to add, as shown in some of the other answers). See notes below.

from numpy import arctan, pi, signbit
from numpy.linalg import norm


def angle_btw(v1, v2):
    u1 = v1 / norm(v1)
    u2 = v2 / norm(v2)

    y = u1 - u2
    x = u1 + u2

    a0 = 2 * arctan(norm(y) / norm(x))

    if (not signbit(a0)) or signbit(pi - a0):
        return a0
    elif signbit(a0):
        return 0.0
    else:
        return pi

This code is adapted from a Julia implementation by Jeffrey Sarnoff (MIT license), in turn based on these notes by Prof. W. Kahan (page 15).

Answer 7

Use some functions from numpy.

import numpy as np

def dot_product_angle(v1,v2):

    if np.linalg.norm(v1) == 0 or np.linalg.norm(v2) == 0:
        print("Zero magnitude vector!")
    else:
        vector_dot_product = np.dot(v1,v2)
        arccos = np.arccos(vector_dot_product / (np.linalg.norm(v1) * np.linalg.norm(v2)))
        angle = np.degrees(arccos)
        return angle
    return 0

Answer 8

Easy way to find angle between two vectors(works for n-dimensional vector),

Python code:

import numpy as np

vector1 = [1,0,0]
vector2 = [0,1,0]

unit_vector1 = vector1 / np.linalg.norm(vector1)
unit_vector2 = vector2 / np.linalg.norm(vector2)

dot_product = np.dot(unit_vector1, unit_vector2)

angle = np.arccos(dot_product) #angle in radian

Answer 9

Building on sgt pepper's great answer and adding support for aligned vectors plus adding a speedup of over 2x using Numba

@njit(cache=True, nogil=True)
def angle(vector1, vector2):
    """ Returns the angle in radians between given vectors"""
    v1_u = unit_vector(vector1)
    v2_u = unit_vector(vector2)
    minor = np.linalg.det(
        np.stack((v1_u[-2:], v2_u[-2:]))
    )
    if minor == 0:
        sign = 1
    else:
        sign = -np.sign(minor)
    dot_p = np.dot(v1_u, v2_u)
    dot_p = min(max(dot_p, -1.0), 1.0)
    return sign * np.arccos(dot_p)

@njit(cache=True, nogil=True)
def unit_vector(vector):
    """ Returns the unit vector of the vector.  """
    return vector / np.linalg.norm(vector)

def test_angle():
    def npf(x):
        return np.array(x, dtype=float)
    assert np.isclose(angle(npf((1, 1)), npf((1,  0))),  pi / 4)
    assert np.isclose(angle(npf((1, 0)), npf((1,  1))), -pi / 4)
    assert np.isclose(angle(npf((0, 1)), npf((1,  0))),  pi / 2)
    assert np.isclose(angle(npf((1, 0)), npf((0,  1))), -pi / 2)
    assert np.isclose(angle(npf((1, 0)), npf((1,  0))),  0)
    assert np.isclose(angle(npf((1, 0)), npf((-1, 0))),  pi)

%%timeit results without Numba

• 359 µs ± 2.86 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

And with

• 151 µs ± 820 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Answer 10

David Wolever's solution is good, but

If you want to have signed angles you have to determine if a given pair is right or left handed (see wiki for further info).

My solution for this is:

def unit_vector(vector):
    """ Returns the unit vector of the vector"""
    return vector / np.linalg.norm(vector)

def angle(vector1, vector2):
    """ Returns the angle in radians between given vectors"""
    v1_u = unit_vector(vector1)
    v2_u = unit_vector(vector2)
    minor = np.linalg.det(
        np.stack((v1_u[-2:], v2_u[-2:]))
    )
    if minor == 0:
        raise NotImplementedError('Too odd vectors =(')
    return np.sign(minor) * np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

It's not perfect because of this NotImplementedError but for my case it works well. This behaviour could be fixed (cause handness is determined for any given pair) but it takes more code that I want and have to write.

Answer 11

For the few who may have (due to SEO complications) ended here trying to calculate the angle between two lines in python, as in (x0, y0), (x1, y1) geometrical lines, there is the below minimal solution (uses the shapely module, but can be easily modified not to):

from shapely.geometry import LineString
import numpy as np

ninety_degrees_rad = 90.0 * np.pi / 180.0

def angle_between(line1, line2):
    coords_1 = line1.coords
    coords_2 = line2.coords

    line1_vertical = (coords_1[1][0] - coords_1[0][0]) == 0.0
    line2_vertical = (coords_2[1][0] - coords_2[0][0]) == 0.0

    # Vertical lines have undefined slope, but we know their angle in rads is = 90° * π/180
    if line1_vertical and line2_vertical:
        # Perpendicular vertical lines
        return 0.0
    if line1_vertical or line2_vertical:
        # 90° - angle of non-vertical line
        non_vertical_line = line2 if line1_vertical else line1
        return abs((90.0 * np.pi / 180.0) - np.arctan(slope(non_vertical_line)))

    m1 = slope(line1)
    m2 = slope(line2)

    return np.arctan((m1 - m2)/(1 + m1*m2))

def slope(line):
    # Assignments made purely for readability. One could opt to just one-line return them
    x0 = line.coords[0][0]
    y0 = line.coords[0][1]
    x1 = line.coords[1][0]
    y1 = line.coords[1][1]
    return (y1 - y0) / (x1 - x0)

And the use would be

>>> line1 = LineString([(0, 0), (0, 1)]) # vertical
>>> line2 = LineString([(0, 0), (1, 0)]) # horizontal
>>> angle_between(line1, line2)
1.5707963267948966
>>> np.degrees(angle_between(line1, line2))
90.0

Answer 12

If you're working with 3D vectors, you can do this concisely using the toolbelt vg. It's a light layer on top of numpy.

import numpy as np
import vg

vec1 = np.array([1, 2, 3])
vec2 = np.array([7, 8, 9])

vg.angle(vec1, vec2)

You can also specify a viewing angle to compute the angle via projection:

vg.angle(vec1, vec2, look=vg.basis.z)

Or compute the signed angle via projection:

vg.signed_angle(vec1, vec2, look=vg.basis.z)

I created the library at my last startup, where it was motivated by uses like this: simple ideas which are verbose or opaque in NumPy.

Answer 13

import math

def dotproduct(v1, v2):
  return sum((a*b) for a, b in zip(v1, v2))

def length(v):
  return math.sqrt(dotproduct(v, v))

def angle(v1, v2):
  return math.acos(dotproduct(v1, v2) / (length(v1) * length(v2)))

Note: this will fail when the vectors have either the same or the opposite direction. The correct implementation is here: https://stackoverflow.com/a/13849249/71522

Answer 14

Note: all of the other answers here will fail if the two vectors have either the same direction (ex, (1, 0, 0), (1, 0, 0)) or opposite directions (ex, (-1, 0, 0), (1, 0, 0)).

Here is a function which will correctly handle these cases:

import numpy as np

def unit_vector(vector):
    """ Returns the unit vector of the vector.  """
    return vector / np.linalg.norm(vector)

def angle_between(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'::

            >>> angle_between((1, 0, 0), (0, 1, 0))
            1.5707963267948966
            >>> angle_between((1, 0, 0), (1, 0, 0))
            0.0
            >>> angle_between((1, 0, 0), (-1, 0, 0))
            3.141592653589793
    """
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

Answer 15

The other possibility is using just numpy and it gives you the interior angle

import numpy as np

p0 = [3.5, 6.7]
p1 = [7.9, 8.4]
p2 = [10.8, 4.8]

''' 
compute angle (in degrees) for p0p1p2 corner
Inputs:
    p0,p1,p2 - points in the form of [x,y]
'''

v0 = np.array(p0) - np.array(p1)
v1 = np.array(p2) - np.array(p1)

angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1))
print np.degrees(angle)

and here is the output:

In [2]: p0, p1, p2 = [3.5, 6.7], [7.9, 8.4], [10.8, 4.8]

In [3]: v0 = np.array(p0) - np.array(p1)

In [4]: v1 = np.array(p2) - np.array(p1)

In [5]: v0
Out[5]: array([-4.4, -1.7])

In [6]: v1
Out[6]: array([ 2.9, -3.6])

In [7]: angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1))

In [8]: angle
Out[8]: 1.8802197318858924

In [9]: np.degrees(angle)
Out[9]: 107.72865519428085

Answer 16

Using numpy and taking care of BandGap's rounding errors:

from numpy.linalg import norm
from numpy import dot
import math

def angle_between(a,b):
  arccosInput = dot(a,b)/norm(a)/norm(b)
  arccosInput = 1.0 if arccosInput > 1.0 else arccosInput
  arccosInput = -1.0 if arccosInput < -1.0 else arccosInput
  return math.acos(arccosInput)

Note, this function will throw an exception if one of the vectors has zero magnitude (divide by 0).