Merge to head

cardiacSPECT/cardiacSPECT.py

@@ -4,8 +4,8 @@ import unittest
 import vtk, qt, ctk, slicer
 from slicer.ScriptedLoadableModule import *
 import logging
-import vtkInterface as vi
 import parseDicom as pd
+import vtkInterface as vi
 #
 # cardiacSPECT
 #

parseDicom.py

@@ -1,286 +0,0 @@
-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

vtkInterface.py

@@ -1,112 +0,0 @@
-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)