C++ - Projet quaternions

ven. 02 avril 2010

Ce projet est basé sur un objet mathématique, le quaternion. On se sert de certains quaternions, ceux qui se trouvent sur la sphere 4D unitaire i.e l'ensemble des quaternions q tel que |q| = 1, en 3D pour représenter des rotations. Votre travail est d'implanter une classe quaternion, une classe matrice3x3, une fonction d'interpolation de quaternions, et une fonction de conversion de quaternion vers matrice.


Définition

Un quaternion est un élément (a, b, c, d) € R4 ; on note H l'ensemble des quaternions. (H, +, x) est un corps, non commutatif.

L'ensemble H possède une structure d'algèbre sur R. C'est en particulier un R-espace vectoriel de dimension 4, dont la famille {1, i, j, k} est une base, donc tout quaternion s'écrit sous la forme d'une combinaison linéraire des éléments de cette base.

Un quaternion s'écrit (avec w, x, y, z réels) :

  • q = w + xi + yj + zk,
  • q = (w, x, y, z), par abus de notation,
  • q = [w, (x, y, z)] avec w la partie réelle, et (x, y, z) la partie imaginaire.

i, j, k suivent les règles ci-dessous :

  • i2 = j2 = k2 = ijk = -1
  • ij = -ji = k
  • jk = -kj = i
  • ki = -ik = j

Seuls les quaternions unitaires permettent de représenter des rotations (de la même manière que les matrices doivent être orthonormales pour représenter une rotation). On appelle quaternion unitaire, un quaternion tel que |q| = 1.

Soit u(x, y, z) un vecteur unitaire de R3, et θ un angle compris entre 0 et 2 π, alors la rotation d'axe u et d'angle θ (dans le sens direct) correspond au quaternion :

q = cos θ/2 + x sin θ/2 i + y sin θ/2 j + z sin θ/2 k

Ce quaternion est encore noté (cos θ/2, sin θ/2 u).

