Slice rotation in viewport when using SetSliceToRASByNTP in python

hey guys,
I have been trying to debug this issue to no avail and I’m hoping someone could help me out or point me in some direction. Please refer to my stackoverflow post, feel free to answer here or there. If solution or discussion occurs here I can post it in stackoverflow.

Thanks!

Determining a smoothly changing curve normal direction is not a trivial task. In recent Slicer versions, we introduced curve node, which can provide curve tangent and normal using Frenet-Serret frame.

Here is a complete example of computing a panoramic X-ray from a cone-beam dental CT. Just a demonstration of how an image can be resampled along a curve - the code is not optimized for performance or quality and distance along curve is not scaled (we simply used all point indices instead of retrieving point indices based on desired distance along curve).

# Get a dental CT scan
import SampleData
volumeNode = SampleData.SampleDataLogic().downloadDentalSurgery()[1]

# Define curve
curveNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLMarkupsCurveNode')
curveNode.CreateDefaultDisplayNodes()
curveNode.GetCurveGenerator().SetNumberOfPointsPerInterpolatingSegment(25) # add more curve points between control points than the default 10
curveNode.AddControlPoint(vtk.vtkVector3d(-45.85526315789473,	-104.59210526315789,	74.67105263157896))
curveNode.AddControlPoint(vtk.vtkVector3d(-50.9078947368421,	-90.06578947368418,	66.4605263157895))
curveNode.AddControlPoint(vtk.vtkVector3d(-62.27631578947368,	-78.06578947368419,	60.7763157894737))
curveNode.AddControlPoint(vtk.vtkVector3d(-71.86705891666716,	-58.04403581456746,	57.84679891116521))
curveNode.AddControlPoint(vtk.vtkVector3d(-74.73084356325877,	-48.67611043794342,	57.00664267528636))
curveNode.AddControlPoint(vtk.vtkVector3d(-88.17105263157895,	-35.75,	55.092105263157904))
curveNode.AddControlPoint(vtk.vtkVector3d(-99.53947368421052,	-35.75,	55.092105263157904))
curveNode.AddControlPoint(vtk.vtkVector3d(-107.75,	-43.96052631578948,	55.092105263157904))
curveNode.AddControlPoint(vtk.vtkVector3d(-112.80263157894736,	-59.118421052631575,	56.355263157894754))
curveNode.AddControlPoint(vtk.vtkVector3d(-115.32894736842104,	-73.01315789473684,	60.144736842105274))
curveNode.AddControlPoint(vtk.vtkVector3d(-125.43421052631578,	-83.74999999999999,	60.7763157894737))
curveNode.AddControlPoint(vtk.vtkVector3d(-132.3815789473684,	-91.96052631578947,	63.934210526315795))
curveNode.AddControlPoint(vtk.vtkVector3d(-137.43421052631578,	-103.96052631578947,	67.72368421052633))

sliceNode = slicer.mrmlScene.GetNodeByID('vtkMRMLSliceNodeRed')
rotationAngleDeg = 180
sliceNode.SetFieldOfView(40,40,1) # zoom in

appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sliceLayerLogic = sliceLogic.GetBackgroundLayer()
reslice = sliceLayerLogic.GetReslice()
reslicedImage = vtk.vtkImageData()

# Straightened volume (useful for example for visualization of curved vessels)
straightenedVolume = slicer.modules.volumes.logic().CloneVolume(volumeNode, 'straightened')

# Capture a number of slices orthogonal to the curve and append them into a volume.
# sliceToWorldTransform = curvePointToWorldTransform * RotateZ(rotationAngleDeg)
curvePointToWorldTransform = vtk.vtkTransform()
sliceToWorldTransform = vtk.vtkTransform()
sliceToWorldTransform.Concatenate(curvePointToWorldTransform)
sliceToWorldTransform.RotateZ(rotationAngleDeg)
sliceNode.SetXYZOrigin(0,0,0)
numberOfPoints = curveNode.GetCurvePointsWorld().GetNumberOfPoints()
append = vtk.vtkImageAppend()

