Le gain le plus important de la méthode matricielle pour calculer les rotations réside dans la grande facilité à gérer des mouvements complexes. Mais nous ne verrons cela que dans les chapitre suivants. Pour le moment, on va se polariser sur les avantages de calculs immédiats.
En l'état, c'est marginalement meilleur que la méthode initiale mais cela permet de gagner quelques multiplications. Moins de 9 muls/point c'est possible ? Et oui, on peut descendre à 6.
Traditionnellement, la multiplication est une opération longue; moins que la division, mais tout de même largement plus que l'addition. Remplacer des multiplications par des additions est donc souvent payant en terme de vitesse de calcul, même si le nombre total d'opérations effectuées est supérieur.
C'est de moins en moins vrai à mesure que les processeurs deviennent plus performants et plus du tout sur les DSP. Mais pour ce qui nous occupe (les 80x86 en entier), chercher à gagner des multiplications peut encore être avantageux.
Une des opérations courante qui souffre de la lenteur des multiplications est le produit de matrices car le produit de 2 matrices NxN demande normalement N^3 multiplications.
En 1967, Winograd a encore surpris son monde (ce n'était pas la 1ère ni la dernière fois !) en trouvant une optimisation permettant de diviser grosso-modo par 2 le nombre de multiplications nécessaires au produit de matrices :
[C] = [A]*[B] , [A] une matrice mxn et [B] une matrice nxs
----- -----
a(i) = \ A(i,2*j)*A(i,2*j-1) b(k) = \ B(2*j-1,k)*B(2*j,k)
/ /
----- -----
1<=j<=n/2 1<=j<=n/2
c(i,k) = 0 si n est pair, A(i,n)*B(n,k) si n est impair.
d(i,k) = c(i,k) - a(i) - b(k)
-----
C(i,k)=d(i,k) + \ [A(i,2*j)+B(2*j-1)]*[A(i,2*j-1)+B(2*j,k)]
/
-----
i<=j<=n/2
Cela demande précisément ceil(n/2)*m*s + floor(n/2)*(m+s) multiplications, soit pas loin de n*m*s/2.
Pour être juste, il faut dire que le nombre d'additions est considérablement augmenté et les mesures de l'époque ne donnaient cette méthode rentable que pour des matrices plus grandes que 20x20.
Mais la voie était ouverte et de nombreuses recherches ont été entreprises pour réduire encore le nombre de multiplications. Les résultats obtenus sont en grande partie très théoriques et ne sont pas utilisables dans notre cas.
Par conséquent, nous ne pouvons guère nous servir que de la méthode de Winograd pour tenter d'optimiser nos rotations. Mais quel rapport avec des rotations justement ?
Et bien, il suffit de voir les points à transformer comme les colonnes d'une matrice. Si on a N points à roter, on peut constituer une matrice [P] de taille Nx3 composée des points et faire le calcul :
[T] = [R]*[P]
Les points rotés se retrouvent comme les colonnes de la matrice [T] (également par conséquent de taille 3xN).
Le nombre de multiplications en utilisant la méthode de Winograd est :
ceil(3/2)*3*N + floor(3/2)*(3+N) = 6*N+N+6 = 7*N+6.
Si le nombre de points est grand, celà tend vers 7 multiplications par point.
Et comme on peut pré-calculer une multiplication parmi les 7, on arrive au résultat de 6 multiplications par point (merci Cédric).
Décortiqué, celà donne ceci :
Précalcul du vecteur A associé à la matrice de rotation :
A(1) = R(1,2)*R(1,1)
A(2) = R(2,2)*R(2,1)
A(3) = R(3,2)*R(3,1)
Précalcul associé à chaque point P :
B = Px*Py
Pour chaque point P calculer T = [R]*P :
Tx = [R(1,2)+Px]*[R(1,1)+Py] - A(1) - B + R(1,3)*Pz
Ty = [R(2,2)+Px]*[R(2,1)+Py] - A(2) - B + R(2,3)*Pz
Tz = [R(3,2)+Px]*[R(3,1)+Py] - A(3) - B + R(3,3)*Pz
soit 6 muls/15 add.
Je suis arrivé là, comme beaucoup de gens je pense, en regardant en détail le merveilleux livre de Donald Knuth "Seminumerical Algorithms".
A ma connaissance, personne n'a publié mieux pour traiter le cas général.
Dans un cas particulier, on peut faire un peu mieux :
si il n'est pas nécessaire de disposer des coordonnées rotées Tx,Ty,Tz mais seulement des rapports Tx/Tz et Ty/Tz correspondants à la projection sur l'écran, on peut gagner une multiplication. En effet, au moins un des 9 éléments de la matrice de rotation n'est pas nul (il y a même au moins un élément non nul par ligne ou par colonne). Il suffit de diviser les 9 éléments de la matrice par cet élément non nul pour forcer un "1" dans la matrice et gagner une multiplication.
Prenons l'exemple de R(3,3). L'application de la méthode de Winograd dans ce cas nous amène à :
Précalcul de la nouvelle matrice de transformation :
R(i,j) = R(i,j)/R(3,3) pour i,j compris entre 1 et 3
Précalcul du vecteur A associé à la matrice de rotation :
A(1) = R(1,2)*R(1,1)
A(2) = R(2,2)*R(2,1)
A(3) = R(3,2)*R(3,1)
Précalcul associé à chaque point P :
B = Px*Py
C = Pz - Px*Py
Pour chaque point P calculer T = [R]*P :
Tx = [R(1,2)+Px]*[R(1,1)+Py] - A(1) - B + R(1,3)*Pz
Ty = [R(2,2)+Px]*[R(2,1)+Py] - A(2) - B + R(2,3)*Pz
Tz = [R(3,2)+Px]*[R(3,1)+Py] - A(3) + C
soit 5 muls/14 add.
On a obtenu le point roté avec un facteur d'échelle (R(3,3) dans l'exemple). De ce fait, les rapports Tx/Tz et Ty/Tz ne sont pas modifiés. Mais celà interdit d'utiliser Tz pour faire du tri en z par exemple. Si on ne travaille qu'avec un objet convexe (ou presque) ou si on utilise le BSP pour déterminer l'ordre d'affichage des facettes, cette idée est peut-être exploitable.
Celui qui voudra se lancer dans le codage de ces optimisations sera certainement heureux de noter que tous les éléments de la matrice de rotation [R] sont inférieurs à 1 en valeur absolue. Il est par conséquent possible de travailler en virgule fixe avec les valeurs Px,Py et Pz « avant la virgule » et les termes R(i,j) ou leur opposé « après la virgule ». Associé à la capacité des processeurs Intel de faire des multiplications 32 bits x 32 bits résultat sur 64 bits, une telle technique est peut-être rentable.
Pendant qu'on est dans les optimisations "originales", il est possible de réduire à 10 le nombre de multiplications nécessaires au précalcul des 9 coefficients de la matrice de rotation. Voici comment :
cox = cos(a)
coy = cos(b)
coz = cos(c)
coxy = cos(a + b)
coxz = cos(a + c)
six = sin(a)
siy = sin(b)
siz = sin(c)
sixy = sin(a + b)
sixz = sin(a + c)
R(1,1) = coy*coz
R(2,1) = coy*siz
R(3,1) = -siy
R(1,2) = siy*six
R(2,2) = cox*coz
R(3,2) = six*coy
R(1,3) = sixy - R(3,2)
R(2,3) = six*coz
R(3,3) = coxy + R(1,2)
t1 = R(1,2)
R(1,2) = R(1,2)*coz - sixz + R(2,3)
t2 = R(1,3)
R(1,3) = R(1,3)*coz - coxz + R(2,2)
R(2,2) = t1*siz + R(2,2)
R(2,3) = t2*siz - R(2,3)
Voilà, 10 muls/12 add à comparer avec 16 muls/4 add de la méthode directe.
C'était ma propre optimisation, mais sur le point de rendre l'article, Cédric (Rixed) Cellier, toujours lui; a mis sur le BBS Eden un petit texte d'un de ses amis , HLI donnant une méthode largement plus efficace :
Tout d'abord, on précalcule ces angles:
A = a + b B = a - b C = a + c D = a - c E = b + c F = b - cEt on pose:
W = Sin(b)Je suppose, bien entendu, que Cos et Sin sont entrées dans une table pour une obtention plus rapide de leurs valeurs.
On précalcule également:
W(a) = Sin(a)*Cos(c) = (Sin(C) + Sin(D))/2 W(b) = Cos(a)*Cos(c) = (Cos(C) + Cos(D))/2 W(c) = Sin(a)*Sin(c) = (Cos(D) - Cos(C))/2 W(d) = Cos(a)*Sin(c) = (Sin(C) - Sin(D))/2Voici alors les coefficients de la matrice de rotation R:
R(1,1) = Cos(b)*Cos(c) = (Cos(E) + Cos(F))/2 R(1,2) = Sin(a)*Sin(b)*Cos(c) - Cos(a)*Sin(c) = W(a)*W - W(d) R(1,3) = Cos(a)*Sin(b)*Cos(c) - Sin(a)*Sin(c) = W(b)*W + W(c) R(2,1) = Cos(b)*Sin(c) = (Sin(E) - Sin(F))/2 R(2,2) = Sin(a)*Sin(b)*Sin(c) - Cos(a)*Cos(c) = W(c)*W + W(b) R(2,3) = Cos(a)*Sin(b)*Sin(c) - Sin(a)*Cos(c) = W(d)*W - W(a) R(3,1) = -W R(3,2) = Sin(a)*Cos(b) = (Sin(A) + Sin(B))/2 R(3,3) = Cos(a)*Cos(b) = (Cos(A) + Cos(B))/2Soit 4 mul/ 18 add/ 1 neg! (et 8 sar! si on tient à en faire à ce niveau là)
Note: quand je dis "add", je ne fais pas la distinction entre "add" est "sub".
Je peux vous dire que, bien codé en assembleur, on peut faire tenir le tout (non pas simultanément, mais en respectant les bons étapes de calcul qui évite de recalculer sottement!) dans les registres 32 bits des 386+. Résultat: sur les 18 add, 14 sont avec des opérandes registres - soit un 1 cycle-machine sur 486 pour chacun!
hum... si cela vous tente, on peut également se débarrasser des 4 mul sous forme d'additions:
On précalcule également:
I = C + b = a + c + b J = C - b = a + c - b K = D + b = a - c + b L = D - b = a - c - b W(a) = (Sin(C) + Sin(D)*2 W(b) = (Cos(C) + Cos(D)*2 W(c) = (Cos(D) - Cos(C)*2 W(d) = (Sin(C) - Sin(D)*2 W(i) = Sin(I) - Sin(J) W(j) = Cos(J) - Cos(I) W(l) = Cos(L) - Cos(K) W(k) = Sin(K) - Sin(L)On remarque que:
W(a)*W - W(d) = ([Sin(b)-1]*Sin(C) + Sin(D)*[1+Sin(b)])/2 = ([Cos(J)-Cos(I)-2*Sin(C)] + [Cos(L)-Cos(K)+4*Sin(D)])/4 = [W(j) + W(l) - 4*W(d)]/4 W(d)*W - W(a) = ([Sin(b)-1]*Sin(C) - Sin(D)*[1-Sin(b)])/2 = ([Cos(J)-Cos(I)-2*Sin(C)] - [Cos(L)-Cos(K)+2*Sin(D)])/4 = [W(j) - W(l) - 4*W(a)]/4 W(b)*W + W(c) = ([Sin(b)-1]*Cos(C) + Cos(D)*[1+Sin(b)])/2 = ([Sin(I)-Sin(J)-2*Cos(C)] + [Sin(K)-Sin(L)+2*Cos(D)])/4 = [W(i) + W(k) + 4*W(c)]/4 W(c)*W + W(b) = ([Sin(b)-1]*Cos(D) - Cos(C)*[1+Sin(b)])/2 = ([Sin(K)-Sin(L)+2*Cos(C)] - [Sin(I)-Sin(J)-2*Cos(D)])/4 = [W(k) - W(i) + 4*W(b)]/4On obtient finalement 0 mul/ 34 add/ 1 neg (8 sar, 4 shl) - ou 0 mul/ 38 add/ 1 neg (8 sar), selon son goût - au lieu des 4 mul/ 18 add/ 1 neg (8 sar) pour calculer les coefficients de la matrice de rotation R!
[...]
Mais quelque soit la réduction du nombre des multiplications, la question de la rentabilité de ces méthodes se pose lorsque l'on est devant son assembleur favori ... La réduction du nombre de multiplications n'est pas tout.
En fait, il est bien plus interressant de réduire carrément le nombre d'objets (points ou vecteurs) à roter. C'est ce que nous allons voir de suite en commençant par le problème des normales.
Même si il existe des techniques pour se passer des normales dans le 1er cas, il est interressant de se demander si on peut se passer complètement de les roter avec les points qui composent l'objet.
On peut répondre par l'affirmative à cette question. Pour décrire pourquoi, on va prendre 2 exemples : le calcul de visibilité lorsque l'observateur est à une distance finie et le calcul de l'éclairement lorsque la source lumineuse est située à l'infini. Toutes les autres situations se déduisent de ces 2 là.
La méthode "traditionnelle" consiste à effectuer la rotation des points de l'objet (et de ce fait le point P qui va nous servir pour le calcul de V) ainsi que celle des normales. Une fois cela effectuée, on calcule et teste le produit scalaire.
L'idée est de procéder en se posant la question suivante : où serait l'observateur si au lieu de roter/translater l'objet, on le déplacait à l'envers ? Son point de vue sur l'objet serait le même que si il restait fixe avec un objet qui tourne/translate mais cela ne fait qu'un seul point à traiter au lieu de tous ceux de l'objet (avec les normales en particulier).
On se le fait mathématiquement ?
Le point P est transformé comme on l'a vu plus haut en [R]*P + D. Le vecteur normal N est transformé par la rotation seule (les vecteurs sont invarients par translation) donc il devient [R]*N. Et on cherche le produit scalaire ps des deux vecteurs, soit :
ps = {([R]*P + D) - O} . {[R]*N}
Or :
De ces 2 points on déduit qu'on peut appliquer à nos 2 vecteurs la rotation inverse de [R] sans changer le résultat de leur produit scalaire. Ce que nous allons faire de suite :
-1
(pour alléger, j'écris [S] = [R]
ps = [S]*{([R]*P + D) - O} . [S] * {[R]*N}
ps = {([S]*[R]*P + [S]*D) - [S]*O} . {[S]*[R]*N}
et comme [S]*[R] = [I], le calcul se simplifie en :
ps = {P + [S]*(D-O)} . N
Ce résultat est exactement ce qu'on espérait, P et N n'ont plus besoin d'être transformés, seul O l'est une fois...
Il reste un détail pratique, comment calculer [R]-1 ??
J'ai écrit plus haut qu'on n'inverserait pas de matrice, alors je sors de mon chapeau une propriété difficilement intuitive, l'inverse d'une matrice de rotation est sa transposée (ses lignes et colonnes sont échangées). Attention, c'est loin d'être une règle générale mais les matrices de rotation ont le bon goût d'être orthogonales ...
De ce fait, pour calculer le produit scalaire nécessaire au calcul de visibilité, il suffit de faire :
Si on veut chasser la multiplication, cela réduit le calcul de la visibilité d'une facette à 2 multiplications/facette seulement.
Deux, car N étant maintenant fixe, il fait partie de la description immuable de la facette, et comme on ne désire que le signe du produit scalaire, rien n'empêche de tester le signe de (x*Nx/Nz + y*Ny/Nz ± z) au lieu de (x*Nx + y*Ny + z*Nz).
Il faut juste prendre garde à ce qu'aucune facette n'ait une normale avec Nz = 0. Cela se règle très facilement en pré-rotant l'objet. Au passage, j'ai factorisé sur z mais rien n'empêche d'en prendre un autre.
Après le cas du test de visibilité, on ne va faire qu'une bouchée du calcul de l'illumination avec une source à l'infini.
Pour éviter de roter la normale avec l'objet, il suffit de procéder comme on vient de la faire, on rote le vecteur d'éclairement "à l'envers". Ce qui fait que notre produit scalaire :
ps = {[R]*N} . L
devient pour les même raisons que dans le paragraphe précédent :
Et on peut évaluer ce produit scalaire en 3 multiplications (je réserve pour une autre fois, quelques idées sur ce calcul...) sans avoir besoin d'autre calcul que la rotation inverse de l'unique vecteur lumineux.
Et après les normales, les points ?
Pour un objet de taille raisonnable et pas trop "biscornu", la moitié des facettes sont visibles et donc la moitié seulement des points doivent être translatés/rotés.
Lorsque j'ai exposé l'idée à Cedric (Rixed) Cellier, il m'a répondu que cela posait une difficulté pratique qui réduisait vraisemblablement à néant tout espoir de gain en temps de calcul : chaque point appartient à plusieurs facettes (en moyenne un peu plus de 5 lorsque ce sont des triangles).
Lors du balayage des facettes, il va falloir tester si les points d'une facette donnée ont déjà été calculés. Et malheureusement, ce test devra être fait plusieurs fois (un peu plus de 5 donc pour des facettes triangulaires).
Il pensait par conséquent que ces tests répétés pouvaient être très pénalisant en temps et qu'il était certainement préférable de calculer tous les points par une routine bien optimisée qui tienne dans le cache du processeur.
Cet argument ne me convainc pas, mais j'ai fait un effort dans son sens : on va tenter de minimiser les tests inutiles.
Pour en illustrer l'idée, on va partir d'un volume simple : un cube. Il sera constitué de triangles. Le cube est une figure géométrique trop simple pour montrer le gain de 2 dans le nombre de transformations. En effet, dans le cas général, 7 points sur 8 sont visibles. Néanmoins, il permet de montrer le principe recherché.
Voici une représentation de notre cube "mis à plat" comme on le ferait avec du carton :
H---------G
| 12 / |
| / |
| / 11 |
|/ |
H---------D---------C--------G---------H
| / | \ | / | \ |
| 1 / | \ 4 | 5 / | \ 8 |
| / | \ | / | \ |
| / 2 | 3 \ |/ 6 | 7 \ |
E---------A---------B--------F---------E
| 9 / |
| / |
| / |
| / 10 |
E---------F
Ce cube contient 12 facettes triangulaires numérotées de 1 à 12 et 8 coins de A à H.
Supposons que les facettes 1,2,3,4,9,10 soient visibles. On va faire tourner à la main l'algorithme suivant :
- Pour toutes les facettes du cube (de 1 à 12)
- tester la visibilité de la facette
- si elle est visible :
- pour tous les points de la facette, tester si ils ont
déjà été rotés/translatés. Si ce n'est pas le cas,
effectuer la transformation désirée.
Au résultat, on va tester les triplets de points suivants lors du balayage des facettes triangulaires : HED , DEA, DAB, DBC, AEB, EBF. Soit pour chacun des 7 points visibles :
H 1 test
E 3 tests
D 4 tests
A 3 tests
B 4 tests
C 1 test
F 1 test
Total : 17 pour 7 points affichés donc 10 tests en trop.
Là, on a balayé les facettes dans un ordre arbitraire sans tenir compte des points dont on a déterminé précédement qu'ils ont été rotés ou non.
Il est possible d'être plus malin, par exemple lorsque l'on passe de la facette n°1 à la n°2 puis à la n°3:
L'idée consiste donc à optimiser le parcours des facettes pour n'avoir à chaque fois qu'un seul point à tester. On passera à chaque fois d'une facette à une de ses 3 voisines et on ne testera que le point unique qui reste à déterminer.
L'algorithme de traitement est alors :
- Pour toutes les facettes du cube (de 1 à 12 dans l'ordre précalculé)
- tester la visibilité de la facette
- si elle est visible :
- si la facette précédente dans le balayage est visible :
tester l'unique point et le transformer si
nécessaire
- sinon :
tester les 3 points de la facette et les
transformer si nécessaire.
avec : on teste forcément les 3 points de la première facette si elle
est visible.
Dans l'exemple du cube, un ordre possible des facettes est : 10, 9, 3, 2, 1, 8, 7, 6, 5, 4, 11, 12. Vous pouvez vérifier qu'on passe bien à chaque fois d'une facette à une de ses voisines et qu'on parcours complètement le cube.
Avec ce chemin optimisé, on teste les points suivants : EBF, A, D, E, H et BCD. Soit 10 tests pour 7 points. Il n'y en a plus que 3 en trop.
Je pense qu'avec un balayage de ce type, on ne doit faire que 0.75 test/point de l'objet en moyenne.
Pratiquement, le balayage n'a rien de compliqué. Il suffit de placer les facettes dans le "bon" ordre une fois pour toutes et de placer en tête des 3 points celui qui sera éventuellement à tester.
Donc rien de spécial à faire dans la phase temps-réel.
En revanche, il faut trouver le "bon" ordre avant. Et ça c'est une autre paire de manches. Je suis convaincu (sans en avoir la preuve formelle) que pour les volumes convexes et ceux qui peuvent se transformer en volumes convexes, un tel chemin existe toujours.
Lorsque l'objet possède un "trou", comme un tore par exemple, c'est nettement moins sûr.
Pour moi, la recherche du chemin est encore une question ouverte... Et le gain me parrait suffisant pour que ça vaille la peine de s'y interresser.
Pour continuer notre petit tour sur les rotations, il nous faut voir le cas de plusieurs rotations en cascade; c'est à dire d'un mouvement plus complexe des objets.