Documentation/Nightly/Developers/SlicerExecutionModel/Python

Python Modules

A Python interpreter has been integrated into Slicer3. More details of the Python language can be found at Python.org. This interpreter may be used to execute Python script plugins. The modules have a few specific requirements to be correctly found and used as plugins. Like the standard executable and shared library plugins, Python plugins need to be self describing. To do this, the script must have a top level variable called XML or provide a toXML procedure. For example:

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>
  <category>Filtering.Denoising</category>
  ...

def toXML():
  return XML;

Details of the XML format are found in the main Execution Model documentation. Rather than construct a command line to pass into Python, Slicer3 directly calls an Execute procedure. It is assumed that the Execute function expects positional arguments first, and any optional arguments are passed in as keyword arguments. For instance, the GradientAnisotropicDiffusion.py module provides an Execute:

def Execute ( inputVolume, outputVolume, conductance=1.0, timeStep=0.0625, iterations=1 ):
    print "Executing Python Demo Application!"
    Slicer = __import__ ( "Slicer" )
    slicer = Slicer.slicer
    in = slicer.MRMLScene.GetNodeByID ( inputVolume );
    out = slicer.MRMLScene.GetNodeByID ( outputVolume );

    filter = slicer.vtkITKGradientAnisotropicDiffusionImageFilter()
    filter.SetConductanceParameter ( conductance )
    filter.SetTimeStep ( timeStep )
    filter.SetNumberOfIterations ( iterations )
    filter.SetInput ( in.GetImageData() )
    filter.Update()
    out.SetAndObserveImageData(filter.GetOutput())
    matrix = slicer.vtkMatrix4x4()
    in.GetIJKToRASMatrix(matrix)
    out.GetIJKToRASMatrix(matrix)
    return

The function first constructs a Slicer object by importing the Slicer module. The Slicer object is the main interface into Slicer3 as a whole. The first two arguments, inputVolume and outputVolume are not proper MRMLVolumes, and must be looked up using the Slicer object (that is, the arguments are unique IDs rather than pointers to the objects). The filter is constructed through the Slicer object, and the parameters are set. After the filter is updated, the output image is put using the SetAndObserveImageData method on the output volume.

ToDo

• Progress functionality
• Casting image arguments to proper MRML volume objects before calling Execute

Python Module collection

Here are a few sample modules that should clarify how to write CLI modules in Python. Some of them are official Python modules, other are here just to provide working examples.

• PythonScript.py
• SurfaceConnectivity.py
• SurfaceICPRegistration.py
• SurfaceToolbox.py
• ExplodeVolumeTransform.py
• GaussianSmooth.py
• GradientAnisotropicDiffusion.py
• ResampleVolume.py
• ResliceAsVolume.py

To create a working Python CLI module, just create a directory, put your .py file there and add the directory path to the path list in Application settings - Module settings.

PythonScript.py

This module allows to run external Python code saved in a .py file on an image or surface (or both), and have the results loaded back in Slicer. It doesn't really do anything itself, but it's handy for:

• writing Python modules, since one can store the actual module code in a file without having to copy it to Slicer3-build/lib/Slicer3/Plugins
• writing quick one-shot filtering operations that are not elegible to become full-fledged modules
• experimenting with Python in Slicer.

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Python Modules</category>
  <title>Python Script</title>
  <description>Run external Python code on a surface or image.</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>

  <parameters>
    <label>Python Script Parameters</label>
    <description>Parameters for Python script</description>
    <file>
      <name>scriptFileName</name>
      <longflag>scriptFileName</longflag>
      <description>File containing Python code to run</description>
      <label>Script file</label>
    </file>
  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>

    <geometry>
      <name>inputSurface</name>
      <longflag>inputSurface</longflag>
      <label>Input Surface</label>
      <channel>input</channel>
      <description>Input surface to be filtered</description>
    </geometry>

    <image>
      <name>inputVolume</name>
      <longflag>inputVolume</longflag>
      <label>Input Volume</label>
      <channel>input</channel>
      <description>Input image to be filtered</description>
    </image>

    <geometry>
      <name>outputSurface</name>
      <longflag>outputSurface</longflag>
      <label>Output Surface</label>
      <channel>output</channel>
      <description>Output filtered surface</description>
    </geometry>

    <image>
     <name>outputVolume</name>
     <longflag>outputVolume</longflag>
     <label>Output Volume</label>
     <channel>output</channel>
     <description>Output filtered volume</description>
   </image>


  </parameters>

