Python 更快的numpy笛卡尔到球坐标转换？
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/4116658/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Faster numpy cartesian to spherical coordinate conversion?
提问by BobC
I have an array of 3 million data points from a 3-axiz accellerometer (XYZ), and I want to add 3 columns to the array containing the equivalent spherical coordinates (r, theta, phi). The following code works, but seems way too slow. How can I do better?
我有一个来自 3 轴加速度计 (XYZ) 的 300 万个数据点的数组，我想向包含等效球坐标 (r, theta, phi) 的数组添加 3 列。以下代码有效，但似乎太慢了。我怎样才能做得更好？
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)))
采纳答案by mtrw
This is similar to Justin Peel's answer, but using just numpyand taking advantage of its built-in vectorization:
这类似于Justin Peel的答案，但仅使用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
Note that, as suggested in the comments, I've changed the definition of elevation anglefrom your original function. On my machine, testing with pts = np.random.rand(3000000, 3), the time went from 76 seconds to 3.3 seconds. I don't have Cython so I wasn't able to compare the timing with that solution.
请注意，正如评论中所建议的，我已经从原始函数中更改了仰角的定义。在我的机器上，使用 进行测试pts = np.random.rand(3000000, 3)，时间从 76 秒变为 3.3 秒。我没有 Cython，所以我无法将时间与该解决方案进行比较。
回答by Justin Peel
Here's a quick Cython code that I wrote up for this:
这是我为此编写的快速 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
It took the time down from 62.4 seconds to 1.22 seconds using 3,000,000 points for me. That's not too shabby. I'm sure there are some other improvements that can be made.
对我来说，使用 3,000,000 点将时间从 62.4 秒缩短到 1.22 秒。那还不算太简陋。我相信还有一些其他的改进可以做。
回答by rth
To complete the previous answers, here is a Numexprimplementation (with a possible fallback to Numpy),
为了完成前面的答案，这里是一个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
For large array sizes, this allows a factor of 2 speed up compared to pure a Numpy implementation, and would be comparable to C or Cython speeds. The present numpy solution (when used with the ceval=evalargument) is also 25% faster than the appendSpherical_npfunction in the @mtrwanswer for large array sizes,
对于大数组大小，与纯 Numpy 实现相比，这允许提高 2 倍的速度，并且与 C 或 Cython 的速度相当。对于大数组大小，当前的 numpy 解决方案（与ceval=eval参数一起使用时）也比@mtrw答案中的appendSpherical_np函数快 25% ，
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
although for smaller sizes, appendSpherical_npis actually faster,
虽然对于较小的尺寸，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
回答by Vincent
! There is an error still in all the code above.. and this is a top Google result.. TLDR: I have tested this with VPython, using atan2 for theta (elev) is wrong, use acos! It is correct for phi (azim). I recommend the sympy1.0 acos function (it does not even complain about acos(z/r) with r = 0 ) .
！上面所有代码中仍然存在错误..这是谷歌的顶级结果.. TLDR：我已经用 VPython 测试过这个，使用 atan2 作为 theta (elev) 是错误的，使用 acos！phi (azim) 是正确的。我推荐使用 sympy1.0 acos 函数（它甚至不会抱怨 acos(z/r) 与 r = 0 ）。
http://mathworld.wolfram.com/SphericalCoordinates.html
http://mathworld.wolfram.com/SphericalCoordinates.html
If we convert that to the physics system (r, theta, phi) = (r, elev, azimuth) we have:
如果我们将其转换为物理系统 (r, theta, phi) = (r, elev, azimuth) 我们有：
r = sqrt(x*x + y*y + z*z)
phi = atan2(y,x)
theta = acos(z,r)
Non optimized but correctcode for right-handed physics system:
右手物理系统未优化但正确的代码：
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]
you can test it yourself with a function like:
您可以使用如下函数自行测试：
test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319]))
some other test data for some quadrants:
某些象限的一些其他测试数据：
[[ 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]]
I used VPython additionally to easily visualize vectors:
我另外使用 VPython 来轻松可视化向量：
test = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)
