3D numpy配列に3引数関数を適用する

Question

3D numpy配列A（2133,3,3）があります。基本的にこれは3つの3Dポイントを持つ2133リストのリストです。さらに私は、3つの3Dポイントを取り、長さ3のa、b、c、x numpyの配列を持つx = f(a, b, c)という1つの3Dポイントを返す関数を持っています。今度は、出力を形状の配列にするためにfにAを適用したい（2133,3）。だから、numpy.array([f(*A[0]),...,f(*A[2132]))のようなものです。

numpy.apply_along_axisとnumpy.vectorizeを試しましたが、成功しませんでした。

def f(a, b, c, r1, r2=None, r3=None): 
    a = np.asarray(a) 
    b = np.asarray(b) 
    c = np.asarray(c) 

    if np.linalg.matrix_rank(np.matrix([a, b, c])) != 3: 
     # raise ValueError('The points are not collinear.') 
     return None 

    a, b, c, = sort_triple(a, b, c) 

    if any(r is None for r in (r2, r3)): 
     r2, r3 = (r1, r1) 

    ex = (b - a)/(np.linalg.norm(b - a)) 
    i = np.dot(ex, c - a) 
    ey = (c - a - i*ex)/(np.linalg.norm(c - a - i*ex)) 
    ez = np.cross(ex, ey) 
    d = np.linalg.norm(b - a) 
    j = np.dot(ey, c - a) 

    x = (pow(r1, 2) - pow(r2, 2) + pow(d, 2))/(2 * d) 
    y = ((pow(r1, 2) - pow(r3, 2) + pow(i, 2) + pow(j, 2))/(2*j)) - ((i/j)*x) 
    z_square = pow(r1, 2) - pow(x, 2) - pow(y, 2) 
    if z_square >= 0: 
     z = np.sqrt(z_square) 
     intersection = a + x * ex + y*ey + z*ez 
     return intersection 

A = np.array([[[131.83, 25.2, 0.52], [131.51, 22.54, 0.52],[133.65, 23.65, 0.52]], [[13.02, 86.98, 0.52], [61.02, 87.12, 0.52],[129.05, 87.32, 0.52]]]) 

r1 = 1.7115 

Answer 1

これは正常に動作します（sort_tripleをコメントアウトした後）：

は私が考える関数fによって与えられ、より正確には

res = [f(*row,r1) for row in A] 
print(res) 

が生産：

[array([ 132.21182324, 23.80481826, 1.43482849]), None] 

1列に（3）配列が生成され、もう1行に何らかの問題があり、None。そのNoneがその並べ替えを削除したかどうかはわかりません。しかし、どのような場合でも、配列の組み合わせとNoneを配列に戻すことは問題になります。 resのすべての項目が一致する配列であれば、stackを2次元配列に戻すことができます。

（このリストの解説と比較して）控えめな速度向上を得る方法があります。しかし、このような複雑な関数では、関数に費やされた時間（2000回と呼ばれます）は、反復メカニズムが費やした時間を支配します。

そして、あなたが第一次元で反復し、他の2（3として配列）を通過しているので、この明示的なループはvectorize、frompyfuncまたはapply_along/over..より使いずっと簡単です。

時間を大幅に節約するには、f()と書いて3dアレイを直接操作する必要があります。

Answer 2

@jdehesaのおかげで、私は@hpauljによって与えられたものに代わる解決策を作り出すことができました。このソリューションが最もエレガントなソリューションであるかどうかはわかりませんが、それはこれまでのところうまくいきました。コメントは高く評価されます。

def sort_triple(a, b, c): 
    pts = np.stack((a, b, c), axis=1) 
    xSorted = pts[np.arange(pts.shape[0])[:, None], np.argsort(pts[:, :, 0])] 
    orientation = np.cross(xSorted[:, 1] - xSorted[:, 0], xSorted[:, 2] - 
          xSorted[:, 0])[:, 2] >= 0 
    xSorted_flipped = np.stack((xSorted[:, 0], xSorted[:, 2], xSorted[:, 1]), 
           axis=1) 
    xSorted = np.where(orientation[:, np.newaxis, np.newaxis], xSorted, 
         xSorted_flipped) 
    return map(np.squeeze, np.split(xSorted, 3, axis=1)) 


def f(A, r1, r2=None, r3=None): 

    a, b, c = map(np.squeeze, np.split(A, 3, axis=1)) 
    a, b, c = sort_triple(a, b, c) 

    if any(r is None for r in (r2, r3)): 
     r2, r3 = (r1, r1) 

    ex = (b - a)/(np.linalg.norm(b - a, axis=1))[:, np.newaxis] 
    i = inner1d(ex, (c - a)) 
    ey = ((c - a - i[:, np.newaxis]*ex)/
      (np.linalg.norm(c - a - i[:, np.newaxis]*ex, axis=1))[:, np.newaxis]) 
    ez = np.cross(ex, ey) 
    d = np.linalg.norm(b - a, axis=1) 
    j = inner1d(ey, c - a) 

    x = (np.square(r1) - np.square(r2) + np.square(d))/(2 * d) 
    y = ((np.square(r1) - np.square(r3) + np.square(i) + np.square(j))/(2*j) - 
     i/j*x) 
    z_square = np.square(r1) - np.square(x) - np.square(y) 
    mask = z_square < 0 
    z_square[mask] *= 0 
    z = np.sqrt(z_square) 
    z[mask] = np.nan 
    intersection = (a + x[:, np.newaxis] * ex + y[:, np.newaxis] * ey + 
        z[:, np.newaxis] * ez) 
    return intersection 

おそらく、各機能のmapの部分がうまく実行される可能性があります。多分np.newaxisの過剰使用も。