Plus rapide numpy cartésienne sphérique conversion de coordonnées?

Question

I possède un réseau de 3 millions de points de données à partir d'un accellerometer 3-Axiz (XYZ), et je veux ajouter 3 colonnes de la matrice contenant les coordonnées sphériques équivalents (r, thêta, phi). Les travaux de code suivant, 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 à Justin Peel » s réponse, 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, les essais 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 moment avec cette solution.

Answer 2

Voici un code Cython rapide que je l'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 pris le temps de 62,4 secondes vers le bas à 1,22 secondes à l'aide des points 3.000.000 pour moi. Ce n'est pas trop mal. Je suis sûr qu'il ya d'autres améliorations qui peuvent être faites.

Answer 3

! Il y a une erreur encore dans tout le code ci-dessus .. et c'est un meilleur résultat Google .. TLDR: Je l'ai testé avec VPython, en utilisant atan2 pour thêta (Asc) est faux, l'utilisation ACOS! Il est correct pour phi (Azim). Je recommande la sympy1.0 ACOS fonction (il ne se plaint même pas 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 correcte Code pour le système de physique à droite:

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)

Answer 4

Pour compléter 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, ce qui permet un facteur de 2 par rapport à la vitesse jusqu'à une mise en œuvre pur Numpy, et serait comparable à une vitesse C ou cython. La présente solution numpy (lorsqu'elle est utilisée avec l'argument ceval=eval) est 25% plus rapide que la fonction appendSpherical_np dans la réponse @mtrw 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 que 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