从 3d 体积数据中提取 40 mm 立方体

Question

我从 dicom 格式创建了一个 3d numpy 数据数组（形状：133 x 512 x 512）。给定结节位置的中心点，我将如何提取尺寸为 40mm x 40mm x 40mm 的 3d 体积。我不确定像素到毫米的转换会如何发生？我在这里附上了数据。结节的位置是 (317, 363, 89)，即 (x,y,z)，其中 z 表示切片数。所以对于这个例子，结节位于第 89 层。 这是sample data。它是 nrrd 格式。在 dicom 标头信息中，切片厚度为 2.5，像素间距为 ['0.703125', '0.703125']。

Answer 1

看起来 NumPy 数组坐标被反转为 (z, y, x) 又名（平面、行、列），有关更多详细信息，请参阅 this scikit-image documentation。此外，让我们将体素间距设为 [2.5, 0.703125, 0.703125]，尽管我无法使用 pynrrd 找到此信息。但是，想象一下，我这些测量结果看起来是对的。这个单位是毫米/像素。

现在，您想将 40 毫米（或者说是 20 毫米的一半）转换为像素。您可以通过 除 除以毫米/像素来做到这一点，因此您得到的盒子大小为：

>>> 20 / np.array([2.5, 0.703125, 0.703125])
array([ 8.        , 28.44444444, 28.44444444])

这意味着框的顶部应该在：

>>> np.array([89, 363, 317]) - [8, 28, 28]
array([ 81, 335, 289])

和底部在：

>>> np.array([89, 363, 317]) + [8, 29, 29]
array([ 97, 392, 346])

（我在另一边设置了 29 以获得中心像素，也是因为 28.4 * 2 近似于 57，而不是 56。）

所以，要得到你的盒子，你会这样做：

>>> box40mm = data[81:97, 335:392, 289:346]

有关 NumPy 索引的更多详细信息，请参阅 NumPy indexing 文档。

将它们放在一个函数中，假设所有输入都是 NumPy 数组：

def get_box(data, centre_pixel, box_size, pixel_spacing):
    box_size_pixels = box_size / pixel_spacing
    box_size_left = np.round(box_size_pixels / 2).astype(int)
    box_size_right = np.round(box_size_pixels - box_size_left).astype(int)
    starts = centre_pixel - box_size_left
    ends = centre_pixel + box_size_right
    slices = tuple(
        slice(start, end)
        for start, end in zip(starts, ends)
    )
    return data[slices]

请注意，此代码不考虑了触及音量的边界。