[Résolu][cartographie / débutant] : Euler / Quaternion (C)

[Moove]

Bonjour,

Dans le cadre de mes études, je dois créer un programme capable de modéliser un terrain en 3 dimensions, grâce au placement de points dans l'espace. Je souhaite utiliser les angles d'Euler et les quaternions pour calculer les coordonnées de mes points, mais après plusieurs jours à regarder des explications théoriques et mathématiques je me rends compte que je suis incapable d'implémenter les fonctions permettant de placer correctement les points.

Je tiens à préciser également que je suis débutant en mathématiques et que je dois réaliser cet exercice en autodidacte. J'ai donc appris les bases de la trigonométrie pour comprendre à quoi servent les cosinus et sinus d'un triangle rectangle, mais j'ignore si cela suffit à comprendre comment fonctionnent les matrices de rotation et tout ce qui s'en suit. Peut-être quelqu'un est-il en mesure de m'orienter vers des sources "accessibles" aux débutants pour que je puisse consolider mes bases.

Actuellement, mon programme parse les maps données dans les sources, en voici un exemple :



Je peux donc manipuler un tableau de structure coordonnée avec x, y et z. Les espaces entre les valeurs incrémentent les x et les sauts de ligne incrémentent les y, les valeurs entre les espaces représentent l'altitude.

Je peux donc afficher une map de la manière suivante :



et je dois désormais faire une projection isométrique pour un résultat comme celui-ci :



La deuxième partie de l'exercice consiste à implémenter la rotation de cette représentation graphique d'où mon désir de comprendre les matrices de rotations, les angles d'Euler et les quaternions.

Selon vous, qu'est ce qu'il serait important de mettre en place pour pouvoir déplacer correctement mes points ? Un étudiant de mon école m'a donné ces 2 documents et prétend que la solution est dedans, mais je ne suis pas capable de faire le tri de ce qui peut m'être utile.

https://en.wikipedia.org/wiki/Conversio ... ler_angles
https://jeux.developpez.com/faq/math/?page=quaternions

J'ai intégré la fonction pour obtenir le quaternion grâce au radian, maintenant j'aimerais utiliser cette information pour placer mes points sur ma fenêtre.

Je vous remercie d'avoir pris le temps de me lire. Si vous avec besoin de plus de précision je me ferais un plaisir de vous les donner.

Moove.

[GaBuZoMeu]

Bonjour,

Peux-tu préciser le type de rotation qu'on te demande de faire ?
S'agit-il d'une rotation d'axe l'axe des , suivie de toujours la même projection isométrique ?

[Moove]

Je cherche à initialiser ma figure avec une rotation de -45° sur l'axe des z (z pointant vers nous) puis de 40° sur l'axe des x. Ensuite je dois pouvoir à l'aide de hooks créer des raccourcis clavier ;
Par exemple : flèche de droite tourne de 1° sur l'axe des z. flèche du haut tourne de 1° sur l'axe des x.

L'exercice demande un résultat nous sommes libre de faire le type de rotation que l'on veut.
Cependant le mieux serait d'obtenir ce genre de résultat :
https://www.youtube.com/watch?v=X3pcTRCgQF4

Si tu veux je peux te montrer le code qui entoure la fonction qui calcule les coordonnée (celle que je cherche à faire) et les paramètres que je veux lui envoyer.

[GaBuZoMeu]

Je suppose que tu tiens à jour la matrice de rotation qui fait passer des coordonnées de ta structure fixe aux coordonnées qui te servent pour l'affichage : . Pour l'affichage tu gardes seulement les deux premières coordonnées .
Il suffit alors de multiplier à gauche la matrice par la matrice de rotation de 1° par rapport à l'axe des ou par celle de la matrice de rotation de 1° par rapport à l'axe des .
Je ne suis pas sûr que les quaternions te soient d'une grande utilité ici.

[Moove]

D'accord, je vais commencer à mettre ça en place avec ce que j'ai comme ressources pour comprendre et si j'ai des questions je reviendrais les poser ici. Merci d'avoir pris le temps de me répondre.

[GaBuZoMeu]

Avec plaisir. N'hésite pas à demander des précisions en cas de besoin (pour l'aspect mathématique, pas pour l'implémentation en C)..

[GaBuZoMeu]

En y réfléchissant, je pense que je t'ai induit en erreur en te conseillant de toujours multiplier la matrice de rotation à gauche.
Je pense (mais est-ce vraiment le cas ?) que tu cherches à visualiser le terrain depuis un point de vue "à l'infini" dans la direction donnée en coordonnées sphériques par la longitude et la colatitude . La matrice de rotation est alors (avec matrice de rotation d'axe et matrice de rotation d'axe ).
Si tu veux agir sur la longitude, tu multiplies par une matrice de rotation d'axe à droite, si tu veux agir sur la latitude tu multiplies par une matrice de rotation d'axe à gauche.
Et toujours, quand tu as calculé , la projection s'obtient en gardant les deux premières coordonnées.

