Last active
July 4, 2021 13:06
-
-
Save fepegar/b5ea2354fd0160ca0ad150f3d149a575 to your computer and use it in GitHub Desktop.
3D Slicer code to create an ACPC transform
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment