Adding resampler

cardiacSPECT/cardiacSPECT.py

@@ -10,6 +10,7 @@ import fileIO
 import slicer
 import numpy as np
 import slicerNetwork
+import resample
 #
 # cardiacSPECT
 #
@@ -97,16 +98,35 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     self.patientId=qt.QLineEdit();
     dataFormLayout.addRow('Patient ID', self.patientId)
 
+    self.refPatientId=qt.QLineEdit();
+    dataFormLayout.addRow('Reference Patient ID', self.refPatientId)
+
+
     patientLoadButton = qt.QPushButton("Load")
     patientLoadButton.toolTip="Load data from DICOM"
     dataFormLayout.addRow("Load Patient",patientLoadButton)
     patientLoadButton.clicked.connect(self.onPatientLoadButtonClicked)
 
-    saveVolumeButton = qt.QPushButton("Save Volume to NRRD")
-    saveVolumeButton.toolTip="Save to NRRD"
-    dataFormLayout.addRow("Save to NRRD",saveVolumeButton)
+    loadSegmentationButton = qt.QPushButton("Load")
+    loadSegmentationButton.toolTip="Load segmentation from server"
+    dataFormLayout.addRow("Load Segmentation",loadSegmentationButton)
+    loadSegmentationButton.clicked.connect(self.onLoadSegmentationButtonClicked)
+
+    saveVolumeButton = qt.QPushButton("Save")
+    saveVolumeButton.toolTip="Save volume to NRRD"
+    dataFormLayout.addRow("Volume",saveVolumeButton)
     saveVolumeButton.clicked.connect(self.onSaveVolumeButtonClicked)
 
+    saveSegmentationButton = qt.QPushButton("Save")
+    saveSegmentationButton.toolTip="Save segmentation to NRRD"
+    dataFormLayout.addRow("Segmentation",saveSegmentationButton)
+    saveSegmentationButton.clicked.connect(self.onSaveSegmentationButtonClicked)
+
+    saveTransformationButton = qt.QPushButton("Save")
+    saveTransformationButton.toolTip="Save transformation to NRRD"
+    dataFormLayout.addRow("Transformation",saveTransformationButton)
+    saveTransformationButton.clicked.connect(self.onSaveTransformationButtonClicked)
+
     transformNodeButton = qt.QPushButton("Transform Node")
     transformNodeButton.toolTip="Transform node with patient based transform"
     dataFormLayout.addRow("Transform Nodes",transformNodeButton)
@@ -338,12 +358,21 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
   def onPatientLoadButtonClicked(self):
       self.logic.loadPatient(self,self.patientId.text)
 
+  def onLoadSegmentationButtonClicked(self):
+      self.logic.loadSegmentation(self.patientId.text)
+
   def onSaveVolumeButtonClicked(self):
       self.logic.storeVolumeNodes(self.patientId.text,
             self.time_frame_select.minimum,self.time_frame_select.maximum)
 
+  def onSaveSegmentationButtonClicked(self):
+      self.logic.storeSegmentation(self.patientId.text)
+
+  def onSaveTransformationButtonClicked(self):
+      self.logic.storeTransformation(self.patientId.text)
+
   def onTransformNodeButtonClicked(self):
-      self.logic.applyTransform(self.patientId.text,
+      self.logic.applyTransform(self.patientId.text, self.refPatientId.text,
           self.time_frame_select.minimum,self.time_frame_select.maximum)
 
 
@@ -399,11 +428,15 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
         'Slicer Python','Data loaded')
 
 
+  def formatToGetfileSpec(self,path):
+      #path=self.net.GetRelativePathFromLabkeyPath(path)
+      path="labkey://"+path
+      return path
+
   def loadPatient(self,widget,patientId):
       print("Loading {}").format(patientId)
       ds=self.net.loadDataset("dinamic_spect/Patients","Imaging")
       for r in ds['rows']:
-        print r
         if r['aliasID']==patientId:
             visit=r
 
@@ -414,23 +447,23 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
               continue
           if abs(r['SequenceNum']-float(visit['nmMaster']))<0.1:
               path=r['_labkeyurl_Series']
-              path=path[:path.rfind('/')]
-              masterPath="labkey://"+path
+              masterPath=path[:path.rfind('/')]
+              #masterPath="labkey://"+path
           if abs(r['SequenceNum']-float(visit['nmCorrData']))<0.1:
               path=r['_labkeyurl_Series']
-              path=path[:path.rfind('/')]
-              nmPath="labkey://"+path
+              nmPath=path[:path.rfind('/')]
+              #nmPath="labkey://"+path
           if abs(r['SequenceNum']-float(visit['ctData']))<0.1:
               path=r['_labkeyurl_Series']
-              path=path[:path.rfind('/')]
-              ctPath="labkey://"+path
+              ctPath=path[:path.rfind('/')]
+              #ctPath="labkey://"+path
 
 
-      self.pd.readMasterDirectory(masterPath)
+      self.pd.readMasterDirectory(self.formatToGetfileSpec(masterPath))
       self.frame_data, self.frame_time, self.frame_origin, \
-          self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(nmPath)
+          self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(self.formatToGetfileSpec(nmPath))
       self.ct_data,self.ct_origin,self.ct_pixel_size, \
-         self.ct_orientation=self.pd.readCTDirectory(ctPath)
+         self.ct_orientation=self.pd.readCTDirectory(self.formatToGetfileSpec(ctPath))
 
       self.ct_orientation=vi.completeOrientation(self.ct_orientation)
       self.frame_orientation=vi.completeOrientation(self.frame_orientation)
@@ -440,6 +473,16 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
 
       widget.time_frame_select.setMaximum(self.frame_data.shape[3]-1)
 
+  def loadSegmentation(self,patientId):
+      project="dinamic_spect/Patients"
+      relativePath=project+'/@files/'+patientId
+      segName="Heart"
+      labkeyFile=relativePath+'/'+segName+'.nrrd'
+      print ("Remote: {}").format(labkeyFile)
+      self.net.loadNode(labkeyFile,'SegmentationFile')
+
+
+
   def addNode(self,nodeName,v, lpsOrigin, pixel_size, lpsOrientation, dataType):
 
        # if dataType=0 it is CT data, which gets propagated to background an is
@@ -607,12 +650,13 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       return segNode.GetSegmentation().GetSegment(segNode.GetSegmentation().GetNthSegmentID(i)).GetName()
 
   def storeNodeRemote(self,relativePath,nodeName):
-      labkeyPath=self.pd.net.GetLabkeyWebdavUrl()+'/'+relativePath
+      labkeyPath=self.pd.net.GetLabkeyPathFromRelativePath(relativePath)
       print ("Remote: {}").format(labkeyPath)
       #checks if exists
       self.pd.net.mkdir(labkeyPath)
 
-      localPath=self.pd.net.GetLocalCacheDirectory()+'/'+relativePath
+      localPath=self.pd.net.GetLocalPathFromRelativePath(relativePath)
+      localPath.replace('/',os.path.sep)
 
       node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
       if node==None:
@@ -625,7 +669,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       if node.__class__.__name__=="vtkMRMLTransformNode":
           suffix=".h5"
       fileName=nodeName+suffix
-      file=localPath+'/'+fileName
+      file=os.path.join(localPath,fileName)
       slicer.util.saveNode(node,file)
       print("Stored to: {}").format(file)
       f=open(file)
@@ -650,29 +694,34 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
           nodeName=patientId+'Volume'+str(it)
           self.storeNodeRemote(relativePath,nodeName)
 
+  def storeSegmentation(self,patientId):
+      project="dinamic_spect/Patients"
+      relativePath=project+'/@files/'+patientId
       segNodeName="Heart"
       self.storeNodeRemote(relativePath,segNodeName)
 
+  def storeInputFunction(self,patientId):
       doubleArrayNodeName=patientId+'_Ventricle'
       self.storeNodeRemote(relativePath,doubleArrayNodeName)
 
+  def storeTransformation(self,patientId):
       transformNodeName=patientId+"_DF"
       self.storeNodeRemote(relativePath,transformNodeName)
 
-  def applyTransform(self, patientId,n1,n2):
+  def applyTransform(self, patientId,refPatientId,n1,n2):
       transformNodeName=patientId+"_DF"
       transformNode=slicer.util.getFirstNodeByName(transformNodeName)
       if transformNode==None:
           print("Node {} not found").format(transformNodeName)
           return
+
       try:
-          if self.transformApplied:
-              print("Transform already applied")
-              return
+          self.resampler.printMe()
+          return
       except AttributeError:
-          print("Not defined yet")
+          print("Initializing resampler")
+          self.resampler=resample.resampleLogic(None)
 
-      self.transformApplied=True
 
       n=n2-n1+1
       for i in range(n):
@@ -681,18 +730,27 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
           node=slicer.util.getFirstNodeByName(nodeName)
           if node==None:
               continue
-          node.ApplyTransform(transformNode.GetTransformToParent())
+          node.SetAndObserveTransformNodeID(transformNode.GetID())
+
+          refNodeName=refPatientId+'Volume'+str(it)
+          refNode=slicer.util.getFirstNodeByName(refNodeName)
+          if refNode!=None:
+              self.resampler.rebinNode(node,refNode)
 
       nodeName=patientId+'CT'
       node=slicer.util.getFirstNodeByName(nodeName)
       if not node==None:
-          node.ApplyTransform(transformNode.GetTransformToParent())
+          node.SetAndObserveTransformNodeID(transformNode.GetID())
+          refNodeName=refPatientId+'CT'
+          refNode=slicer.util.getFirstNodeByName(refNodeName)
+          if refNode!=None:
+              self.resampler.rebinNode(node,refNode)
 
 
 
 
 
-   def calculateInputFunction(self):
+  def calculateInputFunction(self):
        n=self.frame_data.shape[3]
 
        dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode','Ventricle')

cardiacSPECT/parseDicom.py

@@ -56,13 +56,15 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
             if not ok:
                 print "Error accessing path"
                 return []
+            files=[self.net.GetRelativePathFromLabkeyPath(f) for f in files]
 
-            if mypath.find('file://')==0:
-                print("Using local files")
-                localPath=re.sub('file://','',mypath)
-                files = [os.path.join(localPath,f) for f in os.listdir(localPath)
-                    if os.path.isfile(os.path.join(localPath, f))]
-            return files
+        if mypath.find('file://')==0:
+            print("Using local files")
+            localPath=re.sub('file://','',mypath)
+            files = [os.path.join(localPath,f) for f in os.listdir(localPath)
+                if os.path.isfile(os.path.join(localPath, f))]
+
+        return files
 
     def getfile(self,origin,f):
 
@@ -70,9 +72,8 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
             try:
                 #not sure if labkey is available, but try it
                 print("Using labkey")
-                url=self.net.GetHostName()
-                print("Sever:{0}, file:{1}".format(url,f))
-                return [self.net.readFile(str(url),f),1]
+                print("Server:{0}, file:{1}".format(self.net.GetHostName(),f))
+                return [self.net.readFileToBuffer(f),1]
             except:
                 print('Could not access labkey. Exiting')
                 return ['NULL',0]
@@ -267,10 +268,11 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
     def readMasterDirectory(self,mypath):
         self.axisShift=(2,1,0)
 
+        print("Reading master from {}").format(mypath)
         origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
 
-        onlyfiles=self.filelist(mypath)
-        for f in onlyfiles:
+        relativeFiles=self.filelist(mypath)
+        for f in relativeFiles:
             print '{}:'.format(f)
 
             g,ok=self.getfile(origin,f)

cardiacSPECT/resample.py

@@ -0,0 +1,187 @@
+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!')