文章目录

  • 1. 变换到一个新的spacing
    • 1.1 使用Resample()的方法
      • 1.1.1 完整代码
      • 1.1.2 转换前后图像
    • 1.2 使用ResampleImageFilter()的方法
  • 2. 变换到一个新的size
    • 2.1 代码
    • 2.2 转换后图像
  • 3. 变换成和另一个图一致
    • 3.1 代码
  • 4. 相关内容说明
    • 4.1 体素和spacing
    • 4.2 Resample() 函数说明
      • 4.2.1 用法三
      • 4.2.2 用法二
    • 4.3 itk::simple::ResampleImageFilter函数说明

1. 变换到一个新的spacing

1.1 使用Resample()的方法

1.1.1 完整代码

def resample_volume(volume_path, interpolator = sitk.sitkLinear, new_spacing = [0.4, 0.4, 5]):    volume = sitk.ReadImage(volume_path)    original_spacing = volume.GetSpacing()    original_size = volume.GetSize()    new_size = [int(round(osz*ospc/nspc)) for osz,ospc,nspc in zip(original_size, original_spacing, new_spacing)]    resampled_image = sitk.Resample(volume, new_size, sitk.Transform(), interpolator,                         volume.GetOrigin(), new_spacing, volume.GetDirection(), 0,                         volume.GetPixelID())    sitk.WriteImage(resampled_image,"001.nii.gz")

关于sitk.Resample函数的详细解释，参考4.2 Resample() 部分

关于new_size的计算，可以参考4.1 体素和spacing部分内容理解

参考:

  • ITK论坛：Resample volume to specific voxel spacing – SimpleITK
  • SimpleITK官方jupyter：Transforms and Resampling，其中Cell17差不多就是这个方法

1.1.2 转换前后图像

其实看不到很明显的变化。

但是需要说明的是：space变小之后，像素数/分辨率就会变高，所以图像会变大。

  • 还是把图像想象成方格图，每个小立方体是一个像素点，每个小立方体的内部空白就是space（小立方体长宽高就是三个方向的space）
  • 则space减小，则小立方体就变小，则整个图里可以容纳的小立方体就更多，即分辨率更大。
  • 最直观的改变是：图像所占据的存储空间会变大。

1.2 使用ResampleImageFilter()的方法

感谢 Python SimpleITK.sitkNearestNeighbor() Examples，给了一个非常整洁的代码！

def resample_image(itk_image, out_spacing=[1.0, 1.0, 1.0], is_label=False):    original_spacing = itk_image.GetSpacing()    original_size = itk_image.GetSize()    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])))    ]    # 上述也可以直接用下面这句简写     out_size = [int(round(osz*ospc/nspc)) for osz,ospc,nspc in zip(original_size, original_spacing, out_spacing)]    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())    if is_label: # 如果是mask图像，就选择sitkNearestNeighbor这种插值        resample.SetInterpolator(sitk.sitkNearestNeighbor)    else: # 如果是普通图像，就采用sitkBSpline插值法        resample.SetInterpolator(sitk.sitkBSpline)    return resample.Execute(itk_image) # 使用示例resample_rs=resample_image(sitk.ReadImage("./001.nii.gz"))sitk.WriteImage(resample_rs,"../test.nii.gz")

