37

3 軸加速度計 (XYZ) からの 300 万のデータ ポイントの配列があり、同等の球面座標 (r、シータ、ファイ) を含む配列に 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)))
4

5 に答える 5

40

これは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 を持っていないので、タイミングをそのソリューションと比較することができませんでした。

于 2010-11-07T07:54:27.480 に答える
15

このために私が書いた簡単な 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 秒に短縮されました。それはあまりにもみすぼらしいことではありません。他にも改善できる点があると確信しています。

于 2010-11-07T07:10:32.167 に答える
9

前の回答を完了するために、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 ソリューション (引数を使用した場合) は、大きな配列サイズの@mtrw回答の関数ceval=evalよりも 25% 高速です。appendSpherical_np

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
于 2015-05-12T08:37:42.317 に答える