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)

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