关于3Dslicer和SimpleITK各自的origin，direction，spacing，（和size）-Toy模板网

这篇具有很好参考价值的文章主要介绍了关于3Dslicer和SimpleITK各自的origin，direction，spacing，（和size）。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

众所周知3Dslicer和SimpleITK对三维医学图像的某些参数并不一致

。。。

在正文之前首先要介绍IJK坐标系和RAS坐标系

IJK坐标系是对图像本身而言，三维图像是一个一个一个体素组成，IJK坐标就是这些体素的坐标或者索引，IJK一般只有非负整数值。如果以三维矩阵获取图像，那么IJK就是这个矩阵的索引。在下文的GetArrayFromImage()和arrayFromVolume()函数就是以三维矩阵获取图像。

RAS坐标系也叫世界坐标系，是物理空间中的绝对坐标系或者叫大地坐标系，如果不习惯叫RAS的话，可以记住X指向R面，Y指向A面，Z指向S面（我也经常搞混，不如就XYZ等同于RAS）。这个坐标系原点是0，长度是按mm，cm这种现实中的度量。举个例子，图像的位置对于RAS坐标系可能有各种千奇百怪的状态，如下图

3d slicer对比图像,图像处理,python

这也就是为什么会有origin和direction这两个量

简而言之

size是体素图像的[I, J, K]，I长，J宽，K高

origin是体素图像的原点坐标对于RAS坐标系的位置

spacing体素图像每个维度上像素或体素之间的间距，一般单位mm

direction是图像的方向，即图像坐标系相对世界坐标系的角度，一般采用的是方向余弦矩阵

3Dslicer和SimpleITK里这些参数相互转化的代码见文末

在SimpleITK中，图像被read成了<class 'SimpleITK.SimpleITK.Image'>的类，这个类的函数就可以获取origin，direction，spacing，和size

import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)

# size是[I, J, K]形list，I长，J宽，K高
print('image1.GetSize(): ', image1.GetSize())

# origin图像的原点坐标对于RAS坐标系的位置，顺序[R,A,S]
print('image1.GetOrigin() :', image1.GetOrigin())

# spacing每个维度上像素或体素之间的间距，单位mm，顺序[I, J, K]
print('image1.GetSpacing(): ', image1.GetSpacing())

# direction图像的方向，即图像坐标系相对世界坐标系的角度，采用的是方向余弦矩阵按行主序展成一# 维
print('image1.GetDirection(): ', image1.GetDirection())

# 因为是行主序，可以用reshape变成方向余弦矩阵
direction = np.array(list(itkimage.GetDirection())).reshape(3, 3)
print("directionmatrix",direction)

在slicer中，图像的形式是'vtkMRMLScalarVolumeNode'或者'vtkMRMLLabelMapVolumeNode'，这些是slicer特有的node的形式，可以直接获取origin，spacing

# 获取所有'vtkMRMLScalarVolumeNode'的名称（字符串）列表
ModelNames = getNodesByClass('vtkMRMLScalarVolumeNode')
ModelName = ModelNames[0] # 这里你需要确保列表里的ModelNames[？]是上述的类型
# 用getNode方法按名称获取节点
volumeNode = getNode(ModelName.GetName())
# 获取origin
volumeNode.GetOrigin()
# 获取spacing
volumeNode.GetSpacing()
print("volumeNode.GetOrigin()",volumeNode.GetOrigin())
print("volumeNode.GetSpacing()",volumeNode.GetSpacing())

注意注意注意！！！

这里是第一个坑点，slicer里origin的符号和sitk不一样，在I和J是反的，如果我们希望sitk里获取的origin和slicer一样，需要这样翻符号。

import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)

origin = np.array(list(image1.GetOrigin()))
origin[0] = -origin[0]
origin[1] = -origin[1]
print('origin',origin)

spacing在sitk和slicer里顺序和符号都一样。

至于size，我没有找到slicer中有直接从node中获取的方法，但是可以先获取矩阵，然后查看矩阵的shape，这就等同于sitk里的size

import numpy as np
# 获取所有'vtkMRMLScalarVolumeNode'的名称（字符串）列表
ModelNames = getNodesByClass('vtkMRMLScalarVolumeNode')
ModelName = ModelNames[0] # 这里你需要确保列表里的ModelNames[？]是上述的类型
# 用getNode方法按名称获取节点
volumeNode = getNode(ModelName.GetName())
a = arrayFromVolume(volumeNode)
a = np.array(a)
print(a.shape)

注意注意注意！！！

这里是第二个坑点，在3Dslicer和SimpleITK里的这种arrayfrom的函数，返回的都是[K,J,I]的顺序

所以一般都需要transpose一下

# 对sitk
import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)
label = sitk.GetArrayFromImage(image1).transpose(2,1,0)

# 对slicer
volumeNode = getNode(ModelName.GetName())
a = arrayFromVolume(volumeNode)
a = np.array(a)
a.transpose(2,1,0)

到了direction这里，事情就很麻烦了

slicer里有3个函数分别获取I，J，K轴对于RAS坐标系的方向

# 因为这些函数是直接改变list，所以要先定义后使用
I_AXISToRAS=[0.0,0.0,0.0]
J_AXISToRAS=[0.0,0.0,0.0]
K_AXISToRAS=[0.0,0.0,0.0]

volumeNode.GetIToRASDirection(I_AXISToRAS)
volumeNode.GetJToRASDirection(J_AXISToRAS)
volumeNode.GetKToRASDirection(K_AXISToRAS)

