更快的数值笛卡尔坐标到球坐标的转换？

Question

我有一个来自3轴加速度计(XYZ)的300万个数据点的数组，我想向包含等效球坐标(r，θ，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)))

Answer 1

这与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

请注意，正如注释中所建议的，我已经从原始函数更改了仰角的定义。在我的机器上，使用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

对于我来说，使用300万分，它将时间从62.4秒缩短到1.22秒。这还不算太寒酸。我相信还可以做一些其他的改进。

Answer 3

好了！上述所有代码中仍有一个错误。这是谷歌的顶级搜索结果..我已经用VPython测试过了，使用atan2的theta (elev)是错误的，使用acos！这对于phi (azim)是正确的。我推荐symmy1.0 acos函数(它甚至不会抱怨r=0的acos(z/r) )。

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

如果我们将其转换为物理系统(r，θ，φ)= (r，elev，方位角)，我们得到：

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)