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

Question

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

- Peter

---

这是最好的答案，来自@MK83，因为它恰好是数学表达式theta = atan2(u^v, u.v)。即使在u=[0 0]或v=[0 0]的情况下也有涵盖，因为这是atan2产生NaN的唯一时间，在其他答案中，NaN将由/norm(u)或/norm(v)产生。 - PilouPili

Answer 1

注意：如果两个向量具有相同的方向（例如，(1, 0, 0)，(1, 0, 0)）或相反的方向（例如，(-1, 0, 0)，(1, 0, 0)），则这里的所有其他答案均会失败。

以下是一个可以正确处理这些情况的函数：

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

- David Wolever

---

我的numpy（版本==1.12.1）可以直接而安全地使用arccos函数。示例代码：In [140]: np.arccos(np.dot(np.array([1,0,0]),np.array([-1,0,0]) )) Out[140]: 3.1415926535897931In [141]: np.arccos(np.dot(np.array([1,0,0]),np.array([1,0,0]) )) Out[141]: 0.0 - ene

特殊情况是至少一个输入向量为零向量，这在“ unit_vector”中的除法中会出现问题。一种可能性是在此情况下，在此函数中仅返回输入向量。 - kafman

angle_between((0, 0, 0), (0, 1, 0)) 的结果将会是 nan，而不是 90 度。 - FabioSpaghetti

@kafman，0向量的角度在数学上是未定义的。因此，它引发错误是好事。 - user

即使您使用随机单位向量并且它们的乘积为1，@ene仍然可能返回nan。由于浮点精度问题，dot()函数将返回1.0000000000000002，即使是在版本1.16.1中也是如此。 - user

Answer 2

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

注意：当向量具有相同或相反方向时，此方法将失败。 正确的实现在这里：https://dev59.com/g3E85IYBdhLWcg3wVR-g#13849249

- Alex Martelli

---

另外，如果你只需要角度的cos、sin、tan值而不需要角度本身，那么你可以跳过math.acos获取余弦值，并使用叉积获取正弦值。 - mbeckish

鉴于math.sqrt(x)等价于x**0.5，math.pow(x,y)等价于x**y，我很惊讶它们在Python 2.x到3.0的过渡期间幸存了下来。在实践中，我通常将这些数值操作作为计算密集型进程的一部分进行处理，解释器对'**'直接转化为字节码BINARY_POWER的支持，相对于查找'math'，访问其属性'sqrt'，然后是痛苦缓慢的字节码CALL_FUNCTION，可以在没有编码或可读性成本的情况下显着提高速度。 - PaulMcG

就像在numpy的答案中一样：如果舍入误差产生作用，这种方法可能会失败！这可能会发生在平行和反平行的向量中！ - BandGap

注意：如果向量相同（例如，angle((1., 1., 1.), (1., 1., 1.))），则此方法将失败。请参见我的答案，其中包含稍微更正的版本。 - David Wolever

如果你在谈论上面的实现，那么它失败是因为舍入误差，而不是向量平行。 - Pace

Answer 3

使用numpy（强烈推荐）可以这样做：

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

- Olivier Verdier

---

由于存在舍入误差，最后一行可能会导致错误。因此，如果您执行 dot(u,u)/norm(u)**2，则结果为1.0000000002，然后arccos会失败（对于反平行向量也“有效”）。 - BandGap

我已经使用u=[1,1,1]进行了测试。u=[1,1,1,1]正常工作，但是每添加一个维度，返回的值都比1略大或略小... - BandGap

注意：当两个向量的方向相同时或相反时，这将失败（产生nan）。请参考我的答案获取更正确的版本。 - David Wolever

在Neo的评论中添加，最后一行应该是angle = arccos(clip(c, -1, 1))以避免舍入问题。这解决了@DavidWolever的问题。 - Tim Tisdall

对于使用以上代码片段的人：应该将 clip 添加到 numpy 导入列表中。 - Liam Deacon

Answer 4

另一种可能性是仅使用numpy，它会给出内角。

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)

