Page 1 sur 1

Rotations en coordonnées sphériques

MessagePosté: Mardi 23 Décembre 2014, 16:45
par FiReTiTi
Bonjour,

je souhaite effectuer les rotations des éléments d'un tableau 3D dans un repère sphérique.
Pour cela, j'ai créé cette fonction à partir des formules [url=http://fr.wikipedia.org/wiki/Coordonnées_sphériques]données ici[/url] :
Code: Tout sélectionner
public static void SphericalRotations(double[][][] array, double theta, double phi, double[][][] result)
   {
   int x, y, z, px, py, pz ;
   int sizex = result[0][0].length ; // Dimensions du tableau résultat.
   int sizey = result[0].length ;
   int sizez = result.length ;
   int ami = array[0][0].length >> 1 ; // Centre du tableau d'origine
   int amy = array[0].length >> 1 ;
   int amz = array.length >> 1 ;
   int resmx = sizex >> 1 ; // Centre du tableau résultat
   int resmy = sizey >> 1 ;
   int resmz = sizez >> 1 ;
   double cx, cy, cz, st, sp, r ;
   
   
   for (z=0 ; z < sizez ; z++)
      for (y=0 ; y < sizey ; y++)
         for (x=0 ; x < sizex ; x++)
            {
            cx = x - resmx ; // Coordonnées cartésienne centrées.
            cy = y - resmy ;
            cz = z - resmz ;
            
            r = Math.sqrt(cx*cx + cy*cy + cz*cz) ; // Transformation en coordonnées sphériques.
            sp = Math.acos(cz / r) ;
            st = Math.atan2(cy, cx) ;
            
            sp -= phi ; // Rotations en sphérique.
            st -= theta ;
            
            cx = r * Math.sin(sp) * Math.cos(st) ; // Retour en cartésien.
            cy = r * Math.sin(sp) * Math.sin(st) ;
            cz = r * Math.cos(sp) ;
            
            px = ami + (int)(cx+0.5) ; // Vérifications pour éviter les débordements.
            if ( px < 0 || array[0][0].length <= px ) continue ;
            py = amy + (int)(cy+0.5) ;
            if ( py < 0 || array[0].length <= py ) continue ;
            pz = amz + (int)(cz+0.5) ;
            if ( pz < 0 || array.length <= pz ) continue ;
            
            result[z][y][x] = array[pz][py][px] ;
            }
   }

En résumé :
- pour chaque élément de mon tableau
- je calcule les coordonnées sphériques (depuis les coordonnées cartésiennes)
- je fais la rotation en coordonnées sphérique (simple soustraction)
- je repasse en coordonnées cartésiennes
- je copie dans mon résultat la case du tableau de départ correspondante.