</executable>
"""


def Execute (scriptFileName="", inputSurface="", inputVolume="", outputSurface="", outputVolume=""):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene

    if inputSurface:
        inputSurface = scene.GetNodeByID(inputSurface)

    if inputVolume:
        inputVolume = scene.GetNodeByID(inputVolume)

    if outputSurface:
        outputSurface = scene.GetNodeByID(outputSurface)

   if outputVolume:
       outputVolume = scene.GetNodeByID(outputVolume)

    try:
        execfile(scriptFileName)
    except Exception, error:
        slicer.Application.ErrorMessage("Python script error: %s" % error)

    return

Note how the try...except block intercepts Python exceptions generated while running the external script and prints messages out in Slicer's message console.

Now let's see what an external script code could look like. Say we want to take the input surface and decimate it.

inputPolyData = inputSurface.GetPolyData()

decimation = slicer.vtkDecimatePro()
decimation.SetInput(inputSurface.GetPolyData())
decimation.SetTargetReduction(0.99)
decimation.Update()

outputSurface.SetAndObservePolyData(decimation.GetOutput())

In the script (which we can save in a .py file), we don't have to start by importing Slicer module, since it's already done in the PythonScript.py module. We have to assume that inputSurface, outputSurface, inputVolume and outputVolume are already available as MRML nodes (vtkMRMLModelNode and vtkMRMLVolumeNode, respectively), and that the slicer module is also already imported.

Setting of the filtered poly data to the outputSurface is left to the user (see last line).

In general, the syntax used to use VTK classes in a Slicer Python module is very similar to using VTK from Python. There is just one main difference: classes have to be looked for in slicer, not vtk (i.e. slicer.vtkImageData() instead of vtk.vtkImageData())

SurfaceConnectivity.py

This module runs a connectivity filter on an input surface. There are two modes of operation: one is to extract all regions and color them accordingly, the other is to extract the region closest to a fiducial. Optionally, the module will generate a list of fiducials named "Region #" placed on a point of each connected region. It's a nice example on how one can interact with the MRML scene. This is possible because Python modules are run in the main thread, in contrast with executables and shared library command line modules that are run in a separate processing thread.

   XML = """<?xml version="1.0" encoding="utf-8"?>
   <executable>
   
     <category>Python Modules</category>
     <title>Python Surface Connectivity</title>
     <description>
   Identify connected regions of a surface.
   </description>
     <version>1.0</version>
     <documentation-url></documentation-url>
     <license></license>
     <contributor>Luca Antiga and Daniel Blezek</contributor>
   
     <parameters>
       <label>Surface Connectivity Parameters</label>
       <description>Parameters for surface connectivity</description>
   
       <string-enumeration>
         <name>connectivityMode</name>
         <longflag>connectivityMode</longflag>
         <description>Mode of operation of connectivity filter</description>
         <label>Connectivity mode</label>
         <default>AllRegions</default>
         <element>AllRegions</element>
         <element>ClosestToSeed</element>
       </string-enumeration>
   
       <boolean>
         <name>enableOutputFiducials</name>
         <longflag>enableOutputFiducials</longflag>
         <description>Toggle generation of labeled fiducials corresponding to each extracted region</description>
         <label>Enable output fiducials</label>
         <default>false</default>
       </boolean>
   
       <point multiple="false">
         <name>seedPoint</name>
         <longflag>seedPoint</longflag>
         <label>Seed Point</label>
         <description>Fiducial to extract closest surface</description>
       </point>
     </parameters>
   
     <parameters>
       <label>IO</label>
       <description>Input/output parameters</description>
   
       <geometry>
         <name>inputSurface</name>
         <label>Input Surface</label>
         <channel>input</channel>
         <index>0</index>
         <description>Input surface to be filtered</description>
       </geometry>
   
       <geometry>
         <name>outputSurface</name>
         <label>Output Surface</label>
         <channel>output</channel>
         <index>1</index>
         <description>Output filtered surface</description>
       </geometry>
     </parameters>
   
   </executable>
   """
   
   
   def Execute (inputSurface, outputSurface, connectivityMode="AllRegions", enableOutputFiducials=False, seedPoint=[0.0,0.0,0.0]):
   
       Slicer = __import__("Slicer")
       slicer = Slicer.slicer
       scene = slicer.MRMLScene
       inputSurface = scene.GetNodeByID(inputSurface)
       outputSurface = scene.GetNodeByID(outputSurface)
   
       connectivityFilter = slicer.vtkPolyDataConnectivityFilter()
       connectivityFilter.SetInput(inputSurface.GetPolyData())
       if connectivityMode == "AllRegions":
           connectivityFilter.SetExtractionModeToAllRegions()
           connectivityFilter.ColorRegionsOn()
       elif connectivityMode == "ClosestToSeed":
           connectivityFilter.SetExtractionModeToClosestPointRegion()
           connectivityFilter.SetClosestPoint(*seedPoint)
       connectivityFilter.Update()
   
       if enableOutputFiducials:
           fiducialList = scene.CreateNodeByClass("vtkMRMLFiducialListNode")
           fiducialList.DisableModifiedEventOn()
           fiducialList.SetScene(scene)
           thresholdPoints = slicer.vtkThresholdPoints()
           thresholdPoints.SetInput(connectivityFilter.GetOutput())
           numberOfRegions = connectivityFilter.GetNumberOfExtractedRegions()
           for i in range(numberOfRegions):
               lower = i - 0.5
               upper = i + 0.5
               thresholdPoints.ThresholdBetween(lower,upper)
               thresholdPoints.Update()
               if thresholdPoints.GetOutput().GetNumberOfPoints() > 0:
                   point = thresholdPoints.GetOutput().GetPoint(0)
                   fid = fiducialList.AddFiducial()
                   fiducialList.SetNthFiducialXYZ(fid,*point)
                   fiducialList.SetNthFiducialLabelText(fid,"Region %d" % i)
           fiducialList.InvokePendingModifiedEvent()
           fiducialList.DisableModifiedEventOff()
           scene.AddNode(fiducialList)
   
       outputSurface.SetAndObservePolyData(connectivityFilter.GetOutput())
   
       inputSurface.GetDisplayNode().VisibilityOff()
   
       return

SurfaceICPRegistration.py

This module provides surface registration using the Iterative Closest Point algorithm as implemented in VTK. Nothing fancy going on, just nice functionality added to Slicer.

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Python Modules</category>
  <title>Python Surface ICP Registration</title>
  <description>
Performes registration of the input surface onto a target surface using on the Iterative Closest Point algorithm.
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>

  <parameters>
    <label>Surface ICP Registration Parameters</label>
    <description>Parameters for surface registration</description>

    <string-enumeration>
      <name>landmarkTransformMode</name>
      <longflag>landmarkTransformMode</longflag>
      <description>Mode of operation of internal landmark transform</description>
      <label>Landmark transform mode</label>
      <default>RigidBody</default>
      <element>RigidBody</element>
      <element>Similarity</element>
      <element>Affine</element>
    </string-enumeration>

    <string-enumeration>
      <name>meanDistanceMode</name>
      <longflag>meanDistanceMode</longflag>
      <label>Mean distance mode</label>
      <description>Set how the mean distance should be computed, RMS (square root of the average of the sum of squares 
        of the closest point distances) or AbsoluteValue (mean of the sum of absolute values of the closest point distances)</description>
      <default>RMS</default>
      <element>RMS</element>
      <element>AbsoluteValue</element>
    </string-enumeration>

    <integer>
      <name>maximumNumberOfIterations</name>
      <longflag>maximumNumberOfIterations</longflag>
      <label>Maximum number of iterations</label>
      <description>Maximum number of Interative Closest Point iterations</description>
      <default>50</default>
      <constraints>
        <minimum>1</minimum>
        <maximum>1000</maximum>
        <step>1</step>
      </constraints>
    </integer>

    <integer>
      <name>maximumNumberOfLandmarks</name>
      <longflag>maximumNumberOfLandmarks</longflag>
      <label>Maximum number of landmarks</label>
      <description>Maximum number of landmarks on source and target surfaces</description>
      <default>200</default>
      <constraints>
        <minimum>1</minimum>
        <maximum>1000</maximum>
        <step>1</step>
      </constraints>
    </integer>

    <boolean>
      <name>startByMatchingCentroids</name>
      <longflag>startByMatchingCentroids</longflag>
      <label>Start by matching centroids</label>
      <description>Toggle starting registration by first translating the source centroid to the target centroid</description>
      <default>false</default>
    </boolean>

    <boolean>
      <name>checkMeanDistance</name>
      <longflag>checkMeanDistance</longflag>
      <label>Check mean distance</label>
      <description>Force the algorithm to check the mean distance between two iterations</description>
      <default>false</default>
    </boolean>

    <double>
      <name>maximumMeanDistance</name>
      <longflag>maximumMeanDistance</longflag>
      <label>Maximum mean distance</label>
      <description>Maximum distance between two iterations used for determining convergence</description>
      <default>0.01</default>
    </double>

  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>

    <geometry>
      <name>inputSurface</name>
      <label>Input Surface</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input surface to be registered</description>
    </geometry>

    <geometry>
      <name>targetSurface</name>
      <label>Target Surface</label>
      <channel>input</channel>
      <index>1</index>
      <description>Target surface for registration</description>
    </geometry>

    <geometry>
      <name>outputSurface</name>
      <label>Output Surface</label>
      <channel>output</channel>
      <index>2</index>
      <description>Output registered surface</description>
    </geometry>
  </parameters>

</executable>
"""


