Chuyển đổi nhanh hơn để chuyển đổi tọa độ hình cầu?

Question

Tôi có một mảng gồm 3 triệu điểm dữ liệu từ bộ đếm công suất 3 trục (XYZ) và tôi muốn thêm 3 cột vào mảng chứa tọa độ hình cầu tương đương (r, theta, phi). Mã sau hoạt động, nhưng dường như quá chậm. Làm thế nào tôi có thể làm tốt hơn?

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

này cũng tương tự như câu trả lời Justin Peel 's, nhưng chỉ sử dụng numpy và tận dụng lợi thế của nó được xây dựng trong vector hóa:

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 

Lưu ý rằng, theo đề nghị trong các ý kiến, tôi đã thay đổi định nghĩa góc độ cao từ chức năng ban đầu của bạn. Trên máy của tôi, thử nghiệm với pts = np.random.rand(3000000, 3), thời gian đã trôi qua từ 76 giây đến 3,3 giây. Tôi không có Cython nên tôi không thể so sánh thời gian với giải pháp đó.

Answer 2

Dưới đây là một mã Cython nhanh chóng mà tôi đã viết lên cho việc này:

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 

Phải mất thời gian giảm từ 62,4 giây để 1.22 giây sử dụng 3.000.000 điểm đối với tôi. Đó không phải là quá tồi tàn. Tôi chắc rằng có một số cải tiến khác có thể được thực hiện.

Answer 3

Để hoàn thành các câu trả lời trước, đây là một thực hiện Numexpr (với một fallback thể 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 

Đối với kích thước mảng lớn, điều này cho phép hệ số 2 tốc độ lên so với thực hiện NumPy tinh khiết và có thể so sánh với tốc độ C hoặc Cython. Các giải pháp NumPy hiện nay (khi được sử dụng với đối số ceval=eval) cũng là nhanh hơn so với appendSpherical_np chức năng trong câu trả lời @mtrw Đối với kích thước mảng lớn 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 

mặc dù đối với kích thước nhỏ hơn, appendSpherical_np là thực sự nhanh hơn,

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 

Answer 4

! Đã xảy ra lỗi trong tất cả các mã ở trên .. và đây là kết quả Google hàng đầu .. TLDR: Tôi đã thử nghiệm điều này với VPython, sử dụng atan2 cho theta (độ cao) là sai, sử dụng acos! Nó là đúng cho phi (azim). Tôi đề nghị chức năng acos sympy1.0 (nó thậm chí không phàn nàn về acos (z/r) với r = 0).

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

Nếu chúng ta chuyển đổi đó để hệ thống vật lý (r, theta, phi) = (r, elev, góc phương vị) ta có:

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

không tối ưu nhưng đúng mã cho hệ thống vật lý thuận tay phải:

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] 

bạn có thể tự mình thử nghiệm với chức năng như:

test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319])) 

một số dữ liệu thử nghiệm khác đối với một số trục tọa:

[[ 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]] 

tôi đã sử dụng VPython thêm để dễ dàng hình dung vectơ:

test = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)