Dicom deformable image registration

Question

I have dicom files for two CT scans and deformable registration between them. I would like to ask how I can deform one CT scan based on the deformable registration file: so far, I managed to extract the deformation vector field from the dicom file:

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])

I can plot the deformation field on each CT slice as can be seen in the image below. 1

So, the remaining task is how to deform the CT data based on the deformation vectors (dx, dy, dz)

Many thanks

Answer 1

I appreciate this question is almost a year old, but hey - better late than never!

You are actually pretty close to a solution. I think the missing piece of the puzzle is to switch from numpy for image manipulation to SimpleITK. This is an awesome python library available here (and on 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)

I hope this helps you and anyone else who finds this post!