Strona główna » Artykuły » Teoria » Matematyka »

Kwaterniony w praktyce

Wersja do druku

Parametry:
Dodany: 2008-05-20 przez !Reg
Autor: Adam Sawicki
Opis: Dla tych którzy już znają wektory i macierze, artykuł wprowadza kwaterniony jako sposób reprezentowania obrotów i orientacji.
Poziom:
Ocena: 132

Wstęp

Artykuł poświęcony jest w całości kwaternionom. To obiekty matematyczne, które służą programistom grafiki i gier do reprezentowania rotacji (obrotów) i orientacji w przestrzeni 3D. Aby zrozumieć ten artykuł, musisz już znać podstawy matematyki 3D - wektory i macierze. Wiedzę na temat tych zagadnień można zdobyć np. z [1]. Jeśli spełniasz ten warunek, a wcześniej omijałeś temat kwaternionów myśląc że to coś bardzo trudnego, dzięki temu artykułowi (mam nadzieję) przekonasz się, że nie ma w nich żadnej czarnej magii.

W tekście tym pokażę, jak można samodzielnie zaimplementować strukturę kwaterniona i funkcje do operacji na kwaternionach. Programiści używający DirectX mają ułatwione zadanie, gdyż taką strukturę i takie funkcje mają już dostępne w bibliotece D3DX (struktura D3DXQUATERNION) - muszą tylko zrozumieć, jak ich używać. Mimo tego pokażę wewnętrzną implementację wszystkich tych obliczeń. Używam tu języka C++, biblioteki Direct3D i lewoskrętnego układu współrzędnych, ale użytkownicy OpenGL również mogą skorzystać z tego artykułu.

Wprowadzenie

Żeby od razu oswoić się z kwaternionem, zobaczmy jak ten potwór wygląda:

struct QUATERNION
{
  float x, y, z, w;

  QUATERNION() { }
  QUATERNION(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) { }
};

Jest to jak widać struktura zbudowana z czterech liczb, nazywanych x, y, z, w. Mimo podobnej budowy nie należy go jednak mylić z wektorem czterowymiarowym (w D3DX to struktura D3DXVECTOR4) używanym do reprezentowania punktów i wektorów 3D we współrzędnych jednorodnych (ang. homogeneous coordinates). Kwaternion służy do zupełnie czego innego i obowiązują go inne zasady wykonywania działań. Mówi się też czasami, że składowe x, y, z to część wektorowa kwaterniona (choć nie jest to zwyczajny wektor), a w to jego część skalarna.

Czym wobec tego jest kwaternion? Teoretycznie jest to rozszerzenie liczb zespolonych na cztery wymiary i posiada aż trzy pierwiastki urojone. Brzmi strasznie? Owszem, ale na szczęście nie musisz tego rozumieć. Nie musisz nawet znać liczb zespolonych, żeby używać kwaternionów. Takie sprawy zostawmy matematykom. Dla nas kwaternion jest obiektem, który reprezentuje rotację (obrót) lub orientację w przestrzeni 3D. Zobacz rys. 1.

Rys. 1. Imbryczek obrócony do różnej orientacji w przestrzeni 3D.

Czym się różni rotacja (obrót) od orientacji? Tym samym czym translacja (przesunięcie) od pozycji. Pozycja to miejsce w którym się znajdujesz - np. Warszawa (punkt o podanych współrzędnych, np. szerokości i długości geograficznej). Translacja nastąpi, jeśli się gdzieś przemieścisz - np. z Warszawy do Londynu (wektor iluś tam kilometrów na zachód i iluś tam na północ). Podobnie twoja orientacja to ułożenie twojego ciała w danej chwili (np. stoisz, leżysz). Rotacja nastąpi, jeśli się przewrócisz (o 90 stopni wokół osi X :)

Do reprezentowania zarówno pozycji (punktu), jak i translacji (wektora) używamy w programowaniu tej samej struktury - np. D3DXVECTOR3. Podobnie do reprezentowania zarówno orientacji, jak i rotacji używamy kwaterniona. Punkt możemy traktować jako wektor przesuwający od początku układu współrzędnych do swojego końca i na odwrót. Tak więc różnica jest bardzo subtelna i występuje tylko w naszych myślach. Warto jednak zdawać sobie z niej sprawę.

