CTCauto0.6.tar

CTC_auto.py0000664000175100017510000004744113165066156011544 0ustar  mcqamcqa#!/usr/bin/env python

# This is the main file for CTC_auto
# CTC_auto is script that converts DICOM files (RP+RD+RS-CT) to a DOSXYZnrc
# compatible egsphant phantom
# It is a modified version of CTC_ask as published in:
# RO Ottosson and CF Behrens. CTC-ask: a new algorithm for conversion of CT
# numbers to tissue parameters for Monte Carlo dose calculations applying
# DICOM RS knowledge. Phys. Med. Biol. 56 (2011) N1-N12
#
#    Copyright [2016] [Rickard Cronholm] Licensed under the
#    Educational Community License, Version 2.0 (the "License"); you may
#    not use this file except in compliance with the License. You may
#    obtain a copy of the License at
#
#http://www.osedu.org/licenses/ECL-2.0
#
#    Unless required by applicable law or agreed to in writing,
#    software distributed under the License is distributed on an "AS IS"
#    BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
#    or implied. See the License for the specific language governing
#    permissions and limitations under the License.
#
# v0.4
#
# Usage
# python CTC_auto RPRDfile RSfile CTfile fileName
# or
# python CTC_auto DICOMdirectory fileName
#

import sys
import dicom
import numpy as np
import CTCtools
import copy
import os
from scipy import ndimage

# full path to configuration file
confFile = '/home/mcqa/MCQA/samc.conf'


class struct:
    def __init__(self):
        pass

    def init_name(self):
        self.name = []
        self.type = []
        self.logicMatrix = []
        self.materialRamp = []
        self.material = []


