Faster numpy cartesiano de coordenadas esféricas conversión?

Question

I tienen una serie de 3 millones de puntos de datos de un accellerometer 3-axiz (XYZ), y yo quiero añadir 3 columnas a la matriz que contiene las coordenadas esféricas equivalentes (r, theta, PHI). El siguiente código funciona, pero parece demasiado lento. ¿Cómo puedo hacer mejor?

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

Esto es similar a Justin Peel ' s respuesta, pero utilizando sólo numpy y aprovechando su base de vectorización:

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

Tenga en cuenta que, como se sugiere en los comentarios, he cambiado la definición del ángulo de elevación de su función original. En mi máquina, probando con pts = np.random.rand(3000000, 3), el tiempo pasó de 76 segundos a 3,3 segundos. No tengo Cython así que no era capaz de comparar el tiempo con esa solución.

Answer 2

Aquí hay un código Cython rápida que escribí para esto:

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

Se tomó el tiempo de inactividad de 62,4 segundos a 1,22 segundos usando 3.000.000 puntos para mí. Eso no está nada mal. Estoy seguro de que hay algunas otras mejoras que se pueden hacer.

Answer 3

! Hay un error todavía en todo el código anterior .. y este es un buen resultado Google .. TLDR: He probado esto con VPython, utilizando atan2 para theta (elev) es incorrecto, el uso acos! Es correcto para phi (azim). Recomiendo el sympy1.0 acos función (incluso no se quejan de acos (z / r) con r = 0).

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

Si convertimos que al sistema de la física (r, theta, phi) = (r, elev, acimut) tenemos:

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

No optimizado pero correcta código para el sistema de la física diestro:

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]

Puede probar usted mismo con una función como:

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

algunos otros datos de prueba para algunos cuadrantes:

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

utilicé VPython además de visualizar fácilmente vectores:

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

Answer 4

Para completar las respuestas anteriores, aquí hay una Numexpr aplicación (con un posible retorno a 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

Para grandes tamaños de matriz, esto permite un factor de 2 velocidad en comparación con una implementación pura Numpy, y sería comparable a velocidades de C o Cython. La solución numpy presente (cuando se utiliza con el argumento ceval=eval) también es 25% más rápido que la función appendSpherical_np en el @mtrw respuesta para tamaños grandes de la matriz,

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

A pesar de que para los tamaños más pequeños, appendSpherical_np es realmente más rápido,

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