Conversion plus rapide de coordonnées cartésiennes en coordonnées sphériques?
j'ai un tableau de 3 millions de points de données d'un accelléromètre 3-axiz (XYZ), et je veux ajouter 3 colonnes au tableau contenant les coordonnées sphériques équivalentes (r, thêta, phi). Le code suivant fonctionne, mais semble trop lent. Comment puis-je faire de mieux?
import numpy as np
import math as m
def cart2sph(x,y,z):
XsqPlusYsq = x**2 + y**2
r = m.sqrt(XsqPlusYsq + z**2) # r
elev = m.atan2(z,m.sqrt(XsqPlusYsq)) # theta
az = m.atan2(y,x) # phi
return r, elev, az
def cart2sphA(pts):
return np.array([cart2sph(x,y,z) for x,y,z in pts])
def appendSpherical(xyz):
np.hstack((xyz, cart2sphA(xyz)))
4 réponses
c'est similaire à Justin Peel 'S réponse, mais en utilisant juste numpy et en profitant de son vectorisation intégrée:
import numpy as np
def appendSpherical_np(xyz):
ptsnew = np.hstack((xyz, np.zeros(xyz.shape)))
xy = xyz[:,0]**2 + xyz[:,1]**2
ptsnew[:,3] = np.sqrt(xy + xyz[:,2]**2)
ptsnew[:,4] = np.arctan2(np.sqrt(xy), xyz[:,2]) # for elevation angle defined from Z-axis down
#ptsnew[:,4] = np.arctan2(xyz[:,2], np.sqrt(xy)) # for elevation angle defined from XY-plane up
ptsnew[:,5] = np.arctan2(xyz[:,1], xyz[:,0])
return ptsnew
notez que, comme suggéré dans les commentaires, j'ai modifié la définition de l'angle d'élévation de votre fonction originale. Sur ma machine, testant avec pts = np.random.rand(3000000, 3) , le temps est passé de 76 secondes à 3,3 secondes. Je n'ai pas Cython donc je n'ai pas pu comparer le timing avec cette solution.
voici un code Cython rapide que j'ai écrit pour ceci:
cdef extern from "math.h":
long double sqrt(long double xx)
long double atan2(long double a, double b)
import numpy as np
cimport numpy as np
cimport cython
ctypedef np.float64_t DTYPE_t
@cython.boundscheck(False)
@cython.wraparound(False)
def appendSpherical(np.ndarray[DTYPE_t,ndim=2] xyz):
cdef np.ndarray[DTYPE_t,ndim=2] pts = np.empty((xyz.shape[0],6))
cdef long double XsqPlusYsq
for i in xrange(xyz.shape[0]):
pts[i,0] = xyz[i,0]
pts[i,1] = xyz[i,1]
pts[i,2] = xyz[i,2]
XsqPlusYsq = xyz[i,0]**2 + xyz[i,1]**2
pts[i,3] = sqrt(XsqPlusYsq + xyz[i,2]**2)
pts[i,4] = atan2(xyz[i,2],sqrt(XsqPlusYsq))
pts[i,5] = atan2(xyz[i,1],xyz[i,0])
return pts
il a fallu le temps de 62,4 secondes à 1,22 secondes en utilisant 3 000 000 de points pour moi. Ce n'est pas trop mal. Je suis sûr qu'il ya quelques autres améliorations peuvent être apportées.
Pour compléter les réponses précédentes, voici un Numexpr mise en œuvre (avec une possibilité de repli à Numpy),
import numpy as np
from numpy import arctan2, sqrt
import numexpr as ne
def cart2sph(x,y,z, ceval=ne.evaluate):
""" x, y, z : ndarray coordinates
ceval: backend to use:
- eval : pure Numpy
- numexpr.evaluate: Numexpr """
azimuth = ceval('arctan2(y,x)')
xy2 = ceval('x**2 + y**2')
elevation = ceval('arctan2(z, sqrt(xy2))')
r = eval('sqrt(xy2 + z**2)')
return azimuth, elevation, r
pour les grandes tailles de réseaux, cela permet une vitesse de 2 fois supérieure par rapport à la mise en œuvre pure et simple, et serait comparable aux vitesses de C ou Cython. La solution actuelle de numpy (lorsqu'elle est utilisée avec l'argument ceval=eval ) est aussi 25% plus rapide que la fonction appendSpherical_np dans le @mtrw réponse pour la grande taille des matrices,
In [1]: xyz = np.random.rand(3000000,3)
...: x,y,z = xyz.T
In [2]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 397 ms per loop
In [3]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 280 ms per loop
In [4]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 145 ms per loop
bien que pour les plus petites tailles, appendSpherical_np est en fait plus rapide,
In [5]: xyz = np.random.rand(3000,3)
...: x,y,z = xyz.T
In [6]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 206 µs per loop
In [7]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 261 µs per loop
In [8]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 271 µs per loop
! Il y a encore une erreur dans tout le code ci-dessus.. et C'est un résultat haut de Google.. TLDR: J'ai testé cela avec VPython, en utilisant atan2 pour theta (elev) est faux, utiliser acos! Il est correct pour phi (azim). Je recommande le sympy1.0 fonction acos(il ne se plaint même pas d'acos (z/r) avec r = 0 ) .
http://mathworld.wolfram.com/SphericalCoordinates.html
si nous convertissons cela au système de physique (r, thêta, phi) = (r, elev, azimuth) nous avons:
r = sqrt(x*x + y*y + z*z)
phi = atan2(y,x)
theta = acos(z,r)
Non optimisé mais correct code pour droitier système physique:
from sympy import *
def asCartesian(rthetaphi):
#takes list rthetaphi (single coord)
r = rthetaphi[0]
theta = rthetaphi[1]* pi/180 # to radian
phi = rthetaphi[2]* pi/180
x = r * sin( theta ) * cos( phi )
y = r * sin( theta ) * sin( phi )
z = r * cos( theta )
return [x,y,z]
def asSpherical(xyz):
#takes list xyz (single coord)
x = xyz[0]
y = xyz[1]
z = xyz[2]
r = sqrt(x*x + y*y + z*z)
theta = acos(z/r)*180/ pi #to degrees
phi = atan2(y,x)*180/ pi
return [r,theta,phi]
, vous pouvez tester vous-même avec une fonction comme:
test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319]))
autres données d'essai pour certains quadrants:
[[ 0. 0. 0. ]
[-2.13091326 -0.0058279 0.83697319]
[ 1.82172775 1.15959835 1.09232283]
[ 1.47554111 -0.14483833 -1.80804324]
[-1.13940573 -1.45129967 -1.30132008]
[ 0.33530045 -1.47780466 1.6384716 ]
[-0.51094007 1.80408573 -2.12652707]]
j'ai utilisé VPython en plus pour visualiser facilement les vecteurs:
test = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)