Moved py files to cardiacSPECT

cardiacSPECT/parseDicom.py

@@ -0,0 +1,286 @@
+import os
+import sys
+import dicom
+import numpy as np
+import re
+import slicer
+
+#rom os import listdir
+#from os.path import isfile, join
+#onlyfiles = [f for f in listdir(mypath) if isfile(join(mypath, f))]
+#import Tkinter as tk
+#from Tkinter import filedialog
+
+#root = tk.Tk()
+#root.withdraw()
+#file_path = filedialog.askopenfilename()
+class parseDicom:
+  def __init__(self, parent):
+    parent.title = "parse dicom"
+    parent.categories = ["Examples"]
+    parent.dependencies = []
+    parent.contributors = ["Andrej Studen (FMF/JSI)"] # replace with "Firstname Lastname (Org)"
+    parent.helpText = """
+    Parse dynamic SPECT DICOM files
+    """
+    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 = parent
+
+def filelist(mypath):
+#mypath=os.environ['PWD']
+#list files
+    if mypath.find('labkey://')==0:
+        print("Using labkey")
+        labkeyPath=re.sub('labkey://','',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(labkeyPath)
+        print 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(origin,f):
+
+    if origin.find('labkey')==0:
+        try:
+            #not sure if labkey is available, but try it
+            net=slicer.modules.labkeySlicerPythonExtensionWidget.network
+            print("Using labkey")
+            url=slicer.modules.labkeySlicerPythonExtensionWidget.serverURL.text
+            print("Sever:{0}, file:{1}".format(url,f))
+            return [net.readFile(str(url),f),1]
+        except:
+            print('Could not access labkey. Exiting')
+            return ['NULL',0]
+
+    if origin.find('file')==0:
+        print("Using local directory")
+        return [f,1]
+
+    return ['NULL',0]
+
+def read_dynamic_SPECT(mypath):
+    origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
+    onlyfiles=filelist(mypath)
+    for f in onlyfiles:
+        print '{}:'.format(f)
+
+        g,ok=getfile(origin,f)
+        if not(ok):
+                return
+
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("Not a dicom file")
+            continue
+        try:
+            nframe=plan[0x0019,0x10a5].value;
+        except:
+            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]
+
+        print(nframe)
+
+    #nframe now holds for index i total number of frames collected up
+    #to the end of each phase
+
+        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=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:
+
+        g,ok=getfile(origin,f)
+        if not(ok):
+                return
+
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("Not a dicom file")
+            continue
+
+        try:
+            pf=plan[0x0018,0x5020]
+        except:
+            print ("Tag not found")
+            continue
+        try:
+            phase=plan[0x0035,0x1005].value
+            cycle=plan[0x0035,0x1004].value
+        except:
+            print ("Phase/Cycle tag not found")
+            continue
+
+        #convert phase/cycle to frame index
+        off=0
+        if phase > 1:
+            off=nframe[phase-2]
+        ifi=off+cycle-1
+
+    #from values in the master file determine frame time
+    #(as the mid point between starting and ending the frame)
+        frame_time[ifi]=0.5*(frame_start[ifi]+frame_stop[ifi]); #in ms
+
+        print "({},{}) converted to {} at {} for {}".format(
+            phase,cycle,ifi,frame_time[ifi],frame_duration[ifi])
+
+    #play with pixel data
+        if frame_data.shape[0] == 1:
+            sh=plan.pixel_array.shape;
+            sh=list(sh)
+            sh.append(nframe[-1])
+            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
+
+    #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=filelist(mypath)
+    origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
+
+    ct_data = []
+    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)
+
+        g,ok=getfile(origin,f)
+        if not(ok):
+                return
+
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("Not a dicom file")
+            continue
+
+        if plan.Modality != 'CT':
+            print ('Not a CT file')
+            continue
+
+        if re.match("AC",plan.SeriesDescription) == None:
+            print (plan.SeriesDescription)
+            print ('Not a AC file')
+            continue
+        #a slice of pure CT
+        print '.',
+        ct_data.append(plan.pixel_array)
+        ct_idx.append(plan.InstanceNumber)
+        #ct_center.append(plan.ImagePositionPatient)
+
+        pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
+            plan.SliceThickness]
+
+
+        for i in range(0,3):
+            if ct_pixel_size[i] == 0:
+                ct_pixel_size[i] = float(pixel_size_read[i])
+            if abs(ct_pixel_size[i]-pixel_size_read[i]) > 1e-3:
+                print 'Pixel size mismatch {.2f}/{.2f}'.format(ct_pixel_size[i],
+                pixel_size_read[i])
+
+        for i in range(0,2):
+            if ct_center[i] == 0:
+                ct_center[i] = float(plan.ImagePositionPatient[i])
+            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 (!)
+
+        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])
+
+    print
+    nz=len(ct_idx)
+    #not average, again
+    #ct_center[2]/=nz
+    sh=ct_data[-1].shape
+    sh_list=list(sh)
+    sh_list.append(nz)
+    data_array=np.zeros(sh_list)
+
+    for k in range(0,nz):
+        data_array[:,:,ct_idx[k]-1]=ct_data[k]
+
+    return data_array,ct_center,ct_pixel_size,ct_orientation

cardiacSPECT/vtkInterface.py

@@ -0,0 +1,112 @@
+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):
+    parent.title = "vtk Interface"
+    parent.categories = ["Examples"]
+    parent.dependencies = []
+    parent.contributors = ["Andrej Studen (FMF/JSI)"] # replace with "Firstname Lastname (Org)"
+    parent.helpText = """
+    Convert native numpy data structures to vtk
+    """
+    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 = parent
+
+
+def numpyToVTK(numpy_array, shape, data_type=None):
+    v=vtk.vtkImageData()
+    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)