dynamicSPECT/parseDicom.py

import os
import sys
import dicom
import numpy as np
import re
#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 read_dynamic_SPECT(mypath):
#mypath=os.environ['PWD']
#list files
    onlyfiles = [f for f in os.listdir(mypath)
            if os.path.isfile(os.path.join(mypath, f))]

    for f in onlyfiles:
        print '{}:'.format(f)
        try:
            plan = dicom.read_file(os.path.join(mypath,f))
        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) :
            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

        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

    #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_data=np.empty([1,1,1,nframe[-1]])
    for f in onlyfiles:
        try:
            plan = dicom.read_file(os.path.join(mypath,f))
        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

        frame_data[:,:,:,ifi]=plan.pixel_array
    return [frame_data,frame_time]

def read_CT(mypath):
    onlyfiles = [f for f in os.listdir(mypath)
            if os.path.isfile(os.path.join(mypath, f))]

    ct_data = []
    ct_idx = []
    ct_pixel_size = [0,0,0]
    ct_center = [0,0,0]

    for f in onlyfiles:
        print '{}:'.format(f)
        try:
            plan = dicom.read_file(os.path.join(mypath,f))
        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
        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]

    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