Updating vtk interface to allow changing from vtk to itk to numpy image representation; allow cardiac SPECT to use labkey directly

cardiacSPECT/cardiacSPECT.py

@@ -51,7 +51,8 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     # Layout within the sample collapsible button
     dataFormLayout = qt.QFormLayout(dataButton)
 
-    pathGuess="file://"+os.environ['HOME']+"/SPECT"
+    #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)
 

parseDicom.py

@@ -39,9 +39,9 @@ def filelist(mypath):
         #not sure if labkey is available, so try it
         net=slicer.modules.labkeySlicerPythonExtensionWidget.network
         print("Found network")
-        url=slicer.modules.labkeySlicerPythonExtensionWidget.serverURL.text
-        print("Seting url={}".format(url))
-        files=net.listDir(str(url),labkeyPath)
+        #url=slicer.modules.labkeySlicerPythonExtensionWidget.serverURL.text
+        #print("Seting url={}".format(url))
+        files=net.listDir(labkeyPath)
         print files
 
     if mypath.find('file://')==0:

vtkInterface.py

@@ -1,5 +1,10 @@
 import vtk, qt, ctk, slicer
 import numpy as np
+import SimpleITK as sitk
+
+#set of routines to transform images from one form to another, most notably
+#numpy to vtk to itk and all possible combinations inbetween. Keep track of
+#orientation, origin and spacing between transforms
 
 class vtkInterface:
   def __init__(self, parent):
@@ -17,30 +22,91 @@ class vtkInterface:
     """ # replace with organization, grant and thanks.
     self.parent = parent
 
-def numpyToVTK3D(numpy_array, origin, spacing):
+
+def numpyToVTK(numpy_array, shape, data_type=None):
     v=vtk.vtkImageData()
-    v.SetDimensions(
-            numpy_array.shape[0],
-            numpy_array.shape[1],
-            numpy_array.shape[2])
-    v.SetOrigin(origin[0], origin[1], origin[2])
-    v.SetSpacing(spacing[0], spacing[1], spacing[2])
-    scalars = vtk.vtkShortArray()
-    for k in range(0,numpy_array.shape[2]):
-       #print '.',
-       kOffset = k * numpy_array.shape[1]*numpy_array.shape[0]
-       for j in range(0,numpy_array.shape[1]):
-           jOffset = j * numpy_array.shape[0]
-           for i in range(0,numpy_array.shape[0]):
-               offset = i + jOffset + kOffset;
-               scalars.InsertTuple1(offset,numpy_array[i,j,k])
-    #print
-    v.GetPointData().SetScalars(scalars)
+    v.SetDimensions(shape)
+    v.SetOrigin(0,0,0)
+    v.SetSpacing(1,1,1)
+    v.GetPointData().SetScalars(
+        vtk.util.numpy_support.numpy_to_vtk(
+            numpy_array,deep=True, array_type=data_type))
+
     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
+
+
+def ITK2VTK(img):
+    #convert itk to vtk format.
+    #get the array
+    data=sitk.GetArrayFromImage(img)
+    #reverse the shape (don't ask, look at vtk manual if really curios)
+    shape=list(reversed(data.shape))
+    return numpyToVTK(data.ravel(),shape)
+
+def VTK2ITK(v):
+    #convert vtk image to sitk image
+    #convert to numpy first and then go to sitk
+    scalars=v.GetPointData().GetScalars()
+    shape=v.GetDimensions()
+    data=vtk.util.numpy_support.vtk_to_numpy(scalars)
+    #now convert to sitk (notice the little reversal of the shape)
+    return sitk.GetImageFromArray(np.reshape(data,list(reversed(shape))))
+
+def ITKfromNode(nodeName):
+    #use node as data source and generate an itk image
+    node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
+    if node==None:
+        print "Node {0} not available".format(nodeName)
+        return
+
+    img=VTK2ITK(node.GetImageData())
+
+    img.SetOrigin(node.GetOrigin())
+    img.SetSpacing(node.GetSpacing())
+    m=vtk.vtkMatrix4x4()
+    node.GetIJKToRASDirectionMatrix(m)
+    orientation=[0]*9
+    for i in range(0,3):
+        for j in range (0,3):
+            orientation[3*j+i]=m.GetElement(i,j)
+    img.SetDirection(orientation)
+    return img
+
+
+
+def ITKtoNode(img,nodeName):
+    #display itk image and assign it a volume node
+    #useful for displaying outcomes of itk calculations
+    node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
+    if node==None:
+        node=slicer.vtkMRMLScalarVolumeNode()
+        node.SetName(nodeName)
+        slicer.mrmlScene.AddNode(node)
+
+    node.SetAndObserveImageData(ITK2VTK(img))
+
+    #hairy - keep orientation, spacing and origin from node and pass it to itk
+    #for future reference
+    spacing=img.GetSpacing()
+    orientation=img.GetDirection()
+    origin=img.GetOrigin()
+
+    #we should get orientation, spacing and origin from somewhere
+    ijkToRAS = vtk.vtkMatrix4x4()
+    #think how to do this with image orientation
+
+    for i in range(0,3):
+       for j in range(0,3):
+           ijkToRAS.SetElement(i,j,spacing[i]*orientation[3*j+i])
+
+       ijkToRAS.SetElement(i,3,origin[i])
+
+    node.SetIJKToRASMatrix(ijkToRAS)