Contains python functions which should work in the 3D Slicer envrionment to convert between Slicer image registration transformation matrices and FSL/FLIRT image registration transform matrices

Slicer_FSL-FLIRT_Transforms_Conversion.py

import numpy as np
import slicer, vtk


def getVolumeCoordinateTransformMatrixDict(volNode):
    import numpy as np

    def translationMatrix(point):
        trMat = np.eye(4)
        trMat[0:3, 3] = (*point,)
        return trMat

    vtk_IJKDirs = vtk.vtkMatrix4x4()
    volNode.GetIJKToRASDirectionMatrix(vtk_IJKDirs)
    IJKDirs = slicer.util.arrayFromVTKMatrix(vtk_IJKDirs)
    detIJKDirs = np.linalg.det(IJKDirs)
    # Coordinate systems are ijk, hjk, hjkMm, hras0, ras
    # Scaling by voxel size (this is independent of IJKDirs determinant)
    V_hjk_TO_hjkMm = np.diag([*volNode.GetSpacing(), 1])
    if detIJKDirs < 0:
        ### h-axis remains identical to i-axis
        V_ijk_TO_hjk = np.eye(4)
        ### Rotate scaled array coordinate axes to orient them in RAS space
        V_hjkMm_TO_hras0 = IJKDirs
        ### Properly locate the image origin in RAS space
        V_hras0_TO_ras = translationMatrix(volNode.GetOrigin())
    else:
        ### Flip and shift i-axis to become h-axis
        # The h-axis is obtained by flipping the i-axis and moving the origin to
        # opposite the opposite corner of the volume (still along the i/h axis)
        numVoxelsInIDimension = volNode.GetImageData().GetDimensions()[0]
        V_ijk_TO_hjk = translationMatrix([numVoxelsInIDimension - 1, 0, 0]) @ np.diag(
            [-1, 1, 1, 1]
        )
        # The h coordinate is N - 1 - i, where N is the number of voxels in the i
        # dimension of the image volume array.  Voxel indices are numbered 0 to
        # N-1 for both i and h, but in reversed order from each other.
        ### Rotate scaled array coordinate axes to orient them in RAS space
        # The h-axis is flipped relative to i-axis, so that needs to be reflected here
        HJKDirs = IJKDirs @ np.diag([-1, 1, 1, 1])
        V_hjkMm_TO_hras0 = HJKDirs
        ### Properly locate the image origin in RAS space
        original_origin_destination = np.array(volNode.GetOrigin())
        voxel_spacing_h = volNode.GetSpacing()[0]  # h spacing = i spacing
        # The origin shift needs to be constructed so that the original origin ends up
        # at the originally intended location in RAS space. In the V_hras0
        # coordinate system, the original origin has been moved by the HJKDirs
        # matrix, so it is
        original_origin_location_in_hras0 = (
            V_hjkMm_TO_hras0 @ V_hjk_TO_hjkMm @ V_ijk_TO_hjk @ np.array([0, 0, 0, 1])
        )
        needed_origin_shift_hras0_TO_ras = (
            original_origin_destination - original_origin_location_in_hras0[:3]
        )
        # original_origin_destination - HJKDirs[:3, :3] @ np.array([voxel_spacing_h * (numVoxelsInIDimension - 1), 0, 0])
        V_hras0_TO_ras = translationMatrix(needed_origin_shift_hras0_TO_ras)
    # Add some grouped ones needed for easy FLIRT matrix calculations
    V_hjkMm_TO_ras = V_hras0_TO_ras @ V_hjkMm_TO_hras0
    # Store calculated transform matrices into output dict
    transformMatrixDict = {
        "ijk_to_hjk": V_ijk_TO_hjk,
        "hjk_to_hjkMm": V_hjk_TO_hjkMm,
        "hjkMm_to_hras0": V_hjkMm_TO_hras0,
        "hras0_to_ras": V_hras0_TO_ras,
        "hjk_to_ijk": np.linalg.inv(V_ijk_TO_hjk),
        "hjkMm_to_hjk": np.linalg.inv(V_hjk_TO_hjkMm),
        "hras0_to_hjkMm": np.linalg.inv(V_hjkMm_TO_hras0),
        "ras_to_hras0": np.linalg.inv(V_hras0_TO_ras),
        "ras_to_hjkMm": np.linalg.inv(V_hjkMm_TO_ras),
        "hjkMm_to_ras": V_hjkMm_TO_ras,
    }

    return transformMatrixDict


