Schneller numpy kartesian zu kugelförmiger Koordinatenumwandlung?

Question

Ich habe ein Array von 3 Millionen Datenpunkten von einem 3-Achsen-Beschleunigungsmesser (XYZ), und ich möchte 3 Spalten zu dem Array hinzufügen, die die äquivalenten Kugelkoordinaten (r, Theta, Phi) enthält. Der folgende Code funktioniert, scheint aber viel zu langsam zu sein. Wie kann ich es besser machen?

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

Dies ist ähnlich Justin Peel ‚s Antwort, aber mit nur numpy und unter Ausnutzung seiner eingebauten in Vektorisierung:

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 

zu beachten, dass, wie in den Kommentaren vorgeschlagen, ich habe die Definition geändert Höhenwinkel von Ihrer ursprünglichen Funktion. Auf meiner Maschine, die mit pts = np.random.rand(3000000, 3) testete, ging die Zeit von 76 Sekunden bis 3,3 Sekunden. Ich habe Cython nicht, so dass ich das Timing mit dieser Lösung nicht vergleichen konnte.

Answer 2

Hier ist ein schneller Cython Code, den ich für diesen aufschrieb:

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 

Es nahm sich die Zeit nach unten von 62,4 Sekunden auf 1,22 Sekunden mit 3.000.000 Punkten für mich. Das ist nicht zu schäbig. Ich bin mir sicher, dass es noch weitere Verbesserungen gibt.

Answer 3

die vorherigen Antworten zu vervollständigen, ist hier eine Numexpr Implementierung (mit einem möglichen Rückfall auf 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 

Für große Feldgrößen, ermöglicht dies einen Faktor von 2 Geschwindigkeit im Vergleich zur reinen einer Numpy Implementierung , und wäre vergleichbar mit C oder Cython-Geschwindigkeiten. Die vorliegende numpy Lösung (wenn sie mit dem ceval=eval Argument verwendet) ist ebenfalls um 25% schneller als die appendSpherical_np Funktion in der @mtrw Antwort für große Array Größen

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 

obwohl für kleinere Größen appendSpherical_np ist eigentlich schneller,

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

! Es gibt einen Fehler immer noch in allen Code oben .. und das ist ein Top-Google-Ergebnis .. TLDR: Ich habe dies mit VPython getestet, mit atan2 für Theta (elev) ist falsch, verwenden Sie acos! Es ist korrekt für Phi (Azim). Ich empfehle die sympy1.0 acos Funktion (es beschweren sich nicht einmal über acos (z/r) mit r = 0).

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

Wenn wir konvertieren, dass in das Physiksystem (r, theta, phi) = (r, elev, Azimut) haben wir:

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

Nicht optimiert aber richtigen Code für rechtshändig Physik-System:

012.351:

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] 

Sie es selbst mit einer Funktion wie testen

test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319])) 

einige andere Testdaten für einige Quadranten:

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

I VPython verwendet zusätzlich Vektoren leicht sichtbar zu machen:

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