Skip to content

Instantly share code, notes, and snippets.

@lassoan
Created October 22, 2019 13:05
Show Gist options
  • Star 5 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save lassoan/b445c734f118a5fb7643f3fb05f98b07 to your computer and use it in GitHub Desktop.
Save lassoan/b445c734f118a5fb7643f3fb05f98b07 to your computer and use it in GitHub Desktop.
Computing a panoramic X-ray from a cone-beam dental CT. 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)
@lassoan
Copy link
Author

lassoan commented Feb 9, 2023

That's fine, you can execute path/to/Slicer --python-script path/to/someprocessing.py from flask script. You can also have a Slicer instance running and use its findings via the Slicer REST API (then you don't need to wait for Slicer application to start, but you need to wait or start another Slicer instance if there are multiple requests at the range time).

Note that the having the CBCT images are not sufficient for computing panoramic X-ray. You also need the curve among you want to reformat (curve connecting the center of the teeth crowns). There are many AI models for teeth segmentation that you could use to get the curve but I'm not sure if there are any good quality ones available for free.

@AUHS2001
Copy link

AUHS2001 commented Feb 9, 2023

Thanks Again.

But I only want to use my dataset not the SampleData you have used. How can I proceed in this ?

I have 288 dicom images of cbct series, I have generated axial, saggital, coronal and also 3d rendering view from these images, the last thing I want to generate is the panoramic view. I am unable to generate it.

@lassoan
Copy link
Author

lassoan commented Feb 9, 2023

Let's continue this at the Slicer forum - https://discourse.slicer.org - where more people can contribute to and learn from the discussion.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment