¿Número más rápido de conversión de coordenadas cartesianas a esféricas?

Tengo una matriz de 3 millones de puntos de datos de un acelerómetro de 3 ejes (XYZ), y quiero agregar 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 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))) 

Esto es similar a la respuesta de Justin Peel , pero usando solo numpy y aprovechando su vectorización incorporada:

 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 de á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.

Aquí hay un rápido código de 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 

Tomó el tiempo de 62.4 segundos a 1.22 segundos usando 3,000,000 puntos para mí. Eso no es tan lamentable. Estoy seguro de que hay algunas otras mejoras que se pueden hacer.

! 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, usar atan2 para theta (once) está mal, ¡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 eso al sistema de física (r, theta, phi) = (r, elev, azimut) tenemos:

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

Código no optimizado pero correcto para el sistema de 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] 

Puedes probarlo tú 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 para visualizar fácilmente vectores:

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

Para completar las respuestas anteriores, aquí hay una implementación Numexpr (con una posible alternativa 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 tamaños de matriz grandes, esto permite un factor de aceleración de 2 en comparación con una implementación Numpy pura, y sería comparable a las velocidades C o Cython. La solución numpy actual (cuando se usa con el argumento ceval=eval ) también es un 25% más rápida que la función appendSpherical_np en la respuesta @mtrw para tamaños de matriz grandes,

 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