Python 中两个 n 维向量之间的角度

我需要在 Python 中确定两个 n 维向量之间的角度。例如,输入可以是如下两个列表: [1,2,3,4][6,7,8,9]

220004 次浏览

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

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

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

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

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.

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

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.

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)

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

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

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