SimpleITK - 冠状/矢状视图大小问题

Question

我正在尝试使用 SimpleItk 库从 DICOM 格式的 CTA 中提取所有三个视图（轴向、矢状和冠状）。

我可以从给定目录中正确读取该系列：

    ...
    import SimpleITK as sitk
    ...
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(input_dir)
    reader.SetFileNames(dicom_names)
    # Execute the reader
    image = reader.Execute()
    ...

然后，使用 this questions 中所述的 numpy 数组，我可以提取并保存 3 个视图。

    ...
    image_array = sitk.GetArrayFromImage(image)
    ...
    for i in range(image_array.shape[0]):
        output_file_name = axial_out_dir + 'axial_' + str(i) + '.png'
        logging.debug('Saving image to ' + output_file_name)
        imageio.imwrite(output_file_name, convert_img(image_array[i, :, :], axial_min, axial_max), format='png')
    ...

另外两个是通过保存image_array[:, i, :]和image_array[:, :, i]制作的，而convert_img(..)是一个只转换数据类型的函数，所以不会改变任何形状。

但是，冠状和矢状视图被拉伸、旋转并带有宽黑色带（在某些切片中它们非常宽）。

这是 Slicer3d 的屏幕截图：

虽然这是我的代码的输出：

轴向

矢状面

冠状动脉

图像形状为 512x512x1723，这导致轴向 png 为 512x512 像素，冠状和矢状为 512x1723，因此这似乎是正确的。

我应该尝试使用 PermuteAxes 过滤器吗？问题是我找不到任何关于它在 python 中使用的文档（由于文档页面中的 404 而不是其他语言）

还有办法提高对比度吗？我使用了 simpleitk 的 AdaptiveHistogramEqualization 过滤器，但它比 Slicer3D 可视化效果差得多，而且速度很慢。

感谢您的帮助，谢谢！

Answer 1

当您将 SimpleITK 图像转换为 NumPy 数组时，所有像素间距信息都会丢失（如上面的 cmets 所示）。如果您在 SimpleITK 中执行所有操作，它会保留该间距信息。

在 SimpleITK 中使用 python 的数组切片很容易从图像中提取 X、Y 和 Z 的切片：

import SimpleITK as sitk

# a blank test image
img = sitk.Image([100, 101, 102], sitk.sitkUInt8)

# non-uniform spacing, for illustration
img.SetSpacing([1.0, 1.1, 1.2])

# select the 42nd Z slice
zimg = img[:, :, 42]

#select the 0th X slice
ximg = img[0, :, :]

#select the 100th Y slice
yimg = img[:, 100, :]

#print the spacing to show it's retained
print(yimg.GetSpacing())

Answer 2

如果有人需要回答我自己的问题。

鉴于我需要在深度学习框架中使用切片并进行数据增强，我需要以 (1.0, 1.0, 1.0) 的新间距对它们进行重新采样。

用这个函数解决了：

def resample_image(itk_image, out_spacing=(1.0, 1.0, 1.0)):
    """
    Resample itk_image to new out_spacing
    :param itk_image: the input image
    :param out_spacing: the desired spacing
    :return: the resampled image
    """
    # get original spacing and size
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()
    # calculate new size
    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))
    ]
    # instantiate resample filter with properties and execute it
    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resample.SetInterpolator(sitk.sitkNearestNeighbor)
    return resample.Execute(itk_image)

然后按照原始问题中的说明使用numpy arrays 进行保存。

Answer 3

我可能会迟到，但您可以使用 Torchio。我认为适合您的情况的一个很好的解决方案是使用随 TorchIO 一起安装的 CLI 工具：

$ tiohd your_image.nii.gz
ScalarImage(shape: (1, 512, 512, 1723); spacing: (0.50, 0.50, 1.00); orientation: RAS+; memory: 1.7 GiB; dtype: torch.ShortTensor)
$ torchio-transform your_image.nii.gz Resample one_iso.nii.gz
$ tiohd one_iso.nii.gz
ScalarImage(shape: (1, 256, 256, 1723); spacing: (1.00, 1.00, 1.00); orientation: RAS+; memory: 430.8 MiB; dtype: torch.ShortTensor)

这是因为 1 毫米是 Resample 变换的默认目标分辨率。

当然，您也可以使用 TorchIO 的普通 Python 接口来操作图像。

免责声明：我是 TorchIO 的主要开发者。