4D ITK image to vtkMRMLGridTransformNode

Support
,
1

Hi all,

I have a displacement field “field.nii” that I need to convert from RAS to LPS so that I can use it in Slicer. What I’m doing now is:

  1. Read the file with SimpleITK
  2. Modify the data
  3. Save it
  4. Load it with slicer.util.loadTransform()

Is there a way to create a vtkMRMLGridTransformNode and set the vector data or to modify the data of the transform node directly?

My code:

def getRASFieldFromLPSField(self, displacementFieldPath):
    image = sitk.ReadImage(displacementFieldPath)
    arr = sitk.GetArrayFromImage(image)
    arr[..., :2] *= -1  # RAS to LPS

    # Create new image
    newImage = sitk.GetImageFromArray(arr)
    newImage.SetOrigin(image.GetOrigin())
    newImage.SetDirection(image.GetDirection())
    newImage.SetSpacing(image.GetSpacing())

    sitk.WriteImage(newImage, displacementFieldPath)
    transformNode = slicer.util.loadTransform(displacementFieldPath, returnNode=True)[1]
    return transformNode

Best,
Fernando

2

I have a displacement field “field.nii”

Where does that data come from? When Slicer saves transforms to file it is always converted to LPS. When you load a displacement field from nrrd or ITK transform file (hdf) then Slicer automatically converts to LPS. If you process data in CLI modules you typically do everything in LPS.

You can set displacement field in grid transform using SetDisplacementGridData method -VTK: vtkGridTransform Class Reference

Probably the easiest is to load transfrom from .nii and then modify the grid data image in place using numpy.

3

It comes from a registration tool that saves displacement fields in RAS.

What is the method to access the grid data image from the loaded vtkMRMLGridTransformNode?

4

GetDisplacementGridData / SetDisplacementGridData

5

Nowadays the de facto standard in medical image computing research is to store images in LPS coordinate system when they are written to files. That registration tool should be updated accordingly. It may be difficult to update a software to change coordinate systems internally (that’s why Slicer still uses RAS internally), but it should be easy to change data export to write as LPS. Although nrrd, mha, and maybe nifty can store the used coordinate system reference as metadata, as far as I know, ITK only uses that metadata when it reads nrrd files. So, another option is to save as nrrd.

6

I can’t access those methods from Python (Nightly):

>>> _, t = loadTransform('/tmp/transform.nii', returnNode=True)
>>> t

(vtkCommonCorePython.vtkMRMLGridTransformNode)0x12d39d3f8
>>> t.GetDisplacementGridData

Traceback (most recent call last):
  File "<console>", line 1, in <module>
AttributeError: 'vtkCommonCorePython.vtkMRMLGridTransformNode' object has no attribute 'GetDisplacementGridData'

I’ll let the developers know, thanks for the information.

7

You need to get the VTK transform from the transform node: t.GetTransformFromParent().GetDisplacementGrid()

I think the coordinate system situation is particularly bad in .nii files (I don’t know the details, because I avoid .nii as much as possible but see for example this discussion thread). So, it may not be possible to simply use LPS in nifty.

Option A. Everything works very reliably in nrrd, so you could use that file format.

Option B. Try to clean up the coordinate system mess in nifti. It should be started at the lowest level, making sure ITK can read/write nifti files consistently (e.g., make sure the data is converted to LPS upon loading, based on information stored in the header; make sure the used coordinate system is stored in the image header of the written files). Once that is cleaned up, Slicer would probably work correctly.

1 Like
8

I was looking for how to load my own numpy-based transform into a Slicer Grid Transform, and found that this was the closest post on here. After getting something to work, I thought I’d share my example code below.

  • You may have to change your directionmatrix to account for differences in coordinates systems (as discussed in other posts on this topic).
  • The warpfile is a .npy file holding a numpy array with my transform. This had shape (256, 256, 256, 3) for me.
  • All my volumes were 256x256x256, so its hardcoded in, but you should change it according to your needs.

I hope this is helpful for others and will save them the time I spent figuring this out. Also, please feel free to offer corrections if you notice something odd. I’m just beginning to code in Slicer, so there is probably a better way to go about this.

Cheers

import numpy as np
import vtk.util.numpy_support