Les rotations en fonction de theta (autour de l'axe Z) fonctionnent parfaitement, mais j'ai un résultat totalement aberrant pour la rotation avec phi (inclinaison par rapport à l'axe Z).
Est ce que quelqu'un pourrait m'expliquer mon erreur ?

Re: Rotations en coordonnées sphériques

MessagePosté: Dimanche 28 Décembre 2014, 00:19
par Framboise
Bonjour,

Un peu de détails sur le problème encontré seraient bienvenus.
Le résultat faux, des distorsions, un manque de précision,... ?
Quel langage ( C++ ? ) et compilateur ?

Quelques remarques:
Les équations sources utilisées sont exactes en théorie mais sources de problèmes pour une application "pro" numérique.
On a des pertes de précision dans certaines plages trigo ( -> distorsions ).
Des plantages avec les fonctions ACOS, ASIN pour des arguments du style 1.000 000 000 000 000 1 pouvant parfois résulter des limites numériques de précision des calculs.

Il vaut mieux tout concevoir en passant par des coordonnées rect. -> polaire 3D -> rotation -> polaire 3D -> rect. effectivement.

Les fonctions ACOS, ASIN ( et ATAN, ASEC, ACSC ) ne doivent pas apparaitre sauf s'il fallait les utiliser strictement pour implémenter une fonction ATAN2 ( présente dans notre cas, donc sans ce problème ici si elle est bien conçue ).
Je vise le: sp = Math.acos(cz / r) ; incongru. Il faut choisir ensuite le bon quadrant vu que acos est ambigu sur ce point:
sp = Math.acos(cz / r) ;
sp = 2*PI - Math.acos(cz / r) ;
selon les cas.
Je préfère utiliser atan2 qui n'a pas d’ambiguïté et remplace avantageusement acos, pire encore si l'on rencontre acos (1.000 000 000 000 000 1) sans soucis pour atan2 non plus...

Quelques références:

Jean Meeus
Astronomical Algorithms

Numerical Recipes : The Art of Scientific Computing
http://www.numerical-recipes.com/
LA 2ed est dispo sur Internet en freeware par chapitre.

Re: Rotations en coordonnées sphériques

MessagePosté: Lundi 29 Décembre 2014, 00:43
par FiReTiTi
Bonjour, merci pour la réponse.

C'est un code Java avec une JVM 1.7

Pour les résultats aberrants, lorsque j'essaie d'incliner un segment, je me retrouve avec une sorte de V.

Tu écris "Il vaut mieux tout concevoir en passant par des coordonnées rect. -> polaire 3D -> rotation -> polaire 3D -> rect. effectivement.", mais il me semble que c'est ce que j'essaie de faire.

Comment atan2 peut remplacer ados ?

Re: Rotations en coordonnées sphériques

MessagePosté: Lundi 29 Décembre 2014, 11:49
par Framboise
je me retrouve avec une sorte de V

Ce serait bien un problème de quadrant qui n'est pas le bon selon le signe d'un des paramètres comme je le soupçonne plus haut.
Je ne vois rien dans le code qui traite ce problème du bon quadrant, donc selon les cas cela fonctionne ou provoque un déraillage ( le V ).

Pour passer de acos ou asin à atan2 : [ extrait d'un help fortran ]
ATAN2(Y, X) computes the principal value of the argument function of the complex number X + i Y. This function can be used to transform from Cartesian into polar coordinates and allows to determine the angle in the correct quadrant.


Un coup de magie avec Pythagore et l'on a les 2 côtés du triangle rectangle permettant d'appliquer atan2 au lieu de acos.
Cela évite les risques d'un acos ( 1.0....01 ) et l’ambiguïté de quadrant.

mais il me semble que c'est ce que j'essaie de faire.

Oui, je ne fais que confirmer que la voie est bonne [ d'où le 'effectivement' ].

Note: je ne pratique pas Java, mais le C qui en est très proche.

Re: Rotations en coordonnées sphériques

MessagePosté: Lundi 29 Décembre 2014, 12:27
par FiReTiTi
Merci pour ta réponse.
Je viens donc de remplacer
Code: Tout sélectionner
sp = Math.acos(cz / r)

par :
Code: Tout sélectionner
double d = cz / r ;
sp = Math.atan2(Math.sqrt(1-d*d), d) ;


J'ai malheureusement exactement le même "mauvais" résultat :-(
Je continue à regarder pour ces histoires de quadrants.

Re: Rotations en coordonnées sphériques

MessagePosté: Lundi 29 Décembre 2014, 15:21
par Framboise
Il semble que l'on perde en route l'info du bon quadrant.

Je vais regarder cela.

Un aperçu de mon code en C:

Code: Tout sélectionner
void RP3 ( double *x, double *y, double *z )
{ /* Conv. RECT(x,y,z)->POL(rho,a,d) 3-dim */
   RP( x, y );
   RP( x, z );
   return;
}

void RP( double *x, double *y )
{ /* RECT(X,Y) => POL(RHO,THETA) */
   double t;
   t = ANG( *x, *y );
   *x = RHO( *x, *y );
   *y = t;
   return;
}

double RHO( double x, double y )
{ /* Distance (Pythagore) */
   return sqrt( x * x + y * y );
}


avec ANG ( x, y ) correspondant à atan2 ( y, x )

C'est un fragment de programme à usage astro de calculs de planètes.
Les fonctions RP sont la transposition de cette fonction des calculatrices en RPN ( HP41,... ).

Re: Rotations en coordonnées sphériques

MessagePosté: Mardi 30 Décembre 2014, 01:43
par FiReTiTi
Il semblerait bien que ce soit un problème de quadrant.
Si j'ajoute
Code: Tout sélectionner
if ( cx < 0.0 ) sp = 2.0*Math.PI - sp ;
, alors cela règle le problème dans certains cas, mais pas tous.

Re: Rotations en coordonnées sphériques

MessagePosté: Mardi 30 Décembre 2014, 14:02
par Framboise
cx ... ou cz < 0.0 ?

Mais le problème ne devrait pas apparaitre avec atan2.
C'est :
sp = Math.acos(cz / r) ;
initial qui ne va pas.
sp = atan2( sqrt(cx*cx +cy*cy) , cz )
me semble la solution.

J'ai imaginé la figure:
axes Ox, Oy, Oz
Le point M, sa projection M' sur le plan xOy
cx, cy, cz les proj. sur les axes xyz.
sp = angle zOM "distance zénitale"
st = angle xOM'

donc atan2 doit prendre le triangle rect. cz O M' et non pas cz O M.

Sans la figure cela reste abscons.

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 01:36
par FiReTiTi
Merci pour ton aide, mais malheureusement cela ne fonctionne pas :-(

Tu écris : "donc atan2 doit prendre le triangle rect. cz O M' et non pas cz O M."
Mais si cz est la projection de M sur l'axe Z et M' la projection de M sur le plan xOy, alors cet angle est toujours 90°.

En regardant ces formules ( http://en.wikipedia.org/wiki/Inverse_trigonometric_functions ), j'avais plutôt conclu qu'il fallait remplacer
sp = Math.acos(cz / r)
par
double d = cz / r
sp = Math.atan2(Math.sqrt(1.0-d*d), d)
Mais ça ne fonctionne pas mieux :-(

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 02:12
par FiReTiTi
En fait, les trois formules sont identiques (je viens de tester) :

- 1 - double d = cz / r ;
phi = Math.atan2(Math.sqrt(1.0-d*d), d) ;

- 2 - phi = Math.atan2(Math.sqrt(cx*cx+cy*cy), cz) ;

- 3 - phi = Math.acos(cz / r) ;

Sauf que sur Wikipedia, la formule 3 est couplée avec "if (cy<0) phi = 2*PI-phi". Ce problème ne semble donc pas réglé en utilisant atan2.
De toute façon en ajoutant le test ou non, le résultat est exactement le même.

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 10:08
par kojak
bonjour,


FiReTiTi a écrit:Sauf que sur Wikipedia, la formule 3 est couplée avec "if (cy<0) phi = 2*PI-phi".


Non, ce n'est pas l'angle $\varphi$ mais l'angle $\theta$.

Tu as bien $\varphi=\arccos\dfrac{z}{r}$ ce qui te donnera toujours un angle entre $0$ et $\pi$, alors qu'avec l'arctangente, tu auras n angle entre $]-\pi/2,\pi/2[$ donc ça ne colle pas.

Ton $\theta$ vaut bien $\arccos\dfrac{x}{\sqrt{x^2+y^2}$ si $y\geq 0$ car l'arccosinus te sort un angle entre $0$ et $\pi$ et l'autre $2\pi-\arccos\dfrac{x}{\sqrt{x^2+y^2}$ si $y\leq 0$. comme ça tu auras bien un angle entre $\pi$ et $2\pi$.

Pour finir avec ces coordonnées sphériques, il faut faire attention car il y a plusieurs définitions : celles des matheux et celle des physiciens. Là, j'ai pris ta référence du début.

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 11:54
par Framboise
Mais si cz est la projection de M sur l'axe Z et M' la projection de M sur le plan xOy, alors cet angle est toujours 90°.

Le TR cz O M est bien rectangle, mais c'est l'angle cz O M ( non rect. dans le cas général ) qui nous intéresse et qui doit être fourni par atan2. cz O M' est lui un angle rect. effectivement.

Le choix des angle est un peu inhabituel, on utilise plutôt un système longitude - latitude, alors que là, sp est la distance au pôle N.

Ça change un peu des formules habituelles.

Pour info atan( y, x ) est un alias de atan2( y, x ) en C depuis quelques années.

2 - phi = Math.atan2(Math.sqrt(cx*cx+cy*cy), cz) ;

- 3 - phi = Math.acos(cz / r) ;

ne peuvent pas être identique car r = sqrt(cx*cx+cy*cy + cz*cz)

PS: Je ne trouve pas le Ç majuscule sur mon clavier QWERTZ-CH, j'ai bien le ç avec MAJ 4, mais rien en format majuscule. Une astuce autre que alt + pavé num. ?
En QWERTY-US int., c'est facile avec ' puis C, mais ici cela ne fonctionne pas.

Image

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 16:38
par kojak
Framboise a écrit:Le choix des angle est un peu inhabituel, on utilise plutôt un système longitude - latitude, alors que là
Tout dépend de quel point de vue on se place : matheux, physiciens ou géographes.

PS : sur ton image jointe, ton repère n'est pas direct, et ton cx n'est pas correct, car la projection de $M'$ sur les abscisses n'est pas parallèle à l’axe des ordonnées.

Re: Rotations en coordonnées sphériques

MessagePosté: Mercredi 31 Décembre 2014, 16:59
par Framboise
Hello kojak,

Je me place d'un point de vue astronomique ascension droite et déclinaison mais j'ai préféré la terminologie longitude/latitude mieux connue et similaire.
C'est vrai, ma // est un peu de travers...
L'inversion "miroir": heureusement cela ne change pas le calcul.
Je n'ai pas d'excuse :roll: , je ne bois quasiment jamais d'alcool ( -> migraines :evil: ).

Bonnes fêtes à toi et à tout le monde

Re: Rotations en coordonnées sphériques

MessagePosté: Vendredi 02 Janvier 2015, 17:03
par FiReTiTi
Le problème avec la rotation en Phi est en fait que ce n'est pas une rotation complète : un point est éloigné ou rapproché de l'origine Oz.
Donc actuellement, lorsque je fais une rotation suivant Phi, tout les points se rapprochent ou s'éloignent de l'axe, au lieu de tourner autour du point O.
Pour que ce soit une rotation, il faudrait que Phi soit codé sur [0, 2Pi[, mais je ne pense pas que ce soit possible.

Re: Rotations en coordonnées sphériques

MessagePosté: Vendredi 02 Janvier 2015, 18:12
par FiReTiTi
En fait, si c'est possible. Il faut juste que je fasse une rotation afin que le point à tourner soit dans un plan ( xOz par exemple ), une simple rotation autour de Y, puis retour dans le plan d'origine.