dynamicSPECT/cardiacSPECT/resample.py

import slicer
import vtk
import os
from slicer.ScriptedLoadableModule import *

class resample(ScriptedLoadableModule):
  """Uses ScriptedLoadableModule base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """
  def __init__(self, parent):
    ScriptedLoadableModule.__init__(self, parent)
    parent.title = "Resample"
    parent.categories = ["Examples"]
    parent.dependencies = []
    parent.contributors = ["Andrej Studen (FMF/JSI)"] # replace with "Firstname Lastname (Org)"
    parent.helpText = """
    Resample to different shapes
    """
    parent.acknowledgementText = """
    This module was developed within the frame of the ARRS sponsored medical
    physics research programe to investigate quantitative measurements of cardiac
    function using sestamibi-like tracers
    """ # replace with organization, grant and thanks.
    self.parent.helpText += self.getDefaultModuleDocumentationLink()
    self.parent = parent

#
# resampleWidget
#
class resampleWidget(ScriptedLoadableModuleWidget):
  """Uses ScriptedLoadableModuleWidget base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setup(self):
    ScriptedLoadableModuleWidget.setup(self)
    self.logic=resampleLogic(self)


class resampleLogic(ScriptedLoadableModuleLogic):
  """This class should implement all the actual
  computation done by your module.  The interface
  should be such that other python code can import
  this class and make use of the functionality without
  requiring an instance of the Widget.
  Uses ScriptedLoadableModuleLogic base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """
  def __init__(self,parent):
      ScriptedLoadableModuleLogic.__init__(self, parent)

  def printMe(self):
      print("resampleLogic ready")

  def rebinNode(self,node,refNode):
      #refNodeName="2SBMIRVolume19"
      #nodeName="2SMobrVolume19"
      #node=slicer.util.getFirstNodeByName(nodeName)
      #refNode=slicer.util.getFirstNodeByName(refNodeName)
      #transformNodeName="2SMobr_DF"
      #transformNode=slicer.util.getFirstNodeByName(transformNodeName)
      #node.SetAndObserveTransformNodeID(transformNode.GetID())

      print("rebinNode: {} {}").format(node,refNode)
      refImage=refNode.GetImageData()
      n=refImage.GetNumberOfPoints()
      refIJKtoRAS=vtk.vtkMatrix4x4()
      refNode.GetIJKToRASMatrix(refIJKtoRAS)

      nodeRAStoIJK=vtk.vtkMatrix4x4()
      node.GetRASToIJKMatrix(nodeRAStoIJK)
      nodeName=node.GetName()

      coeff=vtk.vtkImageBSplineCoefficients()
      coeff.SetInputData(node.GetImageData())
      coeff.SetBorderMode(vtk.VTK_IMAGE_BORDER_CLAMP)
      #between 3 and 5
      coeff.SetSplineDegree(5)
      coeff.Update()
      #interpolation ready to use

      #this is tought. COpy only links (ie. shallow copy)
      newImage=vtk.vtkImageData()
      newImage.DeepCopy(refNode.GetImageData())
      newScalars=newImage.GetPointData().GetScalars()
      #doesn't set the scalars


      for i in range(0,n):
            refIJK=refImage.GetPoint(i)
            refIJK=list(refIJK)
            refIJK.append(1)
            #shift to world coordinates to match global points
            refPos=refIJKtoRAS.MultiplyPoint(refIJK)
            refPos=refPos[0:3]
            fWorld=[0,0,0]
            #apply potential transformations
            refNode.TransformPointToWorld(refPos,fWorld)

            #now do the opposite on the target node; reuse fPos
            nodePos=[0,0,0]
            node.TransformPointFromWorld(fWorld,nodePos)
            nodePos.append(1)
            nodeIJK=nodeRAStoIJK.MultiplyPoint(nodePos)

            #here we should apply some sort of interpolation
            nodeIJK=nodeIJK[0:3]
            v=coeff.Evaluate(nodeIJK)
            v0=newScalars.GetTuple1(i)
            newScalars.SetTuple1(i,v)
            #print("[{}]: {}->{}").format(i,v0,v)

      node.SetName(nodeName+"_BU")


      newNode=slicer.vtkMRMLScalarVolumeNode()
      newNode.SetName(nodeName)
      newNode.SetOrigin(refNode.GetOrigin())
      newNode.SetSpacing(refNode.GetSpacing())
      newNode.SetIJKToRASMatrix(refIJKtoRAS)

      newNode.SetAndObserveImageData(newImage)
      slicer.mrmlScene.AddNode(newNode)



class resampleTest(ScriptedLoadableModuleTest):
  """
  This is the test case for your scripted module.
  Uses ScriptedLoadableModuleTest base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setUp(self):
    """ Do whatever is needed to reset the state - typically a scene clear will be enough.
    """
    slicer.mrmlScene.Clear(0)
    refNodeName="2SBMIRVolume19"
    nodeName="2SMobrVolume19"
    transformNodeName="2SMobr_DF"
    path="c:\\Users\\studen\\labkeyCache\\dinamic_spect\\Patients\\@files"

    refPath=os.path.join(path,"2SBMIR")
    refPath=os.path.join(refPath,refNodeName+".nrrd")
    slicer.util.loadNodeFromFile(refPath,'VolumeFile')

    transformPath=os.path.join(path,"2SMobr")
    transformPath=os.path.join(transformPath,transformNodeName+".h5")
    slicer.util.loadNodeFromFile(transformPath,'TransformFile')

    nodePath=os.path.join(path,"2SMobr")
    nodePath=os.path.join(nodePath,nodeName+".nrrd")

    slicer.util.loadNodeFromFile(nodePath,'VolumeFile')

    self.node=slicer.util.getFirstNodeByName(nodeName)
    self.refNode=slicer.util.getFirstNodeByName(refNodeName)
    self.transformNode=slicer.util.getFirstNodeByName(transformNodeName)
    self.node.SetAndObserveTransformNodeID(self.transformNode.GetID())

    self.resampler=resampleLogic(None)

  def runTest(self):
    """Run as few or as many tests as needed here.
    """
    self.setUp()
    self.test_resample()

  def test_resample(self):
    """ Ideally you should have several levels of tests.  At the lowest level
    tests should exercise the functionality of the logic with different inputs
    (both valid and invalid).  At higher levels your tests should emulate the
    way the user would interact with your code and confirm that it still works
    the way you intended.
    One of the most important features of the tests is that it should alert other
    developers when their changes will have an impact on the behavior of your
    module.  For example, if a developer removes a feature that you depend on,
    your test should break so they know that the feature is needed.
    """

    self.delayDisplay("Starting the test")
    #
    # first, get some data
    #
    self.resampler.rebinNode(self.node,self.refNode)

    self.delayDisplay('Test passed!')