SimpleITK使用——1. 进行ResampleResize *** 作

SimpleITK使用——1. 进行ResampleResize *** 作,第1张

文章目录
  • 1. 变换到一个新的spacing
    • 1.1 使用Resample()的方法
      • 1.1.1 完整代码
      • 1.1.2 转换前后图像
    • 1.2 使用ResampleImageFilter()的方法
  • 2. 变换到一个新的size
    • 2.1 代码
    • 2.2 转换后图像
  • 3. 相关内容说明
    • 3.0 体素和spacing
    • 3.1 Resample() 函数说明
    • 3.2 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函数的详细解释,参考3.1 Resample() 部分

关于new_size的计算,可以参考3.0 体素和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? #561
  • SimpleITK官方培训文档:images_and_resampling,直接去page 12就是Resampling
  • stack overflow:SimpleITK Resize images
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的写法了,

  • 只需要修改new_space这个变量,其他基本和1.2 使用ResampleImageFilter()的方法一样
2.2 转换后图像

上半部分是原图,下半部分是转换后的

连同标签一起转换的话,如下图

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

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

  • 图像从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.0 体素和spacing


图自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

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

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

3.2 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中看到上述参数介绍

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

欢迎分享,转载请注明来源:内存溢出

原文地址: https://outofmemory.cn/langs/796512.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-05-06
下一篇 2022-05-06

发表评论

登录后才能评论

评论列表(0条)

保存