4D ITK image to vtkMRMLGridTransformNode

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

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 -http://www.vtk.org/doc/nightly/html/classvtkGridTransform.html

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

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?

GetDisplacementGridData / SetDisplacementGridData

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.

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.

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

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
2 Likes