插值/调整 3D 数组

Question

我有一个 3D 数组，其中包含来自 mri 数据集的体素。模型可以沿一个或多个方向拉伸。例如。体素大小 (x,y,z) 可以是 0.5x0.5x2 mm。现在我想将 3D 数组重新采样为一个包含 1、1、1 毫米体素的数组。为此，我需要使 x/y 尺寸更小，z 尺寸更大，并插入体素值。 我的问题是：在 numpy 或 scipy 中是否有一个简单的函数可以对简单的 3d 数组进行这种重新采样？

从 *.nii 文件加载模型：

img = nib.load(sFileName)
array = np.array(img.dataobj).astype("uint8") # 3d array with e.g. 0.5,0.5,2 mm
# TODO resample

Answer 1

nd​​image.zoom

这可能是最好的方法，zoom method 正是为此类任务而设计的。

from scipy.ndimage import zoom
new_array = zoom(array, (0.5, 0.5, 2))

按指定因子更改每个维度的大小。如果数组的原始形状是(40, 50, 60)，那么新的形状将是(20, 25, 120)。

signal.resample_poly

SciPy 有一个large set of methods 用于信号处理。这里最相关的是decimate 和resample_poly。我在下面使用后者

from scipy.signal import resample_poly
factors = [(1, 2), (1, 2), (2, 1)]
for k in range(3):
    array = resample_poly(array, factors[k][0], factors[k][1], axis=k)

因子（必须是整数）是上采样和下采样的。那就是：

• (1, 2) 表示大小除以 2
• (2, 1) 表示大小乘以 2
• (2, 3) 表示先增加 2，然后再减少 3，因此大小乘以 2/3

可能的缺点：该过程在每个维度上独立发生，因此可能无法像 ndimage 方法那样考虑空间结构。

RegularGridInterpolator

这更需要动手操作，但也更费力，而且没有过滤的好处：直接下采样。我们必须为插值器制作一个网格，在每个方向上使用原始步长。插值器创建后，需要在新的网格上进行评估；它的调用方法采用另一种网格格式，用mgrid 准备。

values = np.random.randint(0, 256, size=(40, 50, 60)).astype(np.uint8)  # example

steps = [0.5, 0.5, 2.0]    # original step sizes
x, y, z = [steps[k] * np.arange(array.shape[k]) for k in range(3)]  # original grid
f = RegularGridInterpolator((x, y, z), values)    # interpolator

dx, dy, dz = 1.0, 1.0, 1.0    # new step sizes
new_grid = np.mgrid[0:x[-1]:dx, 0:y[-1]:dy, 0:z[-1]:dz]   # new grid
new_grid = np.moveaxis(new_grid, (0, 1, 2, 3), (3, 0, 1, 2))  # reorder axes for evaluation
new_values = f(new_grid)

缺点：例如，当一个维度减少 2 时，它实际上会每隔一个值下降一次，这就是简单的下采样。理想情况下，在这种情况下应该平均相邻值。在信号处理方面，decimation 中的低通滤波先于下采样。

Answer 2

您可以为此使用TorchIO：

import torchio as tio
image = tio.ScalarImage(sFileName)
resample = tio.Resample(1)
resampled = resample(image)

Answer 3

如果你想直接对MR图像进行重采样操作，可以使用antpy，它会对数据进行重采样，并改变体素间距信息。

    spacing = (1,1,1)
    import ants
    img = ants.image_read(file)
    filin = ants.resample_image(img,spacing,False,1)
    ants.image_write(filin,output)

Answer 4

最好的现代解决方案必须包括使用 Python 的 ITK 库。

假设

我假设您的数据针对每个对象进行了二值化！如果您正在处理颜色甚至灰度，您可能应该使用深度神经网络，但这超出了范围

方法

我会为 ITK 使用 morphological interpolation extension。您需要手动添加每个图像平面之间的间距。假设您的 XYZ 比率为 1:1:7，在每个图像平面之间添加 7 个空平面，并使用新矩阵作为函数的主要参数。

确保遵守库支持的dtype。我将np.uint8 用于我的数组dtype。

实施

from typing import Union

import itk
import numpy as np
from yaspin import yaspin
from yaspin.spinners import Spinners

INTERPOLATION_TEXT = "Interpolating for the first time"


@yaspin(Spinners.bouncingBall, text=INTERPOLATION_TEXT)
def morphological_interpolation(
    labeled_stack: np.ndarray, stack_depth_spacing: Union[float, int]
):
    spacing_array = np.zeros(
        (round(stack_depth_spacing), *labeled_stack.shape[1:]),
        dtype=labeled_stack.dtype,
    )

    spaced_stack = []
    for xy_plane in labeled_stack:
        spaced_stack.append(np.expand_dims(xy_plane, axis=0))
        spaced_stack.append(spacing_array)
    spaced_stack.pop(-1)
    spaced_stack = np.concatenate(spaced_stack)
    spaced_stack.astype(labeled_stack.dtype, copy=False)

    return itk.GetArrayFromImage(
        itk.morphological_contour_interpolator(itk.GetImageFromArray(spaced_stack))
    )

Answer 5

我在torch.nn.functional 中使用grid_sample。当您想在非均匀分布的点上重新采样 3d 数组时，这很方便。

以下是示例代码。您可以通过更改代码中的grid 来指定重采样点。

import numpy as np
import torch
import torch.nn.functional as F

def resize_by_grid_sample(x):

    dx = torch.linspace(-1, 1, 8)
    dy = torch.linspace(-1, 1, 32)
    dz = torch.linspace(-1, 1, 32)
    meshx, meshy, meshz = torch.meshgrid((dx, dy, dz))
    grid = torch.stack((meshx, meshy, meshz), 3)
    grid = grid.unsqueeze(0)

    x = x[np.newaxis, np.newaxis, :, :, :]
    x = torch.tensor(x, requires_grad=False)

    out = F.grid_sample(x, grid, align_corners=True)
    out = out.data.numpy()
    out = np.squeeze(out)

    return out


data = np.random.normal(size=(32,32,32)).astype(np.float32)

# grid_sample
resized_data = resize_by_grid_sample(data)

我参考了代码here。

align_corners=True 和align_corners=False 之间的区别实际上是我们是否在调整大小之前和之后对齐角像素。 nice figure here