print(I_AXISToRAS)
print(J_AXISToRAS)
print(K_AXISToRAS)

返回的是I对R，A，S分别的余弦值，J对R，A，S分别的余弦值，K对R，A，S分别的余弦值

再对比前面sitk得到的direction，会发现有很大出入，有2行负号是相反的

import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)
direction = np.array(list(image1.GetDirection())).reshape(3, 3)
print('direction\n',direction)

>>> direction
[[ 9.94355000e-01 -4.89225598e-08 -1.06104344e-01]
 [-3.40554847e-05  9.99999948e-01 -3.19611396e-04]
 [ 1.06104348e-01  3.21420675e-04  9.94354950e-01]]

I_AXIS=[0.0,0.0,0.0]
J_AXIS=[0.0,0.0,0.0]
K_AXIS=[0.0,0.0,0.0]

volumeNode.GetIToRASDirection(I_AXIS)
volumeNode.GetJToRASDirection(J_AXIS)
volumeNode.GetKToRASDirection(K_AXIS)
>>> print(I_AXIS)
[-0.9943550000701419, 3.4055484670310706e-05, 0.1061043480529044]
>>> print(J_AXIS)
[4.892255980967634e-08, -0.9999999483443722, 0.000321420675237069]
>>> print(K_AXIS)
[0.10610434399770712, 0.0003196113955094529, 0.9943549497203562]

那么如何把这3个向量拼成方向余弦矩阵呢，当时做到这里我直接妈妈生的。。。

SimpleITK官网

3d slicer对比图像,图像处理,python

Slicer官网

3d slicer对比图像,图像处理,python

都没有详细说这里面的元素如何对应

去查了一下方向余弦矩阵的概念

做控制要知道的刚体旋转知识（四）旋转矩阵/方向余弦矩阵 - 知乎 (zhihu.com)

3d slicer对比图像,图像处理,python

这篇文章里面单位向量 I, J, K 表示为世界坐标系（global）的XYZ轴，三个小写的向量 i, j, k是本体坐标系的三个单位向量

其实slicer里获取到的是方向余弦矩阵的列向量

那么slicer里的方向余弦矩阵可以这样得到

# ！！！注意把direction拼成二维！！！
I_AXIS_ras_slicer = [0.0,0.0,0.0]
J_AXIS_ras_slicer = [0.0,0.0,0.0]
K_AXIS_ras_slicer = [0.0,0.0,0.0]

labelmapNode.GetIToRASDirection(I_AXIS_ras_slicer)
labelmapNode.GetJToRASDirection(J_AXIS_ras_slicer)
labelmapNode.GetKToRASDirection(K_AXIS_ras_slicer)

direction = np.array([I_AXIS_ras_slicer,J_AXIS_ras_slicer,K_AXIS_ras_slicer]).reshape(3,3).swapaxes(0, 1)
print('direction',direction)

那么上面的结果中sitk的I轴和K轴负号相反，K轴不变

对于sitk，如果我们希望把得到的direction在slicer里使用，则应该这样计算方向余弦矩阵

# 如果我们希望sitk里获取的origin和slicer一样，需要这样翻符号。

import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)

origin = np.array(list(image1.GetOrigin()))
origin[0] = -origin[0]
origin[1] = -origin[1]
print('origin',origin)

# 如果我们希望把得到的direction在slicer里使用，则应该这样计算方向余弦矩阵
import SimpleITK as sitk
imagepath = "…………/图像名.nii" # 图像路径
image1 = sitk.ReadImage(imagepath)
direction = np.array(list(image1.GetDirection())).reshape(3, 3)
direction[0, 0] = -direction[0, 0]
direction[0, 1] = -direction[0, 1]
direction[1, 0] = -direction[1, 0]
direction[1, 1] = -direction[1, 1]
direction[0, 2] = -direction[0, 2]
direction[1, 2] = -direction[1, 2]

print('direction',direction)

如果我们想把slicer里的结果用于sitk，则应该这样产生origin向量和direction向量

labelmapNode_origin = np.array(labelmapNode.GetOrigin())
labelmapNode_origin[0] = -labelmapNode_origin[0]
labelmapNode_origin[1] = -labelmapNode_origin[1]
# ！！！注意这里sitk要翻负号！！！
labelmapNode_I_direction = [0.0,0.0,0.0]
labelmapNode_J_direction = [0.0,0.0,0.0]
labelmapNode_K_direction = [0.0,0.0,0.0]
labelmapNode.GetIToRASDirection(labelmapNode_I_direction)
labelmapNode.GetJToRASDirection(labelmapNode_J_direction)
labelmapNode.GetKToRASDirection(labelmapNode_K_direction)
# ！！！注意sitk的direction是一维！！！
labelmapNode_direction = np.array([-labelmapNode_I_direction[0],-labelmapNode_J_direction[0],-labelmapNode_K_direction[0],
                                  -labelmapNode_I_direction[1],-labelmapNode_J_direction[1],-labelmapNode_K_direction[1],
                                  labelmapNode_I_direction[2],labelmapNode_J_direction[2],labelmapNode_K_direction[2]])

labelmapNode_spacing = np.array(labelmapNode.GetSpacing())

文章来源地址https://www.toymoban.com/news/detail-808826.html

到了这里，关于关于3Dslicer和SimpleITK各自的origin，direction，spacing，（和size）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！