Pour composer des rotations sous forme de quaternion, on les multiplie par l’opérateur de multiplication des quaternions (qui généralement n'est pas symétrique). Pour inverser une rotation, on inverse le quaternion en multipliant par –1 chaque composante de l’axe.

Ceci est l'explication réduite au plus simple. Pour des compléments, vous êtes libre de chercher sur internet.


Les opérations

Vous devez implanter des opérateurs pour les opérations suivantes. Ces opérateurs seront internes à la classe. Toutes les formules doivent être transformées en C++.

Addition et soustraction

Soient p et q deux quaternions : q = (a, b, c, d) = a + ib + jc + kd et p = (w, x, y, z) = w + ix + jy + kz, les résultats de l'addition et de la soustraction sont les suivants :

  • q + p = (a + w, b + x, c + y, d + z) = (a + w) + i(b + x) + j(c + y) + k(d + z)
  • q - p = (a - w, b - x, c - y, d - z) = (a - w) + i(b - x) + j(c - y) + k(d - z)

Multiplication par un réel

Soit x un réel, et soit q un quaternion, q = (a, b, c, d) = a + ib + jc + kd, le résultat de la multiplication d'un quaternion par un réel est :

  • xq = qx = (xa, xb, xc, xd) = xa + i(xb) + j(xc) + k(xd)

Conjugué et norme

Soit q un quaternion, q = (a, b, c, d) = a + ib + jc + kd, :

  • le conjugué de q est q = (a, -b, -c, -d) = a - ib - jc - kd
  • la norme de q est : |q| = sqrt(a2 + b2 + c2 + d2) = sqrt(qq) avec sqrt() la racine carrée

Multiplication

La multiplication n'est pas commutative. Soit q = (a, b, c, d) et p = (w, x, y, z) deux quaternions, la multiplication à droite de q par p donne :

qp = (a + ib + jc + kd)(w + ix + jy + kz)

= (aw - bx - cy - dz, ax + bw + cz - dy, ay - bz + cw + dx, az + by - cx + dw)

Inverse et division

On a : qq = qq = |q|2

L'inverse de q est : q-1 = q / |q|2

On peut utiliser l'inverse pour définir la division : q / p = qp-1

Sous espaces réels et complexes

(a, 0, 0, 0) = a est un nombre réel ; (a, b, 0, 0) = a + ib, (a, 0, c, 0) = a + jc et (a, 0, 0, d) = a + kd sont des quaternions classiques. Les quaternions purs sont de la forme q = ib + jc + kd.


La classe quaternion

La structure de votre classe quaternion est libre, mais sachez qu'habituellement, on utilise trois floats (x, y, z) pour la partie imaginaire et un quatrième float (w) pour la partie réelle.

Liste des méthodes :

  • constructeur par défaut,
  • constructeur par copie,
  • constructeur avec valeurs d'initialisation des membres,
  • méthode qui normalise le quaternion,
  • opérateur d'addition (+),
  • opérateur de soustraction (-),
  • opérateur de multiplication (*) par un réel (opérateur commutatif) et par un quaternion à gauche (opérateur non commutatif),
  • opérateur de division (/),
  • opérateur unaire de conjugaison (~),
  • opérateur d'affectation,
  • les méthodes d'entrée/sortie dont read et print, explicitées plus loin.

La classe matrice3x3

La structure de la matrice est libre, mais sachez qu'habituellement, on utilise un tableau bidimensionnel (3 par 3) de floats. Nous avons établi que les matrices devront être faites pour se multiplier avec des vecteurs colonnes. Soit m une matrice, v un vecteur 3D, et v' le transformé de v par m alors v' = m * v.

Nous n'exigeons rien de plus pour cette classe, à part les méthodes d'entrée/sortie read et print, explicitées plus loin.


Les fonctions outils

Spherical Linear IntERPolation (SLERP)

Cette méthode permet de passer du quaternion q1 au quaternion q2 de manière linéaire sur une sphère. On veut donc déterminer le chemin d'interpolation le plus court entre les deux quaternions sur la sphere unité des quaternions. Avec une vitesse angulaire constante, la courbe d'interpolation forme un grand arc sur la sphère unité (en termes de géométrie différentielle, ce grand arc est un géodésique ­ correspondant à une ligne droite que l'on projette sur la sphère) qui détermine donc la courbe optimale d'interpolation entre deux rotations.

Cette fonction est donc particulièrement utile pour calculer l'interpolation linéaire entre deux rotations exprimées sous forme de quaternions, par exemple dans le cadre de l'animation d'un objet 3D tel qu'une caméra dans l'animation temps réel dans les jeux vidéo.

Soit q1 et q2 deux quaternions et t un reel compris entre 0 et 1 on a :

q(t) = c1(t) q1 + c2(t) q2

avec c1 et c2 des fonctions réelles de 0 <= t <= 1. L'angle entre q(t) et q1 est tθ et l'angle entre q(t) et q2 est (1-t)θ. En calculant le produit de q(t) et q1 d'une part et de q(t) et q2 d'autre part, il vient (attention aux crochets !) :

Slerp(q1, q2, t) = [q1 sin((1-t) θ) + q2 sin (t θ)] / sin(θ)

Pour le prouver, il faut utiliser le produit scalaire de deux quaternions unitaires, ce qui donne un quaternion unitaire, soit :

q1 . q2 = q1 q2 = cos θ + sin θ u

Appliquons le à q(t) :

q(t) . q1 = q(t) q1 = c1(t) q1 q1 + c2(t) q2 q1

= c1(t) |q1|2 + c2(t) q2 q1

= c1(t) + c2(t) cos θ + c2(t) sin θ u

Remarque : si q1 et q2 sont deux quaternions unitaires, alors q(t) est aussi un quaternion unitaire et nous pouvons alors écrire : Slerp(q1, q2, t) = q1(q1q2)t

Rappel : un quaternion pouvant être assimilé à un axe de rotation et un angle, la mesure de l'angle entre deux quaternions q1 et q2 peut être calculée avec la formule suivante (produit vectoriel de deux vecteurs dans R3) :

sin θ = |p1 ^ p2|

avec p1 et p2 les quaternions purs correspondant à q1 et q2 respectivement.

La fonction à implanter devra prendre en paramètres 2 quaternions, et 1 float représentant le pas de l'interpolation.

Quaternions et rotations dans l'espace

Préambule, représentation matricielle des quaternions : de par la structure d'algèbre de H, la multiplication (à gauche) par un quaternion q donné, soit l'application p -> q . p, est une application linéaire sur H. Si q s'écrit a + bi + cj + dk, cette application a pour matrice, dans la base 1, i , j, k :

a -b -c -d

b a -d c

c d a -b

d -c b a

La conjugaison par un quaternion non nul de norme 1 peut s'interpréter comme une rotation dans l'espace. Soit l'application Sqdéfinie sur H par :

Sq : p -> q . p . q-1 = q . p . q

A un quaternion pur, c'est-à-dire sans partie réel, on fait correspondre le vecteur p(x, y, z) de R3 alors pour transformer ce point p, par la rotation représentée par le quaternion q = [w, (x, y, z)], on réalise l'opération q[0, p]q-1. Sans faire de démonstration, on peut établir le résultat suivant (dans la base canonique {i, j, k}) :

1 - 2y2 - 2z2 2xy - 2zw 2xz + 2yw

Mq = 2xy + 2zw 1 - 2x2 - 2z2 2yz - 2xw

2xz - 2yw 2yz + 2xw 1 - 2x2 - 2y2

La fonction à implanter devra prendre en paramètre un quaternion unitaire et retourner une matrice.


Entrées / sorties

Vous devez implémenter les méthodes qui sortent sur un ofsteam (print) et lisent sur un ifstream (read) un quaternion, présentés sous le format suivant :

q.x q.y q.z q.w

Et une matrice 3x3 sous la forme suivante :

m._11 m._12 m._13

m._21 m._22 m._23

m._31 m._32 m._33

Attention à respecter ce format, ce sera une moulinette qui vous corrigera dans un premier temps.

Vous devrez aussi implémenter l'opérateur binaire de sortie `<<' prenant un ostream à sa gauche ainsi que l'opérateur binaire d'entrée `>>' prennant un istream à sa gauche. Nous restons volontairement imprécis sur la description de ces prototypes, inserez les const et & nécessaires pour rester cohérents.


Programme principal d'interpolation

Vous devez réaliser un programme simple qui prend sur l'entrée standard, des quaternions et donne en sortie des matrices. Le programme prend une liste de quaternions, il interpole les quaternions 2 à 2 en utilisant le Slerp : q0 et q1, puis q1 et q2, etc.

Il prend un certain nombre d'échantillons de ces interpolations (des quaternions), ce nombre étant précisé en entrée. Les matrices en sortie correspondent à la conversion de ces échantillons en matrice 3x3, plus la matrice 3x3 équivalente au dernier quaternion de la liste d'entrée.

Format d'entrée : l'entier n d'échantillons pour l'interpolation entre chaques paires de quaternions puis les quaternions (jusqu'à EOF), que vous lirez grâce à l'opérateur "`>>'", sous la forme :

q0.x q0.y q0.z q0.w

q1.x q1.y q1.z q1.w

q2.x q2.y q2.z q2.w

Format de sortie :

m0._11 m0._12 m0._13

m0._21 m0._22 m0._23

m0._31 m0._32 m0._33

m1._11 m1._12 m1._13

m1._21 m1._22 m1._23

m1._31 m1._32 m1._33

Toutes les valeurs en sortie seront à donner avec les précisions, et tailles par défaut, donc n'utilisez pas les manipulateurs 'setprecision', ou 'setw'.

Pour éviter les ambiguļtés, voici deux exemples :

  1. si le nombre d'échantillons est 1, et qu'il y a 2 quaternions donnés, vous devez sortir exactement les 2 matrices correspondantes aux 2 quaternions donnés,
  2. si le nombre d'échantillons est 2, et que vous avez 3 quaternions en entrée, vous allez retourner 5 matrices, dans l'ordre : celle de q0, celle du quaternion Slerp à 0.5 entre q0 et q1, celle de q1, celle du quaternion Slerp à 0.5 entre q1 et q2 et enfin celle de q2.

Application numérique :

L'entrée :

4
0.182574 0.365148 0.547723 0.730297
0.377964 0.377964 0.377964 0.755929
0.632456 0.316228 0.632456 0.316228

Donne en sortie :

0.133333 -0.666668 0.733333
0.933334 0.333333 0.133333
-0.333333 0.666666 0.666667
0.208262 -0.581009 0.786801
0.926716 0.374462 0.0312234
-0.312769 0.722638 0.616417
0.283293 -0.487891 0.825656
0.911617 0.404352 -0.0738506
-0.297824 0.773603 0.559319
0.357149 -0.388897 0.849237
0.888293 0.422491 -0.1801
-0.288755 0.818694 0.496347
0.428573 -0.285714 0.857141
0.857141 0.428573 -0.285714
-0.285714 0.857141 0.428573
0.301862 -0.265855 0.915533
0.950092 0.163236 -0.265855
-0.0787684 0.950092 0.301862
0.181748 -0.207934 0.96111
0.971739 -0.111745 -0.207934
0.150635 0.971739 0.181748
0.0782489 -0.116781 0.990071
0.920278 -0.373434 -0.116781
0.383362 0.920278 0.0782489
-1.49445e-06 0 1
0.800001 -0.600002 0
0.600001 0.800001 -1.49445e-06

Modalités de rendu

La date de rendu est fixée au ??/??/2005. Vous devez rendre votre miniprojet dans un tarball nommé : login-quat.tgz et placer dans votre répertoire : $HOME avec les droits 444. Vos fichiers devront s'extraire dans le répertoire : ./login-quat/ (attention au point).Votre exécutable doit s'appeler quat.

Votre tarball devra contenir, en plus des modules, un makefile dont les règles suivantes sont nécessaires :

  • all : crée le programme demandé,
  • clean : nettoie tout ce qui a été généré pendant la compilation du programme.

La notation se base sur le passage des tests, le respect du sujet et le respect des principes et conseils énoncés dans le sujet. Les commentaires utiles seront pris en compte dans la notation.