def loadCustomTransform(warpfile):
    #create a blank Grid Transform to load into
    warpnode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLGridTransformNode')
    transformGrid = warpnode.GetTransformFromParent()
    displacementGrid = transformGrid.GetDisplacementGrid()
    displacementGrid.SetDimensions((256,256,256))
    directionmatrix = np.array([[1, 0, 0, 0],[0, 1, 0, 0],[0, 0, 1, 0],[0, 0, 0, 1]])
    directionVTKmatrix = createVTK4x4matrix(directionmatrix)
    transformGrid.SetGridDirectionMatrix(directionVTKmatrix)

    #load a numpy file holding your deformation vector field.  
    #This had shape (256, 256, 256, 3) for me
    warpdata = np.load(warpfile)
    
    #process your numpy file so that it works for Slicer
    warpdata_reshaped = warpdata.reshape((256*256*256, 3))
    vtkwarpdata = vtk.util.numpy_support.numpy_to_vtk(warpdata_reshaped)

    #Plug in your own deformation vector field and Voila!
    #The warpnode now has your custom data
    displacementGrid.GetPointData().SetScalars(vtkwarpdata)
    transformGrid.SetDisplacementGridData(displacementGrid)
    warpnode.SetAndObserveTransformFromParent(transformGrid)
    
def createVTK4x4matrix(numpyarray):
    #a helper function to change numpy arrays into 4x4 VTK matrices
    vtkmatrix = vtk.vtkMatrix4x4()
    for i in np.arange(4):
        for j in np.arange(4):
            vtkmatrix.SetElement(i,j, numpyarray[i,j])
    return vtkmatrix
4 Likes
9

@emckenzi123 I have the exact same use-case.
An external program is providing a numpy_array (LPS coordinate system) that represents a dense vector field. My goal is to load this dense vector field in Slicer for vizualisation purposes (the corresponding 3D image is already loaded).

Best option for me would be to write a .nrrd file and load it using the simple AddData in the GUI:

import nrrd
import numpy as np

# At first, a variable 'numpy_array' is created (LPS coordinate system). 
# It contains a numpy.ndarray with the following signature:
# > numpy_array.shape
#   (512, 512, 256, 3)
# > numpy_array.dtype
#   dtype('float32')
# > numpy_array.flags
#   C_CONTIGUOUS : True
#   F_CONTIGUOUS : False
#   OWNDATA : True
#   WRITEABLE : True
#   ALIGNED : True
#   WRITEBACKIFCOPY : False

header = {
		'type': numpy_array.dtype.__str__(),
		'sizes': np.array(numpy_array.shape), 
		'dimension': len(numpy_array.shape),
		'endian': 'little',
		'encoding': 'raw',
		'space': 'right-anterior-superior',
		'space directions': [[1, 0., 0.], [0., 1, 0.], [0., 0., 1], None],
		'space origin': np.array([0, 0., 0., 0.]),
		'kinds': ['domain', 'domain', 'domain', 'vector'],
		}
nrrd.write("my_file.nrrd", numpy_array, header)

Slicer is not reading this when I try AddData (GUI) and importing my file as a Transform // Volume (i tried both).
Among the logs I get I can find this “ReadData: This is not a nrrd file”.

Anyone (@lassoan if you could take a look would be awesome) with an idea of what I’m doing wrong is more than welcome.
Thanks :sunny:

10

The file cannot be loaded because the dimension of space origin is incorrect. You have 3 spatial dimensions, so space origin has to be set to [0., 0., 0.].

You don’t need to specify type and size (pynrrd sets these based on the input array).

DICOM and most medical image computing software store files in LPS (left-posterior-superior) coordinate system, so I would recommend to use left-posterior-superior in space field. Slicer converts the vector values and the axis orientations to RAS when loading the image, so keeping RAS could make sense, too, but then you risk incompatibility with other software (that ignore the space field and blindly assume LPS).

Complete working example:

import nrrd
import numpy as np

# Use a random numpy array as displacement field
displacement_field = np.random.rand(50, 40, 30, 3)

header = {
    'endian': 'little',
    'encoding': 'raw',
    'space': 'left-posterior-superior',
    'space directions': [[1, 0., 0.], [0., 1, 0.], [0., 0., 1], None],
    'space origin': np.array([0., 0., 0.]),
    'kinds': ['domain', 'domain', 'domain', 'vector'],
    }
nrrd.write("c:/tmp/my_displacement_field.nrrd", displacement_field, header)
1 Like
11

Thank you for the clarification sir, your example works like a charm.
You made my day.