Encontrar el cuaternión que representa la rotación de un vector a otro

105

Tengo dos vectores uy v. ¿Hay alguna forma de encontrar un cuaternión que represente la rotación de u a v?

sdfqwerqaz1
fuente

Respuestas:

115
Quaternion q;
vector a = crossproduct(v1, v2);
q.xyz = a;
q.w = sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + dotproduct(v1, v2);

No olvide normalizar q.

Richard tiene razón en que no hay una rotación única, pero lo anterior debería dar el "arco más corto", que es probablemente lo que necesita.

Polaris878
fuente
30
Tenga en cuenta que esto no se aplica al caso de vectores paralelos (ambos en la misma dirección o apuntando en direcciones opuestas). crossproductno será válido en estos casos, por lo que primero debe verificar dot(v1, v2) > 0.999999y dot(v1, v2) < -0.999999, respectivamente, y devolver un quat de identidad para vectores paralelos, o devolver una rotación de 180 grados (sobre cualquier eje) para vectores opuestos.
sinisterchipmunk
11
Una buena implementación de esto se puede encontrar en el código fuente de ogre3d
João Portela
4
@sinisterchipmunk En realidad, si v1 = v2, el producto cruzado sería (0,0,0) yw sería positivo, lo que se normaliza a la identidad. Según gamedev.net/topic/… debería funcionar bien también para v1 = -v2 y en sus inmediaciones.
jpa
3
¿Cómo ha conseguido alguien que funcione esta técnica? Por un lado, se sqrt((v1.Length ^ 2) * (v2.Length ^ 2))simplifica a v1.Length * v2.Length. No pude obtener ninguna variación de esto para producir resultados sensibles.
Joseph Thomson
2
Si, esto funciona. Ver código fuente . L61 maneja si los vectores están orientados en direcciones opuestas (devuelve PI, de lo contrario, devolvería la identidad según el comentario de @ jpa). L67 maneja vectores paralelos: matemáticamente innecesario, pero más rápido. L72 es la respuesta de Polaris878, asumiendo que ambos vectores son unidades de longitud (evita una raíz cuadrada). Consulte también pruebas unitarias .
sinisterchipmunk
63

Solución vectorial a mitad de camino

Se me ocurrió la solución que creo que Imbrondir estaba tratando de presentar (aunque con un error menor, que probablemente fue la razón por la que sinisterchipmunk tuvo problemas para verificarlo).

Dado que podemos construir un cuaternión que represente una rotación alrededor de un eje así:

q.w == cos(angle / 2)
q.x == sin(angle / 2) * axis.x
q.y == sin(angle / 2) * axis.y
q.z == sin(angle / 2) * axis.z

Y que el punto y el producto cruzado de dos vectores normalizados son:

dot     == cos(theta)
cross.x == sin(theta) * perpendicular.x
cross.y == sin(theta) * perpendicular.y
cross.z == sin(theta) * perpendicular.z

Al ver que se puede lograr una rotación de u a v girando por theta (el ángulo entre los vectores) alrededor del vector perpendicular, parece que podemos construir directamente un cuaternión que represente dicha rotación a partir de los resultados del punto y los productos cruzados. ; sin embargo, tal como está, theta = ángulo / 2 , lo que significa que hacerlo daría como resultado el doble de la rotación deseada.

Una solución es calcular un vector a mitad de camino entre u y v , y usar el punto y el producto cruzado de u y el vector de mitad de camino para construir un cuaternión que represente una rotación del doble del ángulo entre u y el vector de mitad de camino , lo que nos lleva hasta v !

Hay un caso especial, donde u == -v y un vector único a mitad de camino se vuelve imposible de calcular. Esto es de esperar, dadas las infinitas rotaciones del "arco más corto" que pueden llevarnos de u a v , y simplemente debemos rotar 180 grados alrededor de cualquier vector ortogonal a u ( ov ) como nuestra solución de caso especial. Esto se hace tomando el producto cruzado normalizado de u con cualquier otro vector que no sea paralelo a u .

A continuación, se muestra un pseudocódigo (obviamente, en realidad, el caso especial tendría que tener en cuenta las inexactitudes del punto flotante, probablemente al comparar los productos punto con un umbral en lugar de un valor absoluto).

También tenga en cuenta que no hay un caso especial cuando u == v (se produce el cuaternión de identidad; compruébelo usted mismo).

