Á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].
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
defdotproduct(v1, v2):return sum((a*b) for a, b in zip(v1, v2))
deflength(v):return math.sqrt(dotproduct(v, v))
defangle(v1, v2):return math.acos(dotproduct(v1, v2) / (length(v1) * length(v2)))
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
defunit_vector(vector):""" Returns the unit vector of the vector. """return vector / np.linalg.norm(vector)
defangle_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))
¿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.
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
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
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:
defunit_vector(vector):""" Returns the unit vector of the vector"""return vector / np.linalg.norm(vector)
defangle(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.
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)defangle(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 = 1else:
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)defunit_vector(vector):""" Returns the unit vector of the vector. """return vector / np.linalg.norm(vector)
deftest_angle():defnpf(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)
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.0defangle_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° * π/180if line1_vertical and line2_vertical:
# Perpendicular vertical linesreturn0.0if 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))
defslope(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)
Respuestas:
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
fuente
math.sqrt(x)
es equivalentex**0.5
ymath.pow(x,y)
es equivalente ax**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.angle((1., 1., 1.), (1., 1., 1.))
). Vea mi respuesta para una versión un poco más correcta.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))
fuente
np.isnan
lugar 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.arccos
directa 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.0unit_vector
. Una posibilidad es simplemente devolver el vector de entrada en esta función cuando este sea el caso.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
fuente
nan
) cuando la dirección de los dos vectores sea idéntica u opuesta. Vea mi respuesta para una versión más correcta.angle = arccos(clip(c, -1, 1))
evitar problemas de redondeo. Esto resuelve el problema de @DavidWolever.clip
debe agregarse a la lista de importaciones numpy.La otra posibilidad es usar solo
numpy
y te da el ángulo interiorimport 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
fuente
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:
O calcule el ángulo con signo mediante proyección:
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.
fuente
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,
NotImplementedError
pero 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.fuente
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 NumbaY con
fuente
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
fuente
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).
fuente
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 elshapely
mó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
fuente