Tworzenie kwaterniona

Kwaternion tworzy się podając wektor, którego kierunek wskazuje oś obrotu oraz kąt, o jaki chcemy obracać wokół tego wektora (nie muszę chyba przypominać, że kąty podajemy w radianach?). Zobacz rys. 2. Który kierunek jest dodatni - to zależy od skrętności układu współrzędnych. Jeśli układ jest lewoskrętny, weź swoją lewą rękę, wskaż kciukiem kierunek wektora osi obrotu, a kiedy zegniesz pozostałe palce, wskażą one dodatni kierunek obrotu.

Rys. 2. Kwaternion powstaje przez podanie wektora i kąta obrotu wokół niego.

Tych informacji nie wpisujemy jednak do składowych kwaterniona bezpośrednio. Trzeba je zakodować według algorytmu jak na poniższym listingu, obliczając najpierw sinus i cosinus połowy podanego kąta.

void AxisToQuaternion(QUATERNION *Out, const VEC3 &Axis, float Angle)
{
  Angle *= 0.5f;
  float Sin = sinf(Angle);
  Out->x = Sin * Axis.x;
  Out->y = Sin * Axis.y;
  Out->z = Sin * Axis.z;
  Out->w = cosf(Angle);
}

Przypadkiem szczególnym jest obracanie wokół osi X, Y, Z. Algorytm znacznie się wówczas upraszcza i dla optymalizacji warto przygotować osobne funkcje:

void QuaternionRotationX(QUATERNION *Out, float a)
{
  a *= 0.5f;
  Out->x = sinf(a); Out->y = 0.0f; Out->z = 0.0f;
  Out->w = cosf(a);
}

void QuaternionRotationY(QUATERNION *Out, float a)
{
  a *= 0.5f;
  Out->x = 0.0f; Out->y = sinf(a); Out->z = 0.0f;
  Out->w = cosf(a);
}

void QuaternionRotationZ(QUATERNION *Out, float a)
{
  a *= 0.5f;
  Out->x = 0.0f; Out->y = 0.0f; Out->z = sinf(a);
  Out->w = cosf(a);
}

Można też napisać funkcje, które robią rzecz odwrotną - wyciągają z kwaterniona oś obrotu i kąt, z których został utworzony:

void CalcQuaternionAxis(VEC3 *Out, const QUATERNION &Q)
{
  float SinThetaOver25q = 1.0f - Q.w * Q.w;
  if (SinThetaOver25q <= 0.0f)
  {
    Out->x = Out->y = Out->z = 0.0f;
    return;
  }
  float OneOverSinThetaOver2 = 1.0f / sqrtf(SinThetaOver25q);
  Out->x = Q.x * OneOverSinThetaOver2;
  Out->y = Q.y * OneOverSinThetaOver2;
  Out->z = Q.z * OneOverSinThetaOver2;
}

float CalcQuaternionAngle(const QUATERNION &Q)
{
  return safe_acos(Q.w) * 2.0f;
}

// Funkcja pomocnicza
float safe_acos(float x)
{
  if (x <= -1.0f) return PI;
  if (x >= 1.0f) return 0.0f;
  return acosf(x);
}

Mnożenie kwaternionów

Kwaternion reprezentuje obrót, podobnie jak macierz 4x4 może reprezentować pewne przekształcenie (przesunięcie, obrót, skalowanie, rzutowanie itd.). Macierze, jak pewnie wiesz, można mnożyć. W ten sposób powstaje nowa macierz, której zastosowanie spowoduje wykonanie wszystkich przekształceń, z których powstała, w kolejności w jakiej zostały pomnożone ich macierze. Z kwaternionami jest podobnie - je również można mnożyć składając obroty, które one reprezentują, w jeden. Mnożenie kwaternionów, tak jak macierzy, jest łączne, ale nie jest przemienne. To ma sens, ponieważ wykonanie jednego przekształcenia, a potem drugiego to nie to samo co wykonanie najpierw drugiego, a dopiero potem pierwszego.

Oryginalny wzór na mnożenie kwaternionów działa w taki dziwny sposób, że powoduje złożenie reprezentowanych przez nie obrotów od ostatniego do pierwszego (od prawej do lewej). Żeby nie komplikować zanadto sprawy, przedstawiona niżej funkcja na mnożenie dwóch kwaternionów ma zamienione argumenty. Dzięki temu kwaterniony można mnożyć tak jak macierze, a ich obroty złożą się od lewej do prawej.

void Mul(QUATERNION *Out, const QUATERNION &q1, const QUATERNION &q2)
{
  Out->x = q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y;
  Out->y = q1.w*q2.y + q1.y*q2.w + q1.z*q2.x - q1.x*q2.z;
  Out->z = q1.w*q2.z + q1.z*q2.w + q1.x*q2.y - q1.y*q2.x;
  Out->w = q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z;
}

Oto przykład:

QUATERNION Q, Q1, Q2;
// Ten kwaternion obraca o 45 stopni wokół osi Y
QuaternionRotationY(&Q1, PI * 0.25f);
// Ten kwaternion obraca o 180 stopni wokół osi (0.707, 0.707, 0)
AxisToQuaterion(&Q2, VEC3(0.707f, 0.707f, 0.0f), PI);
// Ten kwaternion reprezentuje obrót najpierw o 45 stopni wokół osi Y,
// a potem o 180 stopni wokół osi (0.707, 0.707, 0)
Mul(&Q, Q1, Q2);

Używanie kwaternionów

OK, mamy już kwaternion, ale jak teraz za jego pomocą obrócić jakiś punkt? Prawdę mówiąc istnieje sposób na obracanie punktu (oczywiście wokół pozycji (0, 0, 0)) bezpośrednio za pomocą kwaterniona, ale to nie jest najlepsze rozwiązanie, bo wymaga dużo obliczeń. Dlatego nie pokażę funkcji, która to robi, a jedynie podam wzór. Żeby to wykonać, trzeba potraktować punkt (x, y, z) jako kwaternion (x, y, z, 0), użyć też odwrotności kwaterniona (o odwracaniu kwaternionów będzie mowa później) i policzyć dwa mnożenia kwaternionów:

P_rotated = ( Inverse(Q) * QUATERNION(P.x,P.y,P.z,0) * Q ).xyz;

Dużo lepszym sposobem na użycie kwaterniona jest jego przekształcenie na macierz. Taka macierz będzie reprezentowała ten sam obrót, który reprezentował kwaternion. Mając ją, możesz już swobodnie przekształcać z jej użyciem punkty i wektory mnożąc je przez tą macierz tak samo, jak to się robi z innymi macierzami, zawierającymi np. translację (przesunięcie) czy skalowanie. Możesz też mnożyć tą macierz z innymi łącząc ich przekształcenia. Oto funkcja przekształcająca kwaternion na macierz:

void QuaternionToRotationMatrix(MATRIX *Out, const QUATERNION &q)
{
  float
    xx = q.x * q.x, yy = q.y * q.y, zz = q.z * q.z,
    xy = q.x * q.y, xz = q.x * q.z,
    yz = q.y * q.z, wx = q.w * q.x,
    wy = q.w * q.y, wz = q.w * q.z;

  Out->_11 = 1.0f - 2.0f * ( yy + zz );
  Out->_12 = 2.0f * ( xy + wz );
  Out->_13 = 2.0f * ( xz - wy );
  Out->_14 = 0.0f;
  
  Out->_21 = 2.0f * ( xy - wz );
  Out->_22 = 1.0f - 2.0f * ( xx + zz );
  Out->_23 = 2.0f * ( yz + wx );
  Out->_24 = 0.0f;
  
  Out->_31 = 2.0f * ( xz + wy );
  Out->_32 = 2.0f * ( yz - wx );
  Out->_33 = 1.0f - 2.0f * ( xx + yy );
  Out->_34 = 0.0f;
  
  Out->_41 = 0.0f;
  Out->_42 = 0.0f;
  Out->_43 = 0.0f;
  Out->_44 = 1.0f;
}