def Execute (inputSurface, targetSurface, outputSurface, \
            landmarkTransformMode="RigidBody", meanDistanceMode="RMS", \
            maximumNumberOfIterations=50, maximumNumberOfLandmarks=200, \
            startByMatchingCentroids=False, checkMeanDistance=False, maximumMeanDistance=0.01):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene
    inputSurface = scene.GetNodeByID(inputSurface)
    targetSurface = scene.GetNodeByID(targetSurface)
    outputSurface = scene.GetNodeByID(outputSurface)

    icpTransform = slicer.vtkIterativeClosestPointTransform()
    icpTransform.SetSource(inputSurface.GetPolyData())
    icpTransform.SetTarget(targetSurface.GetPolyData())
    if landmarkTransformMode == "RigidBody":
        icpTransform.GetLandmarkTransform().SetModeToRigidBody()
    elif landmarkTransformMode == "Similarity":
        icpTransform.GetLandmarkTransform().SetModeToSimilarity()
    elif landmarkTransformMode == "Affine":
        icpTransform.GetLandmarkTransform().SetModeToAffine()
    if meanDistanceMode == "RMS":
        icpTransform.SetMeanDistanceModeToRMS()
    elif meanDistanceMode == "AbsoluteValue":
        icpTransform.SetMeanDistanceModeToAbsoluteValue()
    icpTransform.SetMaximumNumberOfIterations(maximumNumberOfIterations)
    icpTransform.SetMaximumNumberOfLandmarks(maximumNumberOfLandmarks)
    icpTransform.SetStartByMatchingCentroids(int(startByMatchingCentroids))
    icpTransform.SetCheckMeanDistance(int(checkMeanDistance))
    icpTransform.SetMaximumMeanDistance(maximumMeanDistance)

    transformFilter = slicer.vtkTransformPolyDataFilter()
    transformFilter.SetInput(inputSurface.GetPolyData())
    transformFilter.SetTransform(icpTransform)        
    transformFilter.Update()
    
    outputSurface.SetAndObservePolyData(transformFilter.GetOutput())

    return

