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].

Broome

This is the best answer is @MK83's as it is exactly the mathematical expression theta = atan2(u^v, u.v). even the case where u=[0 0] or v=[0 0] is covered because this is only time atan2 will produce the NaN in the other answers NaN will be produced by the / norm(u) or / norm(v) – Canvasback

Answer 1

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://mcmap.net/q/180949/-angles-between-two-n-dimensional-vectors-in-python

Stane

Also, if you only need cos, sin, tan of angle, and not the angle itself, then you can skip the math.acos to get cosine, and use cross product to get sine. – Lezley

Given that math.sqrt(x) is equivalent to x**0.5 and math.pow(x,y) is equivalent to x**y, I'm surprised these survived the redundancy axe wielded during the Python 2.x->3.0 transition. In practice, I'm usually doing these kinds of numeric things as part of a larger compute-intensive process, and the interpreter's support for '**' going directly to the bytecode BINARY_POWER, vs. the lookup of 'math', the access to its attribute 'sqrt', and then the painfully slow bytecode CALL_FUNCTION, can make a measurable improvement in speed at no coding or readability cost. – Pleurodynia

As in the answer with numpy: This can fail if the rounding error comes into play! This can happen for parallel and anti-parallel vectors! – Michamichael

Note: this will fail if the vectors are identical (ex, angle((1., 1., 1.), (1., 1., 1.))). See my answer for a slightly more correct version. – Karolekarolina

(err, and "by identical" I mean "have the same or opposite directions") – Karolekarolina

@DavidWolever This tripped me up so I'm going to have to clarify. The only time this is undefined is when the magnitude of one of the vectors is 0. Any other time an the cauchy-schwarz inequality states that the input to acos will be between -1 and 1. – Milks

@Milks uuhh… I don't believe that's true… As the vectors (1., 1., 1.), (1., 1., 1.) obviously have non-zero magnitudes, but the calculation results in nan… – Karolekarolina

norm(1,1,1)*norm(1,1,1) = 3. dot((1,1,1),(1,1,1)) = 3. acos(3/3) = 0. – Milks

If you're talking about the implementation above then it fails because of rounding errors, not because the vectors are parallel. – Milks

Simply use sklearn.metrics.pairwise import cosine_similarity and then acos – Britanybritches

you are genius, – Quintan

The most general answer was given by @MK83. Mathematically the angle is given by theta = atan2(u^v, u.v) is python terms angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1)) – Canvasback

Also, it is best to use the approach given by PilouPili, since acos and asin are not numerically stable for some values. Why does nobody teach this anymore? – Selenodont

Answer 2

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))

Karolekarolina

Wouldn't it be better to use np.isnan instead of the one from the math library? In theory they should be identical, but I'm not quite sure in practice. Either way I'd imagine it would be safer. – Petitionary

The only difference is that np.isnan will do something sensible if the input is an array, which will never be the case here. However, using np.isnan would definitely be cleaner (not sure why I used math.isnan…), so I'll switch that up. – Karolekarolina

My numpy (version==1.12.1) can use arccos directly and safely. : In [140]: np.arccos(np.dot(np.array([1,0,0]),np.array([-1,0,0]) )) Out[140]: 3.1415926535897931 In [141]: np.arccos(np.dot(np.array([1,0,0]),np.array([1,0,0]) )) Out[141]: 0.0 – Chemotherapy

The special case where at least one input vector is the zero vector is omitted, which is problematic for the division in unit_vector. One possibility is to just return the input vector in this function when this is the case. – Lulululuabourg

If you apply this to an array of vectors that are all positive (including zero), the maximum angle between two vectors is 90 degrees. i've tried scaling my data from -1 to 1, but i'm still getting angles strictly between 0 and 90 (instead of 0 and 180) – Cecilycecity

angle_between((0, 0, 0), (0, 1, 0)) will give nan as result, and not 90 – Gudrun

@Lulululuabourg 0-vectors' angle is undefined (in math). So the fact it raises an error is good. – Survival

@Chemotherapy it can still return nan if you use random unit vectors that happen to have a product of 1. Due to floating-point precision dot() would return 1.0000000000000002. This is true even for version 1.16.1. – Survival

I tried this out, but for v1=np.array([0,0,-1]) and v2=np.array([0,2,-4]), the output is 0. However, the angle between these two vectors is clearly not 0. Am I missing something? – Unanimous

@Unanimous I get 0.46364760900080615 as a result. – Ineducation

Answer 3

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

Beeves

The last line can result in an error as I've found out because of rounding errors. Thus if you to dot(u,u)/norm(u)**2 it results in 1.0000000002 and the arccos then fails (also 'works' for antiparallel vectors) – Michamichael

I've tested with u=[1,1,1]. u=[1,1,1,1] works fine but every dimension added returns slightly larger or smaler values than 1... – Michamichael

Note: this will fail (yield nan) when the direction of the two vectors is either identical or opposite. See my answer for a more correct version. – Karolekarolina

adding neo's comment to this, the last line should be angle = arccos(clip(c, -1, 1)) to avoid rounding issues. This solves @DavidWolever 's issue. – Chassidychassin

For the folks using the code snippet above: clip should be added to the list of numpy imports. – Cureton

Answer 4

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

Baxter

This is the best answer as it is exactly the mathematical expression theta = atan2(u^v, u.v). And this never fails ! – Canvasback

This is for 2-D. The OP was asking for n-D – Lexis

I used Right Angle Triangle as >>> p0, p1, p2 = [0, 0], [0, 3], [4, 0] >>> v1 = np.array(p2) - np.array(p1) >>> v0 = np.array(p0) - np.array(p1) >>> angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1)) >>> angle 0.9272952180016122 >>> np.degrees(angle) 53.13010235415598 .. Expected Angles in the triangles as 30,60,90 , not 53 >>> – Eidson

@Eidson a 3,4,5 right triangle is NOT a 30-60-90 triangle. The 30-60-90 triangle has sides 1, sqrt(3) and hypotenuse 2. – Mol

Given X,Y as [0,0],[0,3],[4,0] .. If would form a Triangle of Angles 30,60,90. Since one of the angle is 90, its right angle triangle. I am sorry, I didt understand why 'a 3,4,5 right triangle is NOT a 30-60-90 triangle' ? – Eidson

Below I have posted an answer that vectorizes this solution for any n-dimensional array. – Antoinetteanton

Answer 5

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.

Subsumption

Answer 6

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

Welt

Answer 7

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.

Reexamine

Answer 8

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).

Brader

I tried this out, but for v1=np.array([0,0,-1]) and v2=np.array([0,2,-4]), the output is 0. However, the angle between these two vectors is clearly not 0. Am I missing something? – Unanimous

@Unanimous I tested the values you suggested for v1 and v2, and the function is returning the correct result. Check here: imgur.com/a/EKTnBny Maybe you passed the incorrect variables to the function? – Brader

@Unanimous here's the notebook: gist.github.com/nunofachada/8cbe6ab7f856ae492a7587f3bbdc96f1 – Brader

Answer 9

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

Denominate

Answer 10

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)

Advantageous

Answer 11

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

Savannahsavant

Answer 12

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).

Milks

Answer 13

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}')

Minuet

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Barroom

Answer 14

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.

Mol

Answer 15

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))

Hypodermic

Nice succinct way of putting it, but you might want to fix a couple of typos: (a) second line should be _nv; (b) line 3 you should multiply the result by 2 – Caddoan

Answer 16

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:

Antoinetteanton