MATRICES ET ROTATIONS 3D

L'utilisation des matrices pour effectuer des transformations 3D fait parfois peur car elles ont la réputation d'être difficiles à manier. En réalité, elles permettent de se simplifier la vie et de gérer à peu de frais des mouvements complexes.

• Conventions et notations
• Translation
• Rotation
• Un peu de pratique

CONVENTIONS ET NOTATIONS

Dans ce message, on se limite a des objets indéformables composés de polygones plans. Par conséquent, chaque objet est complètement défini par les sommets des polygones qui le composent.

On parlera des points de l'objet. Pour clarifier les notations, je note les points et vecteurs "à plat" et par des lettres majuscules, comme ceci : P = (0,0,0).

Je note les matrices entre crochets comme ceci : [T].
Le produits de matrices (ou de n'importe-quoi sauf ambiguité) est noté par une *, le produit scalaire par un . et il n'y a pas de produit vectoriel.

En ce qui concerne le sens des angles et l'orientation du repère, j'ai considéré un repère "DROITIER", à savoir que :
- une rotation d'un angle positif fait tourner un point dans le sens inverse des aiguilles d'une montre. C'est le sens traditionnel en mathématiques.
- si on se figure le repère avec les 3 doigts de la main DROITE, le pouce figure l'axe des x, l'ongle vers les x croissants; l'index l'axe des y (l'ongle vers les y croissants) et le majeur l'axe des z (l'ongle vers les z croissants). Ce n'est qu'une question de convention, et il se trouve des gens pour utiliser un repère gaucher.

Il y a en gros 3 façons de représenter une rotation générale :
1 - les matrices de rotation 3x3 :
les matrices 3x3 permettent de représenter les rotations générales. Elles permettent d'optimiser les calculs, c'est pourquoi je ne traiterais que ce sujet.
2 - les matrices généralisées 4x4 : les librairies professionnelles utilisent des matrices 4x4 pour représenter les transformations géométriques. La vitesse de calcul est dégradée mais en contrepartie, il est possible d'inclure la translation (et d'autres transformations comme la projection 3D->2D) dans la matrice.
3 - les quaternions : c'est une extension des complexes.

Les déplacement des objets 3D peuvent se décomposer en rotations et translations. Dans un simulateur de vol par exemple, les maisons sont translatées lorsque l'avion avance et tournent lorsqu'il vire sur l'aile (l'observateur est dans le cockpit). Mais avant de passer aux rotations proprement-dites, on va parler brièvement de la translation.

TRANSLATION

Un petit point sur la translation car ça va nous servir dans la suite : translater un objet c'est le déplacer en ligne droite d'une distance donnée. Comme notre objet ne se déforme pas, pour faire celà on va déplacer tous les points de l'objet de la distance et dans la direction spécifiée. Ceux qui ont fait un peu de géométrie savent que ce déplacement peut être caractérisé par un vecteur : sa direction et sa longueur (sa norme) représentent le déplacement à faire et il n'y a qu'à ajouter ce vecteur à chaque sommet des polygones de l'objet. Pour être mathématiquement correct, on ajoute pas le point P au vecteur D pour le translater mais le vecteur OP au vecteur D. Et l'addition vectorielle c'est l'ajout des composantes des vecteurs.

Par exemple (1, 2, 3) + (4, 5, 6) = (1+4, 2+5, 3+6) = (5, 7, 9).

Chaque composante du vecteur translation représente le déplacement à faire sur les 3 axes :
en x, en y et en z.

Par exemple un cube est représenté par les 8 points :

P1 = (-5, -5, -5)
P2 = (5, -5, -5)
P3 = (5, 5, -5)
P4 = (-5, 5, -5)
P5 = (5, -5, 5)
P6 = (-5, -5, 5)
P7 = (-5, 5, 5)
P8 = (5, 5 , 5)

C'est un cube de 10 de coté et centré sur (0, 0, 0). On va translater le cube de manière à ce que le point P1 se retrouve en (5, 5, 5). Le déplacement en x de P1 sera 10, 10 également en y et en z.

Le vecteur déplacement sera donc (10, 10, 10). Une autre façon de le dire : le vecteur déplacement D est la différence des points d'arrivée et de départ :

D = (5, 5, 5) - (-5, -5, -5) = (10, 10, 10)

Ensuite, il ne reste plus qu'à ajouter D à tous les points du cube, ce qui donne :

A1 = P1 + D = (5, 5, 5)
A2 = P2 + D = (15, 5, 5)
A3 = P3 + D = (15, 15, 5)
A4 = P4 + D = (5, 15, 5)
A5 = P5 + D = (15, 5, 15)
A6 = P6 + D = (5, 5, 15)
A7 = P7 + D = (5, 15, 15)
A8 = P8 + D = (15, 15 , 15)

Plus généralement, pour déterminer ce vecteur translation on cherche 2 points qui vont caractériser le déplacement: un point de départ et un point d'arrivée. Dans l'exemple du cube, on pouvait partir du principe que le centre du cube passe de (0, 0, 0) à (10, 10, 10).

Le vecteur translation est la différence donc :

D = (10, 10, 10) - (0, 0, 0) = (10-0, 10-0, 10-0) = (10, 10, 10).

Et au niveau de la complexité de calcul, la translation d'un point d'un objet c'est 3 additions de précalcul (pour déterminer le vecteur déplacement) et ensuite 3 additions par point translaté.

ROTATION AUTOUR DE (0,0,0)

Un peu de théorie

Pour commencer on va faire la rotation d'un point P en 2D, cela va nous donner les formules de base.
On a donc un point P de coordonnées (x,y). Une autre façon de le représenter est de considérer sa distance d à l'origine O et l'angle b entre PO et x'Ox. Les vieux souvenirs de trigonométrie permettent de passer de la représentation polaire (d,b) à la représentation cartésienne (x,y) : x = d*cos(b) y = d*sin(b)

On va donc faire subir à notre point P une rotation d'un angle a autour de O. La représentation polaire du point une fois tourné est : P' = (d, a+b) car sa distance à l'origine n'a pas changé, il a tourné comme sur un cercle de centre O et de rayon d et l'angle entre P'O et x'Ox a été augmenté de l'angle a.
La représentation cartésienne du point P' est (on reprend les mêmes formules) :

x' = d*cos(a+b)
y' = d*sin(a+b)

Comme on ne connait que x,y et pas d,b il faut s'en débarrasser dans le résultat.
Pour cela on développe :

sin(a+b) et cos(a+b) :
sin(a+b) = sin(a)*cos(b) + sin(b)*cos(a)
cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)

