NumPy 直方图索引

Question

假设我有一个3D直方图或简单的形状为（X，Y，Z）的3D numpy数组

import numpy as np
array = np.random.random((100,100,100))

使用numpy或scipy最好的方法是获取满足球形条件的数组值的索引。

(index_x**2 + index_y**2 + index_z**2) <= radius**2

显然，在后一种情况下，数组的中心是 (0, 0, 0)。通常情况下，条件将会是

((index_x-center_x)**2 + (index_y-center_y)**2 +(index_z-center_z)**2) <= radius**2

问题很容易用Python循环解决，但我需要进行优化。

非常感谢您的帮助。

- Cobry

Answer 1

你可以先使用 ogrid()高效地获取索引，然后再使用 nonzero() 获取符合条件的索引。

可以通过 nonzero()来获取索引：

indexes = numpy.transpose((x**2+y**2+z**2 <= radius**2).nonzero())  # transpose() might be unnecessary: it depends on your needs

索引数组可以通过高效地使用ogrid()函数获得：

x, y, z = numpy.ogrid[:100, :100, :100]

或者，针对您的输入data数组的任意形状：

x, y, z = ogrid[tuple(slice(None, dim) for dim in data.shape)]

- Eric O Lebigot

---

太棒了...我想使用ogrid，但不知道它在多维度中如何工作。非常感谢。 - Cobry

Answer 2

为了让 @EOL 的好方法更加通用，可以在数组的形状内定义一个中心点。

array = np.random.random((100,100,100))
center = (30,10,25)
radius = 5.0
x, y, z = np.ogrid[-center[0]:array.shape[0]-center[0],-center[1] :array.shape[1]-center[1], -center[2]:array.shape[2]-center[2]]
indexes = numpy.transpose((x**2+y**2+z**2 <= radius**2).nonzero())

- Cobry