Calculer la matrice de projection de WebGL

Tout développeur qui s'initie à la 3D avec WebGL est rapidement confronté à une entité mathématique, la matrice de projection :
[ Xe Ye Ze 1 ] x [ 2NR-L 0 0 0 0 2NT-B 0 0 R+LR-L T+BT-B F+NN-F -1 0 0 2FNN-F 0 ]
Les différentes constantes utilisées pour définir le volume de clipping se présentant par exemple ainsi (coupe Y/Z) :
Pyramide tronquée de clipping dans WebGL
Cette matrice peut sembler d'autant plus étrange au développeur qu'elle comporte 4 lignes et 4 colonnes (une matrice 4x4), alors qu'il semblerait que l'on puisse se contenter d'une matrice 3x3 pour réaliser les calculs. D'ailleurs, où est passée l'intervention de Ze au dénominateur dans les calculs ?
Bref, pourquoi cette matrice ? Comment calculer ses coefficients ? Et comment en déduire les coordonnées de la projection d'un point ?
A écouter en lisant l'article...

Pourquoi une matrice 4x4

Soit le point (Xe,Ye,Ze). Le principe des coordonnées homogènes est qu'au terme de la multiplication d'une matrice [XeYeZe1] par une matrice 4x4 de projection, on obtienne une matrice [XcYcZc1] telle que (XcWc,YcWc,ZcWc) constituent les versions normalisées des coordonnées la projection du point (Xp,Yp,Zp) :
-1XcWc1 si LXpR
-1YcWc1 si BYpT
-1ZcWc1 si N-ZeF
Sans aller jusqu'à dire que recourir à des coordonnées homogènes ne soit justifié que par cela, il faut noter que cette matrice présente deux intérêts :
  • Tout d'abord, tous ses éléments sont indépendants des coordonnées du point, puisque la division par une grandeur dépendant de Ze, sans laquelle aucune projection ne serait possible, est repoussée à une étape ultérieure. Dans ces conditions, il sera possible de projeter une série de points après avoir calculé les coefficients de la matrice une seule fois.
  • Ensuite, le résultat de la multiplication par la matrice suivi de la division par Wc sera normalisé. Pourquoi normaliser ? Car cela permettra de savoir rapidement si le point transformé se trouvait dans le volume de clipping. Si au moins une des coordonnées de sa projection ne respecte pas les inégalités correspondantes, le point était en dehors ce volume, ce qui doit entraîner découpage et/ou exclusion des surfaces dont il est un sommet.
Pour rappel, le volume de clipping n'est pas un cube : c'est une pyramide tronquée. En effet, il faut bien noter que les inégalités concernent Xp et Yp, et non Xe et Ye. En imposant des contraintes sur les coordonnées de la projection du point et non du point, on signifie bien que seuls les points dont les projections tiennent dans la surface de projection, et qui se trouvent entre telle et telle profondeur, doivent pouvoir être retenu pour l'affichage en testant de simples inégalités au terme des calculs.
Pourquoi -Ze ? Car dans WebGL, l'axe des affixes est orienté vers l'observateur, qui est positionné à son origine. Le plan de projection est confondu avec le plan de clipping avant (situé à une distance N, exprimée positivement). Le plan de clipping arrière est situé à une distance F de l'observateur (exprimée aussi positivement).
La question est de savoir quels doivent être les éléments de la matrice 4x4 sachant que (on le répète), ils doivent être totalement indépendants des coordonnées (Xe,Ye,Ze).
Pour faciliter la lecture, posons que la matrice des coordonnées [XeYeZe1] est un vecteur ligne. La matrice 4x4 recherchée est :
[ Ax Bx Cx Dx Ay By Cy Dy Az Bz Cz Dz Aw Bw Cw Dw ]
Par exemple, on a donc :
Xc=AxXe+AyYe+AzZe+Aw

Les coefficients pour calculer Xc

En appliquant le théorème de Thalès, on détermine que :
XpXe=N-Ze
Xp=N-ZeXe
Partons du résultat. On souhaite finalement que si... (1) :
LXpR
...on calcule en fait XcWc, version normalisée de Xp, telle que :
-1XcWc1
Donc :
-WcAxXe+AyYe+AzZe+AwWc
Partons de (1). Il se trouve que rapporter une valeur V de l'intervalle [L,R] à l'intervalle [-1,1] est assez simple :
L V R 0 V-L R-L 0 V-LR-L 1 0 2(V-L)R-L 2 -1 2(V-L)R-L-1 1 -1 2V-R-LR-L 1
Donc si V=Xp=N-ZeXe :
-1 2N-ZeXe-R-LR-L 1 -1 2N-Ze(R-L)Xe-R+LR-L 1
A partir de là, si Ze0 (ce qui correspond bien à la situation, comme il est possible de le constater sur la figure au début de cet article, la pyramide tronquée résidant dans les affixes négatives), alors il est possible de multiplier par -Ze0 sans inverser les inégalités :
Ze 2NR-LXe+R+LR-LZe -Ze
C'est une couple d'inégalités de la forme qu'on recherche :
-WcAxXe+AyYe+AzZe+AwWc
A condition donc de poser que Wc=-Ze.
Partant, les calculs peuvent être inscrits dans la matrice :
[ Ax Ay Az Aw ] = [ 2NR-L 0 R+LR-L 0 ]
Sans oublier :
[ Dx Dy Dz Dw ] = [ 0 0 -1 0 ]

Les coefficients pour calculer Yc

Même raisonnement que pour les coefficients pour calculer Xc, ce qui débouche sur :
[ Bx By Bz Bw ] = [ 0 2NT-B T+BT-B 0 ]

