Conversion plus rapide des coordonnées cartésiennes en coordonnées sphériques?

Question

J'ai un tableau de 3 millions de points de données d'un accéléromètre 3-axiz (XYZ), et je veux ajouter 3 colonnes au tableau contenant les coordonnées sphériques équivalentes (r, theta, phi). Le code suivant fonctionne, mais semble trop lent. Comment puis-je faire 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))) 

Answer 1

Ceci est similaire à la réponse de Justin Peel, mais en utilisant seulement numpy et en tirant parti de son intégré vectorisation:

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 l'a suggéré dans les commentaires, j'ai changé la définition de l'angle d'élévation de votre fonction d'origine. 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.

Answer 2

est ici un code Cython rapide que je l'ai écrit pour cela:

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 d'arrêt de 62,4 secondes à 1,22 secondes à l'aide des points 3.000.000 pour moi. Ce n'est pas trop moche. Je suis certain que d'autres améliorations peuvent être apportées.

Answer 3

Pour terminer les réponses précédentes, voici une mise en œuvre Numexpr (avec une solution de repli possible 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 tableau, cela permet un facteur de 2 vitesses en hausse par rapport à pur une mise en œuvre Numpy , et serait comparable aux vitesses C ou Cython. La solution numpy présente (lorsqu'il est utilisé avec l'argument ceval=eval) est 25% plus rapide que la fonction appendSpherical_np dans la @mtrw réponse pour les grandes tailles de tableau,

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 pour les petites tailles, appendSpherical_np est 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 

Answer 4

! Il y a encore une erreur dans tout le code ci-dessus .. et ceci est un top résultat Google. TLDR: J'ai testé cela avec VPython, en utilisant atan2 pour theta (elev) est faux, utiliser acos! C'est correct pour phi (azim). Je recommande la fonction sympy1.0 acos (elle ne se plaint même pas d'acos (z/r) avec r = 0).

http://mathworld.wolfram.com/SphericalCoordinates.html

Si nous convertissons que le système de physique (r, thêta, phi) = (r, Elév, azimut) nous avons:

r = sqrt(x*x + y*y + z*z) 
phi = atan2(y,x) 
theta = acos(z,r) 

non optimisé, mais bon code pour système physique droitier:

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])) 

d'autres données de test pour certains quarts de cercle:

[[ 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]] 

je VPython en plus de visualiser facilement des vecteurs:

test = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)