def main(RPRDfile, RSfile, CTfile, fileName, addStructType=[], addRampName=[]):
    # set some variables
    with open(confFile) as f:
        confCont = f.readlines()
    whatToGet = ['common.RPfilePrefix', 'common.RDfilePrefix',
        'common.RSfilePrefix', 'common.CTfilePrefix', 'common.DICOMfileEnding',
        'CTC_auto.fixedMedDens', 'CTC_auto.extName',
        'CTC_auto.suppName', 'CTC_auto.suppOuter',
        'CTC_auto.suppInner', 'CTC_auto.densRampName',
        'CTC_auto.externalRamp', 'CTC_auto.outsideRamp',
        'CTC_auto.otherwiseRamp', 'CTC_auto.airDens',
        'CTC_auto.setAir', 'CTC_auto.lowerDens',
        'CTC_auto.spaceDelimit', 'CTC_auto.relElecFile']

    # get the variables from confFile
    cv = struct()
    cv = CTCtools.getConfVars(cv, whatToGet, confCont)

    # define additional variables
    replaceType = 'NONE'  # the structure type to be given to replaced types
    replaceList = []  # list of structure types to be replaced
    extFilterSize = 3  # number of voxel to add to external contour
    extFilt = (extFilterSize, extFilterSize, extFilterSize)
    filt = np.ones(extFilt).astype(int)

    # load density ramp
    # densRamp = CTCtools.getStructFromMatFile(cv.densRampName, 'densRamp', 1)
    densRamp = CTCtools.grabData(cv.densRampName, 'densRamp', 1, 1)

    # Use DICOM to obtain data

    # CT
    with open(CTfile) as f:
        ct = f.readlines()
    ct = map(str.strip, ct)
    # get x and y coords, just use first file in list as all will have same grid
    ct_xmesh, ct_ymesh = CTCtools.getDICOMcoords(ct[0], True)  # get DICOMcoords in cm
    # get CT matrix and z coords
    ct_zmesh, ct_mtrx = CTCtools.getCTinfo(ct, True)

    # RD
    with open(RPRDfile) as f:
        rp = f.readline().strip()
        rd = f.readline().strip()
    # get x and y DICOMcoords in cm
    rd_xmesh, rd_ymesh = CTCtools.getDICOMcoords(rd, True)
    # get z coords
    rd_zmesh = CTCtools.getDICOMzCoord(rd, True)

    # get orientation of CT and RD
    orientCT = CTCtools.getOrient(ct[0])
    orientRD = CTCtools.getOrient(rd)

    # inverse rd_*mesh depending on orientation relative ct
    if sum(orientRD[:3]) == -sum(orientCT[:3]):
        rd_xmesh = rd_xmesh[::-1]
    if sum(orientRD[3:]) == -sum(orientCT[3:]):
        rd_ymesh = rd_ymesh[::-1]
    if len(np.unique([np.all(np.diff(ct_zmesh) > 0), np.all(np.diff(rd_zmesh) > 0)])) > 1:
        step = np.unique(np.around(np.diff(ct_zmesh), decimals=3))
        if len(step) == 1:
            step = float(step)
        else:
            step = float(np.average(step))
        rd_zmesh = CTCtools.create1dDICOMcoord(rd_zmesh[0], step,
        len(rd_zmesh), -1)
        rd_zmesh = rd_zmesh[::-1]
    '''
    try:
        int(np.argwhere(ct_zmesh == rd_zmesh[-1]))
    except TypeError:
        step = np.unique(np.around(np.diff(ct_zmesh), decimals=3))
        if len(step) == 1:
            step = float(step)
        else:
            step = float(np.average(step))
        rd_zmesh = CTCtools.create1dDICOMcoord(rd_zmesh[0], step,
        len(rd_zmesh), -1)
        rd_zmesh = rd_zmesh[::-1]
    '''

    # inverse arrays and flip matrix depending on meshes
    if ct_xmesh[-1] < ct_xmesh[0]:
        ct_xmesh = ct_xmesh[::-1]
        rd_xmesh = rd_xmesh[::-1]
        ct_mtrx = ct_mtrx[:, :, ::-1]
    if ct_ymesh[-1] < ct_ymesh[0]:
        ct_ymesh = ct_ymesh[::-1]
        rd_ymesh = rd_ymesh[::-1]
        ct_mtrx = ct_mtrx[:, ::-1, :]
    if ct_zmesh[-1] < ct_zmesh[0]:
        ct_zmesh = ct_zmesh[::-1]
        rd_zmesh = rd_zmesh[::-1]
        ct_mtrx = ct_mtrx[::-1, :, :]

    # write doseGrid
    f = open(''.join([fileName, '.doseGrid']), 'w')
    f.write('{0:12.8f} {1:12.8f} {2:d}\n'.format(rd_xmesh[0], rd_xmesh[-1],
    len(rd_xmesh)))
    f.write('{0:12.8f} {1:12.8f} {2:d}\n'.format(rd_ymesh[0], rd_ymesh[-1],
     len(rd_ymesh)))
    f.write('{0:12.8f} {1:12.8f} {2:d}\n'.format(rd_zmesh[0], rd_zmesh[-1],
    len(rd_zmesh)))
    f.close()

    cts_xmesh = ct_xmesh[:]
    cts_ymesh = ct_ymesh[:]
    cts_zmesh = ct_zmesh[:]
    rds_xmesh = rd_xmesh[:]
    rds_ymesh = rd_ymesh[:]
    rds_zmesh = rd_zmesh[:]

    refROIs = []
    refContSeq = []

    # Read RS file
    with open(RSfile) as f:
        rs = f.readline().strip()
    RS = dicom.read_file(rs)  # open file
    # replace empty strucutre types with specified string
    for i in range(0, len(RS.RTROIObservationsSequence)):
        if len(RS.RTROIObservationsSequence[i].RTROIInterpretedType) == 0:
            RS.RTROIObservationsSequence[i].RTROIInterpretedType = replaceType
        try:
            replaceList.index(RS.RTROIObservationsSequence[i].RTROIInterpretedType)
            RS.RTROIObservationsSequence[i].RTROIInterpretedType = replaceType
        except ValueError:
            pass

    # start reading RS struct types
    allTypesFull = []
    for i in range(0, len(RS.RTROIObservationsSequence)):
        allTypesFull.append(RS.RTROIObservationsSequence[i].RTROIInterpretedType)
    allTypes = list(set(allTypesFull))
    # Check for support structures
    try:
        allTypes.index(cv.suppName)
        supportStructures = True
    except ValueError:
        supportStructures = False

    if supportStructures:
        suppDim = struct()
        suppDim.X = [None] * 2
        suppDim.Y = [None] * 2
        suppDim.Z = [None] * 2
        # get indices of support structures
        suppNr = [i for i, x in enumerate(allTypesFull) if x == cv.suppName]
        # get the ref ROI number(s)
        for elem in suppNr:
            refROIs.append(int(RS.ROIContourSequence[elem].ReferencedROINumber))
        # find the corresponding ContourSequence
        for elem in refROIs:
            refContSeq.append(CTCtools.getCorrContSeq(RS.ROIContourSequence, elem))
        # get the extreme of the support structures
        for elem in refContSeq:
            suppDim = CTCtools.getExtremeOfContour(CTCtools.getContour(RS.ROIContourSequence[elem], ct_zmesh, abs(orientCT[1]), True), suppDim)
        # extend rd_*mesh to encompass support structures
        rds_xmesh = CTCtools.extendMesh(rd_xmesh, suppDim.X)
        rds_ymesh = CTCtools.extendMesh(rd_ymesh, suppDim.Y)
        rds_zmesh = CTCtools.extendMesh(rd_zmesh, suppDim.Z)
        suppDim.X = [rds_xmesh[0], rds_xmesh[-1]]
        suppDim.Y = [rds_ymesh[0], rds_ymesh[-1]]
        suppDim.Z = [rds_zmesh[0], rds_zmesh[-1]]
        cts_xmesh = CTCtools.extendMesh(ct_xmesh, suppDim.X)
        cts_ymesh = CTCtools.extendMesh(ct_ymesh, suppDim.Y)
        cts_zmesh = CTCtools.extendMesh(ct_zmesh, suppDim.Z)
        # pad ct using map_coordinates if cts != ct
        if ct_mtrx.shape != (len(cts_zmesh), len(cts_ymesh), len(cts_xmesh)):
            ct_mtrx = CTCtools.map_coordinates(ct_mtrx, ct_xmesh, ct_ymesh, ct_zmesh, cts_xmesh, cts_ymesh, cts_zmesh, 0)

    # Deinterpolate CT data onto dose grid
    # check if rds_zmesh is beyond cts_zmesh, if so eliminate slices
    rds_zmesh = np.intersect1d(np.around(rds_zmesh, decimals=3), np.around(cts_zmesh, decimals=3), assume_unique=True)
    # density matrix is computed using cubic interpolation
    dens_mtrx = CTCtools.map_coordinates(ct_mtrx, cts_xmesh, cts_ymesh,
    cts_zmesh, rds_xmesh, rds_ymesh, rds_zmesh, 3)
    #ct_mtrx = CTCtools.map_coordinates(ct_mtrx, cts_xmesh, cts_ymesh,
    #cts_zmesh, rds_xmesh, rds_ymesh, rds_zmesh, 0)
    # change: use cubic interpolation also for material
    #ct_mtrx = CTCtools.map_coordinates(ct_mtrx, cts_xmesh, cts_ymesh,
    #cts_zmesh, #rds_xmesh, rds_ymesh, rds_zmesh, 3)
    ct_mtrx = copy.deepcopy(dens_mtrx)
    # set nans to 0
    dens_mtrx = np.nan_to_num(dens_mtrx)
    ct_mtrx = np.nan_to_num(ct_mtrx)

    # Locate external contour
    extNr = [i for i, x in enumerate(allTypesFull)
    if x.startswith(cv.extName)][0]
    refROIs.append(int(RS.ROIContourSequence[extNr].ReferencedROINumber))
    refContSeq.append(CTCtools.getCorrContSeq(RS.ROIContourSequence,
    refROIs[-1]))

    # Check if additional structure types were requested
    if len(addStructType) > 0:
        for i in range(0, len(addStructType)):
            addNr = [j for j, x in enumerate(allTypesFull)
            if x == addStructType[i]]
            for j in range(0, len(addNr)):
                refROIs.append(int(RS.ROIContourSequence[addNr[j]].ReferencedROINumber))
                refContSeq.append(CTCtools.getCorrContSeq(RS.ROIContourSequence, refROIs[-1]))

    # Correlate between ROIContourSequence and RTROIObservationsSequence
    refObsSeq = []
    for elem in refROIs:
        refObsSeq.append(CTCtools.getCorrContSeq(RS.RTROIObservationsSequence,
        elem))

    # create and init structures
    structures = []
    structureShell = struct()

    # get name, types and properties
    names = []
    cnt = 0
    skips = []
    for elem in refObsSeq:
        structure = copy.deepcopy(structureShell)
        # structure.init_name()
        try:
            structure.name = RS.RTROIObservationsSequence[elem].ROIObservationLabel
            names.append(RS.RTROIObservationsSequence[elem].ROIObservationLabel)
        except AttributeError:
            structure.name = ''
            names.append('')
        # names.append(RS.RTROIObservationsSequence[elem].ROIObservationLabel)
        structure.type = RS.RTROIObservationsSequence[elem].RTROIInterpretedType
        try:
            if RS.RTROIObservationsSequence[elem].ROIPhysicalPropertiesSequence[0].ROIPhysicalProperty == 'REL_ELEC_DENSITY':
                structure.RelElecDens = format(float(RS.RTROIObservationsSequence[elem].ROIPhysicalPropertiesSequence[0].ROIPhysicalPropertyValue), '.4f')
        except AttributeError:
            pass
        #if not structure.type == cv.suppName:
        #structure.RelElecDens = 1.0000  # set to 1 for all structures
        # don't append if name starts with Z_Ring
        skipNames = ['z_Tring']
        if not any(structure.name.lower().startswith(item) for item in skipNames):
            structures.append(structure)
        else:
            skips.append(cnt)
        cnt += 1

    # remove the skippers
    [refObsSeq.pop(x) for x in list(reversed(skips))]
    [refContSeq.pop(x) for x in list(reversed(skips))]

    # create and append structure for outside
    structure = copy.deepcopy(structureShell)
    names.append('OUTSIDE')
    structure.name = 'OUTSIDE'
    structure.type = 'OUTSIDE'
    structures.append(structure)

    # get contour sequences and create logicMatrix
    cnt = 0

    for elem in refContSeq:
        # get contour
        structures[cnt].contour = CTCtools.getContour(RS.ROIContourSequence[elem], ct_zmesh, abs(orientCT[1]), True)
        # deInterpolate contour onto dose grid and generate boolean matrix
        structures[cnt].logicMatrix = CTCtools.interpStructToDose(structures[cnt].contour, rds_xmesh, rds_ymesh, rds_zmesh, cts_xmesh, cts_ymesh, cts_zmesh)
        cnt += 1

    # outside
    out = [i for i, x in enumerate(names) if 'OUTSIDE' in x][0]
    structures[out].logicMatrix = np.ones(structures[0].logicMatrix.shape)
    
    # drop empty structures (due to being too small)
    structures = CTCtools.dropEmpty(structures)  

    # drop structures that are not engulfed by external save dedicated types
    saveTypes = ['OUTSIDE', cv.suppName, 'BOLUS']
    structures = CTCtools.dropOutOfExt(structures, cv.extName, saveTypes)

    # expand external contour by convolution
    for cnt in range(0, len(structures)):
        if structures[cnt].type == cv.extName and not hasattr(structures[cnt],'RelElecDens'):  # added && not releELec !!!
            indx = np.where(ndimage.convolve(structures[cnt].logicMatrix,filt,mode='nearest') >= 1)
            structures[cnt].logicMatrix = np.zeros(structures[cnt].logicMatrix.shape).astype(int)
            structures[cnt].logicMatrix[indx] = 1
            break

    # make sure that each voxel only belongs to one structure
    if supportStructures:
        # remove suppInner from suppOuter
        inner = [i for i, x in enumerate(names) if cv.suppInner in x][0]
        outer = [i for i, x in enumerate(names) if cv.suppOuter in x][0]
        structures[outer].logicMatrix = np.where(structures[inner].logicMatrix == 1, 0, structures[outer].logicMatrix)
    sortOrder = [cv.suppName, 'OUTSIDE', cv.extName, 'CTV', True]
    structures = CTCtools.sortStructures(structures, sortOrder)  # sort
    # belong to one
    structures = CTCtools.belongToOne(structures, cv.extName)
    structures = CTCtools.dropEmpty(structures)  # drop empty

    # perform HU corrections
    if cv.setAir.lower() == 'y':
        out = [i for i, x in enumerate(structures) if 'OUTSIDE' in x.name][0]
        # assumption: the density of air is below the breakpoint in the bilinear HU-dense curve
        airHU = (cv.airDens - densRamp[0][1][0]) / densRamp[0][1][1]
        ct_mtrx = np.where(structures[out].logicMatrix == 1, airHU, ct_mtrx)
        dens_mtrx = np.where(structures[out].logicMatrix == 1, airHU, dens_mtrx)

    # compute density matrix
    density = CTCtools.computeDensity(dens_mtrx, densRamp)
    # set minimum density if value is positive
    if cv.lowerDens >= 0:
        density = np.where(density < cv.lowerDens, cv.lowerDens, density)

    # assign material ramps to each structure
    for i in range(0, len(structures)):
        if structures[i].type == cv.extName:
            rampName = cv.externalRamp
        elif structures[i].type == cv.suppName:
            rampName = cv.otherwiseRamp
        elif structures[i].type == 'OUTSIDE':
            rampName = cv.outsideRamp
        else:
            # find correct rampName
            indx = [j for j, x in enumerate(addStructType)
            if x == structures[i].type][0]
            rampName = addRampName[indx]
        print structures[i].name, rampName
        structures[i].ramp = CTCtools.grabData(rampName, 'materialRamp', 0, 2)

    # Build and sort total media list
    medNr, medium = CTCtools.buildGlobalMediaList(structures)

    # compute media matrix
    media = CTCtools.computeMedia(ct_mtrx, structures, medium, medNr)

    # set uniform media and density for structures with defined relElec
    relElec = struct()
    relElec.data = CTCtools.getFromFile(cv.relElecFile, 3)
    relElec.relElec = [x[1] for x in relElec.data]
    relElec.physDens = [x[2] for x in relElec.data]
    relElec.media = [x[0] for x in relElec.data]
    for struc in structures:
        if hasattr(struc, 'RelElecDens'):
            try:
                rIndex = relElec.relElec.index(getattr(struc, 'RelElecDens'))
                # get medNr corresponding to media
                try:
                    medIndx = medium.index(relElec.media[rIndex])
                    media = np.where(struc.logicMatrix == 1, medNr[medIndx],
                    media)
                    density = np.where(struc.logicMatrix == 1,
                    relElec.physDens[rIndex], density)
                except ValueError:  # add if not in list
                    medium.append(relElec.media[rIndex])
                    medNr.append(max(medNr) + 1)
                    media = np.where(struc.logicMatrix == 1, max(medNr),
                    media)
                    density = np.where(struc.logicMatrix == 1,
                    relElec.physDens[rIndex], density)
            except ValueError:
                pass

    # set fixed density for medias listed in cv.fixedMedDens
    if os.path.isfile(cv.fixedMedDens):
        with open(cv.fixedMedDens) as fmd:
            fixedMedDens = fmd.readlines()
        for line in fixedMedDens:
            try:
                fixedMedia = line.split('\t')[0].strip()
                fixedDensity = float(line.split('\t')[1].strip())
                # get medNr corresponding to media
                try:
                    medIndx = medium.index(fixedMedia)
                    density = np.where(media == medNr[medIndx],
                    fixedDensity, density)
                except ValueError:
                    pass
            except ValueError:
                pass
            except IndexError:
                pass

    # rotate phantom and interchange rds_mesh with rds_ymesh if the orientation calls for it
    if orientCT[0] == 0:
        density = density[:, :, ::-1]
        media = media[:, :, ::-1]
        tmp = rds_ymesh[:]
        rds_ymesh = rds_xmesh[:]
        rds_xmesh = tmp[:]
        density = np.transpose(density, (2, 1, 0))
        density = np.rot90(density, 3)
        density = np.transpose(density, (2, 1, 0))
        media = np.transpose(media, (2, 1, 0))
        media = np.rot90(media, 3)
        media = np.transpose(media, (2, 1, 0))

    # write to file
    estepe = [0.25] * len(medium)  # dummy variable for ESTEPE
    if cv.spaceDelimit.lower() == 'y':
        cv.spaceDelimit = True
    else:
        cv.spaceDelimit = False
    print 'Writing egs4phant file'
    CTCtools.writeEgsphant(fileName, rds_xmesh, rds_ymesh, rds_zmesh,
    medium, estepe, media, density, cv.spaceDelimit)


