Более быстрое переоборудование Numpy Cartesian в сферическое преобразование координат?
-
29-09-2019 - |
Вопрос
У меня есть массив из 3 миллионов точек данных из 3-осового аккаурка (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)))
Решение
Это похоже на Джастин ПилОтвет, но используя просто 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, поэтому я не смог сравнить время с этим решением.
Другие советы
Вот быстрый код цинтона, который я написал для этого:
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 секунды, используя для меня 3 000 000 баллов. Это не слишком потертый. Я уверен, что есть некоторые другие улучшения, которые могут быть сделаны.
! Есть ошибка еще во всем коде выше .. И это главный результат Google. TLDR: Я проверил это с vpython, используя ATAN2 для TETA (EVE) неверно, используйте ACOS! Это правильно для Phi (Azim). Я рекомендую функцию ACOS Sympy1.0 (даже не жалуется на ACOS (Z / R) с R = 0).
http://mathworld.wolfram.com/sphericalcoordinates.html
Если мы преобразуем это в физическую систему (r, theta, phi) = (r, lete, azimuth), у нас есть:
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)
Чтобы завершить предыдущие ответы, вот 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 или цинтона. Нынешнее решение 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