参考：

  • 医学影像重采样
  • SimpleITK Github issue：How to resample sitk_image to a new spacing” />
    连同标签一起转换的话，如下图

    可以滑动一下，确认形状位置之类的特征没有发生变化。

    另外，数据大小变化非常明显。

    • 图像从10MB变成了1.61MB
    • 标签从42KB变成了10KB

    参考：

    • ITK study notes (9) simpleITK resampling data to a fixed size，这是一个别人博客的合集摘录
    • CSDN博客：Python SimpleItk库的医学图像重采样（resample）实现代码
    • Github gist：zivy/resample_to_same_size.ipynb
    • stack overflow：Resample DICOM Images to Size and Spacing and align to same Origin
    • ITK论坛：Resample to same Spacing and Size and Align to same Origin

    3. 变换成和另一个图一致

    3.1 代码

    代码：

    """以下代码是将裁剪后的mask和原始mask使用resample方法对齐，由于会直接使用参考图像的outsize，space，direction和origin等信息，因此这里不需要进行多余的设置"""mask=sitk.Resample(crop_mask,                   referenceImage=raw_mask,                   transform =sitk.Transform(),                   interpolator=sitk.sitkNearestNeighbor,                   defaultPixelValue = 0.0,                   outputPixelType = croped_img.GetPixelID())

    4. 相关内容说明

    4.1 体素和spacing


    图自Sparse Voxel Octree GI (构建体素)

    以LOD 0为例，假设这就是一个CT图，则可以看到，其被分成了很多小立方体

    • 每个立方体作为一个基本单位，但是由于每个立方体里有很多数据，为了从中提取出一个有代表性的值，会对每个小立方体里的数据进行一些操作，最常见的就是平均。
    • 即用均值代表这个小立方体中所有数据，而这个均值会作为这个小立方体的数值代表，即体素
    • 如果volume.GetSize()得到的结果是(512, 512, 29)，也就是说X和Y方向上各有512个space，Z方向上有29个space。也就有512x512x29个小立方体。
    • 体素的大小就是spacing的大小，例如：体素的大小是0.97mm 0.97mm 2.5mm。
    • 如果volume.GetSpacing()得到的结果是(0.4,0.4,5)，那就是说X和Y方向上每隔0.4mm就有一个点，Z方向上每隔5mm就有一个点，即每个小立方体的大小是(0.4mm,0.4mm,5mm)。应该称之为长方体
    • 则原始图像的大小其实就是(512*0.4, 512*0.4, 29*5)，带单位的mm。
    • 医学图像是和真实空间的物理尺寸相对应的
    • 则对应的，如果已知新的space，想求出新的size，则就可以使用原始图像的大小 除以 新的space

    部分用语出自：医学影像重采样

    4.2 Resample() 函数说明

    python版本的说明里，可以接受三种参数调用，如下，用到一种写一种，因此不会全写

    Resample ( Image image1,            Transform transform = itk::simple::Transform(),           InterpolatorEnum interpolator = itk::simple::sitkLinear,           double defaultPixelValue = 0.0,           PixelIDValueEnum outputPixelType = itk::simple::sitkUnknown,           bool useNearestNeighborExtrapolator = false);Resample ( Image image1,           Image referenceImage,           Transform transform = itk::simple::Transform(),           InterpolatorEnum interpolator = itk::simple::sitkLinear,           double defaultPixelValue = 0.0,           PixelIDValueEnum outputPixelType = sitkUnknown,           bool useNearestNeighborExtrapolator = false);Resample ( const Image& image1,             VectorUInt32 size,            Transform transform = itk::simple::Transform(),            InterpolatorEnum interpolator = itk::simple::sitkLinear,           VectorDouble outputOrigin = std::vector<double>(3, 0.0),           VectorDouble outputSpacing = std::vector<double>(3, 1.0),           VectorDouble outputDirection = std::vector<double>(),           double defaultPixelValue = 0.0,           PixelIDValueEnum outputPixelType = sitkUnknown,           bool useNearestNeighborExtrapolator = false);

    4.2.1 用法三

    第三种调用方式适用于，已知输出的size或者输出的space时

    Resample ( const Image& image1,             VectorUInt32 size,            Transform transform = itk::simple::Transform(),            InterpolatorEnum interpolator = itk::simple::sitkLinear,           VectorDouble outputOrigin = std::vector<double>(3, 0.0),           VectorDouble outputSpacing = std::vector<double>(3, 1.0),           VectorDouble outputDirection = std::vector<double>(),           double defaultPixelValue = 0.0,           PixelIDValueEnum outputPixelType = sitkUnknown,           bool useNearestNeighborExtrapolator = false);"""image1：要进行resample的图像size：向量，经过resample操作后图像的size,（x,y,z）`image.GetSize() -> (512, 512, 29)`transform： 一般默认就是  `sitk.Transform()`interpolator:插值方式，比如`sitk.sitkLinear`或者`sitk.sitkNearestNeighbor`等outputOrigin：输出图像的起点/原点，一般就和原图一致。`image.GetOrigin()->(-114.1510009765625, -127.27660369873047, -565.5)`outputSpacing：输出图像的spacing，需要自己指定，如果是3d图像，那就是维度为3的向量outputDirection ：一般也是和原图保持一致，例如：`image.GetDirection()->(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)`defaultPixelValue：默认填充的像素值，一般就是0，不需要改动（背景0黑色）outputPixelType：输出像素的数据类型，默认是`sitkUnknown`，也可以是`sitkUInt8`等"""

    其中

    • 关于插值方式，详见： InterpolatorEnum，常见的就是sitkNearestNeighbor ，sitkLinear ，还有sitkBSpline 等
    • 关于方向的表示，根据ITK文档——GetDirection()。对于2D图像，方向会被表示为一个 2 × 2 2\times 2 2×2的矩阵，对于3d图像，方向会被表示为一个 3 × 3 3\times 3 3×3的矩阵。矩阵一般以一维数组作为主要展现形式、
    • 关于像素类型（其实就是值的数据类型），根据GetPixelIDValueAsString()

      接着继续找，可以发现PixelIDValueEnum，其中有一个枚举表格，列除了ITK支持的数据类型，如下。也可以参考其对应的头文件定义 sitkPixelIDValues.h中第93行。

    C++版本的说明，参考： itk::simple::Resample

    4.2.2 用法二

    第二种调用方式适用于，有要对齐的图像时（即要把当前的图像格式转成和另一个图像一致）

    Resample ( Image image1,           Image referenceImage,           Transform transform = itk::simple::Transform(),           InterpolatorEnum interpolator = itk::simple::sitkLinear,           double defaultPixelValue = 0.0,           PixelIDValueEnum outputPixelType = sitkUnknown,           bool useNearestNeighborExtrapolator = false);
    • 其他参数和用法三都一致，去4.2.1 用法三看就行了，
    • referenceImage这个参数就是一个ITK格式的image
    • 这种调用方式会将输入图像的output size, origin, spacing和direction，变得与要参考/对齐的图像referenceImage一致
    • 因此可以看到，相对于用法三，没有了由outsize，space，direction和origin等信息（都使用参考图像的了，不用我们自己设置了）

    详细说明可以参考C++文档，这里

    以及这里

    4.3 itk::simple::ResampleImageFilter函数说明

    python的函数说明没有很多有效内容，但是提到了

    See alsoitk::ResampleImageFilter for the Doxygen on the original ITK class.Definition at line 53 of file sitkResampleImageFilter.h.

    因此去查看C++的文档，主要是sitkResampleImageFilter.h.这个文件，根据sitkResampleImageFilter.h，其中有一部分代码为：

         private:       using MemberFunctionType = Image (Self::*)( const Image& image1 );         friend struct detail::DualExecuteInternalAddressor<MemberFunctionType>;       template <class TImageType1, class TImageType2> Image DualExecuteInternal ( const Image& image1 );       friend struct detail::DualExecuteInternalVectorAddressor<MemberFunctionType>;       template <class TImageType1, class TImageType2> Image DualExecuteInternalVector ( const Image& image1 );         std::unique_ptr<detail::DualMemberFunctionFactory<MemberFunctionType> > m_DualMemberFactory;             /*  */       std::vector<uint32_t>  m_Size{std::vector<uint32_t>(3, 0)};         Transform  m_Transform{itk::simple::Transform()};         /*  */       InterpolatorEnum  m_Interpolator{itk::simple::sitkLinear};         /*  */       std::vector<double>  m_OutputOrigin{std::vector<double>(3, 0.0)};         /*  */       std::vector<double>  m_OutputSpacing{std::vector<double>(3, 1.0)};         /* Passing a zero sized array, defaults to identiy matrix. The size of the array must exactly match the direction matrix for the dimension of the image. */       std::vector<double>  m_OutputDirection{std::vector<double>()};         double  m_DefaultPixelValue{0.0};         PixelIDValueEnum  m_OutputPixelType{itk::simple::sitkUnknown};         bool  m_UseNearestNeighborExtrapolator{false};

    可以看到，这个ResampleImageFilter的私有内容，其实和Resample()很像，都是：

    • Size：输出向量的维度
    • Transform：转换方式
    • Interpolator：插值方法
    • OutputOrigin：输出的原点
    • OutputSpacing：输出spacing
    • OutputDirection：输出方向
    • DefaultPixelValue：默认填充的像素值
    • OutputPixelType：输出的像素类型
    • UseNearestNeighborExtrapolator：是否使用这个东西

    另外，在 itk::simple::ResampleImageFilter文档中，也可以在Member Data Documentation中看到上述参数介绍

    但是因为是私有成员，因此需要通过相应的set修改器函数和get访问器函数来获取或者修改相应的成员变量。
    （真的是要复习C++了呀）