Scripts to automate dynamic cardiacSPECT registration and inputFunction sampling

cardiacRegistration/register.py

@@ -0,0 +1,67 @@
+import slicer
+import cardiacSPECT
+#must have SlicerElastix extension installed
+import Elastix
+import slicerNetwork
+
+mirID='2SBMIR'
+obrID='2SMobr'
+
+configFile=os.path.join(os.path.expanduser('~'),'.cardiacSPECT','register.json')
+
+dataManager=cardiacSPECT.cardiacSPECTLogic(configFile)
+network=slicerNetwork.labkeyURIHandler()
+
+#parseConfig if a different location than default is chosen
+#network.parseConfig(remote.json)
+#network.initRemote()
+
+dataManager.setURIHandler(network)
+
+
+dataManager.loadPatient(mirID)
+dataManager.loadPatient(obrID)
+
+elastix=Elastix.ElastixLogic()
+registrationIndex=12
+
+"""
+    index 12 corresponds to:
+    ['par0015', '3D CT, monomodal', 'lung',
+   'intrapatient; B-spline transformation; several similarity metrics',
+   'Staring (2013), Towards Local Progression Estimation of Pulmonary Emphysema using CT',
+   ['Parameters.Par0015.expA.patient.NC.affine.txt', 'Parameters.Par0015.expA.patient.NC.bspline.txt']]
+ """
+RegistrationPresets_ParameterFilenames=5
+parameterFilenames = elastix.getRegistrationPresets()[registrationIndex][RegistrationPresets_ParameterFilenames]
+
+fixedVolumeNode=slicer.util.getFirstNodeByName(mirID+'CT')
+movingVolumeNode=slicer.util.getFirstNodeByName(obrID+'CT')
+
+transformNodeName=obrID+'_DF'
+transformNodeFile=outPath=os.path.join(os.path.expanduser('~'),'temp',transformNodeName+'.h5')
+ok, transformNode= slicer.util.loadTransform(transformNodeFile,returnNode=True)
+if not ok:
+    transformNode=slicer.vtkMRMLTransformNode()
+    transformNode.SetName(transformNodeName) # a bspline transform node
+    slicer.mrmlScene.AddNode(transformNode)
+
+    elastix.registerVolumes(fixedVolumeNode, movingVolumeNode,
+        parameterFilenames = parameterFilenames,
+        outputVolumeNode = None,
+        outputTransformNode = transformNode,
+        fixedVolumeMaskNode = None,
+        movingVolumeMaskNode = None)
+        #recent elastix version
+        #forceDisplacementFieldOutputTransform = 1)
+
+    slicer.util.saveNode(transformNode,transformNodeFile)
+#already a grid transform node
+#gridTransformNode=slicer.modules.transforms.logic().ConvertToGridTransform(transformNode, fixedVolumeNode)
+
+dataManager.applyTransform(obrID,mirID,0,19)
+dataManager.storeVolumeNodes(obrID,0,19)
+dataManager.storeVolumeNodes(mirID,0,19)
+dataManager.storeTransformation(obrID)
+
+quit()

cardiacRegistration/storeInputFunction.py

@@ -0,0 +1,24 @@
+import cardiacSPECT
+import slicerNetwork
+
+mirID='2SBMIR'
+obrID='2SMobr'
+
+configFile=os.path.join(os.path.expanduser('~'),'.cardiacSPECT','register.json')
+
+dataManager=cardiacSPECT.cardiacSPECTLogic(configFile)
+network=slicerNetwork.labkeyURIHandler()
+
+#parseConfig if a different location than default is chosen
+#network.parseConfig(remote.json)
+#network.initRemote()
+
+dataManager.setURIHandler(network)
+
+dataManager.loadPatientNRRD(mirID)
+dataManager.storeInputFunction(mirID)
+
+dataManager.loadPatientNRRD(obrID)
+dataManager.storeInputFunction(obrID)
+
+quit()

cardiacSPECT/cardiacSPECT.py

@@ -11,6 +11,8 @@ import slicer
 import numpy as np
 import slicerNetwork
 import resample
+import json
+import re
 #
 # cardiacSPECT
 #
@@ -55,7 +57,8 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     except:
         self.network=slicerNetwork.labkeyURIHandler()
 
-    self.logic=cardiacSPECTLogic(self)
+    configFile=os.path.join(os.path.expanduser('~'),'.cardiacSPECT','cardiacSPECT.json')
+    self.logic=cardiacSPECTLogic(configFile)
     self.logic.setURIHandler(self.network)
     self.selectRemote.setMaster(self)
 