// N.B. the arguments are _not_ axis and angle, but rather the
// raw scalar-vector components.
Quaternion(float w, Vector3 xyz);

Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
  // It is important that the inputs are of equal length when
  // calculating the half-way vector.
  u = normalized(u);
  v = normalized(v);

  // Unfortunately, we have to check for when u == -v, as u + v
  // in this case will be (0, 0, 0), which cannot be normalized.
  if (u == -v)
  {
    // 180 degree rotation around any orthogonal vector
    return Quaternion(0, normalized(orthogonal(u)));
  }

  Vector3 half = normalized(u + v);
  return Quaternion(dot(u, half), cross(u, half));
}

La orthogonalfunción devuelve cualquier vector ortogonal al vector dado. Esta implementación utiliza el producto cruzado con el vector de base más ortogonal.

Vector3 orthogonal(Vector3 v)
{
    float x = abs(v.x);
    float y = abs(v.y);
    float z = abs(v.z);

    Vector3 other = x < y ? (x < z ? X_AXIS : Z_AXIS) : (y < z ? Y_AXIS : Z_AXIS);
    return cross(v, other);
}

Solución de cuaternión a medio camino

Esta es en realidad la solución presentada en la respuesta aceptada, y parece ser marginalmente más rápida que la solución vectorial a mitad de camino (~ 20% más rápido según mis mediciones, aunque no confíe en mi palabra). Lo agrego aquí en caso de que otros como yo estén interesados ​​en una explicación.

Esencialmente, en lugar de calcular un cuaternión utilizando un vector de mitad de camino, puede calcular el cuaternión que da como resultado el doble de la rotación requerida (como se detalla en la otra solución) y encontrar el cuaternión a mitad de camino entre eso y cero grados.

Como expliqué antes, el cuaternión para duplicar la rotación requerida es:

q.w   == dot(u, v)
q.xyz == cross(u, v)

Y el cuaternión para la rotación cero es:

q.w   == 1
q.xyz == (0, 0, 0)

Calcular el cuaternión a mitad de camino es simplemente una cuestión de sumar los cuaterniones y normalizar el resultado, al igual que con los vectores. Sin embargo, como también ocurre con los vectores, los cuaterniones deben tener la misma magnitud, de lo contrario el resultado se desviará hacia el cuaternión con la magnitud mayor.

Un cuaternión construida a partir del punto y el producto vectorial de dos vectores tendrá la misma magnitud que los productos: length(u) * length(v). En lugar de dividir los cuatro componentes por este factor, podemos escalar el cuaternión de identidad. Y si se preguntaba por qué la respuesta aceptada aparentemente complica las cosas con el uso sqrt(length(u) ^ 2 * length(v) ^ 2), es porque la longitud al cuadrado de un vector es más rápida de calcular que la longitud, por lo que podemos guardar un sqrtcálculo. El resultado es:

q.w   = dot(u, v) + sqrt(length_2(u) * length_2(v))
q.xyz = cross(u, v)

Y luego normaliza el resultado. El pseudocódigo sigue:

Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
  float k_cos_theta = dot(u, v);
  float k = sqrt(length_2(u) * length_2(v));

  if (k_cos_theta / k == -1)
  {
    // 180 degree rotation around any orthogonal vector
    return Quaternion(0, normalized(orthogonal(u)));
  }

  return normalized(Quaternion(k_cos_theta + k, cross(u, v)));
}
Joseph Thomson
fuente
12
+1: ¡Genial! Esto funcionó como un encanto. Debería ser la respuesta aceptada.
Rekin
1
La sintaxis de Quaternion está activada en algunos ejemplos (Quaternion (xyz, w) y Quaternion (w, xyz)). También parece que en el último bloque de código se mezclan radianes y grados para expresar ángulos (180 vs. k_cos_theta + k).
Guillermo Blasco
1
Quaternion (float, Vector3) es una construcción a partir de un vector escalar, mientras que Quaternion (Vector3, float) es una construcción a partir de un eje-ángulo. Quizás potencialmente confuso, pero creo que es correcto. ¡Corrígeme si sigues pensando que está mal!
Joseph Thomson
¡Funcionó! ¡Gracias! Sin embargo, encontré otro enlace similar y bien explicado para realizar la operación anterior. Pensé que debería compartir para el registro;)
pecador
1
@JosephThomson La solución de cuaternión a medio camino parece provenir de aquí .
legends2k
6

El problema, como se dijo, no está bien definido: no existe una rotación única para un par de vectores dado. Considere el caso, por ejemplo, donde u = <1, 0, 0> y v = <0, 1, 0> . Una rotación de u a v sería una rotación pi / 2 alrededor del eje z. Otra rotación de u a v sería una rotación pi alrededor del vector <1, 1, 0> .

