Ángulos entre dos vectores n-dimensionales en Python

82

Necesito determinar los ángulos entre dos vectores n-dimensionales en Python. Por ejemplo, la entrada puede ser dos listas como las siguientes: [1,2,3,4]y [6,7,8,9].

Pedro
fuente
1
Esta es la mejor respuesta es @ MK83, ya que es exactamente la expresión matemática theta = atan2 (u ^ v, uv). incluso el caso donde u = [0 0] o v = [0 0] está cubierto porque este es el único momento en que atan2 producirá el NaN en las otras respuestas, NaN será producido por / norm (u) o / norm (v)
PilouPili

Respuestas:

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

Nota : esto fallará cuando los vectores tengan la misma dirección o la opuesta. La implementación correcta está aquí: https://stackoverflow.com/a/13849249/71522

Alex Martelli
fuente
2
Además, si solo necesita cos, sin, tan del ángulo, y no el ángulo en sí, puede omitir math.acos para obtener el coseno y usar el producto cruzado para obtener el seno.
mbeckish
9
Dado que math.sqrt(x)es equivalente x**0.5y math.pow(x,y)es equivalente a x**y, me sorprende que hayan sobrevivido al hacha de redundancia ejercida durante la transición de Python 2.x-> 3.0. En la práctica, normalmente hago este tipo de cosas numéricas como parte de un proceso de cálculo intensivo más grande, y el soporte del intérprete para que '**' vaya directamente al código de bytes BINARY_POWER, frente a la búsqueda de 'matemáticas', el acceso a su atributo 'sqrt', y luego el código de bytes dolorosamente lento CALL_FUNCTION, puede hacer una mejora medible en la velocidad sin costo de codificación o legibilidad.
PaulMcG
5
Como en la respuesta con numpy: ¡esto puede fallar si entra en juego el error de redondeo! ¡Esto puede suceder con vectores paralelos y antiparalelos!
BandGap
2
Nota: esto fallará si los vectores son idénticos (por ejemplo, angle((1., 1., 1.), (1., 1., 1.))). Vea mi respuesta para una versión un poco más correcta.
David Wolever
2
Si está hablando de la implementación anterior, falla debido a errores de redondeo, no porque los vectores sean paralelos.
Ritmo
150

Nota : todas las otras respuestas aquí se producirá un error si los dos vectores tienen ya sea la misma dirección (por ejemplo, (1, 0, 0), (1, 0, 0)) o direcciones opuestas (por ejemplo, (-1, 0, 0), (1, 0, 0)).

Aquí hay una función que manejará correctamente estos casos:

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
fuente
¿No sería mejor usar en np.isnanlugar del de la biblioteca de matemáticas? En teoría deberían ser idénticos, pero no estoy muy seguro en la práctica. De cualquier manera, me imagino que sería más seguro.
Enganchado el
2
My numpy (versión == 1.12.1) se puede utilizar de forma arccosdirecta y segura. : En [140]: np.arccos (np.dot (np.array ([1,0,0]), np.array ([- 1,0,0]))) Fuera [140]: 3.1415926535897931 En [ 141]: np.arccos (np.dot (np.array ([1,0,0]), np.array ([1,0,0]))) Fuera [141]: 0.0
ene
2
Se omite el caso especial en el que al menos un vector de entrada es el vector cero, lo que es problemático para la división en unit_vector. Una posibilidad es simplemente devolver el vector de entrada en esta función cuando este sea el caso.
kafman
1
angle_between ((0, 0, 0), (0, 1, 0)) dará nan como resultado, y no 90
FabioSpaghetti
2
El ángulo de @kafman 0-vectors no está definido (en matemáticas). Entonces, el hecho de que genere un error es bueno.
usuario
45

Usando numpy (muy recomendado), harías:

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
fuente
3
La última línea puede resultar en un error como descubrí debido a errores de redondeo. Así, si usted a punto (u, u) / norma (u) ** 2 da lugar a 1,0000000002 y los arccos luego falla (también obras '' para antiparalelas vectores)
banda prohibida
He probado con u = [1,1,1]. u = [1,1,1,1] funciona bien, pero cada dimensión agregada devuelve valores ligeramente mayores o menores que 1 ...
BandGap
3
Nota: esto fallará (cederá nan) cuando la dirección de los dos vectores sea idéntica u opuesta. Vea mi respuesta para una versión más correcta.
David Wolever
2
agregando el comentario de neo a esto, la última línea debería ser angle = arccos(clip(c, -1, 1))evitar problemas de redondeo. Esto resuelve el problema de @DavidWolever.
Tim Tisdall
4
Para las personas que usan el fragmento de código anterior: clipdebe agregarse a la lista de importaciones numpy.
Liam Deacon
26

La otra posibilidad es usar solo numpyy te da el ángulo interior

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)

y aquí está la salida:

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
fuente
6
Esta es la mejor respuesta ya que es exactamente la expresión matemática theta = atan2 (u ^ v, uv). ¡Y esto nunca falla!
PilouPili
1
Esto es para 2-D. El OP estaba pidiendo nD
normanius
3

Si está trabajando con vectores 3D, puede hacerlo de forma concisa utilizando el toolbelt vg . Es una capa ligera encima de numpy.

import numpy as np
import vg

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

vg.angle(vec1, vec2)

También puede especificar un ángulo de visión para calcular el ángulo a través de la proyección:

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

O calcule el ángulo con signo mediante proyección:

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

Creé la biblioteca en mi última puesta en marcha, donde fue motivada por usos como este: ideas simples que son detalladas u opacas en NumPy.

Paulmelnikow
fuente
3

La solución de David Wolever es buena, pero

Si desea tener ángulos firmados , debe determinar si un par dado es diestro o zurdo (consulte la wiki para obtener más información).

Mi solución para esto es:

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

No es perfecto por esto, NotImplementedErrorpero en mi caso funciona bien. Este comportamiento podría arreglarse (porque la manipulación se determina para cualquier par dado) pero se necesita más código del que quiero y tengo que escribir.

sgt pepper
fuente
1

Sobre la base de la gran respuesta de sgt pepper y agregando soporte para vectores alineados, además de agregar una aceleración de más de 2x usando 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 resultados sin Numba

  • 359 µs ± 2,86 µs por bucle (media ± desviación estándar de 7 corridas, 1000 bucles cada una)

Y con

  • 151 µs ± 820 ns por bucle (media ± desviación estándar de 7 corridas, 10000 bucles cada una)
crizCraig
fuente
1

Manera fácil de encontrar el ángulo entre dos vectores (funciona para el vector n-dimensional),

Código 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
fuente
0

Usando numpy y ocupándose de los errores de redondeo de 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)

Tenga en cuenta que esta función producirá una excepción si uno de los vectores tiene magnitud cero (dividir por 0).

Paso
fuente
0

Para los pocos que pueden haber terminado aquí (debido a complicaciones de SEO) tratando de calcular el ángulo entre dos líneas en Python, como en (x0, y0), (x1, y1)las líneas geométricas, existe la siguiente solución mínima (usa el shapelymódulo, pero se puede modificar fácilmente para no hacerlo):

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)

Y el uso seria

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