@@ -67,34 +70,8 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     # Layout within the sample collapsible button
     dataFormLayout = qt.QFormLayout(dataButton)
 
-    #pathGuess="file://"+os.environ['HOME']+"/SPECT"
-    pathGuess="labkey://" + "dinamic_spect/%40files/Dinamika%20test2/SPECT_Dinamika_Rekonstruirano"
-    self.dataPath=qt.QLineEdit(pathGuess)
-    dataFormLayout.addRow("Data location",self.dataPath)
 
 
-    self.remotePath=qt.QLineEdit();
-    dataFormLayout.addRow('Remote Path', self.remotePath)
-    self.remotePath.textChanged.connect(self.onRemotePathTextChanged)
-
-    browseButton = qt.QPushButton("Browse local")
-    browseButton.toolTip="Set file location"
-    dataFormLayout.addRow("Select local",browseButton)
-    browseButton.connect('clicked(bool)',self.onBrowseButtonClicked)
-
-    browseRemoteButton = qt.QPushButton("Browse remote")
-    browseRemoteButton.toolTip="Set remote location"
-    dataFormLayout.addRow("Select remote",browseRemoteButton)
-    browseRemoteButton.connect('clicked(bool)',self.onRemoteBrowseButtonClicked)
-
-
-    dataLoadButton = qt.QPushButton("Load")
-    dataLoadButton.toolTip="Load data from DICOM"
-    dataFormLayout.addRow("Data",dataLoadButton)
-    dataLoadButton.connect('clicked(bool)',self.onDataLoadButtonClicked)
-
-    self.dataLoadButton = dataLoadButton
-
     self.patientId=qt.QLineEdit();
     dataFormLayout.addRow('Patient ID', self.patientId)
 
@@ -107,6 +84,11 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     dataFormLayout.addRow("Patient",patientLoadButton)
     patientLoadButton.clicked.connect(self.onPatientLoadButtonClicked)
 
+    patientLoadNRRDButton = qt.QPushButton("Load NRRD")
+    patientLoadNRRDButton.toolTip="Load data from NRRD"
+    dataFormLayout.addRow("Patient",patientLoadNRRDButton)
+    patientLoadNRRDButton.clicked.connect(self.onPatientLoadNRRDButtonClicked)
+
     loadSegmentationButton = qt.QPushButton("Load")
     loadSegmentationButton.toolTip="Load segmentation from server"
     dataFormLayout.addRow("Segmentation",loadSegmentationButton)
@@ -358,10 +340,15 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
         cvn.SetChartNodeID(cn.GetID())
 
   def onCountSegmentsClicked(self):
-        self.countSegmentsDisplay.setText(self.logic.countSegments())
+      self.countSegmentsDisplay.setText(self.logic.countSegments())
 
   def onPatientLoadButtonClicked(self):
-      self.logic.loadPatient(self,self.patientId.text)
+      self.logic.loadPatient(self.patientId.text)
+      self.time_frame_select.setMaximum(self.logic.frame_data.shape[3]-1)
+
+  def onPatientLoadNRRDButtonClicked(self):
+      self.logic.loadPatientNRRD(self.patientId.text)
+      self.time_frame_select.setMaximum(len(self.logic.frame_time))
 
   def onLoadSegmentationButtonClicked(self):
       self.logic.loadSegmentation(self.patientId.text)
@@ -404,9 +391,23 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
   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 __init__(self,config):
+      ScriptedLoadableModuleLogic.__init__(self)
       self.pd=parseDicom.parseDicomLogic(self)
+      self.resampler=resample.resampleLogic(None)
+
+      fname=config
+      try:
+
+          f=open(fname)
+      except OSError as e:
+          print "Confgiuration error: OS error({0}): {1}".format(e.errno, e.strerror)
+          return
+
+      self.cfg=json.load(f)
+      self.coreRelativePath=self.cfg["project"]+'/'+self.cfg['atFiles']
+
+
 
   def setURIHandler(self,net):
       self.net=net
@@ -436,12 +437,16 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
         'Slicer Python','Data loaded')
 
 
-  def formatToGetfileSpec(self,path):
-      #path=self.net.GetRelativePathFromLabkeyPath(path)
-      path="labkey://"+path
-      return path
+  def getFilespecPath(self,r):
+      if self.cfg["remote"]==True:
+          path=r['_labkeyurl_Series']
+          path=path[:path.rfind('/')]
+          return "labkey://"+path
+      else:
+          path=os.path.join(self.cfg["dicomPath"],r["Study"],r["Series"])
+          return "file://"+path
 
-  def loadPatient(self,widget,patientId):
+  def loadPatient(self,patientId):
       print("Loading {}").format(patientId)
       ds=self.net.loadDataset("dinamic_spect/Patients","Imaging")
       for r in ds['rows']:
@@ -454,24 +459,21 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
           if not r['PatientId']==visit['aliasID']:
               continue
           if abs(r['SequenceNum']-float(visit['nmMaster']))<0.1:
-              path=r['_labkeyurl_Series']
-              masterPath=path[:path.rfind('/')]
+              masterPath=self.getFilespecPath(r)
               #masterPath="labkey://"+path
           if abs(r['SequenceNum']-float(visit['nmCorrData']))<0.1:
-              path=r['_labkeyurl_Series']
-              nmPath=path[:path.rfind('/')]
+              nmPath=self.getFilespecPath(r)
               #nmPath="labkey://"+path
           if abs(r['SequenceNum']-float(visit['ctData']))<0.1:
-              path=r['_labkeyurl_Series']
-              ctPath=path[:path.rfind('/')]
+              ctPath=self.getFilespecPath(r)
               #ctPath="labkey://"+path
 
 
-      self.pd.readMasterDirectory(self.formatToGetfileSpec(masterPath))
+      self.pd.readMasterDirectory(masterPath)
       self.frame_data, self.frame_time, self.frame_origin, \
-          self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(self.formatToGetfileSpec(nmPath))
+          self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(nmPath)
       self.ct_data,self.ct_origin,self.ct_pixel_size, \
-         self.ct_orientation=self.pd.readCTDirectory(self.formatToGetfileSpec(ctPath))
+         self.ct_orientation=self.pd.readCTDirectory(ctPath)
 
       self.ct_orientation=vi.completeOrientation(self.ct_orientation)
       self.frame_orientation=vi.completeOrientation(self.frame_orientation)
@@ -479,16 +481,44 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       self.addCT(patientId)
       self.addFrames(patientId)
 
-      widget.time_frame_select.setMaximum(self.frame_data.shape[3]-1)
+
+  def loadPatientNRRD(self,patientId):
+      print("Loading NRRD {}").format(patientId)
+
+      self.loadDummyInputFunction(patientId)
+      dnsNode=slicer.util.getFirstNodeByName(patientId+'_Dummy')
+      if dnsNode==None:
+          print("Could not find dummy double array node")
+          return
+
+      n=dnsNode.GetSize();
+      self.frame_time=np.zeros(n);
+
+      a=vtk.reference(1)
+      for i in range(0,n):
+          self.loadVolume(patientId,i)
+          self.frame_time[i]=dnsNode.GetValue(i,0,a)
+
+      self.loadCTVolume(patientId)
+      self.loadSegmentation(patientId)
+
+  def loadDummyInputFunction(self,patientId):
+      self.loadNode(patientId,patientId+'_Dummy','DoubleArrayFile','.mcsv')
+
+  def loadVolume(self,patientId,i):
+      self.loadNode(patientId,patientId+'Volume'+str(i),'VolumeFile')
+
+  def loadCTVolume(self,patientId):
+      self.loadNode(patientId,patientId+'CT','VolumeFile')
 
   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')
+      self.loadNode(patientId,'Heart','SegmentationFile')
 
+  def loadNode(self,patientId,fName,type,suffix='.nrrd'):
+      relativePath=self.coreRelativePath+'/'+patientId
+      labkeyFile=relativePath+'/'+fName+suffix
+      print ("Remote: {}").format(labkeyFile)
+      self.net.loadNode(labkeyFile,type,returnNode=True)
 
 
   def addNode(self,nodeName,v, lpsOrigin, pixel_size, lpsOrientation, dataType):
@@ -701,7 +731,8 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       if (node.__class__.__name__=="vtkMRMLTransformNode" or
         node.__class__.__name__=="vtkMRMLGridTransformNode"):
           suffix=".h5"
-      fileName=nodeName+suffix
+
+      fileName=re.sub(r'_RS$',r'',nodeName)+suffix
 
       if not os.path.isdir(localPath):
          os.mkdir(localPath)
@@ -714,14 +745,20 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       self.pd.net.put(remoteFile,f.read())
 
 
+
+
   def storeVolumeNodes(self,patientId,n1,n2):
       #n1=self.time_frame.minimum;
       #n2=self.time_frame.maximum
-      project="dinamic_spect/Patients"
-      relativePath=project+'/@files/'+patientId
+      relativePath=self.coreRelativePath+'/'+patientId
 
       print("Store CT")
       nodeName=patientId+'CT'
