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

Question

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

Answer 1

使用 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 2

使用 numpy 并处理 BandGap 的舍入误差：

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)

请注意，如果其中一个向量的大小为零（除以 0），此函数将引发异常。

Answer 3

import math

ax, ay = input('输入向量 a 的 x 和 y: ').split()

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

bx, by = input('输入向量 b 的 x 和 y: ' ).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 4

OP 要求 n 维 n>2，但很多人最终都会遇到二维问题，所以我想澄清这种特殊情况的最佳解决方案。

所提出的解决方案都使用反余弦函数（MK83 和 faken 除外），因此都非常不准确，并且对于接近 0 或 180 度的角度容易出错。这是因为角度的巨大变化导致这些值下角度的余弦几乎没有变化。

arccos（对于二维情况）的第二个问题是它无法区分正角和负角。因此，(0,1) 和 (1,0) 之间的角度将与 (1,0) 和 (0,1) 之间的角度相同，尽管第一个角度应为 90 度，第二个角度应为 -90 度。

faken 对 OP 多维问题有一个很好的答案，避免使用 arccos，因此在整个范围内都是准确的。

MK83 使用 atan2 解决二维问题，atan2 是为此问题提供的函数。答案范围是 -180 度到 180 度。我这里提出一个只针对二维的解决方案，比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])

这比MK83的解决方案快大约三倍。

Answer 5

给定两个向量 vec(u) 和 vec(v)，则可以证明以下公式是数值上最稳定的解：

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

该公式相对于任何向量的优点of:

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

are:

1. atan 对于小角度来说数值稳定。对于小角度，acos 函数不像 cos(th) ~ 1 - 0.5*th^2 那样，
2. 对于围绕 Pi 的角度，atan 在数值上是稳定的/2。当公式计算半角时，它对应于使用叉积的公式中的角度 Pi。在这些条件下，叉积计算不稳定，因此使用 asin 的公式也不稳定。

人们可以对向量 u 和 v 进行预归一化并使用以下形式：

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

然而，除法实际上会损失数值精度。在这种情况下，这种精度损失将出现在标准化和最终除法中。

因此，在实现方面，我们应该使用经典的 atan2 函数来避免再次除法。因此，我们获得了最佳的数值稳定解。

因此，使用 numpy，实现很简单：

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

这个问题有很多解决方案，但没有一个采用使用numpy的简单表达式能够处理跨任何轴的n维数组，如下所示：

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

---

示例< /strong>

我们在网格上选择点，并在每个点绘制两个向量：一个为黑色，另一个为灰色。然后，我们计算黑色和灰色向量之间的角度，并根据该角度绘制着色的点。

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

输出：

Answer 7

注意：如果两个向量具有相同的方向（例如，(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))

Answer 8

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://stackoverflow.com/a/13849249/71522

Answer 9

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

Answer 10

另一种可能性是仅使用 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

Answer 11

如果您正在使用 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 中冗长或不透明的简单想法。

Answer 12

求两个向量之间角度的简单方法（适用于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

Answer 13

David Wolever 的解决方案很好，但是

如果你想要有符号的角度，你必须确定给定的角度是右手还是左手（参见wiki 了解更多信息）。

我的解决方案是：

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，它并不完美，但对于我的情况来说，它效果很好。这种行为可以修复（因为任何给定的配对都确定了手感），但它需要我想要并且必须编写的更多代码。

Answer 14

获取两个向量之间角度的传统方法（即 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 许可），又基于 < a href="https://people.eecs.berkeley.edu/%7Ewkahan/MathH110/Cross.pdf" rel="noreferrer">这些笔记，作者：W. Kahan 教授（第 15 页）。

Answer 15

对于少数可能（由于 SEO 复杂性）尝试在 python 中计算两条线之间的角度的人，如 (x0, y0), (x1, y1) 几何线，有以下最小解决方案（使用shapely模块，但可以轻松修改不）：

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

Answer 16

以 sgtpepper 的出色答案为基础，添加对对齐向量的支持，并使用 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 µs ± 2.86 µs（7 次运行的平均值 ± 标准偏差，每次 1000 个循环）

，而使用

• 151 µs ± 820 ns > 每个循环（7 次运行的平均值 ± 标准差，每次 10000 次循环）