Schnell numpy kartesisch zu sphärischer Koordinatenumwandlung?

Question

Ich habe ein Array von 3 Millionen Datenpunkte von einem 3-Axiz accellerometer (XYZ), und ich möchte 3 Spalten zur Anordnung hinzufügen, um die äquivalenten sphärischen Koordinaten, die (r, theta, phi). Der folgende Code funktioniert, aber scheint viel zu langsam. Wie kann ich 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 wie Justin Peel ‘ s Antwort, aber nur numpy verwenden und die Vorteile des integrierten in Vektorisierung unter:

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

Beachten Sie, dass, wie in den Kommentaren vorgeschlagen, ich habe änderte die Definition der Elevationswinkel aus Ihrer ursprünglichen Funktion. Auf meinem Rechner mit pts = np.random.rand(3000000, 3) testen, ging die Zeit von 76 Sekunden bis 3,3 Sekunden. Ich habe nicht Cython so dass ich nicht in der Lage war, das Timing, mit dieser Lösung zu vergleichen.

Answer 2

Hier ist eine schnelle Cython Code, dass 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 sicher, es gibt einige andere Verbesserungen, die gemacht werden können.

Answer 3

! Es ist ein Fehler nach wie vor in dem gesamten Code oben .. und das ist ein Top-Google-Ergebnis .. TLDR: Ich habe dies mit VPython getestet, unter Verwendung von atan2 für Theta (elev) ist falsch, Verwendung acos! Es ist richtig, für phi (azim). Ich empfehle die sympy1.0 acos Funktion (es nicht einmal über acos beschwert (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 optimierte, aber richtig Code für Rechtshänder Physik-System:

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]

können 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]]

Ich habe VPython zusätzlich zu leicht Vektoren zu visualisieren:

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

Answer 4

, um 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 bis zu reinen einer Numpy Implementierung und würde zu C oder Cython Geschwindigkeiten vergleichbar sein. Die vorliegende numpy Lösung (wenn sie mit dem ceval=eval Argument verwendet) ist ebenfalls um 25% schneller als die appendSpherical_np Funktion im @mtrw Antwort für große Feldgröß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