# Function chooser
func_arg = {"-main": main}
# run specifc function as called from command line
if __name__ == "__main__":
    if sys.argv[1] == "-main":
        cnt = 6
        addStructType = []
        addRampName = []
        if len(sys.argv) == 4:
            DICOMdir = sys.argv[2]
            fileName = sys.argv[3]
            # generate RPRD, RS, CT file lists
            with open(confFile) as f:
                confCont = f.readlines()
            whatToGet = ['common.RPfilePrefix', 'common.RDfilePrefix',
            'common.RSfilePrefix', 'common.CTfilePrefix',
            'common.DICOMfileEnding']
            affix = struct()
            affix = CTCtools.getConfVars(affix, whatToGet, confCont)
            CTCtools.genRPRD(DICOMdir, affix.RPfilePrefix, affix.RDfilePrefix,
                affix.DICOMfileEnding, fileName)
            CTCtools.genRS(DICOMdir, affix.RSfilePrefix, affix.DICOMfileEnding,
                fileName)
            CTCtools.genCT(DICOMdir, affix.CTfilePrefix, affix.DICOMfileEnding,
                fileName)
            RPRDfile = os.path.sep.join([DICOMdir, 'RPRD.txt'])
            RSfile = os.path.sep.join([DICOMdir, 'RS.txt'])
            CTfile = os.path.sep.join([DICOMdir, 'CT.txt'])
        else:
            RPRDfile = sys.argv[2]
            RSfile = sys.argv[3]
            CTfile = sys.argv[4]
            fileName = sys.argv[5]
        if len(sys.argv) > cnt:
            while len(sys.argv) > cnt:
                addStructType.append(sys.argv[cnt])
                addRampName.append(sys.argv[cnt + 1])
                cnt += 2
        func_arg[sys.argv[1]](RPRDfile, RSfile, CTfile, fileName, addStructType,
             addRampName)
