第29课时:脚本

教学大纲


:one: https://slicer.readthedocs.io/en/latest/developer_guide/script_repository.html#segmentations


:two:


:three:


去除背板

数据:示例数据自动引用


# This script automatically removes table from a CT volume

# Generate input data
################################################

import SampleData

# Get input image (in this example, download a sample data set)
sampleDataLogic = SampleData.SampleDataLogic()
masterVolumeNode = sampleDataLogic.downloadCTACardio()
threshold = -44

# Process
################################################

# Create segmentation
segmentationNode = slicer.vtkMRMLSegmentationNode()
slicer.mrmlScene.AddNode(segmentationNode)
segmentationNode.CreateDefaultDisplayNodes() # only needed for display
segmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(masterVolumeNode)

# Create segment editor to get access to effects
slicer.app.processEvents()
segmentEditorWidget = slicer.qMRMLSegmentEditorWidget()
# To show segment editor widget (useful for debugging): segmentEditorWidget.show()
segmentEditorWidget.setMRMLScene(slicer.mrmlScene)
segmentEditorNode = slicer.vtkMRMLSegmentEditorNode()
slicer.mrmlScene.AddNode(segmentEditorNode)
segmentEditorWidget.setMRMLSegmentEditorNode(segmentEditorNode)
segmentEditorWidget.setSegmentationNode(segmentationNode)
segmentEditorWidget.setMasterVolumeNode(masterVolumeNode)

# Check that required extensions are installed

if not segmentEditorWidget.effectByName("Wrap Solidify"):
    slicer.util.errorDisplay("Please install 'SurfaceWrapSolidify' extension using Extension Manager.")

if not segmentEditorWidget.effectByName("Mask volume"):
    slicer.util.errorDisplay("Please install 'SegmentEditorExtraEffects' extension using Extension Manager.")

# Create object of interest segment by thresholding
slicer.app.processEvents()
volumeScalarRange = masterVolumeNode.GetImageData().GetScalarRange()
objectSegmentID = segmentationNode.GetSegmentation().AddEmptySegment()
segmentEditorNode.SetSelectedSegmentID(objectSegmentID)
segmentEditorWidget.setActiveEffectByName("Threshold")
effect = segmentEditorWidget.activeEffect()
effect.setParameter("MinimumThreshold",str(threshold))
effect.setParameter("MaximumThreshold",str(volumeScalarRange[1]))
effect.self().onApply()

# Find largest object, remove all other regions from the segment
slicer.app.processEvents()
segmentEditorWidget.setActiveEffectByName("Islands")
effect = segmentEditorWidget.activeEffect()
effect.setParameterDefault("Operation", "KEEP_LARGEST_ISLAND")
effect.self().onApply()

# Fill holes in the segment to create a solid region of interest
slicer.app.processEvents()
segmentEditorWidget.setActiveEffectByName("Wrap Solidify")
effect = segmentEditorWidget.activeEffect()
effect.setParameter("region", "outerSurface")
effect.setParameter("outputType", "segment")
effect.setParameter("remeshOversampling", 0.3)  # speed up solidification by lowering resolution
effect.self().onApply()

# Blank out the volume outside the object segment
slicer.app.processEvents()
segmentEditorWidget.setActiveEffectByName('Mask volume')
effect = segmentEditorWidget.activeEffect()
effect.setParameter('FillValue', -1000)
effect.setParameter('Operation', 'FILL_OUTSIDE')
effect.self().onApply()

# Remove temporary nodes and widget
segmentEditorWidget = None
slicer.mrmlScene.RemoveNode(segmentEditorNode)
slicer.mrmlScene.RemoveNode(segmentationNode)

# Show masked volume
maskedVolume = slicer.mrmlScene.GetFirstNodeByName(masterVolumeNode.GetName()+" masked")
slicer.util.setSliceViewerLayers(background=maskedVolume)

回到原点的快捷键:Ctrl + e

shortcut = qt.QShortcut(qt.QKeySequence('Ctrl+e'), slicer.util.mainWindow())
shortcut.connect('activated()',
  lambda: slicer.modules.markups.logic().JumpSlicesToLocation(0,0,0, True))

关闭软件

exit()

将两个平面设置为’P’, ‘P_1’,工具栏中平面工具应用分别建立两个平面,之后在代码栏中粘贴下面代码即可显示角度。

planeNodeNames = ['P', 'P_1']

# Print angles between slice nodes
def ShowAngle(unused1=None, unused2=None):
    planeNormalVectors = []
    for planeNodeName in planeNodeNames:
        planeNode = slicer.util.getFirstNodeByClassByName('vtkMRMLMarkupsPlaneNode', planeNodeName)
        planeNormalVector = [0.0, 0.0, 0.0]
        planeNode.GetNormalWorld(planeNormalVector)
        planeNormalVectors.append(planeNormalVector)
    angleRad = vtk.vtkMath.AngleBetweenVectors(planeNormalVectors[0], planeNormalVectors[1])
    angleDeg = vtk.vtkMath.DegreesFromRadians(angleRad)
    print('Angle between planes {0} and {1} = {2:0.3f}'.format(planeNodeNames[0], planeNodeNames[1], angleDeg))

# Observe plane node changes
for planeNodeName in planeNodeNames:
    planeNode = slicer.util.getFirstNodeByClassByName('vtkMRMLMarkupsPlaneNode', planeNodeName)
    planeNode.AddObserver(slicer.vtkMRMLMarkupsPlaneNode.PointModifiedEvent, ShowAngle)

# Print current angle
ShowAngle()

计算表面积(CC)

CC=slicer.util.getNodesByClass("vtkMRMLMarkupsClosedCurveNode")
for curve in CC:
    surfaceAreaMm2 = slicer.modules.markups.logic().GetClosedCurveSurfaceArea(curve)
    print("Curve {0}: surface area = {1:.2f} mm2".format(curve.GetName(), surfaceAreaMm2))