将三参数函数应用到3D numpy阵列

Question

我有一个形状（2133,3,3）的3D numpy阵列A.基本上这是包含三个3D点的2133个列表的列表。此外，我有一个函数需要三个三维点，并返回一个三维点，x = f(a, b, c)，长度为3的a，b，c，x numpy数组。现在我想将f应用于A，以便输出是一个形状数组（2133，3）。所以像numpy.array([f(*A[0]),...,f(*A[2132]))。

我试过numpy.apply_along_axis和numpy.vectorize没有成功。

更精确的函数f我考虑的是由下式给出：

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）：

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

生产：

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

这看起来像一行产生了（3）数组，另一行出现了某种问题并产生了None。我不知道None是否是由于删除了这种排序。但无论如何，将数组和None的混合重新转换为数组将会成为问题。如果res的所有项目都是匹配的数组，我们可以将其stack返回到2d数组中。

有些方法可以获得适度的速度提升（与此列表理解相比）。但是，像这样一个复杂的函数，函数花费的时间（被称为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。