CTCtools.py0000664000175100017510000006552713165066111011571 0ustar  mcqamcqa# Toolbox for CTC_auto
#
#    Copyright [2016] [Rickard Cronholm] Licensed under the
#    Educational Community License, Version 2.0 (the "License"); you may
#    not use this file except in compliance with the License. You may
#    obtain a copy of the License at
#
#http://www.osedu.org/licenses/ECL-2.0
#
#    Unless required by applicable law or agreed to in writing,
#    software distributed under the License is distributed on an "AS IS"
#    BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
#    or implied. See the License for the specific language governing
#    permissions and limitations under the License.
#

import glob
import scipy.io
import dicom
import numpy as np
import copy
from progressbar import ProgressBar
from scipy import ndimage, ogrid, mgrid, interpolate
from skimage import draw
from scipy.interpolate import InterpolatedUnivariateSpline
import sys


class struct:

    def __init__(self):
        pass


def getConfVars(cv, wtg, content):
    for item in wtg:  # loop over what to get
        for elem in content:  # loop over content
            if elem.startswith(item):  # if found
                name = elem.split()[0].split('.')[-1]
                val = elem.split()[1]
                try:
                    val = float(val)
                except ValueError:
                    pass
                setattr(cv, name, val)  # add to confvariables
    return cv


def genRPRD(dicomDir, RPFilePrefix, RDFilePrefix, DICOMFileEnding, fileName):
    content = []
    # get full paths
    rp = getFromGlob(dicomDir, RPFilePrefix, DICOMFileEnding)
    rd = getFromGlob(dicomDir, RDFilePrefix, DICOMFileEnding)
    # create list
    content.append(rp[0].strip())
    content.append(rd[0].strip())
    # write file
    writeContFile(dicomDir, 'RPRD.txt', content)


def genRS(dicomDir, RSFilePrefix, DICOMFileEnding, fileName):
    content = []
    # get full paths
    rs = getFromGlob(dicomDir, RSFilePrefix, DICOMFileEnding)
    # create list
    content.append(rs[0].strip())
    # write file
    writeContFile(dicomDir, 'RS.txt', content)


def genCT(dicomDir, CTFilePrefix, DICOMFileEnding, fileName):
    content = []
    # get full paths
    ct = getFromGlob(dicomDir, CTFilePrefix, DICOMFileEnding)
    # create list
    for item in ct:
        content.append(item.strip())
    # write file
    writeContFile(dicomDir, 'CT.txt', content)


def getFromGlob(dicomDir, prefix, suffix):
    return glob.glob(''.join([dicomDir, prefix, '*', suffix]))


def writeContFile(dicomDir, fileName, content):
    f = open(''.join([dicomDir, fileName]), 'w')
    for element in content:
        f.write('{0}\n'.format(element))


def getStructFromMatFile(fileName, attr, level):
    data = struct()
    matfile = scipy.io.loadmat(fileName, squeeze_me=True,
    struct_as_record=False)
    if level == 0:
        for field in matfile[attr]._fieldnames:
            val = getattr(matfile[attr], field)
            setattr(data, field, val)

    if level == 1:
        for elem in matfile[attr]:
            for field in elem._fieldnames:
                try:
                    cont = []
                    val = getattr(data, field)
                    cont.append(val)
                    val = getattr(elem, field)
                    cont.append(val)
                    setattr(data, field, cont)
                except AttributeError:
                    val = getattr(elem, field)
                    setattr(data, field, val)

    return data


def getDICOMcoords(fileName, cm):
    dcm = dicom.read_file(fileName)  # read DICOM file

    # get directional cosines
    xDir = map(int, dcm.ImageOrientationPatient[:3])
    yDir = map(int, dcm.ImageOrientationPatient[3:])

    # get xCoord
    if abs(xDir[0]) == 1:
        xCoord = getCoord(sum(xDir), float(dcm.ImagePositionPatient[0]),
        float(dcm.PixelSpacing[0]), dcm.Columns)
    else:
        xCoord = getCoord(sum(xDir), float(dcm.ImagePositionPatient[1]),
        float(dcm.PixelSpacing[1]), dcm.Columns)
    xCoord = np.asarray(xCoord)

    # get yCoord
    if abs(yDir[1]) == 1:
        yCoord = getCoord(sum(yDir), float(dcm.ImagePositionPatient[1]),
        float(dcm.PixelSpacing[1]), dcm.Rows)
    else:
        yCoord = getCoord(sum(yDir), float(dcm.ImagePositionPatient[0]),
        float(dcm.PixelSpacing[0]), dcm.Rows)
    yCoord = np.asarray(yCoord)

    # if cm is 1 convert to cm from DICOM native mm
    if cm:
        xCoord[:] = [x / 10 for x in xCoord]
        yCoord[:] = [x / 10 for x in yCoord]

    return xCoord, yCoord


