Create GridTransform from 3d motion field

llafitte007 · July 26, 2022, 11:26am

Hi,

I am developing a motion estimation module based on an existing library ; the result of the estimate is 3d vector field of type numpy array (shape = (width, height, depth, 3)) plus other fields (origin, spacing, …).

If I save the motion result from my program (nifti format) and open it as Transformation from slicer it works great: the visualization and the registration of motion are great!

My problem comes in the automatic conversion, in the module, of the estimated motion into a GridTransform. Several utils approach the solution (slicer.util.addVolumeFromArray, slicer.util.loadTransform) but nothing that I can’t use for my case. Would there be an equivalent to slicer.util.addVolumeFromArray for GridTransforms (like slicer.util.addGridTransformFromArray) ?

I tried to manually convert motion field and create GridTransform without success:

# found here https://discourse.slicer.org/t/how-to-convert-3d-numpy-array-to-vtk-and-save-the-vtk-file/22327/3
def numpyToVTK(data):
    data_type = vtk.VTK_FLOAT
    shape = data.shape
    flat_data_array = data.flatten()
    vtk_data = vtk.util.numpy_support.numpy_to_vtk(num_array=flat_data_array, deep=True, array_type=data_type)
    img = vtk.vtkImageData()
    img.GetPointData().SetScalars(vtk_data)
    img.SetDimensions(shape[0], shape[1], shape[2])
    print(img)
    return img

motionNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLGridTransformNode")
motionVtk = numpyToVTK(motionData)
motionNode.GetTransformFromParent().SetDisplacementGridData(motionVtk)
motionNode.GetTransformFromParent().GetDisplacementGrid().SetOrigin(motionOrigin)
motionNode.GetTransformFromParent().GetDisplacementGrid().SetSpacing(motionSpacing)

Can anyone help me ? Thanks

pieper · July 26, 2022, 1:11pm

Here are a couple examples of creating grid transforms in python. Yes, a utility function could be helpful so if you want to try abstracting this code a PR would be great.

github.com

Slicer/Slicer/blob/3328b81211cb2e9ae16a0b49097744171c8c71c0/Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py#L621-L648


      
          # create the grid transform node
          gridTransform = slicer.vtkMRMLGridTransformNode()
          gridTransform.SetName('Acquisition Transform')
          slicer.mrmlScene.AddNode(gridTransform)
          
          
# create a grid transform with one vector at the corner of each slice
          # the transform is in the same space and orientation as the volume node
          gridImage = vtk.vtkImageData()
          gridImage.SetOrigin(*volumeNode.GetOrigin())
          gridImage.SetDimensions(2, 2, slices)
          sourceSpacing = volumeNode.GetSpacing()
          gridImage.SetSpacing(sourceSpacing[0] * columns, sourceSpacing[1] * rows, sourceSpacing[2])
          gridImage.AllocateScalars(vtk.VTK_DOUBLE, 3)
          transform = slicer.vtkOrientedGridTransform()
          directionMatrix = vtk.vtkMatrix4x4()
          volumeNode.GetIJKToRASDirectionMatrix(directionMatrix)
          transform.SetGridDirectionMatrix(directionMatrix)
          transform.SetDisplacementGridData(gridImage)
          gridTransform.SetAndObserveTransformToParent(transform)
          volumeNode.SetAndObserveTransformNodeID(gridTransform.GetID())

This file has been truncated. show original

github.com

Slicer/LandmarkRegistration/blob/master/RegistrationLib/ThinPlatePlugin.py#L83-L118


      
          def onExportGrid(self):
            """Converts the current thin plate transform to a grid"""
            state = self.registrationState()
          
          
  # since the transform is ras-to-ras, we find the extreme points
            # in ras space of the fixed (target) volume and fix the unoriented
            # box around it.  Sample the grid transform at the resolution of
            # the fixed volume, which may be a bit overkill but it should aways
            # work without too much loss.
            rasBounds = [0,]*6
            state.fixed.GetRASBounds(rasBounds)
            from math import floor, ceil
            origin = list(map(int,map(floor,rasBounds[::2])))
            maxes = list(map(int,map(ceil,rasBounds[1::2])))
            boundSize = [m - o for m,o in zip(maxes,origin) ]
            spacing = state.fixed.GetSpacing()
            spacing = [max(spacing)*5]*3
            samples = [ceil(int(b / s)) for b,s in zip(boundSize,spacing)]
            extent = [0,]*6
            extent[::2] = [0,]*3

This file has been truncated. show original

Topic		Replies	Views
4D ITK image to vtkMRMLGridTransformNode Support registration , transforms	10	1678	November 13, 2023
Help on 3D Slicer crash when volume node's image data is invalid Development	6	831	April 3, 2021
Displacement Vector Visualisation Support	15	1313	February 8, 2022
Grid transform Set-/GetIJKToRASMatrix Development transforms	1	968	October 30, 2020
Creating Volume from Numpy Support	19	9342	December 8, 2020

Create GridTransform from 3d motion field

Related topics