I created a BSpline-transformation with the module General Registration (BRAINS) with a CT-image as fixed volume and an MR-image as moving volume. Visualization of the BSpline coefficients and conversion to a displacement field is possible with transforms module.
For further investigation on the occurring displacements I would like to get the 3D distortion matrix with displacement vectors for every voxel in mm. Is there an access to the matrix in 3D-Slicer (e.g. to localize and to get the value and direction of the maximum distortion-vector in every layer)? I already tried to export the data to Matlab and OriginPro, but I only get a table with one column.
Furthermore I would like to take a look at the algorithm. In Elastix module there is a way to display the log files during registration. Is it possible to get a detailed log in General Registration (BRAINS) module? I tried to generate a log file report under debugging parameters, but there wasn’t any data in my CSV-file.
It sounds like you are almost there. In the Transforms module there’s a Convert section where you can create a grid transform within a selected reference grid (defined by a volume). The grid transform is the displacement vectors in mm and you can use arrayFromGridTransform to access them for analysis.
For the BRAINSFit logging information you can add some command line arguments to get extra debug info. You can get the starting command line from the Slicer log and then run it extra debug parameters. Best to look at the BRAINSFit source code for to interpret this.
I managed to get the array from the output displacement field (stored as transform node) with “arrayFromGridTransform”. But I have difficulties in interpreting the matrix. I guess the 3 columns stand for the x-, y-, and z-component of each displacement-vector in LPS coordinate system. Which of the values represent the location of the base of each arrow? Aside from that I don’t know how to get the complete Matrix without “…,” as replacement.
Yes, these are the displacement vectors, but they are in RAS not LPS. The “…” is just to shorten the string for printing; all the values are there in the array. The start point of each arrow is defined by the position of the vector in the 3D volumetric grid. The locations of these start points are calculated based on the origin, spacing, and directions like in an image volume, but the API is a little different. The tip of the arrow would be the base plus the displacement vector.
For another project I created these helper functions that I plan to add to slicer.util at some point. You may find them useful for reference to see the relationship between grid transforms and volumes. I believe they are correct but maybe you can test them and let me know if they work for you. addVolumeFromGridTransform is helpful for debugging because you can get a color image (VectorVolume) that shows the space covered by the grid transform and you can mouse over the image to see the displacement values in the DataProbe.
def addGridTransformFromArray(narray, referenceVolume=None, name=None):
"""Create a new grid transform node from content of a numpy array and add it to the scene.
Displacement values are deep-copied, therefore if the numpy array
is modified after calling this method, voxel values in the volume node will not change.
:param narray: numpy array containing grid transform vectors (shape should be [Nk, Nj, Ni, 3], i.e. one displacement vector for slice, row, column location).
:param referenceVolume: a vtkMRMLVolumeNode or subclass to define the directions, origin, and spacing
:param name: grid transform node name
:return: created new grid transform node
Example::
# create an identity grid transform
import numpy
gridTransformNode = slicer.util.addGridTransformFromArray(numpy.zeros((30, 40, 50, 3)))
"""
import slicer
from vtk import vtkMatrix4x4
from vtk.util import numpy_support
gridTransformNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLGridTransformNode")
if name is not None:
gridTransformNode.SetName(name)
gridTransform = gridTransformNode.GetTransformFromParent()
gridImage = vtk.vtkImageData()
gridImage.SetDimensions(tuple(reversed(narray.shape[:3])))
gridType = numpy_support.get_vtk_array_type(narray.dtype)
gridImage.AllocateScalars(gridType, 3)
if referenceVolume is not None:
gridDirectionMatrix = vtk.vtkMatrix4x4()
referenceVolume.GetIJKToRASDirectionMatrix(gridDirectionMatrix)
gridTransform.SetGridDirectionMatrix(gridDirectionMatrix)
gridImage.SetOrigin(referenceVolume.GetOrigin())
gridImage.SetSpacing(referenceVolume.GetSpacing())
gridTransform.SetDisplacementGridData(gridImage)
transformArray = slicer.util.arrayFromGridTransform(gridTransformNode)
transformArray[:] = narray
slicer.util.arrayFromGridTransformModified(gridTransformNode)
return gridTransformNode
def addVolumeFromGridTransform(gridTransformNode, name=None):
"""Create a new vector volume from grid transform node from content.
:param gridTransformNode: source transform
:param name: created volume node name
:return: created new volume
"""
displacements = arrayFromGridTransform(gridTransformNode)
gridTransform = gridTransformNode.GetTransformFromParent()
gridDirectionMatrix = gridTransform.GetGridDirectionMatrix()
displacementGrid = gridTransform.GetDisplacementGrid()
scratchVolume = slicer.vtkMRMLScalarVolumeNode()
scratchVolume.SetIJKToRASDirectionMatrix(gridDirectionMatrix)
scratchVolume.SetSpacing(displacementGrid.GetSpacing())
scratchVolume.SetOrigin(displacementGrid.GetOrigin())
ijkToRAS = vtk.vtkMatrix4x4()
scratchVolume.GetIJKToRASMatrix(ijkToRAS)
return addVolumeFromArray(displacements, ijkToRAS=ijkToRAS, name=name, nodeClassName="vtkMRMLVectorVolumeNode")
