Conversión cartesiana numerada más rápida a esférica?

Question

Tengo una matriz de 3 millones de puntos de datos de un acelerómetro de 3 ejes (XYZ) y deseo agregar 3 columnas al conjunto que contiene las coordenadas esféricas equivalentes (r, theta, phi). El siguiente código funciona, pero parece demasiado lento. ¿Cómo puedo hacerlo 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 la respuesta Justin Peel 's, 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 desde 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 pude comparar el tiempo con esa solución.

Answer 2

Aquí hay un código rápido Cython 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 con 3.000.000 puntos para mí. Eso no es muy viejo. Estoy seguro de que hay otras mejoras que se pueden hacer.

Answer 3

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 pura una implementación Numpy , y sería comparable a las velocidades C o Cython. La solución numpy presente (cuando se utiliza con el argumento ceval=eval) es también 25% más rápido que la función appendSpherical_np en la respuesta @mtrw para grandes tamaños de 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 

aunque para tamaños más pequeños, appendSpherical_np es en realidad 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 

Answer 4

! Todavía hay un error en todo el código anterior ... y este es un resultado superior de Google ... TLDR: He probado esto con VPython, usando atan2 para theta (elev) es incorrecto, use acos! Es correcto para phi (azim). Recomiendo la función sympy1.0 acos (ni siquiera se queja 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 correcto código para sistema de mano derecha física:

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

Solía ​​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)