从narray中提取Block Around Index

Question

我有一个形状为(80, 180, 144, 160, 11)的5D数组(80个尺寸为180*144*160的3d图像，每个图像有11个通道)和一组指向该形状为(n, 4)的数组的索引(即n索引指向我感兴趣的图像和3d像素)。

现在回到问题，我想提取"block "形状为(18, 18, 20)，以每个指数为中心，保留所有通道。这将生成形状为(n, 18, 18, 20, 11)的数组。此外，如果一个索引太接近3d图像的边界，不适合整个块，那么我想要0-pad图像。

我已经设法在每个索引上使用for循环自己做到这一点，但不幸的是，性能相当差(n=100 ~10 s)。我需要为ns在10 000 - 1 000 000的范围内这样做，所以我的解决方案不是一个真正的选择。

我在images中给出图像的尝试和block_indices中的索引:

block_shape = (18, 18, 20)
blocks = np.empty((0,) + block_shape + (11,))
for index in block_indices:
block = np.pad(images[index[0]], ((block_shape[0], block_shape[0]),
(block_shape[1], block_shape[1]),
(block_shape[2], block_shape[2]),
(0, 0)))[index[1]+int(block_shape[0]/2):index[1]+int(3*block_shape[0]/2),
index[2]+int(block_shape[1]/2):index[2]+int(3*block_shape[1]/2),
index[3]+int(block_shape[2]/2):index[3]+int(3*block_shape[2]/2),
...]
blocks = np.append(blocks, block[np.newaxis, ...], axis=0)

我认为这可能可以用切片和花哨的数组索引来快速完成，但我试图无济于事。你有什么建议可以更快地完成这项工作吗?提前感谢!

PS:给出的数字可能会有所不同，但应该能让你大致了解规模。

Answer 1

对于将来想做同样事情的人

我已经设法想出了另一个解决方案，它更快，规模更好。它涉及到"移位"的使用;块矩阵，np.tile，平坦和一些重塑。需要注意的是，块的索引需要在长度为n的1D数组中给出，其中每个索引对应于3d图像的扁平数组中的索引。但是，可以很容易地在这些不同的表示之间进行转换。

为简洁起见，我将只解释该方法的主要概念，然后发布一个工作代码示例，如下:

主要概念:

1. 首先，我们将images数组平整化，使其形状为(80*180*144*160,11)。
现在我们需要意识到，我们所追求的块可以根据可预测的模式从扁平数组中访问，该模式仅根据块的位置移动。
2. 这些元素可以用np.take去掉，只要我们知道索引。
3. 最后，np.take的结果可以被重塑成一个块数组。

工作代码示例:

# Make a 3D-image which enumerates all pixels.
image_pixel_enumeration = np.arange(180*144*160).reshape(180, 144, 160)
# Get the index pattern of a block.
block_shifts = image_pixel_enumeration[:block_shape[0], :block_shape[1], :block_shape[2]].flatten() 
- image_pixel_enumeration[int(block_shape[0]/2), int(block_shape[1]/2), int(block_shape[2]/2)]
# Tile an array with the pattern, one for each block.
block_shifts = np.tile(block_shifts, (len(block_indices), 1))
# Tile an array with the block center indices add to them the pattern.
validation_data_indices = np.tile(block_indices, (np.prod(block_shape), 1)).transpose() + block_shifts
# Take out elements.
validation_data = np.take(x_test.reshape((-1, 11)), validation_data_indices.flatten(), 0, mode='clip')
# Reshape into blocks.
validation_data = validation_data.reshape((-1,) + block_shape + (11,))

此方法(在我的机器上)分别为10,100和1,000个索引花费大约0.1 s, 0.2 s和1.4 s，而旧方法为相同数量的索引花费大约1 s, 16 s和900 s。一个巨大的进步!

p。请注意，此解决方案不能解决块超出原始图像的问题，并且在这些情况下可能会从错误的图像或错误的切片中挑选像素。