Puis on remplace :

x' = d*cos(b) * cos(a) - d*sin(b) * sin(a)
y' = d*sin(b) * cos(a) + d*cos(b) * sin(a)

Et comme :

x = d*cos(b) y = d*sin(b)

Ca nous donne le résultat attendu : Bon, on sait faire tourner un point d'un angle donné autour de l'origine en 2D. Comment ça se passe en 3D ? Et bien on ne fait plus tourner autour d'un point mais d'une droite.

Commençons par les cas simples :

Rotation autour de l'axe z'Oz :
Dans ce cas, la composante z du point reste identique et les composantes x et y se calculent comme on vient de le faire en 2D donc:
Rotation autour de l'axe x'Ox :
Dans ce cas, la composante x du point reste identique et les composantes y et z se calculent comme on vient de le faire en 2D en transposant :
Rotation autour de l'axe y'Oy :
Dans ce cas, la composante y du point reste identique et les composantes x et z se calculent comme on vient de le faire en 2D en transposant :

Arrivé à ce stade, on peut continuer et en déduire les formules nécessaires pour tourner autour de n'importe-quelle droite. La rotation proprement dite n'est pas plus compliquée ensuite.

Mais les demomakers s'arrêtent là. L'avantage est qu'il est plus facile (vraisemblablement) de gérer ainsi les mouvements simples d'objets. De toutes façons nous verrons plus loin comment faire dans ce cas général.
Donc comment faire quand on en reste là ? Et bien on défini un "centre de rotation" de l'objet.
Pour un avion se sera le milieu des ailes par exemple. Puis on fait subir à l'objet 3 rotations à suivre:

1 - autour de la droite parallèle à x'Ox passant par ce "centre"
2 - autour de la droite parallèle à y'Oy passant par ce "centre"
3 - autour de la droite parallèle à z'Oz passant par ce "centre"

L'ordre de ces trois rotations est bien sûr capital et c'est celui que la plupart retiennent. On a donc a spécifier à notre routine de rotation 3 angles et un "centre".
Mais on ne sait pas tourner autour d'une autre droite qu'un axe avec les formules données plus haut ! Qu'à cela ne tienne, on procède en 3 temps :

1 - on translate l'objet de telle sorte que son "centre de rotation" se retrouve à l'origine
2 - on effectue les 3 rotations autour des 3 axes
3 - on translate en sens inverse notre objet roté de manière à ce que son "centre de rotation" se retrouve à sa place initiale. Et le tour est joué !

Un peu de pratique

Un code en C qui réalise la rotation "comme les demomakers" pourrait être :

/* En entrée on a les 3 angles ax,ay,az, les coordonnées du centre Cx,Cy,Cz ainsi que le point à roter Px,Py,Pz.
En sortie on a les coordonnées du point roté dans Rx,Ry,Rz. */
/* Calcul des constantes */
sinx = sin(ax); siny = sin(ay);
sinz = sin(az); cosx = cos(ax);
cosy = cos(ay); cosz = cos(az);
/* 1ère translation pour amener le "centre" à l'origine */
x = Px - Cx;
y = Py - Cy;
z = Pz - Cz;
/* Rotation autour de l'axe x'Ox */
t = y*cosx - z*sinx; z = y*sinx + z*cosx;
y = t;
/* Rotation autour de l'axe y'Oy */
t = x*cosy + z*siny;
z = z*cosy - x*siny;
x = t;
/* Rotation autour de l'axe z'Oz */
t = x*cosz - y*sinz;
y = x*sinz + y*cosz;
x = t;
/* Translation pour ramener le centre à sa place */
Rx = x + Cx; Ry = y + Cy;
Rz = z + Cz;

Celà fait 12 multiplications et 12 additions. Il est facile de gagner 3 multiplications sans changer grand chose à notre routine. Il suffit pour cela d'utiliser une optimisation très courante en traitement du signal : au lieu de faire :

x' = x*cos(a) - y*sin(a)
y' = x*sin(a) + y*cos(a)

qui prend 4 multiplications, on fait :

x'=cos(a)*(x+y)-y*(sin(a) +cos(a))
y'=x*(sin(a) -cos(a))+cos(a)*(x+y)

Ce qui nous ramène à 3 muls/3 add au lieu de 4 muls/2 add.
Il existe une autre optimisation qui donne aussi 3 muls/3 add mais elle est vraisemblablement moins pratique (à voir ?). En tout cas, voici le petit bijou d'Oscar Buneman en 1973 :

t=x-y*tg(a/2)
y'=y+t*sin(a)
x'=t-y'*tg(a/2)

Résumons-nous : on sait faire la rotation des demomakers en 9 muls /15 add (si les 15 additions vous effraient, nous allons voir ci-après comment en gagner 3).
Ca ne demande aucun précalcul si ce n'est gérer des tables de sinus/cosinus (ou sinus/tangente pour l'optimisation d'O. Buneman).
Si on a pas besoin de faire les 3 rotations et que 2 ou une seule suffisent, c'est certainement la meilleure façon de procéder. Mais si le problème est un peu plus complexe; par exemple modéliser la terre qui tourne sur elle-même et autour du soleil dans le même temps, il va nous falloir réaliser 2 fois cette rotation pour chaque point ? Non bien sûr, mais il va nous falloir passer par le calcul matriciel pour composer les rotations.

Grand Merci Pierre Larbier pour son article parue dans le Reporter 9

C'est un petit début mais, si vous aussi, avez des articles, ou des suggestions, je serai ravi, de les mettre sur cette page. Contactez moi !