def getSlicerTransformNodeFromFlirtTransformMatrix_MRMLNodeBased(
    flirtTransformMatrixNumpyArray,
    fixedVolNode,
    movingVolNode,
    outputSlicerTransformNode=None,
):
    """This method creates a Slicer vtkMRMLLinearTransformNode which
    corresponds to an FSL/FLIRT transformation matrix which registers
    two image volumes. The Slicer transformation matrix is not the same
    as the FLIRT transformation matrix because the coordinate systems
    which are being related are not the same: FLIRT relates scaled (and
    possibly reflected) voxel coordinate systems (which I call _hjkMm),
    while Slicer relates RAS spatial coordinate systems (which are
    after the application of IJKToRASMatrix scaling, rotation, and
    origin translations).

    The inputs are the FLIRT registration transform matrix as a 4x4
    numpy array (np.loadtxt(flirtMatFileName) works for that), the
    MRML scalar volume node for the fixed image (fixedVolNode), the
    MRML scalar volume node for the moving image (movingVolNode), and
    optionally, a MRML linear transform node to hold the calculated
    output Slicer transform matrix.  If an output transform node
    is not supplied, a new one is created and returned.

    An example usage could look like
    flirtTransform = np.loadtxt('highres2standard.mat')
    fixed = slicer.util.loadVolume('standard.nii.gz')
    moving = slicer.util.loadvolume('highres.nii.gz')
    slicerTransformNode = getSlicerTransformNodeFromFlirtTransformMatrix_MRMLNodeBased(flirtTransform, fixed, moving)
    moving.SetTransformNodeID(slicerTransformNode.GetID())
    """
    F_trDict = getVolumeCoordinateTransformMatrixDict(fixedVolNode)
    M_trDict = getVolumeCoordinateTransformMatrixDict(movingVolNode)
    # Slicer_Matrix = F_hjkMm_TO_F_ras @ FSL_Matrix @ M_ras_TO_M_hjkMm
    slicerTransformMatrix = (
        F_trDict["hjkMm_to_ras"]
        @ flirtTransformMatrixNumpyArray
        @ M_trDict["ras_to_hjkMm"]
    )
    # Create output transform node if not supplied
    if outputSlicerTransformNode is None:
        ouputSlicerTransformNode = slicer.mrmlScene.AddNewNodeByClass(
            "vtkMRMLLinearTransformNode", "SlicerTransformFromFlirtTransform"
        )
    outputSlicerTransformNode.SetAndObserveMatrixTransformToParent(
        slicer.util.vtkMatrixFromArray(slicerTransformMatrix)
    )
    return ouputSlicerTransformNode


def getFlirtTransformMatrixFromSlicerTransformNode_MRMLNodeBased(
    slicerTransformNode, fixedVolNode, movingVolNode
):
    """This method calculates a FLIRT transform matrix (as 4x4 numpy array)
    which corresponds to a Slicer transform matrix which registers two image
    volumes. The Slicer transformation matrix is not the same
    as the FLIRT transformation matrix because the coordinate systems
    which are being related are not the same: FLIRT relates scaled (and
    possibly reflected) voxel coordinate systems (which I call _hjkMm),
    while Slicer relates RAS spatial coordinate systems (which are
    after the application of IJKToRASMatrix scaling, rotation, and
    origin translations).
    The inputs are the Slicer MRML nodes corresponding to the registration
    (slicerTransformNode), the fixed image (fixedVolNode), and moving
    image (movingVolNode), and the output is a 4x4 numpy array which
    is the affine transform matrix suitable to save in a .mat file
    (np.savetxt('FLIRT_Matrix.mat', FLIRT_Matrix)). An example usage could be

    fixed = slicer.util.loadVolume('standard.nii.gz')
    moving = slicer.util.loadVolume('highres.nii.gz')

    # Perform any of various Slicer registration methods, and put the result
    # into slicerTransformNode
    slicerTransformNode = slicer.util.getNode(...) # fill in node name

    flirt_transform = getFlirtTransformMatrixFromSlicerTransformNode_MRMLNodeBased(slicerTransformNode, fixed, moving)
    np.savetxt('FLIRT_matrix_from_Slicer.mat', flirt_transform)

    # THEN, in an FSL environment:

    flirt -in standard.nii.gz -ref highres.nii.gz -init FLIRT_matrix_from_Slicer.mat -applyxfm -out standard2highres_image.nii.gz

    # This command will apply the registration as calculated in Slicer using the flirt command, and standard2highres_image.nii.gz
    # will contain the highres image registered and resampled to the voxel grid of the standard.nii.gz image.
    """
    F_trDict = getVolumeCoordinateTransformMatrixDict(fixedVolNode)
    M_trDict = getVolumeCoordinateTransformMatrixDict(movingVolNode)
    vtk_slicerTransformMatrix = vtk.vtkMatrix4x4()
    slicerTransformNode.GetMatrixTransformToParent(vtk_slicerTransformMatrix)
    slicerTransformMatrix = slicer.util.arrayFromVTKMatrix(vtk_slicerTransformMatrix)
    # FLIRT_Matrix = F_ras_TO_FhjkMm @ M_ras_TO_F_ras @ M_hjkMm_TO_M_ras
    FLIRT_Matrix = (
        F_trDict["ras_to_hjkMm"] @ slicerTransformMatrix @ M_trDict["hjkMm_to_ras"]
    )
    return FLIRT_Matrix