def getCoord(dircos, position, spacing, elements):
    coord = np.linspace(position, position + (elements-1) * dircos * spacing, num = elements)
    return coord


def getDICOMzCoord(fileName, cm):
    dcm = dicom.read_file(fileName)  # read DICOM file
    if dcm.Modality == 'CT':
        zCoord = dcm.ImagePositionPatient[2]
    else:
        zCoord = np.asarray(dcm.GridFrameOffsetVector)
        zCoord[:] = [x + dcm.ImagePositionPatient[2] for x in zCoord]

    # if cm is 1 convert to cm from DICOM native mm
    if cm:
        zCoord[:] = [x / 10 for x in zCoord]

    return zCoord


def getCTinfo(fileList, cm):
    # define static variables
    ctOffset = 1000
    zCoord = np.zeros(len(fileList))

    # loop through fileList to get individual z coord
    for i in range(0, len(fileList)):
        zCoord[i] = getDICOMzCoord(fileList[i], False)  # get z coord in mm
    if cm:
        zCoord[:] = [x / 10 for x in zCoord]  # convert to cm

    # sort by ascending
    ct_order = np.argsort(zCoord)
    # arrange zCoord
    zCoordArr = np.zeros(len(zCoord))
    for i in range(0, len(ct_order)):
        zCoordArr[i] = zCoord[ct_order[i]]

    # get CT data
    for i in range(0, len(fileList)):
        ct_slice = getCTdata(fileList[ct_order[i]])
        if i == 0:

            ct_mtrx = np.empty((len(zCoord), len(ct_slice), len(ct_slice[0])))
        ct_mtrx[:][:][i] = ct_slice

    # apply ctOffset
    ct_mtrx[:] = [x + ctOffset for x in ct_mtrx]

    return zCoordArr, ct_mtrx


def getCTdata(fileName):
    dcm = dicom.read_file(fileName)  # read DICOM file
    ct_slice = np.ndarray.astype(dcm.pixel_array, 'int')
    ct_slice[:] = [x * int(dcm.RescaleSlope) + int(dcm.RescaleIntercept)
    for x in ct_slice]
    return ct_slice


def getOrient(fileName):
    dcm = dicom.read_file(fileName)  # open file
    return map(int, dcm.ImageOrientationPatient)


def create1dDICOMcoord(start, pix, dim, direction):
    stop = start + direction * pix * (dim - 1)
    return np.linspace(start, stop, dim)


def deInterpCTmatrix(ct_mtrx, ct_xmesh, ct_ymesh, ct_zmesh, rd_xmesh, rd_ymesh, rd_zmesh, method):
    pbar = ProgressBar()
    mtrx = np.empty((len(rd_zmesh), len(rd_ymesh), len(rd_xmesh)))
    grid_x, grid_y = np.mgrid[rd_xmesh[0]:rd_xmesh[-1]:len(rd_xmesh) * 1j,
    rd_ymesh[0]:rd_ymesh[-1]:len(rd_ymesh) * 1j]

    points = np.empty((len(ct_xmesh) * len(ct_ymesh), 2))
    cnt = 0
    for elem in ct_xmesh:
        for item in ct_ymesh:
            points[cnt][0] = elem
            points[cnt][1] = item
            cnt += 1

    for i in pbar(range (0, len(rd_zmesh))):
        sliceNr = np.argwhere(np.around(ct_zmesh, decimals=1) ==
        np.around(rd_zmesh[i], decimals=1))[0][0]
        ct_slice = np.reshape(ct_mtrx[sliceNr][:][:].T, -1)

        mtrx[i][:][:] = interpolate.griddata(points, ct_slice,
        (grid_x, grid_y), method).T

    return mtrx


def getContour(RSData, zmesh, flip, cm):
# create and initialize structures
    contour = [None] * len(zmesh)

    for i in range(0, len(RSData.ContourSequence)):
        rawCont = map(float, RSData.ContourSequence[i].ContourData)
        if cm:
            rawCont[:] = [x / 10 for x in rawCont]  # convert to cm
        xCont = np.zeros(len(rawCont) / 3)
        yCont = np.zeros(len(rawCont) / 3)
        zCont = np.zeros(len(rawCont) / 3)
        for j in range(1, len(xCont) + 1):
            xCont[j - 1] = rawCont[3 * j - 3]
            yCont[j - 1] = rawCont[3 * j - 2]
            zCont[j - 1] = rawCont[3 * j - 1]
        # add the first element to form closed loop
        xCont = np.append(xCont, xCont[0])
        yCont = np.append(yCont, yCont[0])
        zCont = np.append(zCont, zCont[0])

        if flip == 0:
            cont = np.vstack((xCont, yCont, zCont))
        else:
            cont = np.vstack((yCont, xCont, zCont))

        # locate which slice nr it belongs to
        sliceNr = np.argwhere(np.around(zmesh, decimals=1) ==
        np.around(cont[2][0], decimals=1))[0][0]

        # the first data set will implicitly belong to a new slice
        if i == 0 or sliceNr != lastSlice:
            segments = []
        points = cont[:]
        segments.append(points)
        contour[sliceNr] = segments[:]

        lastSlice = sliceNr

    return contour


def getExtremeOfContour(data, extremes):
    # stack data
    x = np.empty(0)
    y = np.empty(0)
    z = np.empty(0)
    for i in range(0, len(data)):
        if data[i] is not None:
            for j in range(0, len(data[i])):
                x = np.hstack((x, data[i][j][0]))
                y = np.hstack((y, data[i][j][1]))
                z = np.hstack((z, data[i][j][2]))

    # compare extremes
    extremes.X = compExtremes(extremes.X, x)
    extremes.Y = compExtremes(extremes.Y, y)
    extremes.Z = compExtremes(extremes.Z, z)

    return extremes


def compExtremes(ext, arr):
    if ext[0] is None:
        ext[0] = min(arr)
    else:
        ext[0] = min(ext[0], min(arr))
    if ext[1] is None:
        ext[1] = max(arr)
    else:
        ext[1] = max(ext[1], max(arr))
    return ext


def getCorrContSeq(seq, lookfor):
    search = 0
    while True:
        if int(seq[search].ReferencedROINumber) == lookfor:
            break
        search += 1

    return search


def extendMesh(mesh, ext):
    meshMin = min(mesh)
    meshMax = max(mesh)
    step = np.unique(np.around(np.diff(mesh), decimals=3))[0]
    low = 0
    high = 0
    if ext[0] < meshMin:
        while meshMin > ext[0]:
            meshMin -= step
            low += 1
    if ext[1] > meshMax:
        while meshMax < ext[1]:
            meshMax += step
            high += 1

    if high > 0:
        highMesh = np.arange(1, high + 1)
        highMesh = highMesh * step + mesh[-1]
        mesh = np.hstack([mesh, highMesh])
    if low > 0:
        lowMesh = np.arange(low, 0, -1)
        lowMesh = lowMesh * -step + mesh[0]
        mesh = np.hstack([lowMesh, mesh])

    return mesh


def interpStructToDose(contour, rd_x, rd_y, rd_z, ct_x, ct_y, ct_z):
    # start by generating a mask based on CT grid
    maskM = np.zeros((len(ct_z), len(ct_y), len(ct_x)))
    xL = np.arange(len(ct_x))
    yL = np.arange(len(ct_y))

    # iterate through contour to generate a mask for each slice
    # note; there may be empty slices as well as slices with multiple segments
    pbar = ProgressBar()
    for j in pbar(range(0, len(contour))):
        if contour[j] is not None:  # skip empty slices
            sliceNr = np.argwhere(np.around(ct_z, decimals=1) ==
            np.around(contour[j][0][2][0], decimals=1))[0][0]
            for i in range(0, len(contour[j])):
                # take care of unbound contours
                points_x = copy.deepcopy(contour[j][i][0][:])  # get x-points
                points_y = copy.deepcopy(contour[j][i][1][:])  # get y-points
                points_x[np.where(points_x < min(ct_x))] = min(ct_x)
                points_x[np.where(points_x > max(ct_x))] = max(ct_x)
                points_y[np.where(points_y < min(ct_y))] = min(ct_y)
                points_y[np.where(points_y > max(ct_y))] = max(ct_y)
                # interpolate to voxel number instead of absolute coord
                points_x = np.interp(points_x, ct_x, xL)
                points_y = np.interp(points_y, ct_y, yL)
                # assign mask based on polygons
                [rr, cc] = draw.polygon(np.asarray(points_y),
                np.asarray(points_x), (len(ct_y), len(ct_x)))
                tempMask = np.zeros((len(ct_y), len(ct_x)))
                tempMask[rr, cc] = -1
                # add for the current slice
                maskM[sliceNr][:][:] = np.add(maskM[sliceNr][:][:],
                tempMask[:][:])
                maskM[sliceNr][:][:] = np.abs(maskM[sliceNr][:][:])

    del tempMask
    maskM = np.where(maskM > 1, 1, maskM)  # remove duplicates

    # interpolate maskM on to the RD grid
    order = 1
    mask = map_coordinates(maskM, ct_x, ct_y, ct_z, rd_x, rd_y, rd_z, order)
    return np.around(mask).astype(int)


def cropCT(mtrx, ct_x, ct_y, ct_z, rd_x, rd_y, rd_z):
    margin = 5
    # x
    ind = getIndx(ct_x, rd_x, margin)
    ct_x = ct_x[ind[0]:ind[1] + 1]
    mtrx = mtrx[:, :, ind[0]:ind[1] + 1]
    # y
    ind = getIndx(ct_y, rd_y, margin)
    ct_y = ct_y[ind[0]:ind[1] + 1]
    mtrx = mtrx[:, ind[0]:ind[1] + 1, :]
    # z
    ind = getIndx(ct_z, rd_z, margin)
    ct_z = ct_z[ind[0]:ind[1] + 1]
    mtrx = mtrx[ind[0]:ind[1] + 1, :, :]

    return mtrx, ct_x, ct_y, ct_z


def getIndx(a, b, margin):
    ind = []
    if a[0] < a[-1]:
        ind.append((np.abs(a - b[0])).argmin() - margin)
        ind.append((np.abs(a - b[-1])).argmin() + margin)
    else:
        ind.append((np.abs(a - b[-1])).argmin() + margin)
        ind.append((np.abs(a - b[0])).argmin() - margin)

    ind = np.asarray(ind)
    ind = ind.clip(min=0)
    ind = ind.clip(max=len(a) - 1)

    return ind


def map_coordinates(mtrx, ct_x, ct_y, ct_z, rd_x, rd_y, rd_z, order):
    mtrxOut = np.zeros((len(rd_z), len(rd_y), len(rd_x)))  # preallocate

    # create sparse grids
    xct, yct = ogrid[min(ct_x):max(ct_x):len(ct_x) * 1j,
    min(ct_y):max(ct_y):len(ct_y) * 1j]

    # interpolate coord arrays in cm to voxel nr
    xL = np.linspace(0, len(ct_x) - 1, len(ct_x))
    yL = np.linspace(0, len(ct_y) - 1, len(ct_y))
    extrapolatorX = InterpolatedUnivariateSpline(ct_x, xL, k=1)
    x = extrapolatorX(rd_x)
    extrapolatorY = InterpolatedUnivariateSpline(ct_y, yL, k=1)
    y = extrapolatorY(rd_y)
    # create mesh grids
    xrd, yrd = mgrid[min(x):max(x):len(x) * 1j, min(y):max(y):len(y) * 1j]
    coords = np.array([yrd, xrd])

    pbar = ProgressBar()
    # map coordinates slice by slice
    cnt = 0
    for i in pbar(range (0, len(rd_z))):
        try:
            sliceNr = np.argwhere(np.around(ct_z,
            decimals=1) == np.around(rd_z[i], decimals=1))[0][0]
            ct_slice = ndimage.map_coordinates(mtrx[sliceNr][:][:],
            coords, order=order).T
        except IndexError:
            ct_slice = np.zeros((len(rd_y), len(rd_x)))
        mtrxOut[cnt][:][:] = ct_slice[:][:]
        cnt += 1

    return mtrxOut


def computeDensity(mtrx, ramp):
    # remove nans
    mtrx = np.nan_to_num(mtrx)
    densMtrx = np.zeros(mtrx.shape)
    # compute density for all but last 'segment'
    for i in range(0, len(ramp) - 1):
        densMtrx = np.where(mtrx <= ramp[i][0],
        mtrx * ramp[i][1][1] + ramp[i][1][0], densMtrx)
    # compute density for last 'segment'
    densMtrx = np.where(mtrx > ramp[-2][0], mtrx * ramp[-1][1][1] +
    ramp[-1][1][0], densMtrx)
    # remove negatives
    densMtrx = np.where(densMtrx < 0., 0., densMtrx)
    return densMtrx


def computeMedia(mtrx, structures, medium, medNr):
    # remove nans
    mtrx = np.nan_to_num(mtrx)
    medMtrx = np.zeros(mtrx.shape)
    # compute media for each structure in a loop
    for elem in structures:
        tempMtrx = np.zeros(mtrx.shape)
        names = []
        bounds = []
        # loop through appending names and bounds
        for i in range(0, len(elem.ramp)):
            names.append(elem.ramp[i][0])
            bounds.append(elem.ramp[i][1])
        # sort according to bounds
        indx = [bounds.index(x) for x in sorted(bounds)]
        names[:] = [names[x] for x in indx]
        bounds[:] = [bounds[x] for x in indx]
        # find global medNr
        globalNr = medNr[medium.index(names[0])]
        tempMtrx = np.where(mtrx <= bounds[0], globalNr, tempMtrx)
        for i in range(1, len(elem.ramp)):
            globalNr = medNr[medium.index(names[i])]
            tempMtrx = np.where(mtrx > bounds[i - 1], globalNr, tempMtrx)

        medMtrx = np.where(elem.logicMatrix == 1, tempMtrx, medMtrx)
    return medMtrx.astype('int')