SurfaceToolbox.py

This module provides a bunch of operations that one typically performs on a surface, namely decimation, smoothing, computation of normals and surface cleaning. The filters are executed in a pipe after they are manually enabled with the Enabled checkbox.

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Python Modules</category>
  <title>Python Surface Toolbox</title>
  <description>
Process a surface using several filters.
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>

  <parameters>
    <label>Surface Decimation</label>
    <description>Parameters for surface decimation</description>
    <boolean>
      <name>decimationEnabled</name>
      <longflag>decimationEnabled</longflag>
      <description>Toggle decimation</description>
      <label>Enabled</label>
      <default>false</default>
    </boolean>
    <double>
      <name>targetReduction</name>
      <longflag>targetReduction</longflag>
      <description>Target percent reduction of the number of polygons</description>
      <label>Target reduction</label>
      <default>0.8</default>
      <constraints>
        <minimum>0</minimum>
        <maximum>1</maximum>
        <step>0</step>
      </constraints>
    </double>
    <boolean>
      <name>boundaryDeletion</name>
      <longflag>boundaryDeletion</longflag>
      <description>Toggle deletion of boundary vertices</description>
      <label>Boundary vertex deletion</label>
      <default>true</default>
    </boolean>
  </parameters>

  <parameters>
    <label>Surface Smoothing</label>
    <description>Parameters for surface smoothing</description>
    <boolean>
      <name>smoothingEnabled</name>
      <longflag>smoothingEnabled</longflag>
      <description>Toggle smoothing</description>
      <label>Enabled</label>
      <default>false</default>
    </boolean>
    <string-enumeration>
      <name>smoothingMethod</name>
      <longflag>smoothingMethod</longflag>
      <description>Selection of smoothing method between Laplace and Taubin</description>
      <label>Smoothing method</label>
      <default>Taubin</default>
      <element>Laplace</element>
      <element>Taubin</element>
    </string-enumeration>
    <boolean>
      <name>boundarySmoothing</name>
      <longflag>boundarySmoothing</longflag>
      <description>Toggle allow smoothing of boundary</description>
      <label>Boundary smoothing</label>
      <default>true</default>
    </boolean>
    <integer>
      <name>laplaceIterations</name>
      <longflag>laplaceIterations</longflag>
      <description>Number of iterations for the Laplace smoothing method</description>
      <label>Laplace number of iterations</label>
      <default>10</default>
    </integer>
    <double>
      <name>laplaceRelaxation</name>
      <longflag>laplaceRelaxation</longflag>
      <description>Relaxation factor value for the Laplace smoothing method</description>
      <label>Laplace relaxation factor</label>
      <default>0.1</default>
    </double>
    <integer>
      <name>taubinIterations</name>
      <longflag>taubinIterations</longflag>
      <description>Number of iterations for the Taubin smoothing method</description>
      <label>Taubin number of iterations</label>
      <default>10</default>
    </integer>
    <double>
      <name>taubinPassband</name>
      <longflag>taubinPassband</longflag>
      <description>Passband value for the Taubin smoothing method</description>
      <label>Taubin passband</label>
      <default>0.1</default>
    </double>
  </parameters>

  <parameters>
    <label>Surface Normals</label>
    <description>Parameters for surface normals</description>
    <boolean>
      <name>normalsEnabled</name>
      <longflag>normalsEnabled</longflag>
      <description>Toggle normals</description>
      <label>Enabled</label>
      <default>false</default>
    </boolean>
    <boolean>
      <name>flipNormals</name>
      <longflag>flipNormals</longflag>
      <description>Toggle flip normals after computing</description>
      <label>Flip normals</label>
      <default>false</default>
    </boolean>
    <boolean>
      <name>splitting</name>
      <longflag>splitting</longflag>
      <description>Toggle splitting surface based on feature angle</description>
      <label>Splitting</label>
      <default>false</default>
    </boolean>
    <double>
      <name>featureAngle</name>
      <longflag>featureAngle</longflag>
      <description>Feature angle value for surface splitting</description>
      <label>Feature angle</label>
      <default>30</default>
      <constraints>
        <minimum>0</minimum>
        <maximum>180</maximum>
        <step>1</step>
      </constraints>
    </double>
  </parameters>

  <parameters>
    <label>Surface Cleaner</label>
    <description>Parameters for surface cleaner</description>
    <boolean>
      <name>cleanerEnabled</name>
      <longflag>cleanerEnabled</longflag>
      <description>Toggle cleaner</description>
      <label>Enabled</label>
      <default>false</default>
    </boolean>
  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>

    <geometry>
      <name>inputSurface</name>
      <label>Input Surface</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input surface to be filtered</description>
    </geometry>

    <geometry>
      <name>outputSurface</name>
      <label>Output Surface</label>
      <channel>output</channel>
      <index>1</index>
      <description>Output filtered surface</description>
    </geometry>
  </parameters>

</executable>
"""


def Execute (inputSurface, outputSurface, \
    decimationEnabled=False, targetReduction=0.8, boundaryDeletion=True, \
    smoothingEnabled=False, smoothingMethod="Taubin", boundarySmoothing=True, \
    laplaceIterations=10, laplaceRelaxation=0.1, taubinIterations=10, taubinPassband=0.1, \
    normalsEnabled=False, flipNormals=False, splitting=False, featureAngle=30.0, \
    cleanerEnabled=False):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene
    inputSurface = scene.GetNodeByID(inputSurface)
    outputSurface = scene.GetNodeByID(outputSurface)

    surface = inputSurface.GetPolyData()

    if decimationEnabled:
        decimation = slicer.vtkDecimatePro()
        decimation.SetInput(surface)
        decimation.SetTargetReduction(targetReduction)
        if boundaryDeletion:
            decimation.BoundaryVertexDeletionOn()
        else:
            decimation.BoundaryVertexDeletionOff()
        decimation.PreserveTopologyOn()
        decimation.Update()
        surface = decimation.GetOutput()

    if smoothingEnabled:
        if smoothingMethod == "Laplace":
            smoothing = slicer.vtkSmoothPolyDataFilter()
            smoothing.SetInput(surface)
            if boundarySmoothing:
                smoothing.BoundarySmoothingOn()
            else:
                smoothing.BoundarySmoothingOff()
            smoothing.SetNumberOfIterations(laplaceIterations)
            smoothing.SetRelaxationFactor(laplaceRelaxation)
            smoothing.Update()
            surface = smoothing.GetOutput()
        elif smoothingMethod == "Taubin":
            smoothing = slicer.vtkWindowedSincPolyDataFilter.New()
            smoothing.SetInput(surface)
            if boundarySmoothing:
                smoothing.BoundarySmoothingOn()
            else:
                smoothing.BoundarySmoothingOff()
            smoothing.SetNumberOfIterations(taubinIterations)
            smoothing.SetPassBand(taubinPassband)
            smoothing.Update()
            surface = smoothing.GetOutput()

    if normalsEnabled:
        normals = slicer.vtkPolyDataNormals()
        normals.SetInput(surface)
        normals.AutoOrientNormalsOn()
        if flipNormals:
            normals.FlipNormalsOn()
        else:
            normals.FlipNormalsOff()
        if splitting:
            normals.SplittingOn()
            normals.SetFeatureAngle(featureAngle)
        else:
            normals.SplittingOff()
        normals.ConsistencyOn()
        normals.Update()
        surface = normals.GetOutput()

    if cleanerEnabled:
        cleaner = slicer.vtkCleanPolyData()
        cleaner.SetInput(surface)
        cleaner.Update()
        surface = cleaner.GetOutput()

    outputSurface.SetAndObservePolyData(surface)

    return

GaussianSmooth.py

This module runs VTK Gaussian smoothing filter on an input image. Note that the Gaussian smoothing directions here are referred to the IJK directions (VTK's vtkImageData is a non-oriented image), but the volume is then placed back in the appropriate RAS reference system. Take this into account when writing CLI modules relying on non-reference system aware libraries.

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Filtering</category>
  <title>Python Gaussian Smoothing</title>
  <description>
Apply Gaussian smoothing on a volume.
This module uses the implementation given in vtkImageGaussianSmooth.
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>
  <acknowledgements>
This work is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149. 
</acknowledgements>

  <parameters>
    <label>Gaussian Smooth Parameters</label>
    <description>Parameters for Gaussian smoothing</description>

    <double-vector>
      <name>standardDeviations</name>
      <longflag>standardDeviations</longflag>
      <description>Standard deviations in pixel units (ordered as IJK coordinates)</description>
      <label>Standard deviations</label>
      <default>1,1,1</default>
    </double-vector>

    <double-vector>
      <name>radiusFactors</name>
      <longflag>radiusFactors</longflag>
      <description>Radius factors (in standard deviation units) for clamping the Gaussian kernel to zero</description>
      <label>Radius factors</label>
      <default>1.5,1.5,1.5</default>
    </double-vector>

    <integer>
      <name>dimensionality</name>
      <longflag>dimensionality</longflag>
      <description>Dimensionality of the Gaussian filtering</description>
      <label>Dimensionality</label>
      <default>3</default>
      <constraints>
        <minimum>2</minimum>
        <maximum>3</maximum>
      </constraints>
    </integer>
  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>
    <image>
      <name>inputVolume</name>
      <label>Input Volume</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input volume to be filtered</description>
    </image>
    <image>
      <name>outputVolume</name>
      <label>Output Volume</label>
      <channel>output</channel>
      <index>1</index>
      <description>Output filtered volume</description>
    </image>
  </parameters>

</executable>
"""


def Execute (inputVolume, outputVolume, standardDeviations=[1.0,1.0,1.0], radiusFactors=[1.5,1.5,1.5], dimensionality=3):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene
    inputVolume = scene.GetNodeByID(inputVolume)
    outputVolume = scene.GetNodeByID(outputVolume)

    filter = slicer.vtkImageGaussianSmooth()
    filter.SetStandardDeviations(*standardDeviations)
    filter.SetRadiusFactors(*radiusFactors)
    filter.SetDimensionality(dimensionality)
    filter.SetInput(inputVolume.GetImageData())
    filter.Update()

    outputVolume.SetAndObserveImageData(filter.GetOutput())
    matrix = slicer.vtkMatrix4x4()
    inputVolume.GetIJKToRASMatrix(matrix)
    outputVolume.SetIJKToRASMatrix(matrix)

    return

GradientAnisotropicDiffusion.py

This module is yet another version of the ubiquitous GradientAnisotropicDiffusion module. The notable thing of this module is that the filter is implemented in ITK, wrapped in a VTK class using the vtkITK library and made available in Python through Tcl. Enough?

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>
  <category>Filtering</category>
  <title>Python Gradient Anisotropic Diffusion</title>
  <description>
Runs gradient anisotropic diffusion on a volume.

Anisotropic diffusion methods reduce noise (or unwanted detail) in
images while preserving specific image features.  For many
applications, there is an assumption that light-dark transitions
(edges) are interesting.  Standard isotropic diffusion methods move
and blur light-dark boundaries.  Anisotropic diffusion methods are
formulated to specifically preserve edges.
The conductance term for this implementation is chosen as a function
of the gradient magnitude of the image at each point, reducing the
strength of diffusion at edge pixels.
The numerical implementation of this equation is similar to that
described in the Perona-Malik paper, but uses a more robust technique
for gradient magnitude estimation and has been generalized to
N-dimensions.

</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Daniel Blezek</contributor>
  <acknowledgements>This command module was derived from the Gradient Anisotropic Diffusion CLI example</acknowledgements>

  <parameters>
    <label>Anisotropic Diffusion Parameters</label>
    <description>Parameters for the anisotropic diffusion algorithm</description>

    <double>
      <name>conductance</name>
      <longflag>--conductance</longflag>
      <description>Conductance</description>
      <label>Conductance</label>
      <default>1</default>
      <constraints>
        <minimum>0</minimum>
        <maximum>10</maximum>
        <step>.01</step>
      </constraints>
    </double>

    <double>
      <name>timeStep</name>
      <longflag>--timeStep</longflag>
      <description>Time Step</description>
      <label>Time Step</label>
      <default>0.0625</default>
      <constraints>
        <minimum>.001</minimum>
        <maximum>1</maximum>
        <step>.001</step>
      </constraints>
    </double>

    <integer>
      <name>numberOfIterations</name>
      <longflag>--iterations</longflag>
      <description>Number of iterations</description>
      <label>Iterations</label>
      <default>1</default>
      <constraints>
        <minimum>1</minimum>
        <maximum>30</maximum>
        <step>1</step>
      </constraints>
    </integer>

  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>
    <image>
      <name>inputVolume</name>
      <label>Input Volume</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input volume to be filtered</description>
    </image>
    <image>
      <name>outputVolume</name>
      <label>Output Volume</label>
      <channel>output</channel>
      <index>1</index>
      <description>Output filtered volume</description>
    </image>
  </parameters>

</executable>
"""
      


def Execute ( inputVolume, outputVolume, conductance=1.0, timeStep=0.0625, iterations=1 ):
    Slicer = __import__ ( "Slicer" );
    slicer = Slicer.slicer
    scene = slicer.MRMLScene
    inputVolume = scene.GetNodeByID(inputVolume)
    outputVolume = scene.GetNodeByID(outputVolume)

    filter = slicer.vtkITKGradientAnisotropicDiffusionImageFilter()
    filter.SetConductanceParameter(conductance)
    filter.SetTimeStep(timeStep)
    filter.SetNumberOfIterations(iterations)
    filter.SetInput(inputVolume.GetImageData())
    filter.Update()
    outputVolume.SetAndObserveImageData(filter.GetOutput())
    filter.SetOutput()
    matrix = slicer.vtkMatrix4x4()
    inputVolume.GetIJKToRASMatrix(matrix)
    outputVolume.SetIJKToRASMatrix(matrix)

    return

ResampleVolume.py

This module again shows how to handle reference systems correctly: here we get a volume in a vtkImageData, resample it (note where image spacing is fetched from and how it is handled), and put it back in the scene.

XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Converters</category>
  <title>Python Resample Volume</title>
  <description>
Resample a volume based on user-specified parameters.
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>
<acknowledgements>
This work is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149.
</acknowledgements>

  <parameters>
    <label>Resample Parameters</label>
    <description>Parameters for resampling</description>

    <string-enumeration>
      <name>interpolation</name>
      <longflag>interpolation</longflag>
      <description>Interpolation mode</description>
      <label>Interpolation Mode</label>
      <default>linear</default>
      <element>nearest neighbor</element>
      <element>linear</element>
      <element>cubic</element>
    </string-enumeration>

    <double-vector>
      <name>outputSpacing</name>
      <longflag>outputSpacing</longflag>
      <description>Spacing of the output image (in IJK axis ordering)</description>
      <label>Output spacing</label>
      <default>1,1,1</default>
    </double-vector>

  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>

    <image>
      <name>inputVolume</name>
      <label>Input Volume</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input volume to be resampled</description>
    </image>

    <image>
      <name>outputVolume</name>
      <label>Output Volume</label>
      <channel>output</channel>
      <index>1</index>
      <description>Output resampled volume</description>
    </image>

  </parameters>

</executable>
"""


def Execute (inputVolume, outputVolume, interpolation="linear", outputSpacing=[1.0,1.0,1.0]):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene

    inputVolume = scene.GetNodeByID(inputVolume)
    outputVolume = scene.GetNodeByID(outputVolume)

    ijkToRASMatrix = slicer.vtkMatrix4x4()
    inputVolume.GetIJKToRASMatrix(ijkToRASMatrix)

    inputSpacing = inputVolume.GetSpacing()
    spacingRatio = [outputSpacing[0]/inputSpacing[0], 
                    outputSpacing[1]/inputSpacing[1], 
                    outputSpacing[2]/inputSpacing[2]]

    reslice = slicer.vtkImageReslice()
    reslice.SetInput(inputVolume.GetImageData())
    if interpolation == "nearest neighbor":
        reslice.SetInterpolationModeToNearestNeighbor()
    elif interpolation == "linear":
        reslice.SetInterpolationModeToLinear()
    elif interpolation == "cubic":
        reslice.SetInterpolationModeToCubic()
    reslice.SetOutputSpacing(*spacingRatio) 
    reslice.Update()

    changeInformation = slicer.vtkImageChangeInformation()
    changeInformation.SetInput(reslice.GetOutput())
    changeInformation.SetOutputOrigin(0.0,0.0,0.0)
    changeInformation.SetOutputSpacing(1.0,1.0,1.0)
    changeInformation.Update()

    outputVolume.SetAndObserveImageData(changeInformation.GetOutput())
    outputVolume.SetIJKToRASMatrix(ijkToRASMatrix)
    outputVolume.SetSpacing(*outputSpacing)

    return

ResliceAsVolume.py

