Difference between revisions of "Slicer3:Python:DemianExamples"

From Slicer Wiki
Jump to: navigation, search
(New page: <pre> XML = """<?xml version="1.0" encoding="utf-8"?> <executable> <category>Demo Scripted Modules</category> <title>Masked median filtering</title> <description> Perform median fil...)
 
m (Protected "Slicer3:Python:DemianExamples" ([edit=autoconfirmed] (indefinite) [move=autoconfirmed] (indefinite)))
 
(49 intermediate revisions by 30 users not shown)
Line 1: Line 1:
 +
__TOC__
 +
 +
=Introduction=
 +
 +
Demian Wassermann developed a set of tutorial slides and examples for using python and numpy in Slicer3.
 +
 +
=Median Filter in a Masked Region=
 +
 
<pre>
 
<pre>
 
XML = """<?xml version="1.0" encoding="utf-8"?>
 
XML = """<?xml version="1.0" encoding="utf-8"?>
Line 154: Line 162:
 
   outputFilteredVolumeNode.Modified()
 
   outputFilteredVolumeNode.Modified()
  
−
#def maskIJKToImageIJKMatrix( inputVolumeNode, inputMaskVolumeNode ):
+
</pre>
−
#
+
 
−
#  InputTransformNode = inputInputVolumeNode.GetParentTransformNode()
+
=K-Medoids Fiber Clustering=
−
# if InputTransformNode==[]:
+
<pre>
−
#   InputTransformNode=""
+
XML = """<?xml version="1.0" encoding="utf-8"?>
−
#  MaskTransformNode = inputMaskVolumeNode.GetParentTransformNode()
+
<executable>
−
# if MaskTransformNode==[]:
+
 
−
#    MaskTransformNode=""
+
  <category>Demo Scripted Modules</category>
−
#
+
  <title>K-Medoids fiber clustering</title>
−
#  matrixMaskToInputIJK = slicer.vtkMatrix4x4()
+
  <description>
−
#
+
Fiber Clustering simple K-Medoids
−
# if ( MaskTransformNode!="" ):
+
</description>
−
#    MaskTransformNode.GetMatrixTransformToNode( InputTransformNode, matrixMaskToInputIJK )
+
  <version>1.0</version>
−
#  elif ( InputTransformNode!="" ):
+
  <documentation-url></documentation-url>
−
#    InputTransformNode.GetMatrixTransformToNode( MaskTransformNode, matrixMaskToInputIJK )
+
  <license></license>
−
#    matrixMaskToInputIJK.Invert()
+
  <contributor>Demian Wassermann</contributor>
−
#
+
 
−
# 
+
  <parameters>
−
# 
+
  <label>IO</label>
−
# 
+
    <description>Input/output parameters</description>
 +
 
 +
    <geometry type = "fiberbundle" >
 +
      <name>inputFiberBundle</name>
 +
      <longflag>inputFiberBundle</longflag>
 +
      <label>Input Fiber Bundle</label>
 +
      <channel>input</channel>
 +
      <description>Input bundle</description>
 +
    </geometry>
 +
 
 +
    <geometry >
 +
      <name>outputFiberBundle</name>
 +
      <longflag>outputFiberBundle</longflag>
 +
      <label>Output Fiber Bundle</label>
 +
      <channel>output</channel>
 +
      <description>Clustered bundle</description>
 +
    </geometry>
 +
 
 +
    <integer>
 +
      <name>numberOfClusters</name>
 +
      <longflag>numberOfClusters</longflag>
 +
      <label>Number of clusters for K-Medoids</label>
 +
      <default>5</default>
 +
      <step>1</step>
 +
      <channel>input</channel>
 +
      <constraints>
 +
        <minimum>2</minimum>
 +
        <maximum>100</maximum>
 +
      </constraints>
 +
    </integer>
 +
 
 +
  </parameters>
 +
  <parameters advanced="true">
 +
    <label>Advanced</label>
 +
    <integer>
 +
      <name>subsampling</name>
 +
      <longflag>subsampling</longflag>
 +
      <label>Number of fiber points to keep</label>
 +
      <default>15</default>
 +
      <step>1</step>
 +
      <channel>input</channel>
 +
      <constraints>
 +
        <minimum>2</minimum>
 +
        <maximum>1000</maximum>
 +
      </constraints>
 +
    </integer>
 +
    <integer>
 +
      <name>minimumFiberLength</name>
 +
      <longflag>minimumFiberLength</longflag>
 +
      <label>minimum fiber length to consider valid</label>
 +
      <default>15</default>
 +
      <step>1</step>
 +
      <channel>input</channel>
 +
      <constraints>
 +
        <minimum>2</minimum>
 +
        <maximum>1000</maximum>
 +
      </constraints>
 +
    </integer>
 +
  </parameters>
 +
 
 +
</executable>
 +
"""
 +
 
 +
from Slicer import slicer
 +
import numpy
 +
 
 +
# Warning, this example needs the package Pycluster
 +
# http://bonsai.ims.u-tokyo.ac.jp/~mdehoon/software/cluster/software.htm#pycluster
 +
import Pycluster
 +
 
 +
 
 +
 
 +
 
 +
def Execute (inputFiberBundle="", outputFiberBundle="", numberOfClusters=2, subsampling=15, minimumFiberLength=15 ):
 +
 
 +
  scene = slicer.MRMLScene
 +
 
 +
  inputFiberBundleNode = scene.GetNodeByID(inputFiberBundle)
 +
  outputFiberBundleNode = scene.GetNodeByID(outputFiberBundle)
 +
 
 +
 
 +
  #Prepare the output fiber bundle and the Arrays for the atlas labeling and clustering
 +
 
 +
  clusters = setupTheOutputNode( inputFiberBundleNode, outputFiberBundleNode )
 +
  clusters_array = clusters.ToArray().squeeze()
 +
 
 +
  #Get the fibers form the Polydata and susbsample them
 +
  fibers, lines = fibers_from_vtkPolyData( inputFiberBundleNode.GetPolyData(), minimumFiberLength )
 +
 
 +
  subsampledFibers = []
 +
 
 +
  for fiber in fibers:
 +
    subsampledFibers.append( fiber[::max( len(fiber)/subsampling, len(fiber) ) ] )
 +
 
 +
 
 +
  #Generate the distance matrix
 +
  distanceMatrix = numpy.zeros( (len(fibers),len(fibers)), dtype=float )
 +
  for i in xrange( len(fibers) ):
 +
    for j in xrange( 0, i):
 +
      distanceMatrix[ i, j ] = dist_hausdorff_min( subsampledFibers[i], subsampledFibers[j] )
 +
      distanceMatrix[ j, i ] = distanceMatrix[ i, j ]
 +
 
 +
 
 +
 
 +
 
 +
  #Perform the clustering
 +
  fiberClusters = renumberLabels(Pycluster.kmedoids( distanceMatrix, numberOfClusters, npass=100 )[0])
 +
  print fiberClusters
 +
  clusters_array[:]=0
 +
 
 +
 
 +
  for i in xrange(len(lines)):
 +
    clusters_array[ lines[i] ] = fiberClusters[i]
 +
 
 +
  clusters.Modified()
 +
 
 +
 
 +
dist2 = lambda i,j : numpy.sqrt(((i-j)**2).sum(j.ndim-1))
 +
dist_hausdorff_asym_mean = lambda i,j: numpy.apply_along_axis( lambda k: dist2(k,j).min(), 1,i).mean()
 +
dist_hausdorff_min = lambda i,j : numpy.min(dist_hausdorff_asym_mean(i,j),dist_hausdorff_asym_mean(j,i))
 +
 
 +
 
 +
 
 +
def fibers_from_vtkPolyData(vtkPolyData, minimumFiberLength):
 +
    #Fibers and Lines are the same thing
 +
 
 +
    lines = vtkPolyData.GetLines().GetData().ToArray().squeeze()
 +
    points = vtkPolyData.GetPoints().GetData().ToArray()
 +
 
 +
    fibersList = []
 +
    linesList = []
 +
    actualLineIndex = 0
 +
    numberOfFibers = vtkPolyData.GetLines().GetNumberOfCells()
 +
    for l in xrange( numberOfFibers ):
 +
      if lines[actualLineIndex]>minimumFiberLength:
 +
        fibersList.append( points[ lines[actualLineIndex+1: actualLineIndex+lines[actualLineIndex]+1] ] )
 +
        linesList.append( lines[actualLineIndex+1: actualLineIndex+lines[actualLineIndex]+1] )
 +
      actualLineIndex += lines[actualLineIndex]+1
 +
 
 +
    return fibersList, linesList
 +
 
 +
def setupTheOutputNode( inputFiberBundleNode, outputFiberBundleNode ):
 +
  if ( outputFiberBundleNode.GetPolyData()==[] ):
 +
    outputFiberBundleNode.SetAndObservePolyData(slicer.vtkPolyData())
 +
 
 +
  outputPolyData = outputFiberBundleNode.GetPolyData()
 +
  outputPolyData.SetPoints( inputFiberBundleNode.GetPolyData().GetPoints() )
 +
  outputPolyData.SetLines( inputFiberBundleNode.GetPolyData().GetLines() )
 +
  outputPolyData.Update()
 +
 
 +
 
 +
  clusters = outputFiberBundleNode.GetPolyData().GetPointData().GetScalars('Cluster')
 +
  if (clusters==[] or clusters.GetNumberOfTuples()!=outputPolyData.GetPoints().GetNumberOfPoints() ):
 +
    clusters = slicer.vtkUnsignedIntArray()
 +
    clusters.SetNumberOfComponents(1)
 +
    clusters.SetNumberOfTuples( outputPolyData.GetPoints().GetNumberOfPoints() )
 +
    clusters.SetName('Cluster')
 +
    outputPolyData.GetPointData().AddArray( clusters )
 +
 
 +
  return clusters
 +
 
 +
def renumberLabels(labelArray):
 +
  newLabeling=[]
 +
  for a in labelArray:
 +
    if not(a in newLabeling):
 +
      newLabeling.append(a)
 +
 
 +
  newLabelArray=labelArray.copy()
 +
  for i in range(len(labelArray)):
 +
    newLabelArray[i]=newLabeling.index(labelArray[i])+1
  
 +
  return newLabelArray
 
</pre>
 
</pre>

Latest revision as of 16:41, 27 September 2009

Home < Slicer3:Python:DemianExamples

Introduction

Demian Wassermann developed a set of tutorial slides and examples for using python and numpy in Slicer3.

Median Filter in a Masked Region

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

  <category>Demo Scripted Modules</category>
  <title>Masked median filtering</title>
  <description>
Perform median filtering over a masked section of an image
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Demian Wassermann</contributor>

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

    <image type = "scalar" >
      <name>inputVolume</name>
      <longflag>inputVolume</longflag>
      <label>Input Image</label>
      <channel>input</channel>
      <description>Input image to be filtered</description>
    </image>

    <integer>
      <name>medianFilterRadius</name>
      <longflag>medianFilterRadius</longflag>
      <label>Radius of the median filter</label>
      <default>2</default>
      <step>1</step>
      <channel>input</channel>
      <constraints>
        <minimum>2</minimum>
        <maximum>100</maximum>
      </constraints>
    </integer>


    <image type="label">
      <name>inputMaskVolume</name>
      <longflag>inputMaskVolume</longflag>
      <label>Input Mask Volume</label>
      <channel>input</channel>
      <description>Input mask to work on it</description>
    </image>

    <integer>
      <name>labelToUse</name>
      <longflag>labelToUse</longflag>
      <label>Label to use for the mask</label>
      <default>1</default>
      <step>1</step>
      <channel>input</channel>
      <constraints>
        <minimum>0</minimum>
        <maximum>255</maximum>
      </constraints>
    </integer>

    <image type = "scalar">
      <name>outputFilteredVolume</name>
      <longflag>outputFilteredVolume</longflag>
      <label>Output Image</label>
      <channel>output</channel>
      <description>Image that was median filtered</description>
    </image>

  </parameters>

</executable>
"""


from Slicer import slicer
from scipy import ndimage
import numpy

def Execute(\
     inputVolume = "",\
     medianFilterRadius = 0,\
     inputMaskVolume = "",\
     labelToUse = 1,\
     outputFilteredVolume = ""\
     ):


#       Set up the slicer environment
  scene = slicer.MRMLScene

#       Get the nodes from the MRML tree

  inputVolumeNode = scene.GetNodeByID( inputVolume )
  inputMaskVolumeNode = scene.GetNodeByID( inputMaskVolume )
  outputFilteredVolumeNode = scene.GetNodeByID( outputFilteredVolume )

#       Set up the output node
  setupTheOutputNode( inputVolumeNode, outputFilteredVolumeNode )

  #maskToImageIJK = maskIJKToImageIJKMatrix( inputVolumeNode, inputMaskVolumeNode )

#       Do what we are here to do
  input_array = inputVolumeNode.GetImageData().ToArray()
  output_array = outputFilteredVolumeNode.GetImageData().ToArray()
  mask_array = inputMaskVolumeNode.GetImageData().ToArray()

  pointsToProcess = numpy.transpose(numpy.where( mask_array == labelToUse ))
  print pointsToProcess

  # Get the bounding box

  minCorner = pointsToProcess.min(0)   # yeah yeah there's a border problem
  maxCorner = pointsToProcess.max(0)+1

  # Perform the filtering
  medianFiltered = ndimage.median_filter(\
      input_array[\
      minCorner[0]:maxCorner[0],\
      minCorner[1]:maxCorner[1],\
      minCorner[2]:maxCorner[2],\
      ], medianFilterRadius )


# Set it to the output node
  output_array[:]=input_array
  output_array[ tuple(pointsToProcess.T) ] = medianFiltered[ tuple((pointsToProcess-minCorner).T) ]
  outputFilteredVolumeNode.Modified()




def setupTheOutputNode( inputVolumeNode, outputFilteredVolumeNode ):

  inputVolumeNode_imageData = inputVolumeNode.GetImageData()
  outputFilteredVolume_ImageData = outputFilteredVolumeNode.GetImageData()

  if not outputFilteredVolume_ImageData:
    outputFilteredVolume_ImageData = slicer.vtkImageData()
    outputFilteredVolumeNode.SetAndObserveImageData( outputFilteredVolume_ImageData )

  dimensions = inputVolumeNode_imageData.GetDimensions()
  outputFilteredVolume_ImageData.SetDimensions( dimensions[0], dimensions[1], dimensions[2] )
  outputFilteredVolume_ImageData.SetScalarType( inputVolumeNode_imageData.GetScalarType() )
  outputFilteredVolume_ImageData.SetOrigin( 0, 0, 0 )
  outputFilteredVolume_ImageData.SetSpacing( 1, 1, 1 )
  outputFilteredVolume_ImageData.AllocateScalars()

  matrix = slicer.vtkMatrix4x4()
  inputVolumeNode.GetIJKToRASMatrix( matrix )
  outputFilteredVolumeNode.SetIJKToRASMatrix( matrix )


  outputFilteredVolumeNode.Modified()

K-Medoids Fiber Clustering

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

  <category>Demo Scripted Modules</category>
  <title>K-Medoids fiber clustering</title>
  <description>
Fiber Clustering simple K-Medoids 
</description>
  <version>1.0</version>
  <documentation-url></documentation-url>
  <license></license>
  <contributor>Demian Wassermann</contributor>

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

    <geometry type = "fiberbundle" >
      <name>inputFiberBundle</name>
      <longflag>inputFiberBundle</longflag>
      <label>Input Fiber Bundle</label>
      <channel>input</channel>
      <description>Input bundle</description>
    </geometry>

    <geometry >
      <name>outputFiberBundle</name>
      <longflag>outputFiberBundle</longflag>
      <label>Output Fiber Bundle</label>
      <channel>output</channel>
      <description>Clustered bundle</description>
    </geometry>

    <integer>
      <name>numberOfClusters</name>
      <longflag>numberOfClusters</longflag>
      <label>Number of clusters for K-Medoids</label>
      <default>5</default>
      <step>1</step>
      <channel>input</channel>
      <constraints>
        <minimum>2</minimum>
        <maximum>100</maximum>
      </constraints>
    </integer>

  </parameters>
  <parameters advanced="true">
    <label>Advanced</label>
    <integer>
      <name>subsampling</name>
      <longflag>subsampling</longflag>
      <label>Number of fiber points to keep</label>
      <default>15</default>
      <step>1</step>
      <channel>input</channel>
      <constraints>
        <minimum>2</minimum>
        <maximum>1000</maximum>
      </constraints>
    </integer>
    <integer>
      <name>minimumFiberLength</name>
      <longflag>minimumFiberLength</longflag>
      <label>minimum fiber length to consider valid</label>
      <default>15</default>
      <step>1</step>
      <channel>input</channel>
      <constraints>
        <minimum>2</minimum>
        <maximum>1000</maximum>
      </constraints>
    </integer>
  </parameters>

</executable>
"""

from Slicer import slicer
import numpy

# Warning, this example needs the package Pycluster 
# http://bonsai.ims.u-tokyo.ac.jp/~mdehoon/software/cluster/software.htm#pycluster
import Pycluster




def Execute (inputFiberBundle="", outputFiberBundle="", numberOfClusters=2, subsampling=15, minimumFiberLength=15 ):
  
  scene = slicer.MRMLScene

  inputFiberBundleNode = scene.GetNodeByID(inputFiberBundle)
  outputFiberBundleNode = scene.GetNodeByID(outputFiberBundle)


  #Prepare the output fiber bundle and the Arrays for the atlas labeling and clustering

  clusters = setupTheOutputNode( inputFiberBundleNode, outputFiberBundleNode )
  clusters_array = clusters.ToArray().squeeze()

  #Get the fibers form the Polydata and susbsample them
  fibers, lines = fibers_from_vtkPolyData( inputFiberBundleNode.GetPolyData(), minimumFiberLength )

  subsampledFibers = []

  for fiber in fibers:
    subsampledFibers.append( fiber[::max( len(fiber)/subsampling, len(fiber) ) ] )


  #Generate the distance matrix
  distanceMatrix = numpy.zeros( (len(fibers),len(fibers)), dtype=float )
  for i in xrange( len(fibers) ):
    for j in xrange( 0, i):
      distanceMatrix[ i, j ] = dist_hausdorff_min( subsampledFibers[i], subsampledFibers[j] )
      distanceMatrix[ j, i ] = distanceMatrix[ i, j ] 




  #Perform the clustering
  fiberClusters = renumberLabels(Pycluster.kmedoids( distanceMatrix, numberOfClusters, npass=100 )[0])
  print fiberClusters 
  clusters_array[:]=0


  for i in xrange(len(lines)):
    clusters_array[ lines[i] ] = fiberClusters[i]

  clusters.Modified()


dist2 = lambda i,j : numpy.sqrt(((i-j)**2).sum(j.ndim-1))
dist_hausdorff_asym_mean = lambda i,j: numpy.apply_along_axis( lambda k: dist2(k,j).min(),  1,i).mean()
dist_hausdorff_min = lambda i,j : numpy.min(dist_hausdorff_asym_mean(i,j),dist_hausdorff_asym_mean(j,i))



def fibers_from_vtkPolyData(vtkPolyData, minimumFiberLength):
    #Fibers and Lines are the same thing

    lines = vtkPolyData.GetLines().GetData().ToArray().squeeze()
    points = vtkPolyData.GetPoints().GetData().ToArray()

    fibersList = []
    linesList = []
    actualLineIndex = 0
    numberOfFibers = vtkPolyData.GetLines().GetNumberOfCells()
    for l in xrange( numberOfFibers ):
      if lines[actualLineIndex]>minimumFiberLength:
        fibersList.append( points[ lines[actualLineIndex+1: actualLineIndex+lines[actualLineIndex]+1] ] )
        linesList.append( lines[actualLineIndex+1: actualLineIndex+lines[actualLineIndex]+1]  )
      actualLineIndex += lines[actualLineIndex]+1

    return fibersList, linesList

def setupTheOutputNode( inputFiberBundleNode, outputFiberBundleNode ):
  if ( outputFiberBundleNode.GetPolyData()==[] ):
    outputFiberBundleNode.SetAndObservePolyData(slicer.vtkPolyData())

  outputPolyData = outputFiberBundleNode.GetPolyData()
  outputPolyData.SetPoints( inputFiberBundleNode.GetPolyData().GetPoints() )
  outputPolyData.SetLines( inputFiberBundleNode.GetPolyData().GetLines() )
  outputPolyData.Update()


  clusters = outputFiberBundleNode.GetPolyData().GetPointData().GetScalars('Cluster')
  if (clusters==[] or clusters.GetNumberOfTuples()!=outputPolyData.GetPoints().GetNumberOfPoints() ):
    clusters = slicer.vtkUnsignedIntArray()
    clusters.SetNumberOfComponents(1)
    clusters.SetNumberOfTuples( outputPolyData.GetPoints().GetNumberOfPoints() )
    clusters.SetName('Cluster')
    outputPolyData.GetPointData().AddArray( clusters )

  return clusters

def renumberLabels(labelArray):
  newLabeling=[]
  for a in labelArray:
    if not(a in newLabeling):
      newLabeling.append(a)

  newLabelArray=labelArray.copy()
  for i in range(len(labelArray)):
    newLabelArray[i]=newLabeling.index(labelArray[i])+1

  return newLabelArray
Retrieved from "https://www.slicer.org/w/index.php?title=Slicer3:Python:DemianExamples&oldid=10785"