Sorting out relative positioning of images + rotations using slicer's IJStoRAS matrix

cardiacSPECT/cardiacSPECT.py

@@ -151,22 +151,39 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
     sys.path.append(mypath)
     import parseDicom as pd
 
-    self.frame_data,self.frame_time=pd.read_dynamic_SPECT(inputDir)
+    self.frame_data, self.frame_time, self.frame_center, \
+        self.frame_pixel_size, self.frame_orientation=pd.read_dynamic_SPECT(inputDir)
 
-    self.ct_data,self.ct_center,self.ct_pixel_size=pd.read_CT(inputDir)
+    self.ct_data,self.ct_center,self.ct_pixel_size, \
+        self.ct_orientation=pd.read_CT(inputDir)
 
+    self.ct_orientation=vi.completeOrientation(self.ct_orientation)
+    self.frame_orientation=vi.completeOrientation(self.frame_orientation)
     #additional message via qt
     qt.QMessageBox.information(
         slicer.util.mainWindow(),
         'Slicer Python','Data loaded')
 
 
-  def addNode(self,nodeName,v):
+  def addNode(self,nodeName,v, orientation):
        #nodeName='testVolume'+str(it)
        newNode=slicer.vtkMRMLScalarVolumeNode()
        newNode.SetName(nodeName)
+
+       #pixel_size=[0,0,0]
+       pixel_size=v.GetSpacing()
+       print(pixel_size)
+       #origin=[0,0,0]
+       origin=v.GetOrigin()
        ijkToRAS = vtk.vtkMatrix4x4()
-       ijkToRAS.Identity()
+       #think how to do this with image orientation
+
+       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,3,origin[i])
+
        newNode.SetIJKToRASMatrix(ijkToRAS)
        newNode.SetAndObserveImageData(v)
        slicer.mrmlScene.AddNode(newNode)
@@ -180,13 +197,16 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
        frame_data=self.frame_data[:,:,:,it];
        nodeName='testVolume'+str(it)
        self.addNode(nodeName,
-        vi.numpyToVTK3D(frame_data,[0,0,0],[1,1,1]))
+         vi.numpyToVTK3D(frame_data,self.frame_center,
+           self.frame_pixel_size),self.frame_orientation)
 
   def addCT(self):
        nodeName='testCT'
        self.addNode(nodeName,
             vi.numpyToVTK3D(self.ct_data,
-                self.ct_center,self.ct_pixel_size))
+                self.ct_center,self.ct_pixel_size),
+            self.ct_orientation)
+
 
 
 class cardiacSPECTTest(ScriptedLoadableModuleTest):

parseDicom.py

@@ -47,11 +47,13 @@ def read_dynamic_SPECT(mypath):
             print ("Tag not found;")
             continue
         if not (type(nframe) is list) :
+            print("nframe not a list")
             continue
 
     #this is the "master" file where data on other files can be had
     #here we found out the duration of the frame and their distribution through
     #phases and cycles
+        print('Found master file')
 
         for i in range(1,len(nframe)):
             nframe[i]+=nframe[i-1]
@@ -64,13 +66,16 @@ def read_dynamic_SPECT(mypath):
         frame_start=plan[0x0019,0x10a7].value
         frame_stop=plan[0x0019,0x10a8].value
         frame_duration=plan[0x0019,0x10a9].value
-
+        break
     #print "rep [{}] start [{}] stop [{}] duration [{}]".format(
     #len(rep),len(rep_start),len(rep_stop),len(rep_duration))
 
 #select AC reconstructed data
-    frame_time=[None]*nframe[-1]
+    frame_time=np.zeros(nframe[-1]);
     frame_data=np.empty([1,1,1,nframe[-1]])
+    center = [0,0,0]
+    pixel_size =[0,0,0]
+    frame_orientation=[0,0,0,0,0,0]
     for f in onlyfiles:
         try:
             plan = dicom.read_file(os.path.join(mypath,f))
@@ -111,8 +116,43 @@ def read_dynamic_SPECT(mypath):
             frame_data=np.empty(sh)
             print "Setting frame_data to",sh
 
