Version that can compute ROI values on the fly

cardiacSPECT/cardiacSPECT.py

@@ -8,6 +8,7 @@ import parseDicom as pd
 import vtkInterface as vi
 import fileIO
 import slicer
+import numpy as np
 #
 # cardiacSPECT
 #
@@ -274,14 +275,19 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
             dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
             a = dn.GetArray()
             a.SetNumberOfTuples(n)
+            dt=0;
+            t0=0;
             for i in range(0,n):
                 vol="testVolume"+str(i)
                 fx=ft[i]
                 fy=self.logic.meanROI(vol,segment)
+                dt=2*ft[i]-t0
+                t0+=dt
+
                 a.SetComponent(i, 0, fx)
-                a.SetComponent(i, 1, fy)
+                a.SetComponent(i, 1, fy/dt)
                 a.SetComponent(i, 2, 0)
-                print("({0:.2f}:{1:.2f})".format(fx,fy))
+                print("[{0} at {1:.2f}:{2:.2f}]".format(vol,fx,fy))
 
             cn.AddArray(segment, dn.GetID())
 
@@ -342,7 +348,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
         'Slicer Python','Data loaded')
 
 
-  def addNode(self,nodeName,v, origin, pixel_size, orientation, dataType):
+  def addNode(self,nodeName,v, lpsOrigin, pixel_size, lpsOrientation, dataType):
 
        # if dataType=0 it is CT data, which gets propagated to background an is
        #used to fit the view field dimensions
@@ -363,12 +369,15 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
        v.SetSpacing([1,1,1])
        ijkToRAS = vtk.vtkMatrix4x4()
        #think how to do this with image orientation
+       rasOrientation=[-lpsOrientation[i] if (i%3 < 2) else lpsOrientation[i]
+        for i in range(0,len(lpsOrientation))]
+       rasOrigin=[-lpsOrigin[i] if (i%3<2) else lpsOrigin[i] for i in range(0,len(lpsOrigin))]
 
        for i in range(0,3):
            for j in range(0,3):
-               ijkToRAS.SetElement(i,j,pixel_size[i]*orientation[3*j+i])
+               ijkToRAS.SetElement(i,j,pixel_size[i]*rasOrientation[3*j+i])
 
-           ijkToRAS.SetElement(i,3,origin[i])
+           ijkToRAS.SetElement(i,3,rasOrigin[i])
 
        newNode.SetIJKToRASMatrix(ijkToRAS)
        newNode.SetAndObserveImageData(v)
@@ -401,6 +410,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
 
     #get the segmentation mask
     fNode=slicer.mrmlScene.GetNodesByClass("vtkMRMLSegmentationNode").GetItemAsObject(0)
+    print "Found segmentation node: {}".format(fNode.GetName())
     segNode=slicer.vtkMRMLSegmentationNode.SafeDownCast(fNode)
 
     #no python bindings for vtkSegmentation
@@ -409,13 +419,12 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
     #    print("No segments available")
     #    return 0
 
+    segment=segNode.GetSegmentation().GetNthSegmentID(0)
     mask = segNode.GetBinaryLabelmapRepresentation(segment)
     if mask==None:
         print("Segment {} not found".format(segment))
         return s
 
-      #segNode.GetSegmentation().GetNthSegmentID(0))
-
     #get mask at (x,y,z)
     #mask.GetPointData().GetScalars().GetTuple1(mask.FindPoint([x,y,z]))
 
@@ -429,6 +438,10 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
     #iterate over volume pixelData
     n=dataImage.GetNumberOfPoints()
 
+    extent=mask.GetExtent()
+    fM=vtk.vtkMatrix4x4()
+    mask.GetWorldToImageMatrix(fM)
+
     for i in range(0,n):
 
       #get global coordinates of point i
@@ -436,10 +449,19 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
       position_ijk=[ix, iy, iz, 1]
       #ras are global coordinates (in mm)
       position_ras=ijkToRas.MultiplyPoint(position_ijk)
