Forma directa de calcular el ángulo en el sentido de las agujas del reloj entre dos vectores

Quiero encontrar el ángulo en el sentido de las agujas del reloj entre dos vectores (bidimensionales o tridimensionales).

La forma clásica con el producto escalar me da el ángulo interno (0-180 grados), y necesito usar algunos si declaraciones para determinar si el resultado es el ángulo que necesito o su complemento.

¿Hay una forma directa de calcular el ángulo en el sentido de las agujas del reloj?

caso 2D

Al igual que el producto punto es proporcional al coseno del ángulo, el determinante es proporcional a su seno. Entonces puedes calcular el ángulo así:

dot = x1*x2 + y1*y2      # Dot product between [x1, y1] and [x2, y2]
det = x1*y2 - y1*x2      # Determinant
angle = atan2(det, dot)  # atan2(y, x) or atan2(sin, cos)

La orientación de este ángulo coincide con la del sistema de coordenadas. en un sistema de coordenadas para zurdoses decir X apuntando a la derecha y y hacia abajo como es común en los gráficos por computadora, esto significará que obtendrá un signo positivo para los ángulos en el sentido de las agujas del reloj. Si la orientación del sistema de coordenadas es matemática con y hacia arriba, obtienes ángulos en sentido contrario a las manecillas del reloj como es la convención en matemáticas. Cambiar el orden de las entradas cambiará el signo, por lo que si no está satisfecho con los signos, simplemente cambie las entradas.

caso 3D

En 3D, dos vectores colocados arbitrariamente definen su propio eje de rotación, perpendicular a ambos. Ese eje de rotación no viene con una orientación fija, lo que significa que tampoco puede fijar de manera única la dirección del ángulo de rotación. Una convención común es dejar que los ángulos sean siempre positivos y orientar el eje de tal manera que encaje en un ángulo positivo. En este caso, el producto escalar de los vectores normalizados es suficiente para calcular los ángulos.

dot = x1*x2 + y1*y2 + z1*z2    # Between [x1, y1, z1] and [x2, y2, z2]
lenSq1 = x1*x1 + y1*y1 + z1*z1
lenSq2 = x2*x2 + y2*y2 + z2*z2
angle = acos(dot/sqrt(lenSq1 * lenSq2))

Tenga en cuenta que algunos comentarios y respuestas alternativas desaconsejan el uso de acos por razones numéricas, en particular si los ángulos a medir son pequeños.

Plano incrustado en 3D

Un caso especial es el caso en el que sus vectores no se colocan arbitrariamente, sino que se encuentran dentro de un plano con un vector normal conocido norte. Entonces el eje de rotación estará en la dirección norte también, y la orientación de norte fijará una orientación para ese eje. En este caso, puede adaptar el cálculo 2D anterior, incluido norte en el determinante para hacer su tamaño 3×3.

dot = x1*x2 + y1*y2 + z1*z2
det = x1*y2*zn + x2*yn*z1 + xn*y1*z2 - z1*y2*xn - z2*yn*x1 - zn*y1*x2
angle = atan2(det, dot)

Una condición para que esto funcione es que el vector normal norte tiene unidad de longitud. Si no, tendrás que normalizarlo.

Como triple producto

Este determinante también podría expresarse como el triple productocomo señaló @Excrubulent en una edición sugerida.

det = n · (v1 × v2)

Esto podría ser más fácil de implementar en algunas API y brinda una perspectiva diferente de lo que está sucediendo aquí: el producto cruzado es proporcional al seno del ángulo y estará perpendicular al plano, por lo tanto, será un múltiplo de norte. Por lo tanto, el producto punto básicamente medirá la longitud de ese vector, pero con el signo correcto adjunto.

Para calcular el ángulo solo necesitas llamar atan2(v1.s_cross(v2), v1.dot(v2)) para el caso bidimensional. Dónde s_cross es el análogo escalar de la producción cruzada (área firmada del paralelogramo).

Para el caso bidimensional, sería producción de cuña.

Para el caso tridimensional, debe definir la rotación en el sentido de las agujas del reloj, porque desde un lado del plano en el sentido de las agujas del reloj hay una dirección, desde el otro lado del plano hay otra dirección =)

Este es el ángulo en sentido contrario a las agujas del reloj, y el ángulo en el sentido de las agujas del reloj es justo el opuesto.

Esta respuesta es la misma que la de MvG, pero la explica de manera diferente (es el resultado de mis esfuerzos para tratar de entender por qué funciona la solución de MvG).

El ángulo antihorario theta de x a ycon respecto al punto de vista de su dado normal n (||n|| = 1), es dado por

atan2( punto(n, cruz(x,y)), punto(x,y) )