for pointIndex in range(numberOfPoints):
    print(pointIndex)
    curvePointToWorldMatrix = vtk.vtkMatrix4x4()
    curveNode.GetCurvePointToWorldTransformAtPointIndex(pointIndex, curvePointToWorldMatrix)
    curvePointToWorldTransform.SetMatrix(curvePointToWorldMatrix)
    sliceToWorldTransform.Update()
    sliceNode.GetSliceToRAS().DeepCopy(sliceToWorldTransform.GetMatrix())
    sliceNode.UpdateMatrices()
    slicer.app.processEvents()
    tempSlice = vtk.vtkImageData()
    tempSlice.DeepCopy(reslice.GetOutput())
    append.AddInputData(tempSlice)

append.SetAppendAxis(2)
append.Update()
straightenedVolume.SetAndObserveImageData(append.GetOutput())

# Create panoramic volume by mean intensity projection along an axis of the straightened volume
import numpy as np
panoramicVolume = slicer.modules.volumes.logic().CloneVolume(straightenedVolume, 'panoramic')
panoramicImageData = panoramicVolume.GetImageData()
straightenedImageData = straightenedVolume.GetImageData()
panoramicImageData.SetDimensions(straightenedImageData.GetDimensions()[2], straightenedImageData.GetDimensions()[1], 1)
panoramicImageData.AllocateScalars(straightenedImageData.GetScalarType(), straightenedImageData.GetNumberOfScalarComponents())
panoramicVolumeArray = slicer.util.arrayFromVolume(panoramicVolume)
straightenedVolumeArray = slicer.util.arrayFromVolume(straightenedVolume)
panoramicVolumeArray[0, :, :] = np.flip(straightenedVolumeArray.mean(2).T)
slicer.util.arrayFromVolumeModified(panoramicVolume)
panoramicVolume.SetSpacing(4.0, 0.5, 0.5) # just approximate spacing (would need to properly compute from FOV and image size)
sliceNode.SetOrientationToAxial()
slicer.util.setSliceViewerLayers(background=panoramicVolume, fit=True)

Computed panoramic X-ray:

1 Like

Hey,
Really cool algorithm. I am trying to implement it but I am running into this error

Traceback (most recent call last):
File “C:/Users/talmazovg/SlicerDev/SlicerPano/SlicerPano/SlicerPano/SlicerPano.py”, line 234, in onCreatePathButtonClicked
test.construct_pano(self.path)
File “C:/Users/talmazovg/SlicerDev/SlicerPano/SlicerPano/SlicerPano/SlicerPano.py”, line 604, in construct_pano
curveNode.CreateDefaultDisplayNodes()
AttributeError: ‘NoneType’ object has no attribute ‘CreateDefaultDisplayNodes’

invoking curveNode = slicer.mrmlScene.AddNewNodeByClass(‘vtkMRMLMarkupsCurveNode’) return NoneType object.

I am new to slicer programming so I apologize, am I not importing a class or something?
Using Slicer 4.10.2
Thanks!

“vtkMRMLMarkupsCurveNode” is only available in the Slicer nightly (4.11) which is why it is not working for you in Slicer 4.10.2

2 Likes

Based on my original question, how is it that the 3D world view retains proper orientation of the plane while the node view is flipped? Please refer to screenshot.
Is there some way I can anchor the plane to be parallel in the axial direction?

You can arbitrarily rotate the slice view content in the image plane or flip it along any axis, the physical position of the slice will still remain the same.

In your script, you computed the normal direction of the curve inconsistently along the curve (you flipped the normal direction along the curve), that’s why you see a flip in the slice view. If the curve is quasi planar then you can use the curve normal computation method implemented that is implemented in endoscopy module as is (that uses the normal of the plane fit to the curve; a reslicing transform is written to the scene that you can use in Volume reslice driver - see this video for step-by-step instructions). For arbitrary curve shapes, you can use Frenet-Serret frame (that is implemented in recent Slicer preview versions in curves, see my script above for a complete working example).