Schnelleres Umschalten von kartesisch auf sphärische Koordinaten?

Ich habe ein array von 3 Millionen Datenpunkte aus einem 3-axiz accellerometer (XYZ), und ich möchte hinzufügen 3 Spalten array mit den äquivalenten sphärischen Koordinaten (r, theta, phi). Der folgende code funktioniert, scheint aber 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)))

InformationsquelleAutor der Frage BobC | 2010-11-07

29

Dies ist ähnlich zu Justin Peel's Antwort, aber mit nur numpy und die Vorteile der integrierten-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
```
Beachten Sie, dass, wie vorgeschlagen, in die Kommentare, ich habe verändert die definition der Neigungswinkel von Ihrer ursprünglichen Funktion. Auf meinem Rechner testen mit pts = np.random.rand(3000000, 3)die Zeit ging von 76 Sekunden auf 3,3 Sekunden. Ich habe nicht Cython so war ich nicht in der Lage zu vergleichen, die timing mit dieser Lösung.

InformationsquelleAutor der Antwort mtrw

Hier ist eine schnelle Cython-code, den ich schrieb:

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

Er nahm die Zeit von 62,4 Sekunden 1,22 Sekunden mit 3,000,000 Punkte für mich. Das ist nicht zu schäbig. Ich bin sicher, es gibt einige andere Verbesserungen, die gemacht werden können.

InformationsquelleAutor der Antwort Justin Peel

Die vorherigen Antworten, hier ist ein Numexpr Umsetzung (mit einer möglichen fallback 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 array-Größen, dies ermöglicht einen Faktor 2 beschleunigen, im Vergleich zur reinen ein Numpy-Implementierung, und wäre vergleichbar mit C oder Cython-Geschwindigkeiten. Die vorliegende numpy-Lösung (bei Verwendung des ceval=eval - argument) ist auch 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 die kleineren Größen, appendSpherical_np ist tatsächlich 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

InformationsquelleAutor der Antwort rth

! Es ist ein Fehler noch in den code oben.. und das ist ein top-Google-Ergebnis..
TLDR:
Ich habe getestet, mit VPython, mit atan2 für theta (elev) falsch ist, verwenden Sie
acos! Es ist korrekt für phi (azim).
Ich empfehle die sympy1.0 acos-Funktion (es muss nicht einmal beschweren sich über acos(z/r) mit r = 0 ) .

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

Wenn wir zu konvertieren, dass auf die Physik-system (r, theta, phi) = (r, Höhe, Azimut) wir haben:

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

Nicht optimiert, sondern richtige 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 testen mit einer Funktion wie:

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

einige andere test-Daten 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 verwendet VPython zusätzlich zu einfaches visualisieren von Vektoren:

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

InformationsquelleAutor der Antwort Vincent

Schreibe einen Kommentar

Du musst angemeldet sein, um einen Kommentar abzugeben.