[GaBuZoMeu]

Peut-être te montrer un code sera plus clair pour toi que les explications ci-dessus. Je bidouille du Python, je ne travaille pas en C, désolé.

Code: Tout sélectionner

import numpy as np
import matplotlib.pyplot as plt

Je commence par créer des données de terrain. La surface est le graphe de la fonction pour entre 0 et 1 ; peu importe. C'est discrétisé sur une grille 10x10.

Code: Tout sélectionner

x= np.arange(-1,1,0.2)
y= np.arange(-1,1,0.2)
x1,y1= np.meshgrid(x,y)
z= np.cos(x1*2* np.pi)*np.exp(-y1**2-x1**2)
T = np.array([[np.array([x[i],y[j],z[i,j]]).reshape(3,1)\
               for j in range(10)] for i in range(10)])

Je définis les matrices de rotation d'axe x et d'axe y, et la matrice de rotation initiale qui va donner le point de vue de départ :

Code: Tout sélectionner

def Rx(t) :
    u=t*np.pi/180
    return np.array([[1,0,0],[0,np.cos(u),-np.sin(u)],\
                     [0,np.sin(u),np.cos(u)]])

def Rz(f) :
    g=f*np.pi/180
    return np.array([[np.cos(g), -np.sin(g),0],\
                    [np.sin(g),np.cos(g),0], [0,0,1]])

Ri=Rx(-40)@Rz(-135)

Je définis la procédure qui prend une donnée de terrain et une matrice de rotation et fournit la vue depuis le point de vue indiqué par la matrice de rotation :

Code: Tout sélectionner

cox = lambda pt : pt[0,0]
coy = lambda pt : pt[1,0]

def proj(Ro,Te) :
    s=Te.shape ; l=s[0] ; c=s[1]
    vueT=np.array([[(Ro@Te[i,j])[:2] for j in range(c)]\
                     for i in range(l)])
    for j in range(c) :
        rg=[vueT[i,j] for i in range(l)]
        X=list(map(cox,rg)) ; Y=list(map(coy,rg))
        plt.plot(X,Y,color="blue")
    for i in range(l) :
        rg=[vueT[i,j] for j in range(c)]
        X=list(map(cox,rg)) ; Y=list(map(coy,rg))
        plt.plot(X,Y,color="red")
    plt.show()

On voit ce que ça donne pour le point de vue initial.

Code: Tout sélectionner

proj(Ri,T)

On tourne un peu (multiplication par une Rz à droite)

Code: Tout sélectionner

R=Ri@Rz(-20)
proj(R,T)

On prend de la hauteur (multiplication par une Rx à gauche)

Code: Tout sélectionner

R=Rx(20)@R
proj(R,T)

[Moove]

Bonjour GaBuZoMeu, merci pour ton retour, je l'ai vu ce matin.