def convMatToPyth(indata):
    outData = []
    materials = []
    mat_ct_up_bound = []
    if type(indata.materials) == type(unicode('A')):
        materials.append(indata.materials)
        mat_ct_up_bound.append(indata.mat_ct_up_bound)
    else:
        for i in range(0, len(indata.materials)):
            materials.append(indata.materials[i])
            mat_ct_up_bound.append(indata.mat_ct_up_bound[i])
    materials = map(str, materials)
    outData.append(map(str.strip, materials))
    outData.append(mat_ct_up_bound)
    return outData


def buildGlobalMediaList(structures):
    # loop over structures to get all names
    medNames = []
    for i in range(0, len(structures)):
        for j in range(0, len(structures[i].ramp)):
            medNames.append(structures[i].ramp[j][0].strip())
    medNames = flatten(medNames)
    medNames = list(set(medNames))

    # sort medium based on average upper ct bound
    summation = [0] * len(medNames)
    average = [0] * len(medNames)
    for i in range(0, len(medNames)):
        cnt = 0
        for j in range(0, len(structures)):
            for k in range(0, len(structures[j].ramp)):
                try:
                    indx = structures[j].ramp[k][0].index(medNames[i])
                    summation[i] += structures[j].ramp[k][1]
                    cnt += 1
                except ValueError:
                    pass
        average[i] = float(summation[i]) / cnt

    sortIndx = np.ndarray.tolist(np.argsort(np.asarray(average)))
    medium = [medNames[x] for x in sortIndx]

    # generate medNr list
    medNr = np.linspace(1, len(medium), len(medium)).astype('int')
    medNr = np.ndarray.tolist(medNr)

    return medNr, medium


def writeEgsphant(fileName, x, y, z, medium, estepe, media, density, spaceDelimit):
    halfThick = False
    xb = createBoundGrid(x)
    yb = createBoundGrid(y)
    zb = createBoundGrid(z)

    if halfThick:
        # "interpolate" to half slice thickness
        zn = np.linspace(z[0], z[-1], 2 * len(z) - 1)
        z = copy.deepcopy(zn)
        zb = createBoundGrid(z)

    f = open(''.join([fileName, '.egs4phant']), 'w')
    # Write number of media and their names
    f.write('{0:d}\n'.format(len(medium)))
    for elem in medium:
        f.write('{}\n'.format(elem))
    for elem in estepe:
        f.write('{:6.2f}'.format(elem))  # write estepe
    f.write('\n')
    # write num elements and their boundaries
    f.write('{0:5d}{1:5d}{2:5d}'.format(len(x), len(y), len(z)))
    # write x bound
    writeBound(f, xb, 6)
    # write y bound
    writeBound(f, yb, 6)
    # write z bound
    writeBound(f, zb, 6)
    f.write('\n')
    # write media matrix
    pbar = ProgressBar()
    if spaceDelimit:
        for i in pbar(range(0, len(z))):
            for j in range(0, len(y)):
                for k in range(0, len(x)):
                    if not halfThick:
                        # for normal thick
                        f.write('{0:d} '.format(media[i][j][k]))
                    else:
                        # for half thick
                        f.write('{0:d} '.format(media[int(i / 2)][j][k]))
                f.write('\n')
            f.write('\n')
    else:
        for i in pbar(range(0, len(z))):
            for j in range(0, len(y)):
                for k in range(0, len(x)):
                    if not halfThick:
                        # for normal thickness
                        f.write('{0:d}'.format(media[i][j][k]))
                    else:
                        # for half thick
                        f.write('{0:d}'.format(media[int(i / 2)][j][k]))
                f.write('\n')
            f.write('\n')
    # write density matrix
    for i in range(0, len(z)):
        for j in range(0, len(y)):
            if not halfThick:
                writeBound(f, density[i][j][:], 5)  # for normal thickness
            else:
                writeBound(f, density[int(i / 2)][j][:], 5)  # for half thick
        f.write('\n')

    # finish of with an emtpy line to satisfy dosxyznrc.mortran
    f.write('\n')  
    f.close()  # close file


def writeBound(f, grid, num):
    for i in range(0, len(grid)):
        if np.mod(i, num) == 0:
            f.write('\n')
        f.write('{0:12.8f} '.format(grid[i]))


def createBoundGrid(grid):
    step = np.unique(np.around(np.diff(grid), decimals=1))[0]
    return np.linspace(grid[0] - step / 2, grid[-1] + step / 2, len(grid) + 1)


def getFromFile(fileName, switch):
    myList = []
    with open(fileName) as f:
        data = f.readlines()

    data = map(str.strip, data)

    if switch == 0:
        for elem in data:
            try:
                myList.append(map(float, elem.split()))
            except ValueError:
                pass
    elif switch == 1:
        il = []
        for elem in data:
            try:
                map(float, elem.split())
                innerList = copy.deepcopy(il)
                innerList.append(float(elem.split()[0]))
                innerList.append(np.asarray(map(float, elem.split()[1:])))
                myList.append(innerList)
            except ValueError:
                pass
    elif switch == 2:
        il = []
        for elem in data:
            try:
                int(elem.split()[1])
                innerList = copy.deepcopy(il)
                innerList.append(elem.split()[0])
                innerList.append(int(elem.split()[1]))
                myList.append(innerList)
            except ValueError:
                pass
    elif switch == 3:
        for elem in data:
            try:
                elems = filter(None, elem.split('\t'))
                myList.append([elems[2].strip(), elems[0].strip(),
                float(elems[1].strip())])
            except ValueError:
                pass
            except IndexError:
                pass

    return myList


def grabData(fileName, attr, level, switch):
    data = getFromFile(fileName, switch)
    if len(data) == 0:
        data = getStructFromMatFile(fileName, attr, level)
        data = getPythonicList(data, attr)

    return data


def getPythonicList(data, attr):
    myList = []
    tempList = []
    il = []
    if attr == 'materialRamp':
        lookFor = ['materials', 'mat_ct_up_bound']

    else:
        lookFor = ['upperCTbound', 'coefficients']

    for i in range(0, len(lookFor)):
        tempList.append(getattr(data, lookFor[i]))

    # reaarrange
    for i in range(0, len(tempList[0])):
        innerList = copy.deepcopy(il)
        try:
            innerList.append(float(tempList[0][i]))
        except ValueError:
            innerList.append(str(tempList[0][i]).strip())
        innerList.append(tempList[1][i])
        myList.append(innerList)

    return myList


def flatten(seq, container=None):
    if container is None:
        container = []
    for s in seq:
        if hasattr(s, '__iter__'):
            flatten(s, container)
        else:
            container.append(s)
    return container