(1) = atan2( ||x|| ||y|| sen(theta), ||x|| ||y|| cos(theta) )

(2) = atan2( sin(theta), cos(theta) )

(3) = ángulo en sentido antihorario entre el eje x y el vector (cos(theta), sin(theta))

(4) = theta

dónde ||x|| denota la magnitud de x.

El paso (1) sigue al señalar que

cruz(x,y) = ||x|| ||y|| pecado (theta) n,

y entonces

punto(n, cruz(x,y))

= punto(n, ||x|| ||y|| sin(theta) n)

= ||x|| ||y|| pecado(theta) punto(n, n)

que es igual

||x|| ||y|| pecado (theta)

si ||n|| = 1.

El paso (2) se deriva de la definición de atan2señalando que atan2(cy, cx) = atan2(y,x)dónde c es un escalar. El paso (3) se deriva de la definición de atan2. El paso (4) se sigue de las definiciones geométricas de cos y sin.

Dado que una de las soluciones más simples y elegantes está oculta en uno de los comentarios, creo que podría ser útil publicarla como una respuesta separada.

acos puede causar imprecisiones para ángulos muy pequeños, por lo que atan2 generalmente se prefiere. Para el caso tridimensional:

dot = x1*x2 + y1*y2 + z1*z2
cross_x = (y1*z2 – z1*y2)
cross_y = (z1*x2 – x1*z2) 
cross_z = (x1*y2 – y1*x2)
det = sqrt(cross_x*cross_x + cross_y*cross_y + cross_z*cross_z)
angle = atan2(det, dot)

El producto escalar (punto) de dos vectores te permite obtener el coseno del ángulo entre ellos.

Para obtener la ‘dirección’ del ángulo, también debe calcular el producto vectorial. Le permitirá verificar (a través de la coordenada z) si el ángulo es en el sentido de las agujas del reloj o no (es decir, si lo extrae de 360 ​​grados o no).

Para un método bidimensional, podrías usar la ley de los cosenos y el método de “dirección”.

Para calcular el ángulo del segmento P3:P1 barriendo en el sentido de las agujas del reloj hasta el segmento P3:P2.

    P1     P2

        P3

    double d = direction(x3, y3, x2, y2, x1, y1);

    // c
    int d1d3 = distanceSqEucl(x1, y1, x3, y3);

    // b
    int d2d3 = distanceSqEucl(x2, y2, x3, y3);

    // a
    int d1d2 = distanceSqEucl(x1, y1, x2, y2);

    //cosine A = (b^2 + c^2 - a^2)/2bc
    double cosA = (d1d3 + d2d3 - d1d2)
        / (2 * Math.sqrt(d1d3 * d2d3));

    double angleA = Math.acos(cosA);

    if (d > 0) {
        angleA = 2.*Math.PI - angleA;
    }

Esto tiene el mismo número de operaciones trascendentales que las sugerencias anteriores y solo una operación de coma flotante más o menos.

Los métodos que utiliza son:

 public int distanceSqEucl(int x1, int y1,
    int x2, int y2) {

    int diffX = x1 - x2;
    int diffY = y1 - y2;
    return (diffX * diffX + diffY * diffY);
}

public int direction(int x1, int y1, int x2, int y2,
    int x3, int y3) {

    int d = ((x2 - x1)*(y3 - y1)) - ((y2 - y1)*(x3 - x1));

    return d;
}

Si por “vía directa” quiere decir evitar la if declaración, entonces no creo que haya una solución realmente general.

Sin embargo, si su problema específico le permitiría perder algo de precisión en la discretización de ángulos y está de acuerdo con perder algo de tiempo en las conversiones de tipo, puede mapear el [-pi,pi] rango permitido de ángulo phi en el rango permitido de algún tipo de entero con signo. Entonces obtendrías la complementariedad gratis. Sin embargo, realmente no usé este truco en la práctica. Lo más probable es que el costo de las conversiones de flotante a entero y de entero a flotante supere cualquier beneficio de la franqueza. Es mejor establecer sus prioridades en la escritura de código autovectorizable o paralelizable cuando este cálculo de ángulo se realiza mucho.

Además, si los detalles de su problema son tales que hay un resultado definitivo más probable para la dirección del ángulo, entonces puede usar las funciones integradas de los compiladores para proporcionar esta información al compilador, de modo que pueda optimizar la bifurcación de manera más eficiente. ej., en caso de CCGeso es __builtin_expect función. Es algo más útil de usar cuando lo envuelves en tal likely y unlikely macros (como en el kernel de Linux):

#define likely(x)      __builtin_expect(!!(x), 1)
#define unlikely(x)    __builtin_expect(!!(x), 0)