SimpleITK使用——1. 进行Resample/Resize操作

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部分内容理解

参考:

1.1.2 转换前后图像

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

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

SimpleITK使用——1. 进行Resample/Resize操作

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")

参考:

2. 变换到一个新的size

其实和上面大同小异,只是这次是通过新的size来计算新的space,如下

2.1 代码

def resample2size(img,new_size=[192,192,176]):
	original_size=img.GetSize()
	original_spacing = img.GetSpacing()
	
	new_space =[int(round(old_size*old_space/new_size)) for old_size,old_space,new_size in zip(original_size,original_spacing,dst_size)]

	resampled_img=sitk.Resample(img,new_size,sitk.Transform(),sitk.sitkLinear,
                           img.GetOrigin(),new_space,img.GetDirection(),0,img.GetPixelID())
    return resampled_img

这里就不再赘述ResampleImageFilter的写法了,

2.2 转换后图像

上半部分是原图,下半部分是转换后的
SimpleITK使用——1. 进行Resample/Resize操作
连同标签一起转换的话,如下图
SimpleITK使用——1. 进行Resample/Resize操作
可以滑动一下,确认形状位置之类的特征没有发生变化。

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

参考:

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

SimpleITK使用——1. 进行Resample/Resize操作
图自
Sparse Voxel Octree GI (构建体素)

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

  • 每个立方体作为一个基本单位,但是由于每个立方体里有很多数据,为了从中提取出一个有代表性的值,会对每个小立方体里的数据进行一些操作,最常见的就是平均。
  • 即用均值代表这个小立方体中所有数据,而这个均值会作为这个小立方体的数值代表,即体素
  • 如果volume.GetSize()得到的结果是(512, 512, 29),也就是说XY方向上各有512个space,Z方向上有29个space。也就有512x512x29个小立方体。
  • 体素的大小就是spacing的大小,例如:体素的大小是0.97mm 0.97mm 2.5mm
  • 如果volume.GetSpacing()得到的结果是(0.4,0.4,5),那就是说XY方向上每隔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`等
"""

其中

C++版本的说明,参考: itk::simple::Resample
SimpleITK使用——1. 进行Resample/Resize操作

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);

详细说明可以参考C++文档,这里
SimpleITK使用——1. 进行Resample/Resize操作
以及这里
SimpleITK使用——1. 进行Resample/Resize操作

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

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

See also
itk::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中看到上述参数介绍
SimpleITK使用——1. 进行Resample/Resize操作
但是因为是私有成员,因此需要通过相应的set修改器函数和get访问器函数来获取或者修改相应的成员变量。
(真的是要复习C++了呀)

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
扎眼的阳光的头像扎眼的阳光普通用户
上一篇 2023年3月25日
下一篇 2023年3月25日

相关推荐