A variation on the same theme, this module takes a volume and a reference volume and resamples the first on the image grid of the second. Recently this same functionality has been provided in a C++ CLI module, so this one has been relegated here. Again, see how this functionality is achieved using non-RAS aware VTK filters by manually shuffling transform matrices.

 XML = """<?xml version="1.0" encoding="utf-8"?>
<executable>

  <category>Converters</category>
  <title>Python Reslice As Volume</title>
  <description>
Reslice a volume based on the geometry of another volume.
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Luca Antiga and Daniel Blezek</contributor>
<acknowledgements>
This work is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149.
</acknowledgements>

  <parameters>
    <label>Reslice Parameters</label>
    <description>Parameters for the reslice filter</description>

    <string-enumeration>
      <name>interpolation</name>
      <longflag>interpolation</longflag>
      <description>Interpolation mode</description>
      <label>Interpolation Mode</label>
      <default>linear</default>
      <element>nearest neighbor</element>
      <element>linear</element>
      <element>cubic</element>
    </string-enumeration>

    <double>
      <name>background</name>
      <longflag>background</longflag>
      <description>Background level for "outside" values</description>
      <label>Background</label>
      <default>0.0</default>
    </double>

    <boolean>
      <name>autoCrop</name>
      <longflag>autoCrop</longflag>
      <description>Toggle cropping output to ensure that all the input data is contained in the output</description>
      <label>Auto Crop</label>
      <default>false</default>
    </boolean>

  </parameters>

  <parameters>
    <label>IO</label>
    <description>Input/output parameters</description>
    <image>
      <name>inputVolume</name>
      <label>Input Volume</label>
      <channel>input</channel>
      <index>0</index>
      <description>Input volume to be resampled</description>
    </image>
    <image>
      <name>referenceVolume</name>
      <label>Reference Volume</label>
      <channel>input</channel>
      <index>1</index>
      <description>Reference volume defining the geometry</description>
    </image>
    <image>
      <name>outputVolume</name>
      <label>Output Volume</label>
      <channel>output</channel>
      <index>2</index>
      <description>Output resampled volume</description>
    </image>
  </parameters>

</executable>
"""


def Execute (inputVolume, referenceVolume, outputVolume, interpolation="linear", background=0.0, autoCrop=False):

    Slicer = __import__("Slicer")
    slicer = Slicer.slicer
    scene = slicer.MRMLScene

    inputVolume = scene.GetNodeByID(inputVolume)
    referenceVolume = scene.GetNodeByID(referenceVolume)
    outputVolume = scene.GetNodeByID(outputVolume)

    inputIJKToRASMatrix = slicer.vtkMatrix4x4()
    inputVolume.GetIJKToRASMatrix(inputIJKToRASMatrix)
    
    referenceIJKToRASMatrix = slicer.vtkMatrix4x4()
    referenceVolume.GetIJKToRASMatrix(referenceIJKToRASMatrix)

    resliceMatrix = slicer.vtkMatrix4x4()
    inputIJKToRASMatrix.Invert()
    resliceMatrix.Multiply4x4(inputIJKToRASMatrix,referenceIJKToRASMatrix,resliceMatrix)

    resliceTransform = slicer.vtkTransform()
    resliceTransform.SetMatrix(resliceMatrix)

    reslice = slicer.vtkImageReslice()
    reslice.SetInput(inputVolume.GetImageData())
    reslice.SetResliceTransform(resliceTransform)
    if interpolation == "nearest neighbor":
        reslice.SetInterpolationModeToNearestNeighbor()
    elif interpolation == "linear":
        reslice.SetInterpolationModeToLinear()
    elif interpolation == "cubic":
        reslice.SetInterpolationModeToCubic()
    reslice.SetBackgroundLevel(background)
    if autoCrop:
        reslice.AutoCropOutputOn()
    else:
        reslice.AutoCropOutputOff()
        reslice.SetOutputExtent(*referenceVolume.GetImageData().GetWholeExtent())
    reslice.Update()
 
    outputIJKToRASMatrix = slicer.vtkMatrix4x4()

    if autoCrop:
        resliceOrigin = reslice.GetOutput().GetOrigin()
        originMatrix = slicer.vtkMatrix4x4()
        originMatrix.Identity()
        originMatrix.SetElement(0,3,resliceOrigin[0])
        originMatrix.SetElement(1,3,resliceOrigin[1])
        originMatrix.SetElement(2,3,resliceOrigin[2])
        outputIJKToRASMatrix.Multiply4x4(referenceIJKToRASMatrix,originMatrix,outputIJKToRASMatrix)
    else:
        outputIJKToRASMatrix.DeepCopy(referenceIJKToRASMatrix)

    changeInformation = slicer.vtkImageChangeInformation()
    changeInformation.SetInput(reslice.GetOutput())
    changeInformation.SetOutputOrigin(0.0,0.0,0.0)
    changeInformation.SetOutputSpacing(1.0,1.0,1.0)
    changeInformation.Update()

    outputVolume.SetAndObserveImageData(changeInformation.GetOutput())
    outputVolume.SetIJKToRASMatrix(outputIJKToRASMatrix)

    return