Richard Dunlap
fuente
1
De hecho, ¿no hay un número infinito de posibles respuestas? ¿Porque después de alinear el vector "desde" con el vector "hasta" todavía puede girar libremente el resultado alrededor de su eje? ¿Sabe qué información adicional se puede utilizar normalmente para limitar esta elección y definir bien el problema?
Doug McClean
5

¿Por qué no representar el vector usando cuaterniones puros? Quizás sea mejor si los normaliza primero.
q 1 = (0 u x u y u z ) '
q 2 = (0 v x v y v z )'
q 1 q rot = q 2
Pre-multiplique con q 1 -1
q rot = q 1 -1 q 2
donde q 1 -1 = q 1 conj / q norma
Esto se puede considerar como "división a la izquierda". La división derecha, que no es lo que quieres es:
q rot, right = q 2 -1 q 1

madratman
fuente
2
Estoy perdido, ¿no se calcula la rotación de q1 a q2 como q_2 = q_rot q_1 q_rot ^ -1?
yota
4

No soy muy bueno en Quaternion. Sin embargo, luché durante horas con esto y no pude hacer que la solución Polaris878 funcionara. Intenté pre-normalizar v1 y v2. Normalizando q. Normalizando q.xyz. Sin embargo, todavía no lo entiendo. El resultado todavía no me dio el resultado correcto.

Al final, aunque encontré una solución que lo hizo. Si ayuda a alguien más, aquí está mi código de trabajo (python):

def diffVectors(v1, v2):
    """ Get rotation Quaternion between 2 vectors """
    v1.normalize(), v2.normalize()
    v = v1+v2
    v.normalize()
    angle = v.dot(v2)
    axis = v.cross(v2)
    return Quaternion( angle, *axis )

Se debe hacer un caso especial si v1 y v2 son paralelos como v1 == v2 o v1 == -v2 (con cierta tolerancia), donde creo que las soluciones deberían ser Quaternion (1, 0,0,0) (sin rotación) o Quaternion (0, * v1) (rotación de 180 grados)

Imbrondir
fuente
Tengo una implementación que funciona, pero esta tuya es más bonita, así que realmente quería que funcionara. Desafortunadamente, falló en todos mis casos de prueba. Todas mis pruebas se parecen a eso quat = diffVectors(v1, v2); assert quat * v1 == v2.
sinisterchipmunk
Es poco probable que esto funcione, ya que angleobtiene su valor de un producto escalar.
sam hocevar
¿Dónde está la función Quaternion ()?
June Wang
3

Algunas de las respuestas no parecen considerar la posibilidad de que el producto cruzado pueda ser 0. A continuación, el fragmento utiliza la representación del eje del ángulo:

//v1, v2 are assumed to be normalized
Vector3 axis = v1.cross(v2);
if (axis == Vector3::Zero())
    axis = up();
else
    axis = axis.normalized();

return toQuaternion(axis, ang);

Se toQuaternionpuede implementar de la siguiente manera:

static Quaternion toQuaternion(const Vector3& axis, float angle)
{
    auto s = std::sin(angle / 2);
    auto u = axis.normalized();
    return Quaternion(std::cos(angle / 2), u.x() * s, u.y() * s, u.z() * s);
}

Si está utilizando la biblioteca Eigen, también puede hacer:

Quaternion::FromTwoVectors(from, to)
Shital Shah
fuente
toQuaternion(axis, ang)-> olvidó especificar qué esang
Maksym Ganenko
El segundo parámetro es el angleque forma parte de la representación eje-ángulo del cuaternión, medido en radianes.
Shital Shah
Se le pidió que hiciera que el cuaternión rotara de un vector a otro. No tienes ángulo, primero debes calcularlo. Tu respuesta debe contener el cálculo del ángulo. ¡Salud!
Maksym Ganenko
¿Esto es C ++? ¿Qué es ux ()?
June Wang
Sí, esto es C ++. u es el tipo de vector de la biblioteca Eigen (si está usando una).
Shital Shah
2

Desde el punto de vista del algoritmo, la solución más rápida se ve en pseudocódigo

 Quaternion shortest_arc(const vector3& v1, const vector3& v2 ) 
 {
     // input vectors NOT unit
     Quaternion q( cross(v1, v2), dot(v1, v2) );
     // reducing to half angle
     q.w += q.magnitude(); // 4 multiplication instead of 6 and more numerical stable

     // handling close to 180 degree case
     //... code skipped 

        return q.normalized(); // normalize if you need UNIT quaternion
 }

Asegúrese de que necesita cuaterniones de unidad (normalmente, es necesario para la interpolación).

NOTA: Los cuaterniones sin unidad se pueden usar con algunas operaciones más rápido que la unidad.

menor lógica
fuente