如何为SimpleITK图像各切片正确设置Patient Position（Origin）标签？

数据背景 #

我有一个尺寸为(256,256,14)的多切片体数据存储在.vtk文件中，文件头信息如下：

# vtk DataFile Version 3.0    
vtk output
BINARY
DATASET STRUCTURED_POINTS
DIMENSIONS 256 256 14
SPACING 1.25 1.25 8
ORIGIN 111.81 -205.959 104.612
POINT_DATA 917504
SCALARS scalars float
LOOKUP_TABLE default

其中ORIGIN对应切片0的DICOM标签(0020,0032)（Patient Position）。所有14个切片的DICOM标签(0020,0037)（Patient Orientation）取值相同，但每个切片的Patient Position值各不相同。

重现代码 #

我将该.vtk文件读取为SimpleITK图像，调用TransformPhysicalPointToIndex()将世界坐标转换为图像坐标：

import SimpleITK as sitk
import numpy as np

path = "multi_slice_volume.vtk"

reader = sitk.ImageFileReader()
reader.SetFileName(path)
multi_slice_vlm = reader.Execute()

world_coord = np.array([[72.424, -171.130, 125.821]])
image_coord = multi_slice_vlm.TransformPhysicalPointToIndex(world_coord)

已知该世界坐标对应的正确图像坐标应为(5,5,7)（属于切片7），但实际得到的结果是(13, 26, 4)。

问题原因 #

经排查，坐标转换错误的原因有两点：

• a) SimpleITK默认使用了Patient Orientation值(1,0,0,0,1,0,0,0,1)，与数据实际的Orientation不匹配
• b) 转换时仅使用了切片0的Patient Position值，未对应各切片的实际位置

已解决的部分 #

我已经通过以下代码修正了问题a：

multi_slice_vlm.SetDirection((0.43, -0.58, 0.0, 0.85, 0.014, 0.0, -0.29, -0.81, 1.0))

未解决的问题 #

针对问题b，我尝试为每个切片单独设置对应的Patient Position作为Origin：

multi_slice_vlm[:, :, 0].SetOrigin([111.809579849243, -205.95926731824, 104.612316131591]) # slice 0
multi_slice_vlm[:, :, 1].SetOrigin( [106.316583633422, -201.75588673353, 108.632091522216]) # slice 1
    (...)
multi_slice_vlm[:, :, 13].SetOrigin( [100.823587417602, -197.55250614881, 112.651866912841]) # slice 13

但最终得到的图像坐标仍然不正确。

内容的提问来源于stack exchange，提问作者lbrandao