Hi,
I am playing around with Airlab (a nice toolkit and research environment for autograd/pytorch based deformable registration). With Airlab it is possible to non-linearly register two volumes, and save the computed deformation field as a nrrd file, under the hood it uses sitk.WriteImage() to do that. When loading the deformation field, I realized that I need to change the IJKToRASMatrix from sth like diag([-1,-1,1,1]) to diag([1,1,1,1]), to appear at the “right location”, and I need to invert the x-/y-components of the deformation field (by multiplying with -1).
The only way I am right now able to change the affine matrix of the GridTransform, is to load the grid transform via loadVolume(), manipulate its affine transform via GetIJKToRASMatrix() and SetIJKToRASMatrix() (also, change the vector x/y components via arrayFromVolume() and updateVolumeFromArray()), then write it back to disk via saveNode(), and then load it back into Slicer via loadTransform().
My question now is whether there is a way to avoid handling the grid transform as a volume entirely. E.g. when I load the transform via loadTransform(), the vtkMRMLGridTransformNode class does not have the methods GetIJKToRASMatrix and SetIJKToRASMatrix. Is there an equivalent? Alternatively: Is there a way to conveniently convert the grid transform to a volume and back, without the “detour” of having to read/write a volume first?
Files are always in LPS coordinate systems, while Slicer uses RAS internally. If you transfer the volume and displacement field via files then everything should work correctly (ITK will work in LPS, Slicer works in RAS, Slicer knows that it has to transform displacement field vectors from LPS to RAS on load). If you don’t use compression then file reading/writing time should be negligible compared to registration time. You can run your processing script as a Python scripted CLI module, which naturally does all transfers via files and it also has the advantage that it runs the registration in the background (so that the application GUI is not blocked during execution) and any error or exceeding available memory cannot crash the application (because the script runs in a separate process).
You can load the transform using loadTransform. Both the image geometry (origin, spacing, axis directions) and vectors inside the volume are expected to be in LPS.
Since a transform node can contain arbitrarily complex composition of various transforms, the node does not have direct GetIJKToRASMatrix method. If you know that you have a simple transform in the node then you can get the ToParent or FromParent transform as the specified type and get direction matrix from the returned vtkOrientedGridTransform. To convert a displacement field between LPS<->RAS it is of course not enough to change the image geometry but each vector must be transformed, too (as shown in the C++ code above).