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

            #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 = [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]
    ct_center[2]=1e30
    ct_orientation=[0,0,0,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
        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