当前位置： 文江博客 话题详情

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

发布于 2024-10-01 16:44:29 字数 540 浏览 1 评论 0 原文

我有来自 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)))
原文

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?

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)))

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（7）

话少情深 2024-10-08 16:44:29

这类似于 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，所以我无法将时间与该解决方案进行比较。

This is similar to Justin Peel's answer, but using just numpy and taking advantage of its built-in vectorization:

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 angle from 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.

回复 原文
马蹄踏│碎落叶 2024-10-08 16:44:29

以下是我为此编写的快速 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

我使用 3,000,000 个点将时间从 62.4 秒缩短到 1.22 秒。那还不算太寒酸。我确信还可以进行一些其他改进。

Here's a quick Cython code that I wrote up for this:

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.

回复 原文
少年亿悲伤 2024-10-08 16:44:29

！上面的所有代码中仍然存在错误..这是 Google 的顶级结果..
总而言之：
我已经用 VPython 对此进行了测试，使用 atan2 表示 theta (elev) 是错误的，请使用
阿科斯！对于 phi (azim) 来说这是正确的。
我推荐 sympy1.0 acos 函数（它甚至不会抱怨 acos(z/r) 和 r = 0 ）。

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

如果我们将其转换为物理系统(r, theta, phi) = (r, elev, 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)

! 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 ) .

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

If we convert that to the physics system (r, theta, phi) = (r, elev, azimuth) we have:

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

Non optimized but correct code 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:

test   = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)
回复 原文
森罗 2024-10-08 16:44:29

为了完成前面的答案，这里有一个 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

对于大型数组大小，与纯 Numpy 实现相比，这允许速度提高 2 倍，并且与 C 或 Cython 速度相当。当前的 numpy 解决方案（与 ceval=eval 参数一起使用时）也比 appendSpherical_np 函数快 25%。 com/a/4116899/1791279">@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

To complete the previous answers, here is a Numexpr implementation (with a possible fallback to 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=eval argument) is also 25% faster than the appendSpherical_np function in the @mtrw answer for large array sizes,

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_np is actually faster,

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
回复 原文
暖风昔人 2024-10-08 16:44:29

Octave 有一些内置的坐标转换功能，可以使用 oct2py 包访问这些功能，以将笛卡尔坐标中的 numpy 数组转换为球面或极坐标（以及反之）：

from oct2py import octave
xyz = np.random.rand(3000000,3)
%timeit thetaphir = octave.cart2sph(xyz)

724 ms ± 206 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Octave has some built-in functionality for coordinate transformations that can be accessed with the package oct2py to convert numpy arrays in Cartesian coordinates to spherical or polar coordinates (and back):

from oct2py import octave
xyz = np.random.rand(3000000,3)
%timeit thetaphir = octave.cart2sph(xyz)

724 ms ± 206 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
回复 原文
一向肩并 2024-10-08 16:44:29

我的代码基于 mrtw 的答案 - 我添加了轴支持、反向函数并将其基于三元组形状：

def from_xyz(xyz, axis=-1):
    x, y, z = np.moveaxis(xyz, axis, 0)

    lea = np.empty_like(xyz)

    pre_selector = ((slice(None),) * lea.ndim)[:axis]

    xy_sq = x ** 2 + y ** 2
    lea[(*pre_selector, 0)] = np.sqrt(xy_sq + z ** 2)
    lea[(*pre_selector, 1)] = np.arctan2(np.sqrt(xy_sq), z)
    lea[(*pre_selector, 2)] = np.arctan2(y, x)

    return lea


def to_xyz(lea, axis=-1):
    l, e, a = np.moveaxis(lea, axis, 0)

    xyz = np.empty_like(lea)

    pre_selector = ((slice(None),) * xyz.ndim)[:axis]

    xyz[(*pre_selector, 0)] = l * np.sin(e) * np.cos(a)
    xyz[(*pre_selector, 1)] = l * np.sin(e) * np.sin(a)
    xyz[(*pre_selector, 2)] = l * np.cos(e)

    return xyz

I base my code on mrtw's answer - I added axis support, the reverse function and base it off a 3-tuple shape:

def from_xyz(xyz, axis=-1):
    x, y, z = np.moveaxis(xyz, axis, 0)

    lea = np.empty_like(xyz)

    pre_selector = ((slice(None),) * lea.ndim)[:axis]

    xy_sq = x ** 2 + y ** 2
    lea[(*pre_selector, 0)] = np.sqrt(xy_sq + z ** 2)
    lea[(*pre_selector, 1)] = np.arctan2(np.sqrt(xy_sq), z)
    lea[(*pre_selector, 2)] = np.arctan2(y, x)

    return lea


def to_xyz(lea, axis=-1):
    l, e, a = np.moveaxis(lea, axis, 0)

    xyz = np.empty_like(lea)

    pre_selector = ((slice(None),) * xyz.ndim)[:axis]

    xyz[(*pre_selector, 0)] = l * np.sin(e) * np.cos(a)
    xyz[(*pre_selector, 1)] = l * np.sin(e) * np.sin(a)
    xyz[(*pre_selector, 2)] = l * np.cos(e)

    return xyz
回复 原文
尾戒 2024-10-08 16:44:29

您可以使用hyperspherical包。它适用于任何维度。

import numpy as np
from hyperspherical import cartesian2spherical, spherical2cartesian

xyz = np.random.rand(3000000,3)
%timeit r_theta_phi = cartesian2spherical(xyz)

时间：

257 ms ± 7.78 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

You can use hyperspherical package. It works for any dimension.

import numpy as np
from hyperspherical import cartesian2spherical, spherical2cartesian

xyz = np.random.rand(3000000,3)
%timeit r_theta_phi = cartesian2spherical(xyz)

The time:

257 ms ± 7.78 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
回复 原文
~没有更多了~

关于作者

旧情别恋

暂无简介

0 文章
0 评论
21 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

lioqio

文章 0 评论 0

Single

文章 0 评论 0

禾厶谷欠

文章 0 评论 0

alipaysp_2zg8elfGgC

文章 0 评论 0

qq_N6d4X7

文章 0 评论 0

放低过去

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文