3D Slicer code to create an ACPC transform

acpc.py

import numpy as np
import SampleData


def getSampleVolume():
    # Download MRHead
    sampleDataLogic = SampleData.SampleDataLogic()
    volumeNode = sampleDataLogic.downloadMRHead()
    return volumeNode


def getMatrixToACPC(ac, pc, ih):
    # Anteroposterior axis
    pcAc = ac - pc
    yAxis = pcAc / np.linalg.norm(pcAc)

    # Lateral axis
    acIhDir = ih - ac
    xAxis = np.cross(yAxis, acIhDir)
    xAxis /= np.linalg.norm(xAxis)

    # Rostrocaudal axis
    zAxis = np.cross(xAxis, yAxis)

    # Rotation matrix
    rotation = np.vstack([xAxis, yAxis, zAxis])

    # AC in rotated space
    translation = -np.dot(rotation, ac)

    # Build homogeneous matrix
    matrix = np.eye(4)
    matrix[:3, :3] = rotation
    matrix[:3, 3] = translation

    return matrix


def getTransformNodeFromNumpyMatrix(matrix, name=None):
    # Create VTK matrix object
    vtkMatrix = vtk.vtkMatrix4x4()
    for row in range(4):
        for col in range(4):
            vtkMatrix.SetElement(row, col, matrix[row, col])

    # Create MRML transform node
    transformNode = slicer.mrmlScene.AddNewNodeByClass(
        'vtkMRMLLinearTransformNode')
    if name is not None:
        transformNode.SetName(name)
    transformNode.SetAndObserveMatrixTransformToParent(vtkMatrix)
    return transformNode


# I defined these on MRHead
ac = np.array([-0.0641399910762672, 17.61291529545006, 5.009494772024041])
pc = np.array([-0.5843405105866637, -10.4779127581112, 3.044020167647495])
ih = np.array([-1.1045410300970602, 1.746799450383008, 45.04402016764749])

volumeNode = getSampleVolume()
matrix = getMatrixToACPC(ac, pc, ih)
transformNode = getTransformNodeFromNumpyMatrix(matrix, name='World to ACPC')

# Create markups node with AC, PC and an interhemispheric point
markupsNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLMarkupsFiducialNode')
markupsNode.SetName('AC-PC-IH')
markupsNode.AddFiducialFromArray(ac, 'AC')
markupsNode.AddFiducialFromArray(pc, 'PC')
markupsNode.AddFiducialFromArray(ih, 'IH')

# Apply transform to volume node and markups node
volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())
markupsNode.SetAndObserveTransformNodeID(transformNode.GetID())

# Fit image to slices
applicationLogic = slicer.app.applicationLogic()
applicationLogic.FitSliceToAll()

# Center views on AC
markupsLogic = slicer.modules.markups.logic()
acIndex = 0
markupsLogic.JumpSlicesToNthPointInMarkup(markupsNode.GetID(), acIndex)