Ускоренное преобразование декартовой в сферическую координату?

Question

У меня есть массив из 3 миллионов точек данных от 3-axiz accellerometer (XYZ), и я хочу добавить 3 столбца в массив, содержащий эквивалентные сферические координаты (r, theta, phi). Следующий код работает, но кажется слишком медленным. Как я могу сделать лучше?

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

Это похоже на ответ Justin Peel «s, но с использованием только numpy и воспользовавшись встроенной в векторизации:

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 

Обратите внимание, что, как это было предложено в комментариях, я изменил определение угла места от вашего оригинального функция. На моей машине, тестирование с pts = np.random.rand(3000000, 3), время шло от 76 секунд до 3,3 секунды. У меня нет Cython, поэтому я не смог сравнить время с этим решением.

Answer 2

Вот быстрый код Cython, что я написал для этого:

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 

Потребовалось время, по сравнению с 62,4 секунды до 1,22 секунды, используя 3000000 очки для меня. Это не слишком потрепано. Я уверен, что есть некоторые другие улучшения, которые могут быть сделаны.

Answer 3

Чтобы завершить предыдущие ответы, вот Numexpr реализация (с возможностью отката к 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 

Для больших размеров массивов, это позволяет в 2 раза скорости до по сравнению с чистой реализации Numpy , и будет сопоставим с скоростями C или Cython. Присутствует NumPy раствор (при использовании с ceval=eval аргумента) также на 25% быстрее, чем appendSpherical_np функции в @mtrw ответ для больших размеров массива,

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 

хотя для меньших размеров, appendSpherical_np на самом деле быстрее,

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

! Ошибка во всем вышеприведенном коде .. и это лучший результат Google. TLDR: Я тестировал это с помощью VPython, используя atan2 для theta (elev) неправильно, используйте acos! Это правильно для phi (азим). Я рекомендую функцию sympy1.0 acos (он даже не жалуется на acos (z/r) с r = 0).

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

Если мы преобразуем, что в систему физики (г, тета, фи) = (г, отм, азимут) мы имеем:

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

неоптимизированного но правильных код правша система физики:

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] 

вы можете проверить это самостоятельно с помощью функции, как:

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

некоторые другие данные испытаний для некоторых квадрантов:

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

Я VPython дополнительно легко визуализировать векторы:

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