Volume reslice with tracking

VolumeResliceWithTracking.py

# Take a slice out of a volume

import SampleData

# Input parameters

observerTag = None
outputSpacing = [2.0, 2.0, 2.0]  # Millimeters/pixel
outputExtent = [0, 99, 0, 99, 0, 0] # First and last pixel indices along each axis

'''
Extent cannot be negative, because some Slicer components assume that extent starts at zero
outputExtent = [-50, 49, -50, 49, 0, 0]
Reslice and image data content will be correct, but Slicer will not display the negative extent of the slice.
outputNode.ShiftImageDataExtentToZeroStart()
would fix the problem once, but it cannot be used repeatedly with reslice transform modifications.
'''

# Get an input volume

inputVolumeNode = slicer.mrmlScene.GetFirstNodeByName("MRHead")
if inputVolumeNode is None:
  inputVolumeNode = SampleData.SampleDataLogic().downloadMRHead()

# Compute reslice transformation

volumeToIjkMatrix = vtk.vtkMatrix4x4()
inputVolumeNode.GetRASToIJKMatrix(volumeToIjkMatrix)

sliceToRasTransform = vtk.vtkTransform()
sliceToRasTransform.Translate(-100, -100, 0.0)
sliceToRasTransform.RotateX(30)
sliceToRasTransform.RotateY(30)

sliceToIjkTransform = vtk.vtkTransform()
sliceToIjkTransform.Concatenate(volumeToIjkMatrix)
sliceToIjkTransform.Concatenate(sliceToRasTransform)

# Use MRML node to modify transform in GUI

sliceToRasNode = slicer.mrmlScene.GetFirstNodeByName("SliceToRas")
if sliceToRasNode is None:
  sliceToRasNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLinearTransformNode", "SliceToRas")

sliceToRasNode.SetAndObserveTransformToParent(sliceToRasTransform)

# Create illustration for current transform

createModelsLogic = slicer.modules.createmodels.logic()
sliceCoordinateModel = createModelsLogic.CreateCoordinate(30, 2)
sliceCoordinateModel.SetAndObserveTransformNodeID(sliceToRasNode.GetID())
modelDisplay = sliceCoordinateModel.GetDisplayNode()
modelDisplay.SetSliceIntersectionVisibility(True)

# Run reslice to produce output image

reslice = vtk.vtkImageReslice()
reslice.SetInputData(inputVolumeNode.GetImageData())
reslice.SetResliceTransform(sliceToIjkTransform)
reslice.SetInterpolationModeToLinear()
reslice.SetOutputOrigin(0.0, 0.0, 0.0)  # Must keep zero so transform overlays slice with volume
reslice.SetOutputSpacing(outputSpacing)
reslice.SetOutputDimensionality(2)
reslice.SetOutputExtent(outputExtent)
reslice.SetBackgroundLevel(0)
reslice.Update()
reslice.GetOutput().SetSpacing(1,1,1)  # Spacing will be set on MRML node to let Slicer know

# To allow re-run of this script, try to reuse exisiting node before creating a new one

outputNode = slicer.mrmlScene.GetFirstNodeByName("OutputVolume")
if outputNode is None:
  outputNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", "OutputVolume")

outputImageData = reslice.GetOutput()
outputNode.SetAndObserveImageData(reslice.GetOutput())
outputNode.SetSpacing(outputSpacing)
outputNode.CreateDefaultDisplayNodes()

# Transform output image so it is aligned with original volume

outputNode.SetAndObserveTransformNodeID(sliceToRasNode.GetID())

# Show and follow output image in red slice view

redSliceWidget = slicer.app.layoutManager().sliceWidget("Red")
redSliceWidget.sliceController().setSliceVisible(True)
redSliceWidget.sliceLogic().GetSliceCompositeNode().SetBackgroundVolumeID(outputNode.GetID())
driver = slicer.modules.volumereslicedriver.logic()
redView = slicer.util.getNode('vtkMRMLSliceNodeRed')
driver.SetModeForSlice(driver.MODE_TRANSVERSE, redView)
driver.SetDriverForSlice(outputNode.GetID(), redView)

# Callback function for transform updates

def updateReslice(event, caller):
  try:
    slicer.app.pauseRender()
    
    inputTransformId = inputVolumeNode.GetTransformNodeID()
    if inputTransformId is not None:
      inputTransformNode = slicer.mrmlScene.GetNodeByID(inputTransformId)
      rasToVolumeMatrix = vtk.vtkMatrix4x4()
      inputTransformNode.GetMatrixTransformFromWorld(rasToVolumeMatrix)
      sliceToIjkTransform.Identity()
      sliceToIjkTransform.Concatenate(volumeToIjkMatrix)
      sliceToIjkTransform.Concatenate(rasToVolumeMatrix)
      sliceToIjkTransform.Concatenate(sliceToRasTransform)
    
    reslice.Update()
    reslice.GetOutput().SetSpacing(1,1,1)
  finally:
    slicer.app.resumeRender()


if observerTag is not None:
  sliceToRasNode.RemoveObserver(observerTag)
  observerTag = None

observerTag = sliceToRasNode.AddObserver(slicer.vtkMRMLTransformNode.TransformModifiedEvent, updateReslice)