QUATERNIONY
===========
1. Co jsou quaterniony?
--------------------------
Quaterniony jsou zobecnením komplexnich čísel ve 3D. Pan Hamilton někdy
v 17. století přišel na to, že je zapotřebí zavést jednu reálnou a _tři_
imaginární složky, aby byla správně vytvořena tato algebra tvořící grupu.
Quaternion je uspořádaná čtveřice (w, x, y, z). Častěji však bývá
zapisován v úspornějším formátu (w, v), kde w je skalár a v je 3D vektor.
Pro počítačovou grafiku je však nejvýznamnější tzv. jednotkový quaterion |q|=1,
který se dá zapsat jako q=(cos(O), sin(O)*v), kde |v|=1. Tento
quaternion totiž popisuje rotaci/orientaci ve 3D. Jednotkový vektor v
je osou rotace a O je _dvojnásobek_ úhlu otočeni kolem této osy.
Navíc také platí vztah log(q)=O*v (nebo opačně exp(o*v)=q,
v literature se tomuto vztahu říka 'exponencial map').
Největší použití mají quaterniony při interpolaci orientace. Máme-li definovaný
dvě orientace M1 a M2 (obvykle popsané transformačními maticemi), nalézt
interpolovanou orientaci Mo "ležící" na trajektorii spojující M1 a M2,
neni triviální úkol. Snadnou úvahou poznáme, že možných interpolací může být
nekonečně mnoho, stejně jako nekonečně mnoho je možných trajektorií.
Obdobně jako při hledání interpolace v kartézském systému, i na povrchu hyperkoule
se snažíme najít tu nejkratší trajektorii spojující obě orientace.
Quaterniony jsou k řešení tohoto úkolu velmi mocný nástroj. Sférická intepolace
dvou quaternionů totiž poskytuje nejkratší(!) spojnici dvou quaternionů na povrchu
hyperkoule. V důsledku to znamená, že pohyb po takto vytvořené trajektorii zajišťuje
minimální úhlové rychlosti (myslím, že matematik by napsal, že 2-norma diferencí
je minimální...). V aplikaci se to pak pozitivně projeví tak, že (například) rotující
kamera se nám nebude zbytečne otáčet jako čamrda kolem svých os.
Obvyklý postup je pak následující:
1. transformačni matice se převedou na quaterniony
2. provede se sférická interpolace
slerp(q1,q2,t) = q0*sin(O*(1-t))/sin(O) + q1*sin(O*t)/sin(O),
kde plati cos(o)=q1*q2
3. interpolovaný quaternion se převede na transformační matici
Více o této problematice určitě najdete na webu. Hledejte například
Shoemake, K., Quaternions
Dam, E., Koch, M., Lillholm, M., Quaternions, Interpolation and Animation
Další poznatky:
Násobení quaternionů není komutativní a chová se (a má stejné výsledky) jako násobeni
adekvátních matic. Quaterniony se z pohledu počítačové grafiky řadí mezi jednu z metod
parametrizace orientace v prostoru. K tomu stačí u quaternionu čtyři parametry.
Srovnejte s maticí, která jich potřebuje devět! Ješte lépe jsou na tom exponenciální
mapy, které potřebují jenom tři! Úspora paměti je zde evidentní. Zajimavým problémem
je tvorba aproximačních/interpolačních křivek na povrchu hyperkoule.
///////////////////////////////////////////////////////////////////////////
Převzato od autora Defora z diskusního fóra o OpenGL na www.builder.cz
http://forum.builder.cz/read.php?f=124&i=205&t=159
Ostatní informace jsem získal ze stránky
http://www.cs.ualberta.ca/~andreas/math/matrfaq_latest.html
ze které jsem část věnovanou quaternionům volně přeložil
///////////////////////////////////////////////////////////////////////////
2. Jak quaterniony použít pro 3D animaci?
-----------------------------------------------
Pro rotaci pomocí matic musíme mít pro každou osu x,y,z rotační matici.
Vynásobením těchto tří matic nám vyjde konečná rotační matice. Quaternion
má výhodu v tom, že všechny osy jsou obsaženy v jednom quaternionu.
Quaternion pak převedeme rovnou na konečnou rotační matici.
Protože osa rotace je vlastně směrový vektor,
můžeme ji vypočítat pomocí vektorových algoritmů, nebo pomocí
sférických souřadnic (délka (longitude)/ šířka (latitude)).
Quaternions má další výhodu v tom, že můžeme quaterniony mezi sebou interpolovat.
Tím docílíme plynulé a bezchybné rotace.
3. Jak vypočítat konjugaci quaternionu ?
------------------------------------------------------
Toho můžeme dosáhnout otočením polarity ,nebo negací
vektorové části quaternionu:
--
Qr = ( Qr.scalar, -Qr.vector )
----------------------------------------------------------------
quaternion_conjugate( QUAT *qr, QUAT *qa )
{
qr -> qw = qa -> qw;
qr -> qx = -qa -> qx;
qr -> qy = -qa -> qy;
qr -> qz = -qa -> qz;
}
--------------------------------------
4. Jak vypočítat obrácený (opačný) quaternion?
----------------------------------------------------
Toto je protihodnota konjugace quaternionu, jestliže quaternion je normalizován.
Jinak se velikost otočení vypočítá podle vzorce 1/|q|.
5. Jak vypočítat velikost quaternionu?
------------------------------------------------------
Velikost quaternionu se vypočítá vynásobením quaternionu s jeho
konjugací:
------------
/ --
|Qr| = \/ Qr.Qr
Vzorec můžeme napsat takto:
-------------------------------------------------------------------
QFLOAT quaternion_magnitude( QUAT *qa )
{
return( sqrt(qa->qw*qa->qw+
qa->qx*qa->qx+ qa->qy*qa->qy+qa->qz*qa->qz) )
}
-------------------------------------------------------------------
6. Jak normalizovat quaternion?
-------------------------------------
Quaternion se normalizuje stejně jako vektor. Nejdříve se vypočítá velikost
quaternionu. Potom celý quaternion dělíme velikostí quaternionu:
-------------------------------------------------------------
QFLOAT mag = quaternion_magnitude(Quaternion);
Qr = Qr / mag;
-------------------------------------------------------------
Quternion po normalizaci nepřesáhne hodnotu 1.0
7. Jak spojit (násobit) dva quaterniony dohromady?
------------------------------------------------
Jsou stanoveny dva quaterniony Q1 a Q2. Cílem je vypočítat
složenou rotaci Qr:
Qr = Q1.Q2
Toto je dosažitelné pomocí výrazu:
Qr = Q1.Q2 = ( w1.w2 - v1.v2, w1.v2 + w2.v1 + v1 x v2 )
kde
v1 = (x,y,z) of Q1
w1 = (w) of Q1
v2 = (x,y,z) of Q2
w2 = (w) of Q2
_._ je vektor společného bodu
_x_ je kolmý vektor na oba vektory
(dot and cross products).
Vzorec můžeme napsat takto:
---------------------------------------------------
quaternion_multiply( QUAT *qr, QUAT *qa, QUAT *qb )
{
qr.scalar = qa->scalar * qb->scalar - v3_dot( &qa->vector, &qb->vector );
v3_cross( &va, &qa->vector, &qb->vector );
v3_scalef( &vb, &qa->vector, &qb->scalar );
v3_scalef( &vc, &qb->vector, &qa->scalar );
v3_add( &va, &va, &vb );
v3_add( &qr->vector, &va, &vc );
quaternion_normalise( qr );
}
---------------------------------------------------
Optimalizace pak vypadá takto:
w = w1w2 - x1x2 - y1y2 - z1z2
x = w1x2 + x1w2 + y1z2 - z1y2
y = w1y2 + y1w2 + z1x2 - x1z2
z = w1z2 + z1w2 + x1y2 - y1x2
8. Jak převést quaternion na rotační matici ?
--------------------------------------------------------
Předpokládejmet že quaternion je vytvořený ve formě:
Q = |X Y Z W|
Potom můžeme quaternion konvertorovat do 4x4 matice rotace
použitím následujícího výrazu (Varování: musíme mít na zřeteli to,
že jde o OpenGL matici):
¦ 2 2 ¦
¦ 1 - (2Y + 2Z ) 2XY + 2ZW 2XZ - 2YW ¦
¦ ¦
¦ 2 2 ¦
M = ¦ 2XY - 2ZW 1 - (2X + 2Z ) 2YZ + 2XW ¦
¦ ¦
¦ 2 2 ¦
¦ 2XZ + 2YW 2YZ - 2XW 1 - (2X + 2Y ) ¦
¦ ¦
Jelikož jde o matici 4x4 , můžeme přidat spodní řadu a pravý krajní sloup..
Vzorec můžeme napsat takto:
----------------
xx = X * X;
xy = X * Y;
xz = X * Z;
xw = X * W;
yy = Y * Y;
yz = Y * Z;
yw = Y * W;
zz = Z * Z;
zw = Z * W;
mat[0] = 1 - 2 * ( yy + zz );
mat[1] = 2 * ( xy - zw );
mat[2] = 2 * ( xz + yw );
mat[4] = 2 * ( xy + zw );
mat[5] = 1 - 2 * ( xx + zz );
mat[6] = 2 * ( yz - xw );
mat[8] = 2 * ( xz - yw );
mat[9] = 2 * ( yz + xw );
mat[10] = 1 - 2 * ( xx + yy );
mat[3] = mat[7] = mat[11 = mat[12] = mat[13] = mat[14] = 0;
mat[15] = 1;
Výsledná matice má následující pozici:
¦ mat[0] mat[4] mat[ 8] mat[12] ¦
M = ¦ mat[1] mat[5] mat[ 9] mat[13] ¦
¦ mat[2] mat[6] mat[10] mat[14] ¦
¦ mat[3] mat[7] mat[11] mat[15] ¦
----------------
9. Jak převést rotační matici na quaternion ?
--------------------------------------------------------
Otáčení může být převedeno zpět do quaternion pomocí
následujícího algoritmu:
Nejdříve vypočítáme matici T podle rovnice:
2 2 2
T = 4 - 4x - 4y - 4z
2 2 2
= 4( 1 -x - y - z )
= 1 + mat[0] + mat[5] + mat[10]
Jestliže je matice větší než nula, tak
počítáme dál:
Test if ( T > 0.00000001 ) to avoid large distortions!
S = sqrt(T) * 2;
X = ( mat[9] - mat[6] ) / S;
Y = ( mat[2] - mat[8] ) / S;
Z = ( mat[4] - mat[1] ) / S;
W = 0.25 * S;
Jestliže je matice rovna nule, tak vypočítáme,
který prvek úhlopříčky má největší hodnotu.
Závisí to na následujícím výpočtu:
if ( mat[0] > mat[5] && mat[0] > mat[10] ) { // Column 0:
S = sqrt( 1.0 + mat[0] - mat[5] - mat[10] ) * 2;
X = 0.25 * S;
Y = (mat[4] + mat[1] ) / S;
Z = (mat[2] + mat[8] ) / S;
W = (mat[9] - mat[6] ) / S;
} else if ( mat[5] > mat[10] ) { // Column 1:
S = sqrt( 1.0 + mat[5] - mat[0] - mat[10] ) * 2;
X = (mat[4] + mat[1] ) / S;
Y = 0.25 * S;
Z = (mat[9] + mat[6] ) / S;
W = (mat[2] - mat[8] ) / S;
} else { // Column 2:
S = sqrt( 1.0 + mat[10] - mat[0] - mat[5] ) * 2;
X = (mat[2] + mat[8] ) / S;
Y = (mat[9] + mat[6] ) / S;
Z = 0.25 * S;
W = (mat[4] - mat[1] ) / S;
}
Quaternion je definován:
Q = | X Y Z W |
10. Jak převést rotační úhel osy do quaternionu?
----------------------------------------------------------------
Následující algoritmus toto provádí
------------------------------------------------
sin_a = sin( angle / 2 );
cos_a = cos( angle / 2 );
X = axis -> x * sin_a;
Y = axis -> y * sin_a;
Z = axis -> z * sin_a;
W = cos_a;
quaternion_normalise( |X,Y,Z,W| );
------------------------------------------------
Je nutné výsledný quaternionu normalizovat, aby nedošlo k dělení nuly.
11. Jak konvertovat quaternion na rotační úhel osy ?
----------------------------------------------------------------
Quaternion může být převeden na rotační úhel osy
použitím následujícího algoritmu:
---------------------------------------------------
quaternion_normalise( |X,Y,Z,W| );
cos_a = W;
angle = acos( cos_a ) * 2;
sin_a = sqrt( 1.0 - cos_a * cos_a );
if ( fabs( sin_a ) < 0.0005 ) sin_a = 1;
axis -> x = X / sin_a;
axis -> y = Y / sin_a;
axis -> z = Z / sin_a;
---------------------------------------------------
12. Jak převést úhly sférické rotace (rotace po povrchu hyperkoule) na quaternion?
----------------------------------------------------------------
Rotační osa by měla být definována použitím sférických kordinací
šířkou (latitude), délkou (longitude) a úhlem rotace.
Takto vypadá algoritmus:
-----------------------
sin_a = sin( angle / 2 )
cos_a = cos( angle / 2 )
sin_lat = sin( latitude )
cos_lat = cos( latitude )
sin_long = sin( longitude )
cos_long = cos( longitude )
X = sin_a * cos_lat * sin_long
Y = sin_a * sin_lat
Z = sin_a * sin_lat * cos_long
W = cos_a
-----------------------
13. Jak převést quaternion na sférickou (kulovou) rotaci úhlů?
----------------------------------------------------------------
Quaternion může být převeden na sférické koordinace pomocí
tohot výpočtu:
-----------------------
cos_a = W;
sin_a = sqrt( 1.0 - cos_a * cos_a );
angle = acos( cos_a ) * 2;
if ( fabs( sin_angle ) < 0.0005 ) sin_a = 1;
tx = X / sin_a;
ty = Y / sin_a;
tz = Z / sin_a;
latitude = -asin( ty );
if ( tx * tx + tz * tz < 0.0005 )
longitude = 0;
else
longitude = atan2( tx, tz );
if ( longitude < 0 )
longitude += 360.0;
-----------------------
14. Jak převádět rotační úhly zadaných ve stupních na quaternion?
-------------------------------------------------------------------
Převod rotací úhlů (zadaných ve stupních) na quaterion může být vykonán pomocí
součinu dvou quaternionů. Všechny rotace úhlu jsou převedeny
na úhly os obou quaternionů, které odpovídají Euklidovským osám.
Úhly všech os jsou převedeny na quaterniony a vynásobeny.
Tím nám vyjde konečný quaternion.
Následující kód toto demostruje:
---------------------------------------------
quaternion_from_euler( QUATERNION *q, VFLOAT ax, VFLOAT ay, VFLOAT az )
{
VECTOR3 vx = { 1, 0, 0 }, vy = { 0, 1, 0 }, vz = { 0, 0, 1 };
QUATERNION qx, qy, qz, qt;
quaternion_from_axisangle( qx, &vx, rx );
quaternion_from_axisangle( qy, &vy, ry );
quaternion_from_axisangle( qz, &vz, rz );
quaternion_multiply( &qt, &qx, &qy );
quaternion_multiply( &q, &qt, &qz );
}
---------------------------------------------
Více informací naleznete na stránce:
http://vered.rose.utoronto.ca/people/david_dir/GEMS/GEMS.html
//Pitch->X axis, Yaw->Y axis, Roll->Z axis
Quaternion::Quaternion(float fPitch, float fYaw, float fRoll)
{
const float fSinPitch(sin(fPitch*0.5F));
const float fCosPitch(cos(fPitch*0.5F));
const float fSinYaw(sin(fYaw*0.5F));
const float fCosYaw(cos(fYaw*0.5F));
const float fSinRoll(sin(fRoll*0.5F));
const float fCosRoll(cos(fRoll*0.5F));
const float fCosPitchCosYaw(fCosPitch*fCosYaw);
const float fSinPitchSinYaw(fSinPitch*fSinYaw);
X = fSinRoll * fCosPitchCosYaw - fCosRoll * fSinPitchSinYaw;
Y = fCosRoll * fSinPitch * fCosYaw + fSinRoll * fCosPitch * fSinYaw;
Z = fCosRoll * fCosPitch * fSinYaw - fSinRoll * fSinPitch * fCosYaw;
W = fCosRoll * fCosPitchCosYaw + fSinRoll * fSinPitchSinYaw;
}
---------------------------------------------
15. Jak použít quaterniony pro výpočet lineární interpolace mezi maticemí?
-------------------------------------------------------------------------------
Pro mnoho aplikací, používající animace, je nutné vypočítat interpolace
mezi dvěmi rotačními pozicemi daného objektu. Tyto pozice mohou být
použity pro snímky animace, nebo pro inverzní kinematiku.
Pro použití této metody musíme znát dvě rotační matice a
požadovaný cíl interpolace mezi nimi. Dvě matice si označíme
(MS a MF).
Interpolace rotační matice je
vytvořena použitím rovnice průhlednosti a parametrem T, který má
velikost od 0.0 do 1.0.
T=0, interpolační matice je rovna počáteční matici.
T=1, interpolační matice je rovna konečné matici.
Interpolovaná rotační matice (MI) je takto specifikována:
MI = F( MS, MF, T )
kde F je funkce průhlednoti.
První krok. Interpolace mezi dvěmi maticemi je
rozhodující pro rotační matici, která bude převedena z MS do MF.
Toho dosáhneme následujícím výrazem:
-1
T = Ms . Mf
kde Ms je počáteční matice,
Mf je konečná matice,
a T je střední matice.
Druhý krok. Převod matice na rotační úhel osy.
Toho dosáhneme převedením matice na quaternion
a nakonec převedením quaternionu na požadovaný rotační úhel osy.
Použijeme matici 4x4 pro následující výpočet:
m4_transpose( mt, ms ); /* Otočení */
m4_mult( ms, mt, mb ); /* Rotační matice */
m4_to_axisangle( ms, axis, angle ); /* Rotation osa/úhel */
for ( t = 0; t < 1.0; t += 0.05 )
{
m4_from_axisangle( mi, axis, angle * t ); /* Final interpolation */
... whatever ...
}
kde t je interpolační faktor v rozmezí od 0.0 do 1.0
16. Jak použít quaterniony pro kubickou interpolaci mezi maticemi?
------------------------------------------------------------------------------
Některé aplikace nemohou vhodným způsobem použít lineární
interpolaci pro animční účely. Proto tu máme jako alternativu kubickou
interpolaci.
Pro kubickou interpolaci potřebujeme znát čtyři rotační matice.
Všechny matice jsou potom převedeny
přes quaterniony na sférické rotační úhly (délka (longitude), šířka (latitude)
a rotační úhel).
Úhly jsou potom násobeny základní maticí pro Kardinální křivky.
Tato interpolace matice může potom být použita pro určení
střední sférické rotace úhlů.
Jedny interpolované koordinace známe. Sířku (latitude), délku(longitude) a
rotační úhel. Interpolovaná rotace matice může potom být vytvořena
přes konverzi do quaternionů.
Použijeme matici 4x4 pro následující výpočet:
----------------------------------------------------------------------
for ( n = 0; n < 4; n++ )
m4_to_spherical( mat[n], &v_sph[n] ); /* Sférické koordinace */
m4_multspline( m_cardinal, v_sph, v_interp ); /* Interpolační vektor */
...
v3_cubic( v_pos, v_interp, t ); /* Interpolace */
m4_from_spherical( m_rot, v_pos ); /* Konečná matice */
----------------------------------------------------------------------
17. Jak použít quaternions pro rotaci vektoru?
------------------------------------------------------------------------------
Je to velmi elegantní cesta rotace vektoru za použité quaternionu,
protože je výpočet velmi rychlý (qr je rotační quaternion):
-1
v' = qr * v * qr
Takto můžeme jednoduše vytvářet velmi rychlé transformační
a rotační matice.