Documentation/4.6/ScriptRepository

From Slicer Wiki
Jump to: navigation, search
Home < Documentation < 4.6 < ScriptRepository


For the latest Slicer documentation, visit the read-the-docs.



Community-contributed modules

Usage: save the .py file to a directory, add the directory to the additional module paths in the Slicer application settings (choose in the menu: Edit / Application settings, click Modules, click >> next to Additional module paths, click Add, and choose the .py file's location).

Filters

DICOM

Informatics

  • MarkupsInfo.py: Compute the total length between all the points of a markup list.
  • LineProfile.py: Compute intensity profile in a volume along a line.

Community-contributed examples

Usage: Copy-paste the shown code lines or linked .py file contents into Python console in Slicer.

Capture

  • Get a MRML node in the scene based on the node name and call methods of that object. For the MRHead sample data:
 vol=slicer.util.getNode('MR*')
 vol.GetImageData().GetDimensions()
  • Capture the full Slicer screen and save it into a file
 img = qt.QPixmap.grabWidget(slicer.util.mainWindow()).toImage()
 img.save('c:/tmp/test.png')

Launching Slicer

  • How to open an .mrb file with Slicer at the command line?
 Slicer.exe --python-code "slicer.util.loadScene( 'f:/2013-08-23-Scene.mrb' )"
  • How to run a script in the Slicer environment in batch mode (without showing any graphical user interface)?
 Slicer.exe --python-code "doSomething; doSomethingElse; etc." --testing --no-splash --no-main-window

DICOM

  • How to access tags of DICOM images imported into Slicer? For example, to print the first patient's first study's first series' "0020,0032" field:
 db=slicer.dicomDatabase
 patientList=db.patients()
 studyList=db.studiesForPatient(patientList[0])
 seriesList=db.seriesForStudy(studyList[0])
 fileList=db.filesForSeries(seriesList[0])
 print db.fileValue(fileList[0],'0020,0032')
  • How to access tag of a volume loaded from DICOM? For example, get the patient position stored in a volume:
 volumeName='2: ENT IMRT'
 n=slicer.util.getNode(volumeName)
 instUids=n.GetAttribute('DICOM.instanceUIDs').split()
 filename=slicer.dicomDatabase.fileForInstance(instUids[0])
 print slicer.dicomDatabase.fileValue(filename,'0018,5100')
  • How to access tag of an item in the Subject Hierachy tree? For example, get the content time tag of a structure set:
 rtStructName = '3: RTSTRUCT: PROS'
 rtStructNode = slicer.util.getNode(rtStructName)
 rtStructSubjectHierarchyNode = slicer.vtkMRMLSubjectHierarchyNode.GetAssociatedSubjectHierarchyNode(rtStructNode)
 ctSliceInstanceUids = rtStructSubjectHierarchyNode.GetAttribute('DICOM.ReferencedInstanceUIDs').split()
 filename = slicer.dicomDatabase.fileForInstance(ctSliceInstanceUids[0])
 print slicer.dicomDatabase.fileValue(filename,'0008,0033')

Toolbar functions

  • How to turn on slice intersections in the crosshair menu on the toolbar:
 viewNodes = slicer.mrmlScene.GetNodesByClass('vtkMRMLSliceCompositeNode')
 viewNodes.UnRegister(slicer.mrmlScene)
 viewNodes.InitTraversal()
 viewNode = viewNodes.GetNextItemAsObject()
 while viewNode:
   viewNode.SetSliceIntersectionVisibility(1)
   viewNode = viewNodes.GetNextItemAsObject()

How to find similar functions? For this one I searched for "slice intersections" text in the whole slicer source code, found that the function is implemented in Base\QTGUI\qSlicerViewersToolBar.cxx, then translated the qSlicerViewersToolBarPrivate::setSliceIntersectionVisible(bool visible) method to Python.

Manipulating objects in the slice viewer

  • How to define/edit a circular region of interest in a slice viewer?

Drop two markup points on a slice view and copy-paste the code below into the Python console. After this, as you move the markups you’ll see a circle following the markups. 

 # Update the sphere from the fiducial points
 def UpdateSphere(param1, param2):  
   centerPointCoord = [0.0, 0.0, 0.0]
   markups.GetNthFiducialPosition(0,centerPointCoord)
   circumferencePointCoord = [0.0, 0.0, 0.0]
   markups.GetNthFiducialPosition(1,circumferencePointCoord)
   sphere.SetCenter(centerPointCoord)
   radius=math.sqrt((centerPointCoord[0]-circumferencePointCoord[0])**2+(centerPointCoord[1]-circumferencePointCoord[1])**2+(centerPointCoord[2]-circumferencePointCoord[2])**2)
   sphere.SetRadius(radius)
   sphere.SetPhiResolution(30)
   sphere.SetThetaResolution(30)
   sphere.Update()
 
 # Get a reference to the markup
 markups=slicer.util.getNode('F')
 # Create the sphere that will intersect the slice viewer
 sphere = vtk.vtkSphereSource()
 # Initial positioning of the sphere
 UpdateSphere(0,0)
 # Create model node and add to scene
 model = slicer.vtkMRMLModelNode()
 model.SetAndObservePolyData(sphere.GetOutput())
 modelDisplay = slicer.vtkMRMLModelDisplayNode()
 modelDisplay.SetSliceIntersectionVisibility(True) # Show in slice view
 modelDisplay.SetVisibility(False) # Hide in 3D view
 slicer.mrmlScene.AddNode(modelDisplay)
 model.SetAndObserveDisplayNodeID(modelDisplay.GetID())
 modelDisplay.SetInputPolyData(model.GetPolyData())
 slicer.mrmlScene.AddNode(model) 
 # Call UpdateSphere whenever the fiducials are changed
 markups.AddObserver("ModifiedEvent", UpdateSphere, 2)

Add a texture mapped plane to the scene as a model

Note that model textures are not exposed in the GUI and are not saved in the scene 

# use dummy image data here
e = vtk.vtkImageEllipsoidSource()

scene = slicer.mrmlScene

# Create model node
model = slicer.vtkMRMLModelNode()
model.SetScene(scene)
model.SetName(scene.GenerateUniqueName("2DImageModel"))

planeSource = vtk.vtkPlaneSource()
model.SetAndObservePolyData(planeSource.GetOutput())

# Create display node
modelDisplay = slicer.vtkMRMLModelDisplayNode()
modelDisplay.SetColor(1,1,0) # yellow
modelDisplay.SetBackfaceCulling(0)
modelDisplay.SetScene(scene)
scene.AddNode(modelDisplay)
model.SetAndObserveDisplayNodeID(modelDisplay.GetID())

# Add to scene
modelDisplay.SetAndObserveTextureImageData(e.GetOutput())
scene.AddNode(model) 


transform = slicer.vtkMRMLLinearTransformNode()
scene.AddNode(transform) 
model.SetAndObserveTransformNodeID(transform.GetID())

vTransform = vtk.vtkTransform()
vTransform.Scale(50,50,50)
vTransform.RotateX(30)
transform.SetAndObserveMatrixTransformToParent(vTransform.GetMatrix())


Export a model to Blender, including color

plyFilePath = "/tmp/fibers.ply"

lineDisplayNode = getNode("*LineDisplay*")

tuber = vtk.vtkTubeFilter()
tuber.SetInput(lineDisplayNode.GetOutputPolyData())

tubes = tuber.GetOutput()
tubes.Update()
scalars = tubes.GetPointData().GetArray(0)
scalars.SetName("scalars")

triangles = vtk.vtkTriangleFilter()
triangles.SetInput(tubes)

colorNode = lineDisplayNode.GetColorNode()
lookupTable = vtk.vtkLookupTable()
lookupTable.DeepCopy(colorNode.GetLookupTable())
lookupTable.SetTableRange(0,1)

plyWriter = vtk.vtkPLYWriter()
plyWriter.SetInput(triangles.GetOutput())
plyWriter.SetLookupTable(lookupTable)
plyWriter.SetArrayName("scalars")

plyWriter.SetFileName(plyFilePath)
plyWriter.Write()

Clone a volume

This example shows how to clone the MRHead sample volume, including its pixel data and display settings. 

sourceVolumeNode = slicer.util.getNode('MRHead')
volumesLogic = slicer.modules.volumes.logic()
clonedVolumeNode = volumesLogic.CloneVolume(slicer.mrmlScene, sourceVolumeNode, 'Cloned volume')

Create a new volume

This example shows how to create a new empty volume. 

imageSize=[512, 512, 512]
imageSpacing=[1.0, 1.0, 1.0]
voxelType=vtk.VTK_UNSIGNED_CHAR
# Create an empty image volume
imageData=vtk.vtkImageData()
imageData.SetDimensions(imageSize)
imageData.AllocateScalars(voxelType, 1)
thresholder=vtk.vtkImageThreshold()
thresholder.SetInputData(imageData)
thresholder.SetInValue(0)
thresholder.SetOutValue(0)
# Create volume node
volumeNode=slicer.vtkMRMLScalarVolumeNode()
volumeNode.SetSpacing(imageSpacing)
volumeNode.SetImageDataConnection(thresholder.GetOutputPort())
# Add volume to scene
slicer.mrmlScene.AddNode(volumeNode)
displayNode=slicer.vtkMRMLScalarVolumeDisplayNode()
slicer.mrmlScene.AddNode(displayNode)
colorNode = slicer.util.getNode('Grey')
displayNode.SetAndObserveColorNodeID(colorNode.GetID())
volumeNode.SetAndObserveDisplayNodeID(displayNode.GetID())
volumeNode.CreateDefaultStorageNode()

Modify voxels in a volume

This example shows how to change voxels values of the MRHead sample volume. The values will be computed by function f(r,a,s,) = (r-10)*(r-10)+(a+15)*(a+15)+s*s. 

volumeNode=slicer.util.getNode('MRHead')
ijkToRas = vtk.vtkMatrix4x4()
volumeNode.GetIJKToRASMatrix(ijkToRas)
imageData=volumeNode.GetImageData()
extent = imageData.GetExtent()
for k in xrange(extent[4], extent[5]+1):
  for j in xrange(extent[2], extent[3]+1):
    for i in xrange(extent[0], extent[1]+1):
      position_Ijk=[i, j, k, 1]
      position_Ras=ijkToRas.MultiplyPoint(position_Ijk)
      r=position_Ras[0]
      a=position_Ras[1]
      s=position_Ras[2]      
      functionValue=(r-10)*(r-10)+(a+15)*(a+15)+s*s
      imageData.SetScalarComponentFromDouble(i,j,k,0,functionValue)
imageData.SetScalarComponentFromFloat(distortionVectorPosition_Ijk[0], distortionVectorPosition_Ijk[1], distortionVectorPosition_Ijk[2], 0, fillValue)
imageData.Modified()

Access values in a DTI tensor volume

This example shows how to access individual tensors at the voxel level.

First load your DWI volume and estimate tensors to produce a DTI volume called ‘Output DTI Volume’

Then open the python window: View->Python interactor

Use this command to access tensors through numpy: 

tensors = array('Output DTI Volume')

Type the following code into the Python window to access all tensor components using vtk commands: 

volumeNode=slicer.util.getNode('Output DTI Volume')
imageData=volumeNode.GetImageData()
tensors = imageData.GetPointData().GetTensors()
extent = imageData.GetExtent()
idx = 0
for k in xrange(extent[4], extent[5]+1):
  for j in xrange(extent[2], extent[3]+1):
    for i in xrange(extent[0], extent[1]+1):
      tensors.GetTuple9(idx)
      idx += 1

Change window/level (brightness/contrast) or colormap of a volume

This example shows how to change window/level of the MRHead sample volume. 

volumeNode = getNode('MRHead')
displayNode = volumeNode.GetDisplayNode()
displayNode.AutoWindowLevelOff()
displayNode.SetWindow(50)
displayNode.SetLevel(100)

Change color mapping from grayscale to rainbow: 

displayNode.SetAndObserveColorNodeID('vtkMRMLColorTableNodeRainbow')

Manipulate a Slice View

lm = slicer.app.layoutManager()
red = lm.sliceWidget('Red')
redLogic = red.sliceLogic()
# Print current slice offset position
print redLogic.GetSliceOffset()
# Change slice position
redLogic.SetSliceOffset(20)

Save a series of images from a Slice View

Save the following into a file such as '/tmp/record.py' and then in the slicer python console type "execfile('/tmp/record.py')" 

layoutName = 'Green'
imagePathPattern = '/tmp/image-%03d.png'
steps = 10

widget = slicer.app.layoutManager().sliceWidget(layoutName)
view = widget.sliceView()
logic = widget.sliceLogic()
bounds = [0,]*6
logic.GetSliceBounds(bounds)

for step in range(steps):
    offset = bounds[4] + step/(1.*steps) * (bounds[5]-bounds[4])
    logic.SetSliceOffset(offset)
    view.forceRender()
    image = qt.QPixmap.grabWidget(view).toImage()
    image.save(imagePathPattern % step)

Show a volume in the Slice Views

volumeNode = slicer.util.getNode('YourVolumeNode')
applicationLogic = slicer.app.applicationLogic()
selectionNode = applicationLogic.GetSelectionNode()
selectionNode.SetSecondaryVolumeID(volumeNode.GetID())
applicationLogic.PropagateForegroundVolumeSelection(0)

or 

n =  slicer.util.getNode('YourVolumeNode')
for color in ['Red', 'Yellow', 'Green']:
    slicer.app.layoutManager().sliceWidget(color).sliceLogic().GetSliceCompositeNode().SetForegroundVolumeID(n.GetID())

Change opacity of foreground volume in the Slice Views

lm = slicer.app.layoutManager()
sliceLogic = lm.sliceWidget('Red').sliceLogic()
compositeNode = sliceLogic.GetSliceCompositeNode()
compositeNode.SetForegroundOpacity(0.4)

Center the 3D View on the Scene

layoutManager = slicer.app.layoutManager()
threeDWidget = layoutManager.threeDWidget(0)
threeDView = threeDWidget.threeDView()
threeDView.resetFocalPoint()

Display text in a 3D view or slice view

The easiest way to show information overlaid on a viewer is to use corner annotations. 

view=slicer.app.layoutManager().threeDWidget(0).threeDView()
# Set text to "Something"
view.cornerAnnotation().SetText(vtk.vtkCornerAnnotation.UpperRight,"Something")
# Set color to red
view.cornerAnnotation().GetTextProperty().SetColor(1,0,0)
# Update the view
view.forceRender()

Turning off interpolation

You can turn off interpolation for newly loaded volumes with this script from Steve Pieper. 

def NoInterpolate(caller,event):
  for node in slicer.util.getNodes('*').values():
    if node.IsA('vtkMRMLScalarVolumeDisplayNode'):
      node.SetInterpolate(0)
	
slicer.mrmlScene.AddObserver(slicer.mrmlScene.NodeAddedEvent, NoInterpolate)

The below link explains how to put this in your startup script.

http://www.na-mic.org/Wiki/index.php/AHM2012-Slicer-Python#Refining_the_code_and_UI_with_slicerrc


Customize viewer layout

Show a custom layout of a 3D view on top of the red slice view: 

customLayout = ("<layout type=\"vertical\" split=\"true\" >"
  " <item>"
  "  <view class=\"vtkMRMLViewNode\" singletontag=\"1\">"
  "    <property name=\"viewlabel\" action=\"default\">1</property>"
  "  </view>"
  " </item>"
  " <item>"
  "  <view class=\"vtkMRMLSliceNode\" singletontag=\"Red\">"
  "   <property name=\"orientation\" action=\"default\">Axial</property>"
  "   <property name=\"viewlabel\" action=\"default\">R</property>"
  "   <property name=\"viewcolor\" action=\"default\">#F34A33</property>"
  "  </view>"
  " </item>"
  "</layout>")
  
customLayoutId=501

layoutManager = slicer.app.layoutManager()
layoutManager.layoutLogic().GetLayoutNode().AddLayoutDescription(customLayoutId, customLayout)                                         
layoutManager.setLayout(customLayoutId)

See description of standard layouts (that can be used as examples) here: https://github.com/Slicer/Slicer/blob/master/Libs/MRML/Logic/vtkMRMLLayoutLogic.cxx

Running an ITK filter in Python using SimpleITK

Open the "Sample Data" module and download "MR Head", then paste the following snippet in Python interactor: 

import SimpleITK as sitk
import sitkUtils
inputImage = sitkUtils.PullFromSlicer('MRHead')
filter = sitk.SignedMaurerDistanceMapImageFilter()
outputImage = filter.Execute(inputImage)
sitkUtils.PushToSlicer(outputImage,'outputImage')

More information:

Get current mouse coordinates in a slice view

You can get 3D (RAS) coordinates of the current mouse cursor from the crosshair singleton node as shown in the example below: 

def onMouseMoved(observer,eventid):  
  ras=[0,0,0]
  crosshairNode.GetCursorPositionRAS(ras)
  print(ras)

crosshairNode=slicer.util.getNode('Crosshair') 
crosshairNode.AddObserver(slicer.vtkMRMLCrosshairNode.CursorPositionModifiedEvent, onMouseMoved)

Thick slab reconstruction and maximum/minimum intensity volume projections

Set up 'red' slice viewer to show thick slab reconstructed from 3 slices: 

sliceNode = slicer.mrmlScene.GetNodeByID('vtkMRMLSliceNodeRed')
appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sliceLayerLogic = sliceLogic.GetBackgroundLayer()
reslice = sliceLayerLogic.GetReslice()
reslice.SetSlabModeToMean()
reslice.SetSlabNumberOfSlices(10) # mean of 10 slices will computed
reslice.SetSlabSliceSpacingFraction(0.3) # spacing between each slice is 0.3 pixel (total 10 * 0.3 = 3 pixel neighborhood)
sliceNode.Modified()

Set up 'red' slice viewer to show maximum intensity projection (MIP): 

sliceNode = slicer.mrmlScene.GetNodeByID('vtkMRMLSliceNodeRed')
appLogic = slicer.app.applicationLogic()
sliceLogic = appLogic.GetSliceLogic(sliceNode)
sliceLayerLogic = sliceLogic.GetBackgroundLayer()
reslice = sliceLayerLogic.GetReslice()
reslice.SetSlabModeToMax()
reslice.SetSlabNumberOfSlices(600) # use a large number of slices (600) to cover the entire volume
reslice.SetSlabSliceSpacingFraction(0.5) # spacing between slices are 0.5 pixel (supersampling is useful to reduce interpolation artifacts)
sliceNode.Modified()

The projected image is available in a vtkImageData object by calling reslice.GetOutput().

Change default file type for nodes (that have never been saved yet)

Default node can be specified that will be used as a basis of all new storage nodes. This can be used for setting default file extension. For example, change file format to STL for model nodes: 

msn=slicer.vtkMRMLModelStorageNode()
msn.SetDefaultWriteFileExtension('stl')
slicer.mrmlScene.AddDefaultNode(msn)

Change file type for saving for all volumes (with already existing storage nodes)

requiredFileExtension = '.nia'
originalFileExtension = '.nrrd'
volumeNodes = slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
volumeNodes.UnRegister(slicer.mrmlScene)
volumeNodes.InitTraversal()
volumeNode = volumeNodes.GetNextItemAsObject()
while volumeNode:
  volumeStorageNode = volumeNode.GetStorageNode()
  if not volumeStorageNode:
    volumeStorageNode = volumeNode.CreateDefaultStorageNode()
    slicer.mrmlScene.AddNode(volumeStorageNode)
    volumeStorageNode.UnRegister(None)
    volumeNode.SetAndObserveStorageNodeID(volumeStorageNode.GetID())
    volumeStorageNode.SetFileName(volumeNode.GetName()+requiredFileExtension)
  else:
    volumeStorageNode.SetFileName(volumeStorageNode.GetFileName().replace(originalFileExtension,requiredFileExtension))
  volumeNode = volumeNodes.GetNextItemAsObject()

Segmentations

For all operations accessing or manupilating the internals of a segmentation will need this import! 

import vtkSegmentationCorePython as vtkSegmentationCore

Get a segment

segmentation = segmentationNode.GetSegmentation()
segment = segmentation.GetSegment(segmentID)

Get a representation of a segment

# Get representation from a single segment. If it does not exist, it will return None
segment.GetRepresentation(vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName()) # This gets the binary labelmap, but same idea for all others

# Get representation for a single segment. Convert temporarily for that particular segment if needed. Applies parent transforms by default (if not desired, another argument needs to be added to the end: false)
slicer.vtkSlicerSegmentationsModuleLogic.GetSegmentBinaryLabelmapRepresentation(segmentationNode, segmentID, outputOrientedImageData) # Get labelmap
slicer.vtkSlicerSegmentationsModuleLogic.GetSegmentRepresentation(segmentationNode, segmentID, vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationClosedSurfaceRepresentationName(), outputPolyData) # Any representation

Convert using default path and conversion parameters

segmentation.CreateRepresentation(vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName()) # This creates binary labelmap, but same idea for all others

Convert using custom path or conversion parameters

# Custom path
TODO

# Custom conversion parameter
referenceGeometry=vtkSegmentationCore.vtkSegmentationConverter.SerializeImageGeometry(referenceImageData)
segmentation.SetConversionParameter(vtkSegmentationCore.vtkSegmentationConverter.GetReferenceImageGeometryParameterName(), referenceGeometry)
Retrieved from "https://www.slicer.org/w/index.php?title=Documentation/4.6/ScriptRepository&oldid=49896"