dynamicSPECT/cardiacSPECT/parseDicom.py

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))
        ok, files=net.listRemoteDir(labkeyPath)
        if not ok:
            print "Error accessing path"
            return []

    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):
    axisShift=(2,1,0)

    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):
                continue

        try:
            plan = dicom.read_file(g)
        except:
            print ("Not a dicom file")
            continue

        try:
            pf=plan[0x0018,0x5020]
        except:
            print ("ProcessingFunction 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=np.transpose(plan.pixel_array,axisShift).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
        print "Stored center at ({0},{1},{2})".format(center[0],center[1],center[2])
        print "Read   center at ({0},{1},{2})".format(center_read[0],center_read[1],center_read[2])
        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]=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], \
    #    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_z = []
    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

        #this doesn't work in 2019 data version
        #if re.match("AC",plan.SeriesDescription) == None:
        #    print (plan.SeriesDescription)
        #    print ('Not a AC file')
        #    continue
        try:
            iType=plan.ImageType
        except:
            print "Image type not found"
            continue;

        if iType[3].find("SPI")<0:
            print "Not a spiral image"
            continue;


        #a slice of pure CT
        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],
            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 (!) why??

        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):
        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