Display displacement field with SimpleITK

Beginner Questions
1

How could I display a displacement field (in physical space) with arrows on an image slice with SimpleITK ?

I tried using quiver in matplotlib, but I had troubles with the physical to voxel space conversions.

Is there any example available somewhere ?

2

Hello @Arthur_Masson,

SimpleITK is not the tool for this, it has minimal visualization related functionality.

I highly recommend saving the transformation to file and using 3D slicer to do the visualization. Slicer has a very nice module for visualization of transformations.

2 Likes
3

Thanks @zivy ! But I am looking for a method to generate those visualization images programmatically ; to generate a report comparing many different displacement fields.

4

Hi @Arthur_Masson,

Got it. For direct visualization, possibly try the mplot3D quiver or look into using vtk.

Another possibility, which is indirect, would be to visualize the Jacobian determinant image of the displacement field (contraction<1, expansion>1, no change in volume=1). The relevant filters are DisplacementFieldJacobianDeterminantFilter and for color, map the result through the ScalarToRGBColormapImageFilter.

1 Like
5

You should be able to script Slicer to generate renderings.

1 Like
6

Thanks zivy and dchen!

I ended up using the following script:

flair1_path = data_path + 'flair.nii.gz'
df_path = data_path + 'df.nii.gz'

image = sitk.ReadImage(flair1_path)
image_data = sitk.GetArrayFromImage(image)
slice_index = image_data.shape[0]//2
image_slice = image_data[slice_index,:,:][::-1,::-1]

spacing = image.GetSpacing()
size = image.GetSize()

# # extent = (left, right, bottom, top)
extent = (0, size[0] * spacing[0], size[1] * spacing[1], 0)

plt.figure(figsize=(20,20))

t=plt.imshow(image_slice,extent=extent,interpolation=None)

t.set_cmap("gray")

df = sitk.ReadImage(df_path)

df_data = sitk.GetArrayFromImage(df)

slice_index = df_data.shape[0]//2

df_slice = df_data[slice_index,:,:,:]

orientation = df.GetDirection()
spacing = df.GetSpacing()

orientation = np.array(orientation).reshape((3,3))
voxel_to_physical = ( spacing * np.identity(3) ) @ orientation
physical_to_voxel = np.linalg.inv(voxel_to_physical)

# for i in range(df_slice.shape[0]):
#     for j in range(df_slice.shape[1]):
#         p = physical_to_voxel @ df_slice[i,j]
#         df_slice[i,j] = p

x = df_slice[:,:,0][::-1,::-1]
y = -df_slice[:,:,1][::-1,::-1]

coordsX = np.arange(0, size[0] * spacing[0], size[0] * spacing[0] / float(image_data.shape[2]) )
coordsY = np.arange(0, size[1] * spacing[1], size[1] * spacing[1] / float(image_data.shape[1]) )

ic(image_slice.shape)
ic(df_slice.shape)
ic(image_data.shape[2])
ic(image_data.shape[1])
ic(x.shape)
ic(coordsX.shape)
ic(coordsY.shape)


ic(size[0] * spacing[0])
ic(size[0] * spacing[0] / float(image_data.shape[2]))

coordsX, coordsY = np.meshgrid(coordsX, coordsY)

M = np.sqrt(x*x+y*y)

qq=plt.quiver(coordsX, coordsY, x, y, M, cmap=plt.cm.jet, units='x', scale=1)

plt.axis('off')
plt.show()

I was confused by the physical to voxel space conversion ; and the numpy to simpleITK / matlab ordering conversion.

I am surprised that, when generating the displacement field, demons.SetUseImageSpacing(True) does not seem to make any difference (the generated field is the same with demons.SetUseImageSpacing(True) and demons.SetUseImageSpacing(False)), and that I got better results without converting to voxel space.

7

I asked a new question about the FastSymmetricForcesDemonsRegistrationFilter.SetUseImageSpacing() problem.