dynamicSPECT/resample/resample.py

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

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 = ["LabKey"]
    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)

    datasetCollapsibleButton = ctk.ctkCollapsibleButton()
    datasetCollapsibleButton.text = "Dataset"
    self.layout.addWidget(datasetCollapsibleButton)
    # Layout within the dummy collapsible button
    datasetFormLayout = qt.QFormLayout(datasetCollapsibleButton)

    self.transformNode=qt.QLineEdit("NodeToTransform")
    datasetFormLayout.addRow("TransformedNode:",self.transformNode)
    self.referenceNode=qt.QLineEdit("ReferenceNode")
    datasetFormLayout.addRow("ReferenceNode:",self.referenceNode)

    self.transformButton=qt.QPushButton("Transform")
    self.transformButton.clicked.connect(self.onTransformButtonClicked)
    datasetFormLayout.addRow("Volume:",self.transformButton)

    self.transformSegmentationButton=qt.QPushButton("Transform")
    self.transformSegmentationButton.clicked.connect(self.onTransformSegmentationButtonClicked)
    datasetFormLayout.addRow("Segmentation:",self.transformSegmentationButton)


  def onTransformButtonClicked(self):
    node=slicer.util.getFirstNodeByName(self.transformNode.text)
    if node==None:
        print("Transform node [{}] not found").format(self.transformNode.text)
        return
    refNode=slicer.util.getFirstNodeByName(self.referenceNode.text)
    if refNode==None:
        print("Reference node [{}] not found").format(self.referenceNode.text)
        return

    self.logic.rebinNode(node,refNode)

  def onTransformSegmentationButtonClicked(self):

      segNodes=slicer.util.getNodesByClass("vtkMRMLSegmentationNode")
      segNode=None
      for s in segNodes:
          print ("SegmentationNode: {}").format(s.GetName())
          if s.GetName()==self.transformNode.text:
              segNode=s
              break

      if segNode==None:
          print("Segmentation node [{}] not found").format(self.transformNode.text)
          return

      refNode=slicer.util.getFirstNodeByName(self.referenceNode.text)
      if refNode==None:
          print("Reference node [{}] not found").format(self.referenceNode.text)
          return

      self.logic.rebinSegmentation(segNode,refNode)



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)
      try:
          fhome=os.environ["HOME"]
      except:
             #in windows, the variable is called HOMEPATH
          fhome=os.environ['HOMEDRIVE']+os.environ['HOMEPATH']
      self.baseLog=os.path.join(fhome,".resample")
      if not os.path.isdir(self.baseLog):
          os.mkdir(self.baseLog)

  def printMe(self):
      print("resampleLogic ready")
      print("Log: {}").format(self.baseLog)


  def cast(self,newImage,originalImage):
      if newImage.GetPointData().GetScalars().GetDataType()==originalImage.GetPointData().GetScalars().GetDataType():
          return newImage
      outputType=originalImage.GetPointData().GetScalars().__class__.__name__
      shifter=vtk.vtkImageShiftScale()
      shifter.SetInputData(newImage)
      if outputType=="vtkUnsignedShortArray":
          shifter.SetOutputScalarTypeToUnsignedShort()
      if outputType=="vtkShortArray":
          shifter.SetOutputScalarTypeToShort()
      shifter.SetScale(1)
      shifter.SetShift(0)
      shifter.Update()
      return shifter.GetOutput()



  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())

      log=open(os.path.join(self.baseLog,"rebinNode.log"),"w")


      log.write(("rebinNode: volume: {} ref: {}\n").format(node.GetName(),refNode.GetName()))
      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 tough. COpy only links (ie. shallow copy)
      newImage=vtk.vtkImageData()
      newImage.DeepCopy(refNode.GetImageData())
      newImage=self.cast(newImage,node.GetImageData())
      newScalars=newImage.GetPointData().GetScalars()
      #doesn't set the scalars
      log.write(("Iterating: {} points\n").format(n))


      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)
            if i%10000==0:
                log.write(("[{}]: {}->{}\n").format(i,v0,v))

      #node.SetName(nodeName+"_BU")


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

      newNode.SetAndObserveImageData(newImage)
      slicer.mrmlScene.AddNode(newNode)
      log.write(("Adding node {}\n").format(newNode.GetName()))
      log.close()
      return newNode


  def inMask(self,binaryRep,fpos):

       local=[0,0,0]

       segNode=binaryRep['node']
       segNode.TransformPointFromWorld(fpos,local)

       mask=binaryRep['mask']
       maskWorldToImageMatrix=vtk.vtkMatrix4x4()
       mask.GetWorldToImageMatrix(maskWorldToImageMatrix)
       local.append(1)
       maskIJK=maskWorldToImageMatrix.MultiplyPoint(local)
       #mask IJK is in image coordinates. However, binaryLabelMap is a truncated
       #version of vtkImageData, so more work is required
       maskIJK=maskIJK[0:3]#skip last (dummy) coordinate
       maskIJK=[r*s for r,s in zip(maskIJK,mask.GetSpacing())]
       maskIJK=[r+c for r,c in zip(maskIJK,mask.GetOrigin())]

       maskI=mask.FindPoint(maskIJK)
       try:
           return mask.GetPointData().GetScalars().GetTuple1(maskI)
       except:
           return 0

  def rebinSegment(self,refNode,binaryRep):

       refIJKtoRAS=vtk.vtkMatrix4x4()
       refNode.GetIJKToRASMatrix(refIJKtoRAS)
       refImage=refNode.GetImageData()

       #create new node for each segment
       newImage=vtk.vtkImageData()
       newImage.DeepCopy(refNode.GetImageData())
       n=newImage.GetNumberOfPoints()
       newScalars=newImage.GetPointData().GetScalars()

       segNode=binaryRep['node']
       mask=binaryRep['mask']

       for j in range(0,n):
           refIJK=refImage.GetPoint(j)
           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)

           v=self.inMask(binaryRep,fWorld)
           #print("[{}]  Setting ({}) to: {}\n").format(j,fWorld,v)
           newScalars.SetTuple1(j,v)

       newNode=slicer.vtkMRMLScalarVolumeNode()
       newNode.SetName(segNode.GetName()+'_'+binaryRep['segId'])
       newNode.SetOrigin(refNode.GetOrigin())
       newNode.SetSpacing(refNode.GetSpacing())
       newNode.SetIJKToRASMatrix(refIJKtoRAS)

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


  def rebinSegmentation(self,segNode,refNode):

       log=open(os.path.join(self.baseLog,"rebinSegmentation.log"),"w")

       log.write(("rebinNode: {} {}\n").format(segNode.GetName(),refNode.GetName()))

       nSeg=segNode.GetSegmentation().GetNumberOfSegments()
       ## DEBUG:
       #nSeg=1
       #n=1000

       for i in range(0,nSeg):
           #segID
           segID=segNode.GetSegmentation().GetNthSegmentID(i)
           log.write(("Parsing segment {}").format(segNode.GetSegmentation.GetNthSegment(i).GetName()))
           binaryRep={'node':segNode,
                      'mask':segNode.GetBinaryLabelmapRepresentation(segID)}
           newNode=self.rebinSegment(refNode,binaryRep)
           log.write(("Adding node: {}").format(newNode.GetName()))

       log.close()

  def rebinSegmentation1(self,segNode,refNode):


     logfile="C:\\Users\\studen\\labkeyCache\\log\\resample.log"
     print("rebinNode: {} {}\n").format(segNode.GetName(),refNode.GetName())
     refImage=refNode.GetImageData()
     refIJKtoRAS=vtk.vtkMatrix4x4()
     refNode.GetIJKToRASMatrix(refIJKtoRAS)
     refRAStoIJK=vtk.vtkMatrix4x4()
     refNode.GetRASToIJKMatrix(refRAStoIJK)

     nSeg=segNode.GetSegmentation().GetNumberOfSegments()
     ## DEBUG:
     nSeg=1

     for i in range(0,nSeg):
         #segID
         segID=segNode.GetSegmentation().GetNthSegmentID(i)
         binaryRep={'node':segNode,
                    'mask': segNode.GetBinaryLabelmapRepresentation(segID)}

         mask=binaryRep['mask']
         #create new node for each segment
         newImage=vtk.vtkImageData()
         newImage.DeepCopy(refNode.GetImageData())
         newScalars=newImage.GetPointData().GetScalars()
         refN=newImage.GetNumberOfPoints()
         for k in range(0,refN):
             newScalars.SetTuple1(k,0)

         maskN=binaryRep['mask'].GetNumberOfPoints()
         maskScalars=mask.GetPointData().GetScalars()
         maskImageToWorldMatrix=vtk.vtkMatrix4x4()
         binaryRep['mask'].GetImageToWorldMatrix(maskImageToWorldMatrix)

         for j in range(0,maskN):

             if maskScalars.GetTuple1(j)==0:
                 continue

             maskIJK=mask.GetPoint(j)
             maskIJK=[r-c for r,c in zip(maskIJK,mask.GetOrigin())]
             maskIJK=[r/s for r,s in zip(maskIJK,mask.GetSpacing())]
             maskIJK.append(1)
             maskPos=maskImageToWorldMatrix.MultiplyPoint(maskIJK)
             maskPos=maskPos[0:3]
             fWorld=[0,0,0]
             #apply segmentation transformation
             segNode.TransformPointToWorld(maskPos,fWorld)

             refPos=[0,0,0]

             #apply potential reference transformations
             refNode.TransformPointFromWorld(fWorld,refPos)
             refPos.append(1)
             refIJK=refRAStoIJK.MultiplyPoint(refPos)
             refIJK=refIJK[0:3]
             i1=newImage.FindPoint(refIJK)
             if i1<0:
                 continue
             if i1<refN:
                 newScalars.SetTuple1(i1,1)

         newNode=slicer.vtkMRMLScalarVolumeNode()
         newNode.SetName(segNode.GetName()+'_'+segNode.GetSegmentation().GetNthSegmentID(i))
         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!')