Znalezienie kwaternionu reprezentującego obrót z jednego wektora do drugiego

Mam dwa wektory u i v. Czy istnieje sposób na znalezienie kwaternionu reprezentującego obrót od u do v?

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

Nie zapomnij o normalizacji q.

Richard ma rację co do tego, że nie ma unikalnej rotacji, ale powyższe powinno dać „najkrótszy łuk”, co prawdopodobnie jest tym, czego potrzebujesz.

Rozwiązanie wektorowe w połowie drogi

Wymyśliłem rozwiązanie, które moim zdaniem Imbrondir próbował zaprezentować (aczkolwiek z drobnym błędem, który prawdopodobnie był powodem, dla którego sinisterchipmunk miał problem z weryfikacją).

Zakładając, że możemy skonstruować kwaternion reprezentujący obrót wokół osi w następujący sposób:

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

I że kropka i iloczyn poprzeczny dwóch znormalizowanych wektorów to:

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

Widząc, że obrót od u do v można osiągnąć obracając o theta (kąt między wektorami) wokół prostopadłego wektora, wygląda na to, że możemy bezpośrednio skonstruować kwaternion reprezentujący taki obrót z wyników kropki i iloczynów poprzecznych ; jednakże w obecnej postaci theta = kąt / 2 , co oznacza, że ​​spowodowałoby to dwukrotność pożądanego obrotu.

Jednym z rozwiązań jest obliczenie wektora w połowie drogi między u i v i użycie kropki i iloczynu poprzecznego u oraz wektora połowy drogi do skonstruowania kwaternionu reprezentującego obrót o dwukrotności kąta między u a wektorem w połowie drogi , co prowadzi nas aż do v !

Istnieje specjalny przypadek, w którym u == -v i unikatowy wektor połówkowy staje się niemożliwy do obliczenia. Jest to oczekiwane, biorąc pod uwagę nieskończenie wiele obrotów „najkrótszego łuku”, które mogą prowadzić nas od u do v , i musimy po prostu obrócić o 180 stopni wokół dowolnego wektora prostopadłego do u (lub v ), co jest naszym rozwiązaniem w przypadku szczególnym. Odbywa się to poprzez znormalizowany iloczynu U z dowolnym innym wektorem nie równolegle do u .

Następuje pseudokod (oczywiście w rzeczywistości szczególny przypadek musiałby uwzględniać niedokładności zmiennoprzecinkowe - prawdopodobnie poprzez sprawdzanie iloczynów skalarnych względem jakiegoś progu, a nie wartości bezwzględnej).

Zauważ również, że nie ma specjalnego przypadku, gdy u == v (tworzony jest kwaternion tożsamości - sprawdź i przekonaj się sam).

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

orthogonalZwraca dowolny wektor prostopadły do danego wektora. Ta implementacja używa iloczynu krzyżowego z najbardziej ortogonalnym wektorem bazowym.

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

Rozwiązanie Half-Way Quaternion

W rzeczywistości jest to rozwiązanie przedstawione w przyjętej odpowiedzi i wydaje się być nieznacznie szybsze niż rozwiązanie wektorowe w połowie drogi (~ 20% szybciej według moich pomiarów, choć nie wierz mi na słowo). Dodam to tutaj na wypadek, gdyby inni tacy jak ja byli zainteresowani wyjaśnieniem.

Zasadniczo, zamiast obliczać kwaternion przy użyciu wektora połówkowego, można obliczyć kwaternion, która daje dwukrotność wymaganego obrotu (jak opisano w drugim rozwiązaniu) i znaleźć kwaternion w połowie odległości między tym a zerem stopni.

Jak wyjaśniłem wcześniej, kwaternion dla podwójnej wymaganej rotacji wynosi:

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

A kwaternion dla rotacji zerowej to:

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

Obliczanie kwaternionów w połowie drogi polega po prostu na sumowaniu kwaternionów i normalizacji wyniku, podobnie jak w przypadku wektorów. Jednak, podobnie jak w przypadku wektorów, kwaterniony muszą mieć tę samą wielkość, w przeciwnym razie wynik będzie pochylony w kierunku kwaternionu o większej wielkości.

Kwaternion zbudowana z punktów i dwóch wektorów będzie miał taką samą wielkość jak w przypadku produktów: length(u) * length(v). Zamiast dzielić wszystkie cztery składniki przez ten czynnik, możemy zamiast tego skalować w górę kwaternion tożsamości. A jeśli zastanawiałeś się, dlaczego zaakceptowana odpowiedź pozornie komplikuje sprawę przy użyciu sqrt(length(u) ^ 2 * length(v) ^ 2), to dlatego, że długość kwadratu wektora jest szybsza do obliczenia niż długość, więc możemy zapisać jedno sqrtobliczenie. Wynik to:

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

A następnie znormalizuj wynik. Pseudokod następujący:

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

Przedstawiony problem nie jest dobrze zdefiniowany: nie ma unikalnej rotacji dla danej pary wektorów. Rozważmy na przykład przypadek, w którym u = <1, 0, 0> i v = <0, 1, 0> . Jeden obrót od u do v byłby obrotem pi / 2 wokół osi z. Kolejny obrót od u do v oznaczałby obrót pi wokół wektora <1, 1, 0> .

Dlaczego nie przedstawić wektora za pomocą czystych kwaternionów? Może lepiej, jeśli najpierw je znormalizujesz.
q ₁ = (0 u _x u _y u _z ) '
q ₂ = (0 v _x v _y v _z )'
q ₁ q _rot = q ₂
Pomnóż wstępnie przez q ₁ ^-1
q _rot = q ₁ ^-1 q ₂
gdzie q ₁ ^-1 = q ₁ ^spój / q _norm
Można to traktować jako „lewy podział”. Dzielenie od prawej, które nie jest tym, czego chcesz, to:
q _{rot, right} = q ₂ ^-1 q ₁

Nie jestem dobry w Quaternion. Jednak walczyłem z tym godzinami i nie mogłem sprawić, by rozwiązanie Polaris878 działało. Próbowałem przed normalizacją v1 i v2. Normalizowanie q. Normalizowanie q.xyz. Jednak nadal nie rozumiem. Wynik nadal nie dał mi właściwego wyniku.

W końcu jednak znalazłem rozwiązanie, które to zrobiło. Jeśli to pomoże komuś innemu, oto mój działający (w Pythonie) kod:

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 )

Należy zrobić specjalny przypadek, jeśli v1 i v2 są równoległe, jak v1 == v2 lub v1 == -v2 (z pewną tolerancją), gdzie uważam, że rozwiązaniami powinny być Quaternion (1, 0,0,0) (bez rotacji) lub Quaternion (0, * v1) (obrót o 180 stopni)

Wydaje się, że niektóre odpowiedzi nie uwzględniają możliwości, że iloczyn poprzeczny może wynosić 0. Poniższy fragment wykorzystuje reprezentację kąta-osi:

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

toQuaternionMogą być realizowane w następujący sposób:

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

Jeśli korzystasz z biblioteki Eigen, możesz również po prostu:

Quaternion::FromTwoVectors(from, to)

Z punktu widzenia algorytmu najszybsze rozwiązanie wygląda w pseudokodzie

 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
 }

Upewnij się, że potrzebujesz kwaternionów jednostkowych (zwykle jest to wymagane do interpolacji).

UWAGA: Quaternions inne niż jednostki mogą być używane z niektórymi operacjami szybciej niż jednostka.