+      fpos=[int(np.round(x)) for x in fM.MultiplyPoint(position_ras)]
+
+      outOfRange=False
+      for k in range(0,3):
+          if fpos[k]<extent[2*k] or fpos[k]>extent[2*k+1]:
+              outOfRange=True
+              break;
+
+      if outOfRange:
+          continue
 
       #find point in mask with the same global coordinates
-      fr=[position_ras[0],position_ras[1],position_ras[2]]
-      maskValue=mask.GetPointData().GetScalars().GetTuple1(mask.FindPoint(fr))
+      maskValue=mask.GetPointData().GetScalars().GetTuple1(mask.ComputePointId(fpos[0:3]))
 
       if maskValue == 0:
           continue

cardiacSPECT/parseDicom.py

@@ -75,6 +75,8 @@ def getfile(origin,f):
     return ['NULL',0]
 
 def read_dynamic_SPECT(mypath):
+    axisShift=(2,1,0)
+
     origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
     onlyfiles=filelist(mypath)
     for f in onlyfiles:
@@ -139,8 +141,9 @@ def read_dynamic_SPECT(mypath):
         try:
             pf=plan[0x0018,0x5020]
         except:
-            print ("Tag not found")
+            print ("ProcessingFunction not found")
             continue
+
         try:
             phase=plan[0x0035,0x1005].value
             cycle=plan[0x0035,0x1004].value
@@ -163,7 +166,7 @@ def read_dynamic_SPECT(mypath):
 
     #play with pixel data
         if frame_data.shape[0] == 1:
-            sh=plan.pixel_array.shape;
+            sh=np.transpose(plan.pixel_array,axisShift).shape;
             sh=list(sh)
             sh.append(nframe[-1])
             frame_data=np.empty(sh)
@@ -201,7 +204,7 @@ def read_dynamic_SPECT(mypath):
 
 
 
-        frame_data[:,:,:,ifi]=plan.pixel_array
+        frame_data[:,:,:,ifi]=np.transpose(plan.pixel_array,axisShift)
 
     #print('Orientation: ({0:.2f},{1:.2f},{2:.2f}),({3:.2f},{4:.2f},{5:.2f})').format( \
     #    frame_orientation[0],frame_orientation[1],frame_orientation[2], \
@@ -215,6 +218,7 @@ def read_CT(mypath):
 
     ct_data = []
     ct_idx = []
+    ct_z = []
     ct_pixel_size = [0,0,0]
     ct_center = [0,0,0]
     ct_center[2]=1e30
@@ -256,6 +260,7 @@ def read_CT(mypath):
         print '.',
         ct_data.append(plan.pixel_array)
         ct_idx.append(plan.InstanceNumber)
+        ct_z.append(plan.ImagePositionPatient[2])
         #ct_center.append(plan.ImagePositionPatient)
 
         pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
@@ -275,7 +280,7 @@ def read_CT(mypath):
             if abs(ct_center[i]-plan.ImagePositionPatient[i]) > 1e-3:
                         print 'Image center mismatch {.2f}/{.2f}'.format(ct_center[i],
                         plan.ImagePositionPatient[i])
-        #not average, but minimum (!)
+        #not average, but minimum (!) why??
 
         if plan.ImagePositionPatient[2]<ct_center[2]:
             ct_center[2]=plan.ImagePositionPatient[2]
@@ -297,6 +302,7 @@ def read_CT(mypath):
     data_array=np.zeros(sh_list)
 
     for k in range(0,nz):
-        data_array[:,:,ct_idx[k]-1]=ct_data[k]
+        kp=int(np.round((ct_z[k]-ct_center[2])/ct_pixel_size[2]))
+        data_array[:,:,kp]=np.transpose(ct_data[k])
 
     return data_array,ct_center,ct_pixel_size,ct_orientation