Jeżeli nie potrzebujesz macierzy, a jedynie chcesz przekształcić przez kwaternion pojedynczy punkt, nadal opłaca się zrobić to za pośrednictwem macierzy, ale wtedy można te algorytmy zoptymalizować spisując do jednej funkcji:

void QuaternionTransform(VEC3 *Out, const VEC3 &p, const QUATERNION &q)
{
  float
    xx = q.x * q.x, yy = q.y * q.y, zz = q.z * q.z,
    xy = q.x * q.y, xz = q.x * q.z,
    yz = q.y * q.z, wx = q.w * q.x,
    wy = q.w * q.y, wz = q.w * q.z;

  Out->x =
    (1.0f - 2.0f * ( yy + zz )) * p.x +
    (2.0f * ( xy - wz )) * p.y  +
    (2.0f * ( xz + wy )) * p.z;
  Out->y =
    (2.0f * ( xy + wz )) * p.x +
    (1.0f - 2.0f * ( xx + zz )) * p.y +
    (2.0f * ( yz - wx )) *p.z;
  Out->z =
    (2.0f * ( xz - wy )) * p.x +
    (2.0f * ( yz + wx )) * p.y +
    (1.0f - 2.0f * ( xx + yy )) * p.z;
}

Inne reprezentacje obrotów

Rotacje i orientacje w 3D można reprezentować na trzy sposoby:

1. Kąty Eulera, nazywane Yaw, Pitch, Roll i dla układu współrzędnych jak w DirectX oznaczające obroty wokół osi, odpowiednio, Y, X, Z. Patrz rys. 3.
2. Macierz rotacji, czyli taka macierz 3x3, do której wpisany jest sam obrót, bez innych przekształceń.
3. Kwaternion - o nich tutaj mówimy.

Rys. 3. Kąty Eulera i ich interpretacja.

Można swobodnie przechodzić między tymi reprezentacjami, o ile masz pod ręką odpowiednie funkcje. Funkcję przekształcającą kwaternion na macierz już widziałeś. Algorytm odwrotny - znajdujący kwaternion na podstawie macierzy - znajdziesz w [2] lub [3]. Przekształcenia między kątami Eulera a macierzą są poza zakresem tego artykułu, ale na pewno je znasz - przynajmniej w jedną stronę - bo macierz obrotu tworzysz właśnie podając te kąty (funkcja D3DXMatrixRotationYawPitchRoll). Poniżej przedstawione są funkcje przekształcające między kwaternionami a kątami Eulera, każda w dwóch odmianach.

// Konwersja kątów Eulera na kwaternion,
// który obraca ze wsp. obiektu do świata,
// czyli object -> interial space
void EulerAnglesToQuaternionO2I(QUATERNION *Out,
  float Yaw, float Pitch, float Roll)
{
  Yaw *= 0.5f; Pitch *= 0.5f; Roll *= 0.5f;

  float sy = sinf(Yaw), cy = cosf(Yaw);
  float sp = sinf(Pitch), cp = cosf(Pitch);
  float sr = sinf(Roll), cr = cosf(Roll);

  Out->x = cy*sp*cr + sy*cp*sr;
  Out->y = sy*cp*cr - cy*sp*sr;
  Out->z = cy*cp*sr - sy*sp*cr;
  Out->w = cy*cp*cr + sy*sp*sr;
}

// Konwersja kątów Eulera na kwaternion,
// który obraca ze wsp. świata do obiektu,
// czyli interial -> object space
void EulerAnglesToQuaternionI2O(QUATERNION *Out,
  float Yaw, float Pitch, float Roll)
{
  Yaw *= 0.5f; Pitch *= 0.5f; Roll *= 0.5f;

  float sy = sinf(Yaw), cy = cosf(Yaw);
  float sp = sinf(Pitch), cp = cosf(Pitch);
  float sr = sinf(Roll), cr = cosf(Roll);

  Out->x = -cy*sp*cr - sy*cp*sr;
  Out->y =  cy*sp*sr - sy*cp*cr;
  Out->z =  sy*sp*cr - cy*cp*sr;
  Out->w =  cy*cp*cr + sy*sp*sr;
}