Les coefficients pour calculer Zc

On sait désormais que Wc=-Ze et on n'y veut rien changer, car ce serait remettre en cause les calculs précédents. On souhaite que si... :
N-ZeF
...alors :
-1ZcWc1
Comme Xc et Yc, Zc prend la forme d'une combinaison linéaire de Xe, Ye et Ze :
Zc=CxXe+CyYe+CzZe+Cw
On ne va surtout pas s'emmerder en rendant Zc dépendant de Xe et de Ye, car on ne voit pas du tout ce que cela pourrait apporter (Ze est totalement indépendant de Xe et Ye !). On va donc d'emblée poser que :
Cx=0
Cy=0
Bref :
Zc=CzZe+Cw
Ainsi, le problème se résume à déterminer Cz et Cw, indépendants de Xe, Ye et Ze, et tels que si -Ze est dans [N,F], alors :
-1CzZe+CwWc1
Sachant que Wc=-Ze, on peut donc développer :
-1 -Cz-CwZe 1 Cz-1 -CwZe Cz+1
A ce stade, il faut faire un choix. En effet, puisqu'il s'agit de diviser par Cw, faut-il considérer qu'il s'agit d'un nombre positif ou négatif ? Dans le premier cas, les inégalités sont conservées ; dans le second, elles sont inversées.
Je n'ai pas trouvé de solution plus rationnelle que de poursuivre les calculs dans l'un et l'autre cas, pour retenir la solution qui convient. C'est très bête. Si quelqu'un trouve mieux, ce qui est définitivement possible, qu'il me le fasse savoir en passant par la page Contact.
Un test d'une matrice de transformation dans Excel à récupérer ici permet de produire ces deux courbes rendant compte de l'évolution de ZcWc quand Ze évolue entre -N et -F :
Normalisation de l'affixe dans la matrice de projection 3D
Normalisation de l'affixe dans la matrice de projection 3D
C'est donc la solution Cw0 qu'il faut retenir. En divisant par Cw, les inégalités sont donc inversées. Les calculs se poursuivent ainsi :
Cz-1Cw -1Ze Cz+1Cw CwCz-1 -Ze CwCz+1
Or ces inégalités doivent se constater quand N-ZeF. On peut donc poser que :
F=CwCz+1
N=CwCz-1
C'est un système de deux égalités à deux inconnues, qui se résout en injectant la première dans la seconde :
N=F(Cz+1)Cz-1
N(Cz-1)=F(Cz+1)
Cz(N-F)=F+N
Cz=F+NN-F
Donc en injectant ce résultat dans la première (ou la seconde, c'est pareil) :
F=CwF+NN-F+1
Cw=F(F+N)N-F+F
Cw=F(F+N)+F(N-F)N-F
Cw=2FNN-F
C'est fini !
[ Cx Cy Cz Cw ] = [ 0 0 F+NN-F 2FNN-F ]
Il faut noter que dans WebGL, la matrice des coordonnées est un vecteur colonne et non ligne. En conséquence, il faut procéder à une transposition :
[ 2NR-L 0 R+LR-L 0 0 2NT-B T+BT-B 0 0 0 F+NN-F 2FNN-F 0 0 -1 0 ] x [ Xe Ye Ze 1 ] = [ Xc Yc Zc Wc ] (Xp=XcWc,Yp=YcWc,Zp=ZcWc)
Comme on peut le constater, les calculs sont guidés par le souci de trouver des coefficients dans la matrice qui sont indépendants de Xe, Ye et Ze.

Pour en savoir plus

Cela faisait fort longtemps que je voulais tirer au clair la manière dont la matrice de projection 4x4 est construite. En 1999, j'ai publié un ouvrage Direct3D, 3D temps réel sous Windows. A l'époque, j'avais cherché des explications complètes sur la manière de calculer cette matrice pour les reproduire dans les annexes. En vain !
En effet, j'avais pu constater que le calcul de la matrice de projection 4x4 était un de ces sujets que les auteurs présentent sans aller jusqu'au bout, parce que finalement, ils n'ont pas compris - ou ne se sont pas donné la peine de comprendre -, ce dont ils parlent. On tombe parfois sur ce genre de difficulté en lisant un ouvrage technique, et c'est extrêmement frustrant. Par exemple, qu'on me trouve un livre de statistiques relativement abordable qui explique pourquoi on calcule l'écart-type comme il se calcule. Tout le monde semble prendre la formule pour une vérité divine, alors qu'elle a bien été inventée par quelqu'un ! Pénible...
Enfin, pour revenir à cette matrice, je suis tout récemment tombé sur une explication qui répond parfaitement à mes attentes, ici puis ici, sur l'excellent site ScratchPixel 2.0 que je ne saurais trop recommander pour la qualité des explications qu'il contient.
Après avoir survolé les explications données sur ce site, je me suis remis à l'ouvrage en adoptant un angle d'attaque moins rigoureux mais plus intuitif pour moi, aboutissant à cette explication. Il n'est jamais trop tard pour régler ses comptes avec soi-même. Par ailleurs, sa rédaction m'aura donné l'occasion de découvrir MathML pour présenter les équations.
Sur ce point, si quelqu'un sait comment coder correctement le signe négatif pour une variable, c'est-à-dire éviter d'utiliser <mo>-</mo> qui génère un gros signe moins et non un petit (comme il est possible de le faire pour un nombre : utiliser <mn>-1</mn> qui produit un beau -1 plutôt que d'utiliser <mo>-</mo><mn>1</mn> qui produit un affreux -1), qu'il me le fasse savoir en passant par la page Contact. Firefox ne reconnaît pas <minus/> pour l'heure.
Calculer la matrice de projection de WebGL