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

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

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.

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

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

¿Por qué no representar el vector usando cuaterniones puros? Quizás sea mejor si los normaliza primero.
q ₁ = (0 u _x u _y u _z ) '
q ₂ = (0 v _x v _y v _z )'
q ₁ q _rot = q ₂
Pre-multiplique con q ₁ ^-1
q _rot = q ₁ ^-1 q ₂
donde q ₁ ^-1 = q ₁ ^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 ₂ ^-1 q ₁

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)

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)

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.