Il va me falloir un peu de temps pour comprendre le code car il y a beaucoup de fonctionnalités que je vais devoir implémenter en C (Python a l'air très puissant). En tout cas je commence à comprendre le fonctionnement des rotations sur les différents axes.

J'ai lu de la documentation qui parlais d'utiliser les matrice4x4 mais est-ce vraiment utile ?

Pour les matrices de rotation sur les différents axes j'ai trouvé ça mais il semble que le signe des sin(θ) est inversé par rapport à ton exemple, est ce dû au système de coordonné de l'axe y qui est décroissant de haut en bas sur un ordinateur ?


Je reviens vers toi quand j'aurais approfondi le sujet mais j'arrive a initialiser ma map sous différents angles.
Sinon oui, je cherche a faire une rotation autour d'un point central ((0, 0) dans ton exemple) avec des tailles de map plus ou moins grande.

[GaBuZoMeu]

Les matrices de rotation de ton scan sont des matrices de rotation d'angle $-\theta$. Bon, ce n'est pas grave, l'important est de s'y retrouver pour faire ce qu'on veut.
La seule chose que j'utilise et qui n'est pas dans le python de base, mais dans numpy, est la multiplication de matrice (le @ du code). En fait on a juste besoin de la multiplication de matrice 3x3 et de la multiplication d'une telle matrice par un vecteur colonne de taille 3. Tu peux très facilement coder ça en C (il y a aussi sans doute des bibliothèques toutes faites).

Je ne vois pas l'intérêt de travailler avec des matrices 4x4. Si on veut faire des déplacements, oui. Mais ici on a juste des rotations.

[GaBuZoMeu]

Pour récapituler : pour chaque nouvel affichage, tu as à multiplier deux matrices 3x3 : la vieille matrice de rotation R et une matrice de rotation d'axe x ou une matrice de rotation d'axe z pour obtenir la nouvelle matrice R. En fait tu n'as pas besoin de coder la multiplication de matrices en général, il suffit de coder la multiplication par ces deux matrices là, qui est nettement plus simple.

Ensuite il te faut multiplier la nouvelle matrice R par chaque vecteur (x,y,z) de ta discrétisation du terrain (la même matrice R sert pour tous les points de la discrétisation, bien sûr). Là encore, vu en plus que tu ne gardes que les deux premières coordonnées du résultat pour l'affichage, tu peux écrire une petite procédure simple

[Sylviel]

Il y a de nombreuses librairies, extrêmement performante (sans doute parmi les plus optimisées au monde) de calcul linéaire, y compris en C++.

Une liste peut se trouver ici.

Je pense que celle qui est la plus couramment utilisée, en bas niveau, c'est BLAS (en Fortran, wrappé dans de très nombreux langages). Quelques infos complémentaires : https://stackoverflow.com/questions/178 ... -and-atlas

Après, vu la facilité des calculs en jeux, n'importe quelle bibliothèque de calcul scientifique sera plus que suffisante.

[Moove]

Merci Sylviel pour les ressources j'aimerais effectivement pouvoir utiliser des bibliothèques aussi efficaces mais je suis limité par l'énoncé :


Merci pour ces précisions GaBuZoMeu, je commence à mieux comprendre.
Je peux donc facilement coder moi même ces fonctions grâce à tes conseils.
Voilà ce que ça donne, seulement il semble manquer des lignes lors de l'affichage mais ce n'est probablement pas dû aux matrices, je le corrigerais demain.
Le code est volontairement aéré pour respecter la norme imposée pour le premier semestre de l'école.

Angles en radian pour initialisation :

Code: Tout sélectionner

t_map   *ft_init_map(char *path)
{
...
   ret->rx = 0.837758;  // 48°
   ret->ry = 0;
   ret->rz = -0.585398; // -30°
...
   return (ret);

Et pour les matrice avec les mêmes calcules que dans ton exemple.
Cependant je pense ne pas avoir utiliser cette ligne que je ne comprend pas très bien.

Code: Tout sélectionner

z= np.cos(x1*2* np.pi)*np.exp(-y1**2-x1**2)

Code: Tout sélectionner

double   **ft_matrix_rx(t_map *map)
{
   double   **ret;

   ret = ft_alloc_matrix();
   if (!ret)
      exit (EXIT_FAILURE);
   ret[0][0] = 1;
   ret[0][1] = 0;
   ret[0][2] = 0;
   ret[1][0] = 0;
   ret[1][1] = cos(map->rx);
   ret[1][2] = -sin(map->rx);
   ret[2][0] = 0;
   ret[2][1] = sin(map->rx);
   ret[2][2] = cos(map->rx);
   return (ret);
}

double   **ft_matrix_rz(t_map *map)
{
   double   **ret;

   ret = ft_alloc_matrix();
   if (!ret)
      exit (EXIT_FAILURE);
   ret[0][0] = cos(map->rz);
   ret[0][1] = -sin(map->rz);
   ret[0][2] = 0;
   ret[1][0] = sin(map->rz);
   ret[1][1] = cos(map->rz);
   ret[1][2] = 0;
   ret[2][0] = 0;
   ret[2][1] = 0;
   ret[2][2] = 1;
   return (ret);
}

double   **ft_multiply_matrix(double **rx, double **rz)
{
   double   **ret;
   int      i;
   int      j;
   int      k;

   ret = ft_alloc_matrix();
   if (!ret)
      exit (EXIT_FAILURE);
   i = -1;
   while (++i < 3)
   {
      j = -1;
      while (++j < 3)
      {
         k = -1;
         while (++k < 3)
            ret[i][j] += rx[i][k] * rz[k][j];
      }
   }
   return (ret);
}

void   ft_calculate_point(t_map *map, t_point **s)
{
   int      i;
   int      j;
   double   **rx;
   double   **rz;
   double   **m3;

   rx = ft_matrix_rx(map);
   rz = ft_matrix_rz(map);
   m3 = ft_multiply_matrix(rx, rz);
   j = 0;
   while (j < map->hgt)
   {
      i = 0;
      while (i < map->wdt)
      {
         map->init[j][i].x = s[j][i].x * m3[0][0] + s[j][i].y * m3[1][0] +  s[j][i].z * m3[2][0];
         map->init[j][i].y = s[j][i].x * m3[0][1] + s[j][i].y * m3[1][1] +  s[j][i].z * m3[2][1];
         map->init[j][i].z = s[j][i].x * m3[0][2] + s[j][i].y * m3[1][2] +  s[j][i].z * m3[2][2];
         map->init[j][i].color = 0x00FF0000;
         i++;
      }
      j++;
   }
}

Et le résultat :

[GaBuZoMeu]

Pas mal, mais il y a quelque chose qui cloche et ça se voit. Regarde les angles de ta visualisation, aux coins de la grille : ils sont droits ! Or tu vois bien sur l'exemple que tu as mis au début qu'ils ne doivent pas être droits.

Le problème vient de la façon dont tu fais les produits matrice de rotation par vecteur.

Par exemple
s[j][i].x * m3[0][0] + s[j][i].y * m3[1][0] + s[j][i].z * m3[2][0]
ne va pas. Ça devrait être
s[j][i].x * m3[0][0] + s[j][i].y * m3[0][1] + s[j][i].z * m3[0][2]
Pareille modification sur les deux autres lignes

Normalement, on fait d'abord la rotation d'axe z (on s'occupe de la longitude), puis on fait la rotation d'axe x (on s'occupe de la colatitude).
Tandis que toi avec ta multiplication matrice-vecteur mal faite, tu fais le contraire : d'abord rotation d'axe x (les angles droits de la grille restent droits, on garde les horizontales et les verticales de la grille dans la projection), puis rotation d'axe z qui ne change pas les angles dans la projection.

En fait tu multiplies par l'inverse de la matrice de rotation que tu as calculée (qui est aussi la transposée de cette matrice).

En faisant la correction, tu devras faire gaffe aux angles : change-les de signe par rapport à ce que tu as fait (toujours cette histoire d'inverse). Tu devras prendre -48° et 30°.

Fais un retour pour voir si ça va mieux.

PS. La ligne que tu as sautée dans mon code n'a aucune importance. Elle servait juste à fabriquer le modèle de terrain que j'ai utilisé (l'espèce de vague que l'on voit).

[Moove]

Parfait, ça a très bien marché :


J'ai -39° pour rx et 45° pour rz et on est identique à l'exemple.
Il me reste à implémenter les bonus :
-Incluez une projection supplémentaire (ex : parallèle ou conique) !
-Implémentez le zoom avant et arrière.
-Implémentez la translation (déplacement).
-Implémentez la rotation (faire pivoter votre rendu).
-Ajoutez un bonus supplémentaire de votre choix.

Je pense pouvoir y parvenir avec la base que j'ai maintenant.

J'aimerais seulement savoir, si je veux rajouter une rotation sur l'axe des y pour faire la dernière rotation je rajoute une matrice ry et je la multiplie avec la matrice R mise à jour ? A quoi ressemblerait-elle ?
Merci !

[GaBuZoMeu]

La projection isométrique proprement dite est une projection parallèlement au vecteur (1,1,1), c'est-à dire correspond à des angles rx de -55° (plus exactement ) et rz de 45°.
Avec d'autre angles, tu as des projections parallèles qui ne sont pas isométriques.

Avec les rotations effectuées jusqu'à présent (d'axe z puis d'axe x), tu pointes toujours sur le centre de la grille et ton point de vue reste "horizontal", autrement dit tu as du lacet et du tangage, mais pas de roulis.

[Moove]

Je me suis rendu compte que la projection fonctionnait avec une erreur dans le code :

Code: Tout sélectionner

         map->init[j][i].x = s[j][i].x * m3[0][0] + s[j][i].y * m3[0][1] + s[j][i].z * m3[0][2];
         map->init[j][i].y = s[j][i].x * m3[1][0] + s[j][i].y * m3[1][1] + s[j][i].z * m3[2][1]; // <-ici
         map->init[j][i].z = s[j][i].x * m3[2][0] + s[j][i].y * m3[2][1] + s[j][i].z * m3[2][2];

Maintenant si je corrige l'erreur j'obtient ça :


Pour l'axe des y j'ai implémenté ça et ça semble fonctionner mais je ne comprends pas comment multiplier les 3 :

Code: Tout sélectionner

double   **ft_matrix_ry(t_map *map)
{
   double   **ret;

   ret = ft_alloc_matrix();
   if (!ret)
      exit (EXIT_FAILURE);
   ret[0][0] = cos(map->ry);
   ret[0][1] = 0;
   ret[0][2] = sin(map->ry);
   ret[1][0] = 0;
   ret[1][1] = 1;
   ret[1][2] = 0;
   ret[2][0] = -sin(map->ry);
   ret[2][1] = 0;
   ret[2][2] = cos(map->ry);
   return (ret);
}

J'ai trouvé ceci et j'ai essayé la fonction ci-dessous mais sans grand succès:

Code: Tout sélectionner

double   **ft_matrix_euler(t_map *map)
{
   double   **ret;

   ret = ft_alloc_matrix();
   if (!ret)
      exit (EXIT_FAILURE);
   ret[0][0] = cos(map->rz) * cos(map->ry);
   ret[0][1] = (cos(map->rz) * sin(map->ry) * sin(map->rx)) - (sin(map->rz) * cos(map->rx));
   ret[0][2] = (cos(map->rz) * sin(map->ry) * cos(map->rx)) + (sin(map->rz) * sin(map->rx));
   ret[1][0] = sin(map->rz) * cos(map->ry);
   ret[1][1] = (sin(map->rz) * sin(map->ry) * sin(map->rx)) - (cos(map->rz) * cos(map->rx));
   ret[1][2] = (sin(map->rz) * sin(map->ry) * cos(map->rx)) - (cos(map->rz) * sin(map->rx));
   ret[2][0] = -sin(map->ry);
   ret[2][1] = cos(map->ry) * sin(map->rx);
   ret[2][2] = cos(map->ry) * cos(map->rx);
   return (ret);
}

[GaBuZoMeu]

Ne vois-tu pas, en regardant le résultat, qu'il y a un demi-tout en trop (ou en pas assez, ce qui revient au même) ?
Ensuite, que cherches-tu à faire ? Limite-toi à des choses que tu maîtrises, plutôt que d'essayer des trucs au petit bonheur la chance.
Tu parlais de pouvoir faire bouger le point de vue en appuyant sur des flèches haut, bas, gauche, droite. Je t'ai indiqué comment mettre à jour la matrice de rotation pour produire cela. L'as-tu fait ?

[Moove]

Oui je le vois et les seuls solutions que j'ai trouvé pour y remédier sont soit d'initialiser l'angle de rotation sur l'axe des x à 55° au lieu de -55°, soit en initialisant la coordonnée des z en négatif.

Je cherchais un moyen de faire la rotation sur les 3 axes (mais ce n'est pas obligatoire).
Et effectivement c'est tout à fait hasardeux (désolé ) mais j'essaye de chercher un peu par moi même, je veux pas passer pour un assister.

Sinon j'ai bien mis en place la mise à jour de la matrice de rotation mais il arrive que le 42 change de sens et qu'on le lise à l'envers, comme si les coordonnée des z c'était inversé. Est ce qu'il faudrait fixer des limites à la rotation pour modifier ce problème ?

Code: Tout sélectionner

int   key_hook(int keycode, t_running *run)
{
...
   if (keycode == 0)
      run->map->r = ft_multiply_matrix(run->map->ri, ft_matrix_rz(0.0174533));
   if (keycode == 1)
      run->map->r = ft_multiply_matrix(ft_matrix_rx(0.0174533), run->map->ri);
   if (keycode == 2)
      run->map->r = ft_multiply_matrix(run->map->ri, ft_matrix_rz(-0.0174533));
   if (keycode == 13)
      run->map->r = ft_multiply_matrix(ft_matrix_rx(-0.0174533), run->map->ri);
   if (keycode == 0 || keycode == 1 || keycode == 2 || keycode == 13)
      run->map->ri = run->map->r;
   return (0);
}

J'ai aussi essayé d'augmenter et diminuer la valeur du radian passer en paramètre initial et de recalculer par rapport aux points de coordonnées initiales pour pouvoir facilement passer la rotation en négatif ou positif lorsqu'on arrive prêt de 3,14 ou -3,14 radian (Je sais pas si je suis claire), mais pareil c'est quelque chose que j'ai lu quelque pars et je ne maitrise pas le sujet.

[GaBuZoMeu]

Je sais pas si je suis claire

Pas trop.

les seuls solutions que j'ai trouvé pour y remédier sont soit d'initialiser l'angle de rotation sur l'axe des x à 55° au lieu de -55°,

Oui, c'est raisonnable compte tenu de l'orientation du repère. Vu que l'axe des va vers la droite, l'axe des vers le bas et l'axe des vers nous, Une rotation d'axe l'axe des d'angle entre 0 et va emmener le point (0,0,1) sur le point ( : qui se projette en : vers le haut de la figure, donc, puisque le est négatif.