dynamicSPECT/cardiacSPECT/parseDicom.py

import os
import sys
import dicom
import numpy as np
import re
import slicer
from slicer.ScriptedLoadableModule import *

#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(ScriptedLoadableModule):
  def __init__(self, parent):
    ScriptedLoadableModule.__init__(self, parent)
    parent.title = "parseDicom"
    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

class parseDicomWidget(ScriptedLoadableModuleWidget):
    def setup(self):
        ScriptedLoadableModuleWidget.setup(self)
        self.logic=parseDicomLogic(self)

class parseDicomLogic(ScriptedLoadableModuleLogic):

    def __init__(self,parent):
        ScriptedLoadableModuleLogic.__init__(self, parent)

    def setURIHandler(self,net):
        self.net=net

    def filelist(self,mypath):
        if mypath.find('labkey://')==0:
            print("Using labkey")
            labkeyPath=re.sub('labkey://','',mypath)
            #not sure if labkey is available, so try it
            #url=slicer.modules.labkeySlicerPythonExtensionWidget.serverURL.text
            #print("Seting url={}".format(url))
            ok, files=self.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(self,origin,f):

        if origin.find('labkey')==0:
            try:
                #not sure if labkey is available, but try it
                print("Using labkey")
                url=self.net.GetHostName()
                print("Sever:{0}, file:{1}".format(url,f))
                return [self.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(self,mypath):
        axisShift=(2,1,0)

        origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
        onlyfiles=self.filelist(mypath)
        for f in onlyfiles:
            print '{}:'.format(f)

            g,ok=self.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=self.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(self,mypath):
        onlyfiles=self.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=self.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])

            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