Extracting 2D axial slice from 3D volume using current slice view

Support
1

Hi!

Currently I am working on developing a custom extension.

One of the functionalities that I want to implement is a user can scroll in axial view and select the axial slice as the slice-of-interest. Then I want to use this slice a 2D nrrd/mha file in my algorithm.

I looked at the similar topics and implemented suggested ways. The problem is that my original volume has size (512, 512 ,225) but extracted slice has size (932, 548, 1). I want it to be (512, 512, 1). How can I extract with correct shape?

Below you may find the code I used and images [1].

Thank you!

sliceNodeID = "vtkMRMLSliceNodeRed"

# Get image data from slice view
sliceNode = slicer.mrmlScene.GetNodeByID(sliceNodeID)
appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sliceLayerLogic = sliceLogic.GetBackgroundLayer()
reslice = sliceLayerLogic.GetReslice()
reslicedImage = vtk.vtkImageData()
reslicedImage.DeepCopy(reslice.GetOutput())

# Create new volume node using resliced image
volumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode")
volumeNode.SetIJKToRASMatrix(sliceNode.GetXYToRAS())
volumeNode.SetAndObserveImageData(reslicedImage)
volumeNode.CreateDefaultDisplayNodes()
volumeNode.CreateDefaultStorageNode()
slicer.util.exportNode(volumeNode, <SAVE_PATH>)

Screenshot 2023-08-27 155125
Screenshot 2023-08-27 1551251920×826 61.9 KB

  1. Soler, L., A. Hostettler, V. Agnus, A. Charnoz, J. Fasquel, J. Moreau, A. Osswald, M. Bouhadjar, and J. Marescaux. “3D image reconstruction for comparison of algorithm database: A patient specific anatomical and medical image database.” IRCAD, Strasbourg, France, Tech. Rep (2010)
2

Assuming the data was axial (like this CT) you can just access the volume node as an array and get array[k] where k is the slice offset transformed to IJK space. I f you look in the script repository you can find all the pieces pretty easily.

1 Like
3

Hi Steve,

Thank you very much for your answer!

Previously I got slice offset but could not find a way to obtain the slice itself by using that value. With your hint I got the idea.

I am not sure if it is an elegant solution but I am leaving my code as a reference:

# Get image data from  current axial (red) slice view
sliceNodeID = "vtkMRMLSliceNodeRed"
sliceNode = slicer.mrmlScene.GetNodeByID(sliceNodeID)
appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sourceVolumeNode = sliceLogic.GetBackgroundLayer().GetVolumeNode()

# Get current axial slice in RAS:
sliceOffset = sliceLogic.GetSliceOffset()

# Convert RAS to IJK:
volumeRasToIjk = vtk.vtkMatrix4x4()
sourceVolumeNode.GetRASToIJKMatrix(volumeRasToIjk)
point_VolumeRas = [0, 0, sliceOffset, 1]
point_Ijk = [0, 0, 0, 1]
volumeRasToIjk.MultiplyPoint(point_VolumeRas, point_Ijk)
point_Ijk = [ int(round(c)) for c in point_Ijk[0:3] ]

# Get slice of interest as a volume:
imageArray = slicer.util.arrayFromVolume(sourceVolumeNode)
sliceOfInterest = imageArray[point_Ijk[-1], :, :]
volumeOfInterestNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLScalarVolumeNode')
slicer.util.updateVolumeFromArray(volumeOfInterestNode, sliceOfInterest)

#  Copy image origin, spacing and direction:
ijkToRas = vtk.vtkMatrix4x4()
sourceVolumeNode.GetIJKToRASMatrix(ijkToRas)
volumeOfInterestNode.SetIJKToRASMatrix(ijkToRas)

# Export:
slicer.util.saveNode(volumeOfInterestNode,  r"C:\extracted_slice.nrrd")

Best regards!

2 Likes
4

New edit

Much shorter version if anyone is interested:

sliceNodeID = "vtkMRMLSliceNodeRed"
sliceNode = slicer.mrmlScene.GetNodeByID(sliceNodeID)
appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sliceOffset = sliceLogic.GetSliceOffset()

k = sliceLogic.GetSliceIndexFromOffset(sliceOffset) - 1 # -1 since 1-based