Slicer3:Python:DemianExamples

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