通过 3d 数组沿任意线切割一定宽度的盒子

Question

我有一个大的 (600,600,600) numpy 数组，其中填充了我的数据。现在我想从中提取区域，在通过框的任意线周围具有给定宽度。

对于这条线，我在单独的 numpy 数组中拥有每个点的 x、y 和 z 坐标。所以假设这条线在数据框中有 35 个点，那么 x、y 和 z 数组的长度也都是 35。我可以通过使用这样的索引来提取沿线本身的点

extraction = data[z,y,x]

现在理想情况下，我想通过执行以下操作在其周围提取一个框

extraction = data[z-3:z+3,y-3:y+3,z-3:z+3]

但是因为 x、y 和 z 是数组，所以这是不可能的。我能想到的唯一方法是通过每个点的for循环，所以

extraction = np.array([])
for i in range(len(x)):
    extraction = np.append(extraction,data[z[i]-3:z[i]+3,y[i]-3:y[i]+3,z[i]-3:z[i]+3])

然后再对提取数组进行整形。但是，这非常慢，并且我想防止这个 for 循环中的每个切片之间存在一些重叠。

有没有一种简单的方法可以不使用 for 循环直接执行此操作？

编辑： 让我通过我想出的另一个想法来重新表述这个问题，这个想法也很慢。我有一条线穿过数据立方体。我有一个包含定义线的所有点的 x、y 和 z 坐标列表（坐标是数据立方体数组中的索引）。 例如，这些列表如下所示：

 x_line: [345 345 345 345 342 342 342 342 342 342 342 342 342 342 342 342]
 y_line: [540 540 540 543 543 543 543 546 546 546 549 549 549 549 552 552]
 z_line: [84 84 84 87 87 87 87 87 90 90 90 90 90 93 93 93]

如您所见，其中一些坐标是相同的，因为这些线是在不同的坐标中定义的，然后根据数据框的分辨率进行分箱。 现在我想屏蔽数据立方体中距离大于 3 个单元格的所有单元格。 对于沿线的单个点 (x_line[i], y_line[i], z_line[i])，这相对容易。我为数据立方体中的坐标创建了一个网格网格，然后创建了一个零掩码数组并将所有内容都满足条件为1：

data = np.random.rand(600,600,600)
x_box,y_box,z_box = np.meshgrid(n.arange(600),n.arange(600),n.arange(600))
mask = np.zeros(np.shape(data))

for i in range(len(x_line)):
    distance = np.sqrt((x_box-x_line[i])**2 + (y_box-y_line[i])**2 + (z_box-z_line[i])**2)
    mask[distance <= 3] = 1.

extraction = data[mask == 1.]

这样做的好处是掩码数组消除了重复提取的问题。但是，meshgrid 和距离计算都非常慢。那么是否可以直接在整条线上进行距离计算，而不必在每个线点上进行 for 循环，以便它直接掩盖距离任何线点 3 个单元格范围内的所有单元格?

Answer 1

这个怎么样？

# .shape = (N,)
x, y, z = ...

# offsets in [-3, 3), .shape = (6, 6, 6)
xo, yo, zo = np.indices((6, 6, 6)) - 3

# box indices, .shape = (6, 6, 6, N)
xb, yb, zb = x + xo[...,np.newaxis], y + yo[...,np.newaxis], z + zo[...,np.newaxis]

# .shape = (6, 6, 6, N)
extractions = data[xb, yb, zb]

这会提取一系列 6x6x6 的立方体，每个立方体都以 x、y 和 z 中的坐标为“中心”

这将产生重复的坐标，并在靠近边界的情况下失败

如果您将 xyz 保留在一个数组中，这会变得不那么冗长，并且您可以删除重复项：

# .shape = (N,3)
xyz = ...

# offsets in [-3, 3), .shape = (6, 6, 6, 3)
xyz_offset = np.moveaxis(np.indices((6, 6, 6)) - 3, 0, -1)

# box indices, .shape = (6, 6, 6, N, 3)
xyz_box = xyz + xyz_offset[...,np.newaxis,:]

if remove_duplicates:
    # shape (M, 3)
    xyz_box = xyz_box.reshape(-1, 3)
    xyz_box = np.unique(xyz_box, axis=0)

xb, yb, zb = xyz_box
extractions = data[xb, yb, zb]