// Konwersja kwaterniona który obraca ze wsp. obiektu do świata,
// czyli object -> interial space, na kąty Eulera
void QuaternionO2IToEulerAngles(float *Yaw, float *Pitch, float *Roll,
  const QUATERNION &q)
{
  float sp = -2.0f * (q.y*q.z - q.w*q.x);

  if (sp == 1.0f)
  {
    *Pitch = PI_2 * sp;
    *Yaw = atan2f(-q.x*q.z + q.w*q.y, 0.5f - q.y*q.y - q.z*q.z);
    *Roll = 0.0f;
  }
  else
  {
    *Pitch = asinf(sp);
    *Yaw = atan2f(q.x*q.z + q.w*q.y, 0.5f - q.x*q.x - q.y*q.y);
    *Roll = atan2f(q.x*q.y + q.w*q.z, 0.5f - q.x*q.x - q.z*q.z);
  }
}

// Konwersja kwaterniona który obraca ze wsp. świata do obiektu,
// czyli interial -> object space, na kąty Eulera
void QuaternionI2OToEulerAngles(float *Yaw, float *Pitch, float *Roll,
  const QUATERNION &q)
{
  float sp = -2.0f * (q.y*q.z + q.w*q.x);

  if (sp == 1.0f)
  {
    *Pitch = PI_2 * sp;
    *Yaw = atan2f(-q.x*q.z - q.w*q.y, 0.5f - q.y*q.y - q.z*q.z);
    *Roll = 0.0f;
  }
  else
  {
    *Pitch = asinf(sp);
    *Yaw = atan2f(q.x*q.z - q.w*q.y, 0.5f - q.x*q.x - q.y*q.y);
    *Roll = atan2f(q.x*q.y - q.w*q.z, 0.5f - q.x*q.x - q.z*q.z);
  }
}

Operacje na kwaternionach

Z kwaternionami można też zrobić kilka innych rzeczy. Na początek warto wiedzieć, że istnieje kwaternion identycznościowy (inaczej jednostkowy) - taki, który reprezentuje rotację pustą, czyli brak jakiegokolwiek obrotu. Przekształcenie przez niego punktu daje ten sam punkt. Jest to analogiczne do macierzy identycznościowej, która odpowiada pustemu przekształceniu (brak przekształcenia). Kwaternion identycznościowy wygląda tak:

const QUATERNION QUATERNION_IDENTITY = QUATERNION(0.0f, 0.0f, 0.0f, 1.0f);

Pierwszą z operacji, którą chciałbym pokazać, jest normalizacja. Kwaternion, tak jak wektor pokazujący jakiś kierunek, można znormalizować, żeby miał długość równą jeden. Niektóre operacje na kwaternionach mogą wymagać, aby był znormalizowany, bo inaczej będą zwracać błędne wyniki (niestety nie jestem w stanie powiedzieć dokładnie, które to są operacje :)

Dlaczego to jest w ogóle potrzebne? Z reprezentowaniem obrotów w 3D jest trochę tak jak z równaniem prostej na płaszczyźnie i równaniem płaszczyzny w przestrzeni. Trzy liczby reprezentujące obrót w 3D (kąty Eulera) to za mało - istnieją przypadki szczególne, których nie da się za ich pomocą zapisać, bo tracimy jeden stopień swobody (zjawisko gimbal lock). Jednak już cztery liczby (kwaternion) to za dużo - istnieje nieskończenie wiele kwaternionów reprezentujących ten sam obrót, ale posiadających różną długość. Dlatego warto wprowadzić normalizację kwaternionów. Macierz 3x3 (9 liczb) to już stanowczo za dużo, dlatego poprawną macierzą rotacji jest tylko macierz ortogonalna (której wszystkie wiersze jak i kolumny są wektorami, które mają długość 1 i są do siebie wzajemnie prostopadłe).

Dość tej teorii, oto funkcja do normalizowania kwaternionów. Dzieli ona po prostu wszystkie składowe kwaterniona przez jego długość.

void Normalize(QUATERNION *InOut)
{
  float n = 1.0f / sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
  InOut->x *= n;
  InOut->y *= n;
  InOut->z *= n;
  InOut->w *= n;
}