def getHUfromDens(dens, densRamp):
    # iteratively compare against boundaries in densRamp
    cnt = 0
    while cnt < len(densRamp):
        myHU = (dens - densRamp[cnt][1][0]) / densRamp[cnt][1][1]
        if myHU <= densRamp[cnt][0]:
            break
        cnt += 1

    return myHU


def checkEngulf(child, parent):
    if isEngulfed(parent.logicMatrix, child.logicMatrix):
        return parent, child
    else:
        #print '{0:s} is engulfed by {1:s}, switching order'.format(child.name, parent.name)
        return child, parent


def isEngulfed(A, B):
    A = np.reshape(A, -1)
    B = np.reshape(B, -1)
    # return True if all points in A exists in B
    indx = np.where(A == 1)
    check = np.unique(B[indx] == 1)
    if len(check) > 1 or len(check) == 0:
        check = False
    else:
        check = True

    return check


def covered(A, B):
    A = np.reshape(A, -1)
    B = np.reshape(B, -1)
    C = A + B
    return float(len(np.where(C == 2)[0])) / float(len(np.where(A == 1)[0]))


def sortStructures(sIn, order):
    sOut = []
    for item in order:
        if isinstance(item, basestring):
            match = [i for i, x in enumerate(sIn) if x.type == item]
        else:
            match = [i for i, x in enumerate(sIn) if hasattr(x, 'RelElecDens')]
        # append
        for elem in match:
            sOut.append(sIn[elem])
        # remove from sIn
        match.reverse()
        [sIn.pop(x) for x in match]

    # append the rest
    for elem in sIn:
        sOut.append(elem)

    return sOut


def belongToOne(structures, extName):
    for i in range(0, len(structures) - 1):
        #print structures[i].name, '\n----------------'
        for j in range(i + 1, len(structures)):
            #printnamevol(structures[i])
            #printnamevol(structures[j])
            parent = copy.deepcopy(structures[i])
            child = copy.deepcopy(structures[j])
            if parent.type == extName or child.type == extName or parent.name == 'OUTSIDE':
                 parent, child = checkEngulf(child, parent)
            indx = []
            indx = np.where(parent.logicMatrix == child.logicMatrix)
            numPts = []
            numPts = np.where(parent.logicMatrix[indx] == 1)
            numPts = np.asarray(numPts)
            if len(np.ravel(numPts)) > 0:
                child.logicMatrix = np.where(parent.logicMatrix == child.logicMatrix, 0 , child.logicMatrix)
                #print 'removing {2:5d} mutual points in {1:s} and {0:s} from {0:s}'.format(child.name, parent.name, len(np.ravel(numPts)))
                if parent.name == structures[i].name:
                    structures[i] = copy.deepcopy(parent)
                    structures[j] = copy.deepcopy(child)
                else:
                    structures[j] = copy.deepcopy(parent)
                    structures[i] = copy.deepcopy(child)
            #printnamevol(structures[i])
            #printnamevol(structures[j])
    return structures


def dropEmpty(sIn):
    sOut = []
    for elem in sIn:
        if np.sum(elem.logicMatrix) > 0:
            sOut.append(elem)
    return sOut


def printnamevol(s):
    print 'Structure: {0:s} has est volume {1:5d}'.format(s.name, int(np.sum(s.logicMatrix)))



def checkSaveTypes(toCheck, myList):
    try:
        myList.index(toCheck)
        return True
    except ValueError:
        return False


def dropOutOfExt(sIn, extName, saveTypes):
    sOut = []
    ext = [i for i, x in enumerate(sIn) if x.type == extName][0]
    external = sIn[ext]
    for elem in sIn:
        if not elem.name == external.name and not checkSaveTypes(elem.type, saveTypes):
            # check if it is engulfed
            #print '{0:s} covered to {2:.2f}% by {1:s}'.format(elem.name, external.name, covered(elem.logicMatrix, external.logicMatrix)*100.) 
            if covered(elem.logicMatrix, external.logicMatrix) > .95:
                sOut.append(elem)
                #print '{0:s} appended'.format(elem.name)
            else:
                print '{0:s} dropped'.format(elem.name)
        else:
            sOut.append(elem)
            #print '{0:s} appended'.format(elem.name)
    return sOut
SUS.fixedMedDens0000664000175100017510000000012113047337410012435 0ustar  mcqamcqaLook up table for media awith fixed densities
AU700ICRU	19.32
TI700ICRU	4.430
SUS.relElec0000664000175100017510000000074313055533423011464 0ustar  mcqamcqaRelElecDens to MassDens lookup table
RelElecDens	MassDens
1.0000		1.0000	H2O700ICRU
1.0020		1.0000	hu3
0.5000		0.4799	H2O700ICRU
0.0500		0.0275	Rohacell51A
0.7000		0.7109	170C700ICRU
0.0670		0.0500	AIR521ICRU
0.2120		0.2000	LUNG521ICRU
0.4050		0.4000	LUNG521ICRU
0.5980		0.6000	LUNG521ICRU
0.7910		0.8000	hu3
0.8880		0.9000	hu3
0.9990		0.9999	hu3
1.1763		1.2000	hu5
1.3189		1.4000	hu9
1.4614		1.6000	hu11
1.5330		1.7000	hu14
1.1053		1.1000	hu3
1.2478		1.3000	hu7
cirsLund.dat0000664000175100017510000000013713164412274011770 0ustar  mcqamcqaUpper CT limit	Coefficients
988.9	-0.01939081	0.00103581
1867.4  0.35173975      0.00067349
exactIGRT.dat0000664000175100017510000000013413164412350011767 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
Rohacell51A	100
170C700ICRU     1000
justAir700.dat0000664000175100017510000000010213164412370012042 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
AIR700OUT	5000
justWater700.dat0000664000175100017510000000010313164412401012405 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
H2O700ICRU	5000
onlyWater.dat0000664000175100017510000000010313164412411012153 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
H2O700ICRU	5000
21mediaNoCouch+AU.dat0000664000175100017510000000045213164412312013200 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
AIR521ICRU	50
LUNG521ICRU	587
hu3		1015
hu4		1129
hu5		1200
hu6		1300
hu7		1400
hu8		1500
hu9		1600
hu10		1700
hu11		1800
hu12		1900
hu13		2000
hu14		2100
hu15		2200
hu16		2300
hu17		2400
hu18		2500
hu19		2600
hu20		2700
hu21		2800
AU700ICRU	5000
21media+TIT.dat0000664000175100017510000000045213164412325012060 0ustar  mcqamcqaData for media segmentation
MediaName	upperCTbound
AIR521ICRU	50
LUNG521ICRU	587
hu3		1015
hu4		1129
hu5		1200
hu6		1300
hu7		1400
hu8		1500
hu9		1600
hu10		1700
hu11		1800
hu12		1900
hu13		2000
hu14		2100
hu15		2200
hu16		2300
hu17		2400
hu18		2500
hu19		2600
hu20		2700
hu21		2800
TI700ICRU	5000