这里是输出结果：

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

- MK83

---

这是最佳答案，因为它完全符合数学表达式 theta = atan2(u^v, u.v)。而且这种方法从不失败！ - PilouPili

这是针对二维的。原帖是在询问n维度。 - normanius

我使用了直角三角形作为示例： 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 预期的三角形角度应为30、60、90，而不是53。 - user2458922

@user2458922，一个3、4、5的直角三角形并不是一个30-60-90的三角形。30-60-90的三角形的边长分别为1、sqrt(3)和斜边为2。 - Roobie Nuby

@user2458922，一个3、4、5的直角三角形并不是一个30-60-90的三角形。30-60-90的三角形的边长分别为1、sqrt(3)和斜边为2。 - undefined

给定X，Y为[0,0]，[0,3]，[4,0]..如果形成一个角度为30、60、90的三角形。由于其中一个角度是90度，所以它是一个直角三角形。很抱歉，我不明白为什么“一个3,4,5的直角三角形不是一个30-60-90的三角形”？ - user2458922

Answer 5

寻找两个向量间夹角的简便方法（适用于n维向量）,

Python代码:

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

- Kevin Patel

Answer 6

如果你正在处理 3D 向量，你可以使用工具库 vg 来简洁地完成。它是建立在 numpy 上的轻量级扩展。

import numpy as np
import vg

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

vg.angle(vec1, vec2)

你还可以指定一个视角，通过投影计算角度：

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

或者通过投影计算有符号角度：

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

我在上一家创业公司创建了这个库，它的应用场景包括：将NumPy中冗长或难以理解的简单概念变得更加清晰易懂。

- paulmelnikow

Answer 7

David Wolever的解决方案非常好，但是：

如果你想要有带符号的角度，你需要确定给定的一对向量是左手系还是右手系（详见维基）。我的解决方案是：

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

虽然存在 NotImplementedError，但对于我的情况来说它已经很好用了。这种行为可以通过编写更多的代码来修复（因为手性是针对任何给定的一对确定的），但我并不想写那么多。

- sgt pepper

Answer 8

获取两个向量之间的角度的传统方法（即arccos(dot（u，v）/（norm（u）* norm（v））)，如其他答案所述）在几种极端情况下存在数值不稳定性。以下代码适用于n维和所有极端情况（它不检查长度为零的向量，但可以像其他答案中所示那样轻松添加）。请参见下面的注释。

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

这段代码是根据 Jeffrey Sarnoff 的 Julia 实现（使用 MIT 许可证）进行改编的，而 Jeffrey Sarnoff 的实现则基于 W. Kahan 教授的这些笔记（第15页）。

- faken

---

我尝试了一下，但是对于v1=np.array([0,0,-1])和v2=np.array([0,2,-4])，输出结果为0。然而，这两个向量之间的角度显然不是0。我有什么遗漏吗？ - zkytony

@zkytony，我测试了你建议的v1和v2的值，并且该函数返回了正确的结果。请参考这里：https://imgur.com/a/EKTnBny 也许你将错误的变量传递给了该函数？ - faken

@zkytony 这是笔记本链接：https://gist.github.com/nunofachada/8cbe6ab7f856ae492a7587f3bbdc96f1 - faken

Answer 9

对于那些可能因为SEO问题而到达此处，想要在Python中计算两条几何线段之间的角度，例如(x0, y0), (x1, y1)，下面是一个最小的解决方案（使用了shapely模块，但可以轻松修改为不使用）：

Original Answer翻译成"最初的回答"

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)

最初的回答

使用方法将是

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

- Julio Cezar Silva

Answer 10

在 Sgt Pepper 的优秀回答的基础上，增加了对齐向量的支持，并使用 Numba 实现了超过 2 倍的加速。

@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没有使用Numba的结果

• 平均每次循环耗时359微秒，标准偏差为2.86微秒（7次运行，每次运行1000次循环）

使用Numba的结果

• 平均每次循环耗时151微秒，标准偏差为820纳秒（7次运行，每次运行10000次循环）

- crizCraig