Kwaternion możemy zanegować (ang. negation). O dziwo taki zanegowany kwaternion reprezentuje ten sam obrót! Tak więc każdy obrót może być reprezentowany przez dwa różne znormalizowane kwaterniony, przeciwne do siebie. Oto funkcja do negowania kwaternionów. Polega ona, jak widać, na zanegowaniu wszystkich jego składowych.

void Negate(QUATERNION *q)
{
  q->x = -q->x;
  q->y = -q->y;
  q->z = -q->z;
  q->w = -q->w;
}

Jak wobec tego otrzymać kwaternion, który reprezentuje obrót przeciwny do danego? To jest z kolei kwaternion sprzężony (ang. conjugate), wyliczany poprzez zanegowanie części wektorowej kwaterniona (tak jakbyśmy zmienili kierunek obrotu na przeciwny).

void Conjugate(QUATERNION *q)
{
  q->x = -q->x;
  q->y = -q->y;
  q->z = -q->z;
}

Istnieje też kwaternion odwrotny (ang. inverse), ale w przypadku kwaternionów znormalizowanych jest on równy sprzężonemu. Jeszcze inne operacje na kwaternionach to wyliczanie długości, logarytm (log), funkcja wykładnicza (exp) czy potęgowanie (pow). O tych trzech ostatnich nie będziemy mówili, bo nie są nam tu potrzebne.

void Inverse(QUATERNION *Out, const QUATERNION &q)
{
  float NormRc = 1.0f / sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
  Out->w =   q.w * NormRc; 
  Out->x = - q.x * NormRc; 
  Out->y = - q.y * NormRc; 
  Out->z = - q.z * NormRc; 
}

float Length(const QUATERNION &q)
{
  return sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
}

Różnica kwaternionów to kwaternion, który opisuje obrót od orientacji kwaterniona pierwszego do drugiego, tak jak odjęcie dwóch punktów daje wektor przedstawiający przesunięcie między nimi. Wyliczenie różnicy kwaternionów wymaga pomnożenia jednego z nich przez odwrotność drugiego.

void QuaternionDiff(QUATERNION *Out, const QUATERNION &a, const QUATERNION &b)
{
  QUATERNION a_inv;
  Inverse(&a_inv, a);
  Mul(Out, a_inv, b);
}

Iloczyn skalarny (ang. Dot Product) dwóch kwaternionów działa podobnie jak iloczyn skalarny wektorów. Jest liczbą w zakresie [-1; 1], a jego wartość bezwzględna jest tym większa, im bardziej "podobne" do siebie są podane kwaterniony. W przypadku iloczynu skalarnego kwaternionów, inaczej niż przy wektorach, należy zawsze porównywać wartość bezwzględną wyniku, bo policzenie iloczynu skalarnego z wektorem przeciwnym do jakiegoś zwróci przeciwną liczbę, choć obydwa reprezentują ten sam obrót.

float Dot(const QUATERNION &q1, const QUATERNION &q2)
{
  return q1.x*q2.x + q1.y*q2.y + q1.z*q2.z + q1.w*q2.w;
}

Zalety kwaternionów

Jakie korzyści daje używanie kwaternionów, w porównaniu z innymi sposobami na opisywanie obrotów, takimi jak kąty Eulera? Podstawowe zalety są dwie. Pierwsza to brak zjawiska gimbal lock. Kiedy posługujesz się kątami Eulera, może dojść do przypadku szczególnego, kiedy jeden z kątów wyniesie 90 stopni i jedna oś zrówna się z inną. Tracimy wtedy jeden stopień swobody. To wynika z faktu, że budując macierz na podstawie kątów Eulera obroty zawsze składane są w określonej kolejności. Najczęściej to jest Yaw - Pitch - Roll, czyli wokół osi Y, X, Z (a w przypadku przekształcenia z przestrzeni lokalnej modelu do globalnej świata - odwrotnie, czyli Z, X, Y). Próbuje to zilustrować rys. 4. Kwaterniony nie mają tego ograniczenia.

Rys. 4. Obracanie za pomocą kątów Eulera. (Źródło: Wikipedia)

