SimpleITK的IJK坐标和RAS坐标，长度，角度转化并将结果用于3DSlicer

在上一篇文章里提到SimpleITK和Slicer的R轴和A轴很有可能相反，那么必须在负号上做一些修改才能保证SimpleITK的计算结果在Slicer里能正确使用

IJK坐标和RAS坐标转化ras2ijk2ras.py

import SimpleITK as sitk
import numpy as np

segpath = "……/文件名.nii.gz" # label路径
itkimage = sitk.ReadImage(segpath)


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

spacing = np.array(list(itkimage.GetSpacing()))
print('spacing',spacing)

print("itkimage.GetDirection()",itkimage.GetDirection())
direction = np.array(list(itkimage.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)

def RAS_to_IJK(ras, origin, direction, spacing):
    ras = ras - origin
    ras = np.matmul(np.linalg.inv(direction), ras)
    ijk = ras / spacing

    return ijk


def IJK_to_RAS(ijk, origin, direction, spacing):
    ras = ijk * spacing
    ras = np.matmul(direction, ras)
    ras = ras + origin

    return ras
if __name__ == "__main__":
    
    print("RAS to ijk:")
    ras = np.array(
        [-1.1243574619293213, 23.990598678588867, 63.46666717529297]
    )
    print("input RAS : ", ras)
    ijk = RAS_to_IJK(
        ras, origin=origin, direction=direction, spacing=spacing)
    print("output ijk : ", ijk)

    print("ijk to RAS:")
    ijk = np.array(
        [320., 320.,  72.]
    )
    print("input ijk : ", ijk)
    ras = IJK_to_RAS(
        ijk, origin=origin, direction=direction, spacing=spacing)
    print("output RAS : ", ras)

参考文章slicer软件中RAS转换为像素坐标方法

RAS转IJK一共3步

1. 减原点偏移
2. 左乘方向余弦矩阵的逆
3. 除以体素长度

IJK转RAS是上面的逆过程

1. 乘体素长度
2. 左乘方向余弦矩阵
3. 加原点偏移

基于这个原理，还可以写出在sitk中计算RAS长度和角度的函数（角的单位用度）

ijk2ras_len_deg.py

import SimpleITK as sitk
import numpy as np

segpath = "……/文件名.nii.gz" # label路径
itkimage = sitk.ReadImage(segpath)

# Get origin, spacing, and direction
origin = np.array(list(itkimage.GetOrigin()))
origin[0] = -origin[0]
origin[1] = -origin[1]
print('origin: ',origin)
spacing = np.array(list(itkimage.GetSpacing()))
print('spacing :',spacing)
direction = np.array(list(itkimage.GetDirection())).reshape(3, 3)
print()
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]

# Calculate orientation of IJK axes in RAS coordinate system
IToRASDirection = direction[:, 0]  # direction vector of the I-axis
JToRASDirection = direction[:, 1]  # direction vector of the J-axis
KToRASDirection = direction[:, 2]  # direction vector of the K-axis
print("direction: ", direction)
print("Orientation of I-axis in RAS coordinate system:", IToRASDirection)
print("Orientation of J-axis in RAS coordinate system:", JToRASDirection)
print("Orientation of K-axis in RAS coordinate system:", KToRASDirection)

####计算路径ras长度
def IJK_to_RAS_length(ijk_ray_start, ijk_ray_end, origin, direction, spacing):
    ras_ray_start = ijk_ray_start * spacing
    ras_ray_start = np.matmul(direction, ras_ray_start)
    ras_ray_start = ras_ray_start + origin

    ras_ray_end = ijk_ray_end * spacing
    ras_ray_end = np.matmul(direction, ras_ray_end)
    ras_ray_end = ras_ray_end + origin

    ras_length = np.linalg.norm(ras_ray_start - ras_ray_end)
    return ras_length
####计算路径ras长度

####计算路径ras角度
def IJKp_to_IJKaxesinRAS_angle(ijk_ray_start, ijk_ray_end, origin, direction, spacing):
    # 先将IJK点坐标转为RAS点坐标
    ras_ray_start = ijk_ray_start * spacing
    ras_ray_start = np.matmul(direction, ras_ray_start)
    ras_ray_start = ras_ray_start + origin

    ras_ray_end = ijk_ray_end * spacing
    ras_ray_end = np.matmul(direction, ras_ray_end)
    ras_ray_end = ras_ray_end + origin
    # 先将IJK点坐标转为RAS点坐标

    L_origin0thPoint = ras_ray_end
    L_origin1thPoint = np.array(ras_ray_start)
    L_GETNODE_minus_GETNODE = L_origin1thPoint - L_origin0thPoint
    L_originArray = np.array(L_GETNODE_minus_GETNODE)
    L_originLength = np.linalg.norm(L_originArray)
    ####
    # 用图像中的IJK轴
    I_AXIS = np.array(direction[:, 0])
    J_AXIS = np.array(direction[:, 1])
    K_AXIS = np.array(direction[:, 2])
    I_AXISLength = np.linalg.norm(I_AXIS)
    J_AXISLength = np.linalg.norm(J_AXIS)
    K_AXISLength = np.linalg.norm(K_AXIS)

    # 计算a和b的点积（相应位置相乘再把结果相加）
    I_AXIS_dot_route_originArray = I_AXIS.dot(L_GETNODE_minus_GETNODE)
    J_AXIS_dot_route_originArray = J_AXIS.dot(L_GETNODE_minus_GETNODE)
    K_AXIS_dot_route_originArray = K_AXIS.dot(L_GETNODE_minus_GETNODE)
    # 使用余弦定理计算phi的弧度值
    route_originThetaRad = np.arccos(I_AXIS_dot_route_originArray / (I_AXISLength * L_originLength))
    route_originTheta = np.rad2deg(route_originThetaRad)
    route_originThetaTRad = np.arccos(J_AXIS_dot_route_originArray / (J_AXISLength * L_originLength))
    route_originThetaT = np.rad2deg(route_originThetaTRad)
    route_originPhiRad = np.arccos(K_AXIS_dot_route_originArray / (K_AXISLength * L_originLength))
    route_originPhi = np.rad2deg(route_originPhiRad)

    if (route_originPhi < 90):
        route_angle_XOY = 90 - route_originPhi  # XOY是水平面
    else:
        route_angle_XOY = 90 - (180 - route_originPhi)
    if (route_originThetaT < 90):
        route_angle_XOZ = 90 - route_originThetaT
    else:
        route_angle_XOZ = 90 - (180 - route_originThetaT)
    if (route_originTheta < 90):
        route_angle_YOZ = 90 - route_originTheta
    else:
        route_angle_YOZ = 90 - (180 - route_originTheta)

    return route_angle_XOY, route_angle_XOZ, route_angle_YOZ
####计算路径ras角度
if __name__ == "__main__":
    ijk_ray_start = [160., 160. , 72.]
    ijk_ray_end = [320., 320.,  72.]
    print("input ijk_ray_start : ", ijk_ray_start)
    print("input ijk_ray_end : ", ijk_ray_end)
    ang = IJKp_to_IJKaxesinRAS_angle(ijk_ray_start, ijk_ray_end, origin=origin, direction=direction, spacing=spacing)
    print("output IJKp_to_IJKaxesinRAS_angle : ", ang)
    length = IJK_to_RAS_length(ijk_ray_start, ijk_ray_end, origin=origin, direction=direction, spacing=spacing)
    print("output IJKp_to_IJKaxesinRAS_angle : ", length)