+            #check & update pixel size
+        pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
+                        plan.SliceThickness]
+
+        for i in range(0,3):
+            if pixel_size[i] == 0:
+                pixel_size[i] = float(pixel_size_read[i])
+            if abs(pixel_size[i]-pixel_size_read[i]) > 1e-3:
+                print 'Pixel size mismatch {.2f}/{.2f}'.format(pixel_size[i],
+                pixel_size_read[i])
+
+        center_read=plan.DetectorInformationSequence[0].ImagePositionPatient
+        for i in range(0,3):
+            if center[i] == 0:
+                center[i] = float(center_read[i])
+            if abs(center[i]-center_read[i]) > 1e-3:
+                        print 'Image center mismatch {.2f}/{.2f}'.format(center[i],
+                        center_read[i])
+
+        frame_orientation_read=plan.DetectorInformationSequence[0].ImageOrientationPatient
+        for i in range(0,6):
+            if frame_orientation[i] == 0:
+                frame_orientation[i] = float(frame_orientation_read[i])
+            if abs(frame_orientation[i]-frame_orientation_read[i]) > 1e-3:
+                        print 'Image orientation mismatch {.2f}/{.2f}'.format(
+                        frame_rotation[i], frame_orientation_read[i])
+
+
+
+
         frame_data[:,:,:,ifi]=plan.pixel_array
-    return [frame_data,frame_time]
+
+    #print('Orientation: ({0:.2f},{1:.2f},{2:.2f}),({3:.2f},{4:.2f},{5:.2f})').format( \
+    #    frame_orientation[0],frame_orientation[1],frame_orientation[2], \
+    #    frame_orientation[3],frame_orientation[4],frame_orientation[5])
+
+    return [frame_data,frame_time,center,pixel_size,frame_orientation]
 
 def read_CT(mypath):
     onlyfiles = [f for f in os.listdir(mypath)
@@ -122,7 +162,8 @@ def read_CT(mypath):
     ct_idx = []
     ct_pixel_size = [0,0,0]
     ct_center = [0,0,0]
-
+    ct_center[2]=1e30
+    ct_orientation=[0,0,0,0,0,0]
     for f in onlyfiles:
         print '{}:'.format(f)
         try:
@@ -162,9 +203,18 @@ def read_CT(mypath):
                         print 'Image center mismatch {.2f}/{.2f}'.format(ct_center[i],
                         plan.ImagePositionPatient[i])
         #not average, but minimum (!)
+
         if plan.ImagePositionPatient[2]<ct_center[2]:
             ct_center[2]=plan.ImagePositionPatient[2]
 
+        for i in range(0,6):
+            if ct_orientation[i] == 0:
+                ct_orientation[i] = float(plan.ImageOrientationPatient[i])
+            if abs(ct_orientation[i]-plan.ImageOrientationPatient[i]) > 1e-3:
+                print 'Image orientation mismatch {0:.2f}/{1:.2f}'.format(ct_orientation[i],
+                plan.ImageOrientationPatient[i])
+
+
     nz=len(ct_idx)
     #not average, again
     #ct_center[2]/=nz
@@ -176,4 +226,4 @@ def read_CT(mypath):
     for k in range(0,nz):
         data_array[:,:,ct_idx[k]-1]=ct_data[k]
 
-    return data_array,ct_center,ct_pixel_size
+    return data_array,ct_center,ct_pixel_size,ct_orientation

vtkInterface.py

@@ -35,3 +35,10 @@ def numpyToVTK3D(numpy_array, origin, spacing):
                scalars.InsertTuple1(offset,numpy_array[i,j,k])
     v.GetPointData().SetScalars(scalars)
     return v
+
+def completeOrientation(orientation):
+    o=orientation
+    o.append(o[1]*o[5]-o[2]*o[4])#0,3
+    o.append(o[2]*o[3]-o[0]*o[5])#1,4
+    o.append(o[0]*o[4]-o[1]*o[3])#2,5
+    return o