Dicom可变形图像配准

Question

我有两次CT扫描的dicom文件和它们之间的可变形注册。我想问一下，如何根据可变形的注册文件进行一次CT扫描:到目前为止，我设法从dicom文件中提取了变形矢量场：

import pydicom, numpy as np 
from struct import unpack
ds = pydicom.read_file(fn) # fn: dicom deformation file.
v = ds.DeformableRegistrationSequence[0].DeformableRegistrationGridSequence[0].VectorGridData
values = unpack(f"<{len(v) // 4}f", v)
dim = B.GridDimensions
nx, ny, nz  = (dim[0],dim[1],dim[2]) # number of voxels in x,y,z direction 
# exctract the dx,dy,dz as per https://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.20.3.html#sect_C.20.3.1.3
dx3d = values[0::3]
dy3d = values[1::3]
dz3d = values[2::3]
dx = np.reshape(dx3d,[nz,ny,nx])
dy = np.reshape(dy3d,[nz,ny,nx])
dz = np.reshape(dz3d,[nz,ny,nx])

我可以在每个CT切片上绘制变形场，如下图所示。1

因此，剩下的任务是如何基于变形向量(dx，dy，dz)对CT数据进行变形。

非常感谢

Answer 1

我很感激这个问题差不多有一年了，但是嘿--迟到总比没有好！

你实际上离解决方案很近。我认为拼图中缺少的部分是从numpy切换到SimpleITK图像处理。这是一个非常棒的python库，可用的这里 (以及pip上的)。

import pydicom
import SimpleITK as sitk

from struct import unpack

# read DICOM deformable registration file
ds = pydicom.read_file(fn)

# extract the deformation values
dvf_vals_raw = ds.DeformableRegistrationSequence[0].DeformableRegistrationGridSequence[0].VectorGridData

# convert to numbers
dvf_vals = unpack(f"<{len(dvf_vals_raw) // 4}f", dvf_vals_raw)

# we need more image info to properly define the transform
reg = ds["DeformableRegistrationSequence"][0]["DeformableRegistrationGridSequence"][0]

dvf_grid_size = reg["GridDimensions"].value
dvf_grid_res = reg["GridResolution"].value
dvf_grid_origin = reg["ImagePositionPatient"].value

dvf_arr = np.reshape(dvf_vals, dvf_grid_size[::-1]+[3,])
dvf_img = sitk.GetImageFromArray(dvf_arr)
dvf_img.SetSpacing(dvf_grid_res)
dvf_img.SetOrigin(dvf_grid_origin)

# convert this to a transform
tfm = sitk.DisplacementFieldTransform(dvf_img)

# now we can finally apply the transform
# NB: image can be list of str for DICOM filenames, or NIfTI, etc.
img = sitk.ReadImage(image)

# NB: 2 -> linear interpolation, -1000 -> default (background) value
img_def = sitk.Transform(img, tfm, 2, -1000)

我希望这能帮助你和其他找到这篇文章的人！