Difference between revisions of "Slicer3:Python:DemianExamples"
From Slicer Wiki
Line 1: | Line 1: | ||
__TOC__ | __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> | ||
+ | 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() | ||
+ | |||
+ | </pre> | ||
=K-Medoids Fiber Clustering= | =K-Medoids Fiber Clustering= |
Revision as of 16:40, 27 September 2009
Home < Slicer3:Python:DemianExamplesIntroduction
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