выполнение расчетов и сравнений с использованием индексов в массивах

Question

У меня есть трехмерный массив вокселей, т. Е. Где индексы каждой точки соответствуют ее положению в трехмерном пространстве.

Есть ряд вычислений, которые я хотел бы выполнить, зная, удовлетворяют ли точки удовлетворениюразличные геометрические условия.Это означает преобразование индексов в векторы путем умножения на базис, а затем вычисления точечных и перекрестных произведений, норм и т. Д. Я ищу достаточно быстрый способ сделать это, поскольку мои усилия пока кажутся очень медленными.

Если у меня есть произвольная основа a , b , c :

basis = np.array([[a1, a2, a3],
                  [b1, b2, b3],
                  [c1, c2, c3]])

, где a1, a2, a3 - x, компоненты y и z a , а также для b и c .Я могу вычислить декартову координату p = (x, y, z) каждого вокселя следующим образом:

for i in range(vox.shape[0]):
    for j in range(vox.shape[1]):
        for k in range(vox.shape[2]):
            p = np.dot(basis, np.array([i, j, k]))

, где «vox» - трехмерный массив вокселей.Если, например, я хочу вычислить скалярное произведение каждого из этих векторов с одним другим (декартовым) вектором (например, q = np.array([qx, qy, qz])) и сохранить этот индекс, если результат удовлетворяет заданному условию (здесь больше 0,0), яможет сделать что-то вроде этого (в том же цикле, как указано выше):

            if np.dot(p, q) > 0.0:
                desired_vox_indices.append([i, j, k])

Проблема в том, что это очень медленно.Могу ли я сделать это более питоническим способом или с помощью более обалденных инструментов?Я понимаю, что на этом этапе я даже не обращаюсь к значениям массива vox.

РЕДАКТИРОВАТЬ: попытка перекрестного произведения на основе ответа Дивакара

# Get q_x (matrix multiplication version of cross product)
q_x = np.array([[0.0, -q[2], q[1]],
                [q[2], 0.0, -q[0]],
                [-q[1], q[0], 0.0]])

# transpose (as per cross product definition) and matrix multiply with basis
u = np.matmul(q_x.T, basis)

# Get open range arrays
m,n,r = vox.shape
I,J,K = np.ogrid[:m,:n,:r]

# writing everything explicitly, since I am unsure how ogrid objects behave
pxq = np.empty(3)
pxq[0] = u[0,0]*I + u[0,1]*J + u[0,2]*K
pxq[1] = u[1,0]*I + u[1,1]*J + u[1,2]*K
pxq[2] = u[2,0]*I + u[2,1]*J + u[2,2]*K

, который может совпадать спросто:

pxq = np.dot(u, np.array([I, J, K]))

но я не уверен ...

Answer 1

Мы будем строить суммы, используя массивы диапазонов для масштабирования, не генерируя сразу все индексы.Мы сделаем это в три этапа, соответствующих трем вложенным циклам.Идея очень похожа на рассмотренную в this answer to - Python vectorizing nested for loops.Это будет эффективным для памяти и, следовательно, хорошим результатом.Затем мы сравним суммы с порогом 0.0 и используем np.argwhere для получения соответствующих индексов.Таким образом, у нас было бы решение, например, так:

# Get q scaled version
s = q.dot(basis)

# Get open range arrays and scale and sum-reduce s array
m,n,r = vox.shape
I,J,K = np.ogrid[:m,:n,:r]
sums = s[0]*I + s[1]*J + s[2]*K

# Finally compare the sums against threshold amd get corresponding indices
out = np.argwhere(sums > 0.0)

Времена для большого массива vox -

# Setup
In [371]: np.random.seed(0)
     ...: basis = np.random.randn(3,3)
     ...: vox = np.random.randn(100,100,100)
     ...: q = np.random.randn(3)

# Original soln
In [372]: %%timeit
     ...: desired_vox_indices = []
     ...: for i in range(vox.shape[0]):
     ...:     for j in range(vox.shape[1]):
     ...:         for k in range(vox.shape[2]):
     ...:             p = np.dot(basis, np.array([i, j, k]))            
     ...:             if np.dot(p, q) > 0.0:
     ...:                 desired_vox_indices.append([i, j, k])
1 loop, best of 3: 2.13 s per loop

# @jdehesa's soln
In [373]: %%timeit
     ...: ind = np.indices(vox.shape).reshape(3, -1)
     ...: p = basis.dot(ind)
     ...: d = q.dot(p)
     ...: desired_vox_indices = ind[:, d > 0.0].T
10 loops, best of 3: 35.9 ms per loop

# From this post
In [374]: %%timeit
     ...: s = q.dot(basis)
     ...: m,n,r = vox.shape
     ...: I,J,K = np.ogrid[:m,:n,:r]
     ...: sums = s[0]*I + s[1]*J + s[2]*K
     ...: out = np.argwhere(sums > 0.0)
100 loops, best of 3: 7.56 ms per loop

Answer 2

Вы можете сделать это в векторизованном виде, как это:

# Array of indices for your voxel data
ind = np.indices(vox.shape).reshape(3, -1)
# Multiply the basis times each coordinate
p = basis @ ind
# Dot product of each result with vector
d = q @ p
# Select coordinates where criteria is met
desired_vox_indices = ind[:, d > 0.0].T