我有一个二维数组a，包含一组 100 个 x、y、z 坐标：
[[ 0.81 0.23 0.52]
[ 0.63 0.45 0.13]
...
[ 0.51 0.41 0.65]]
我想创建一个 3D 二值图像b，在 x、y、z 维度上各有 101 个像素，坐标范围在 0.00 到 1.00 之间。
由a 定义的位置的像素应取值为 1，所有其他像素的值应为 0。
我可以使用b = np.zeros((101,101,101)) 创建一个形状正确的零数组，但是如何分配坐标并对其进行切片以使用a 创建零数组？
【问题讨论】：
标签： python numpy indexing binary ndimage
首先，从安全地将浮点数舍入为整数开始。在上下文中，请参阅this 问题。
a_indices = np.rint(a * 100).astype(int)
接下来，将b 中的那些索引分配给1。但要小心使用普通的list 而不是数组，否则会触发index arrays 的使用。似乎这种方法的性能与替代方法相当（感谢@Divakar！:-)
b[list(a_indices.T)] = 1
我创建了一个大小为 10 而不是 100，尺寸为 2 而不是 3 的小示例来说明：
>>> a = np.array([[0.8, 0.2], [0.6, 0.4], [0.5, 0.6]])
>>> a_indices = np.rint(a * 10).astype(int)
>>> b = np.zeros((10, 10))
>>> b[list(a_indices.T)] = 1
>>> print(b)
[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
[ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
【讨论】：
你可以这样做 -
# Get the XYZ indices
idx = np.round(100 * a).astype(int)
# Initialize o/p array
b = np.zeros((101,101,101))
# Assign into o/p array based on linear index equivalents from indices array
np.put(b,np.ravel_multi_index(idx.T,b.shape),1)
作业部分的运行时 -
让我们使用更大的网格来进行计时。
In [82]: # Setup input and get indices array
...: a = np.random.randint(0,401,(100000,3))/400.0
...: idx = np.round(400 * a).astype(int)
...:
In [83]: b = np.zeros((401,401,401))
In [84]: %timeit b[list(idx.T)] = 1 #@Praveen soln
The slowest run took 42.16 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 6.28 ms per loop
In [85]: b = np.zeros((401,401,401))
In [86]: %timeit np.put(b,np.ravel_multi_index(idx.T,b.shape),1) # From this post
The slowest run took 45.34 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 5.71 ms per loop
In [87]: b = np.zeros((401,401,401))
In [88]: %timeit b[idx[:,0],idx[:,1],idx[:,2]] = 1 #Subscripted indexing
The slowest run took 40.48 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 6.38 ms per loop
【讨论】：
(a * 100).astype(int) 非常危险！你可能会得到错误的索引！试试：a = np.arange(101) / 100、b = (100 * a).astype(int)、(b == np.arange(101, dtype=int)).all()，你会得到False！例如，由于浮点精度问题，数字 28 在 b 中出现了两次。
np.put 方法很有趣。我会对这之间的时序分析感兴趣，并且只使用list(a_indices.T)。