Na gimbal lock nie należy jednak patrzeć jak na zjawisko zawsze negatywne. Zależnie od sytuacji korzystniejsze może być używanie kątów Eulera zamiast kwaternionów. Na przykład kiedy ustawiasz postacie i przedmioty na terenie rozciągającym się w płaszczyźnie poziomej (jak World of Warcraft), opisywanie ich orientacji przez kąty Eulera będzie OK. Jeśli jednak piszesz grę osadzoną w kosmosie w pełnym 3D (jak Nexus The Jupiter Incident), wtedy lepsze będą kwaterniony. Podobnie, do opisywania orientacji kamery typu FPP najlepsze są kąty Eulera. Zjawisko gimbal lock jest tu czymś normalnym i pożądanym. Jednak jeśli kamera lata swobodnie w przestrzeni pokazując sekwencję filmową, lepsze może być użycie kwaternionów.

Druga olbrzymia zaleta kwaternionów to możliwość interpolacji. Interpolacja, mówiąc w skrócie, to znajdowanie wartości pośrednich między podanymi. Interpolacja między dwoma punktami da w wyniku ruch po odcinku łączącym te punkty. Interpolacja między kolorami da płynne przejście między nimi, czyli gradient. Dzięki kwaternionom w analogiczny sposób możemy interpolować orientacje, znajdując pośrednie między dwiema podanymi. Ilustruje to rys. 5.

Rys. 5. Ilustracja zasady działania interpolacji kwaternionów.

Prawdziwa interpolacja kwaternionów odbywa się oczywiście w 3D. Dla parametru t równego 0, obiekt obracany przez kwaternion będzie zorientowany zgodnie z kwaternionem pierwszym - q0. Dla parametru t = 1 obiekt będzie zorientowany zgodnie z kwaternionem drugim - q1. Dla wartości pośrednich t otrzymamy pośrednie orientacje. Algorytm interpolacji "automagicznie" znajdzie najkrótszą drogę między tymi dwiema orientacjami (tak jakby po powierzchni sfery) i po niej będzie poruszał obiekt, żeby zapewnić jego płynny obrót wraz ze zmianą parametru t. Niesamowite, prawda? :) Ta operacja nosi nazwę "Slerp" (od Spherical Linear Interpolation) i wygląda tak:

void Slerp(QUATERNION *Out, const QUATERNION &q0, const QUATERNION &q1, float t)
{
  float cosOmega = Dot(q0, q1);
  QUATERNION new_q1 = q1;
  if (cosOmega < 0.0f)
  {
    new_q1.x = -new_q1.x;
    new_q1.y = -new_q1.y;
    new_q1.z = -new_q1.z;
    new_q1.w = -new_q1.w;
    cosOmega = -cosOmega;
  }

  float k0, k1;
  if (cosOmega > 0.9999f)
  {
    k0 = 1.0f - t;
    k1 = t;
  }
  else
  {
    float sinOmega = sqrtf(1.0f - cosOmega*cosOmega);
    float omega = atan2f(sinOmega, cosOmega);
    float oneOverSinOmega = 1.0f / sinOmega;
    k0 = sinf((1.0f - t) * omega) * oneOverSinOmega;
    k1 = sinf(t*omega) * oneOverSinOmega;
  }
  Out->x = q0.x*k0 + new_q1.x*k1;
  Out->y = q0.y*k0 + new_q1.y*k1;
  Out->z = q0.z*k0 + new_q1.z*k1;
  Out->w = q0.w*k0 + new_q1.w*k1;
}

Istnieje też interpolacja między czterema kwaternionami "Squad", ale nie będziemy o niej mówić.

Bibliografia

1. Fletcher Dunn, Ian Parberry, 3D Math Primer for Graphics and Game Development, Wordware Publishing, 2002.
2. Ken Shoemake, Quaternion Calculus and Fast Animation, SIGGRAPH, 1987.
3. J.M.P. van Waveren, From Quaternion to Matrix and Back, 2005, http://cache-www.intel.com/cd/00/00/29/37/293748_293748.pdf.

Adam Sawicki
regedit.gamedev.pl
20.05.2008

Komentarze (14)