+
+      #prefer resampled
+      testNode=slicer.util.getFirstNodeByName(nodeName+"_RS")
+      if testNode:
+          nodeName=nodeName+"_RS"
       self.storeNodeRemote(relativePath,nodeName)
 
       print("Storing NM from {} to {}").format(n1,n2)
@@ -729,41 +766,48 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       for i in range(n):
           it=i+n1
           nodeName=patientId+'Volume'+str(it)
+
+          #prefer resampled
+          testNode=slicer.util.getFirstNodeByName(nodeName+"_RS")
+          if testNode:
+              nodeName=nodeName+"_RS"
+
           self.storeNodeRemote(relativePath,nodeName)
 
+      self.storeDummyInputFunction(patientId)
+
   def storeSegmentation(self,patientId):
-      project="dinamic_spect/Patients"
-      relativePath=project+'/@files/'+patientId
+      relativePath=self.coreRelativePath+'/'+patientId
       segNodeName="Heart"
       self.storeNodeRemote(relativePath,segNodeName)
 
   def storeInputFunction(self,patientId):
       self.calculateInputFunction(patientId)
-      project="dinamic_spect/Patients"
-      relativePath=project+'/@files/'+patientId
+      relativePath=self.coreRelativePath+'/'+patientId
       doubleArrayNodeName=patientId+'_Ventricle'
       self.storeNodeRemote(relativePath,doubleArrayNodeName)
 
+  def storeDummyInputFunction(self,patientId):
+      self.calculateDummyInputFunction(patientId)
+      relativePath=self.coreRelativePath+'/'+patientId
+      doubleArrayNodeName=patientId+'_Dummy'
+      self.storeNodeRemote(relativePath,doubleArrayNodeName)
+
   def storeTransformation(self,patientId):
-      project="dinamic_spect/Patients"
-      relativePath=project+'/@files/'+patientId
+      relativePath=self.coreRelativePath+'/'+patientId
       transformNodeName=patientId+"_DF"
       self.storeNodeRemote(relativePath,transformNodeName)
 
   def applyTransform(self, patientId,refPatientId,n1,n2):
+      if patientId == refPatientId:
+          print("Transform [{}] and reference [{}] are the same".format(patientId, refPatientId))
+          return
       transformNodeName=patientId+"_DF"
       transformNode=slicer.util.getFirstNodeByName(transformNodeName)
       if transformNode==None:
-          print("Node {} not found").format(transformNodeName)
+          print("Transform node [{}] not found").format(transformNodeName)
           return
 
-      try:
-          self.resampler.printMe()
-          return
-      except AttributeError:
-          print("Initializing resampler")
-          self.resampler=resample.resampleLogic(None)
-
 
       n=n2-n1+1
       for i in range(n):
@@ -794,7 +838,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
 
 
   def calculateInputFunction(self,patientId):
-       n=self.frame_data.shape[3]
+       n=len(self.frame_time)
 
        dnsNodeName=patientId+'_Ventricle'
        dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
@@ -843,6 +887,31 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
            dn.SetValue(i, 2, 0)
            print("[{0} at {1:.2f}:{2:.2f}]".format(vol,fx,fy))
 
+  def calculateDummyInputFunction(self,patientId):
+
+        n=self.frame_data.shape[3]
+
+        dnsNodeName=patientId+'_Dummy'
+        dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
+        if dns.GetNumberOfItems() == 0:
+            dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
+            dn.SetName(dnsNodeName)
+        else:
+            dn = dns.GetItemAsObject(0)
+
+
+
+        dn.SetSize(n)
+
+
+        ft=self.frame_time
+        for i in range(0,n):
+            fx=ft[i]
+
+            dn.SetValue(i, 0, fx)
+            dn.SetValue(i, 1, 0)
+            dn.SetValue(i, 2, 0)
+
 
 
 class cardiacSPECTTest(ScriptedLoadableModuleTest):

cardiacSPECT/parseDicom.py

@@ -73,14 +73,16 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
                 #not sure if labkey is available, but try it
                 print("Using labkey")
                 print("Server:{0}, file:{1}".format(self.net.GetHostName(),f))
-                return [self.net.readFileToBuffer(f),1]
+                localPath=self.net.DownloadFileToCache(f)
+                #return [self.net.readFileToBuffer(f),1]
+                return [open(localPath,'rb'),1]
             except:
                 print('Could not access labkey. Exiting')
                 return ['NULL',0]
 
         if origin.find('file')==0:
             print("Using local directory")
-            return [f,1]
+            return [open(f,'rb'),1]
 
         return ['NULL',0]