limfomiPET/pythonScripts/runStat.py

import statUtils
import os
import radiomics
import SimpleITK
import sys
import json
import numpy

def main(parFile='../templates/statistics.json'):


    featureExtractor=radiomics.featureextractor.RadiomicsFeatureExtractor
 
    setup=statUtils.loadSetup(parFile)
    rFile='radiomics.json'

    forceLiver=setup.get('forceLiver',False)
    forceSUVMax=setup.get('forceSUVMax',False)
    forceLiver15=setup.get('forceLiver15',False)
    forceGlobal=setup.get('forceGlobal',False)
    forceSUV4=setup.get('forceSUV4',False)
    forceSUV4=setup.get('forceSUV4',False)
    forceSUVMaxLesions=setup.get('forceSUVMaxLesions',False)
    forceSUVMaxMetastases=setup.get('forceSUVMaxMetastases',False)

    petField=setup.get('petField','petResampled')

    doUpload=setup.get('doUpload',True)

    #update threshold values if needed
    with open(rFile,'w') as f:
        f.write(json.dumps(setup['radiomics']))
    setup['db'],setup['fb']=statUtils.connectDB('onko-nix')
    users=statUtils.getUsers(setup['db'],setup['project'])
    p=setup['project']
    qFilter=[]
    try:
        vList=';'.join(setup['participants'])    
        qFilter.append({'variable':'ParticipantId','value':vList,'oper':'in'})
    except KeyError:
        pass

    try:
        vList=';'.join(setup['visits'])
        qFilter.append({'variable':'visitCode','value':vList,'oper':'in'})
    except KeyError:
        pass

    ds=setup['db'].selectRows(p,'study',setup['imagingDataset'],qFilter)
    if not os.path.isdir(setup['localDir']):
        os.mkdir(setup['localDir'])
    #select just the first row; debugging
    rows=ds['rows']
    setup['values']=['COM','MTV','TLG','SUVmean','SUVmax','voxelCount','SUVSD']
    setup['featureExtractor']=featureExtractor(rFile)

    #make sure we calculate Variance
    firstOrder=setup['radiomics']['featureClass']['firstorder']
    if 'Variance' not in firstOrder:
        firstOrder.append('Variance')

    n=setup.get('n',-1)
    if n>0:
        rows=rows[:n]

    for r in rows:

        setup['SUVdataset']='SUVanalysis'
        globalDone=not forceGlobal and checkData(setup,r)

        #check if we have to do calculation
        setup['SUVdataset']='SUVanalysis_liver'
        liverDone=not forceLiver and checkData(setup,r)

        setup['SUVdataset']='SUVanalysis_SUVmax'
        suvMaxDone=not forceSUVMax and checkData(setup,r) 

        setup['SUVdataset']='SUVanalysis_SUVmaxMetastases'
        suvMaxMetastasesDone=not forceSUVMaxMetastases and checkData(setup,r)

        setup['SUVdataset']='SUVanalysis_SUVmaxLesions'
        suvMaxLesionsDone=not forceSUVMaxLesions and checkData(setup,r)

        setup['SUVdataset']='SUVanalysis_liver1p5'
        liver1p5Done=not forceLiver15 and checkData(setup,r)

        setup['SUVdataset']='SUVanalysis_SUV4'
        suv4Done=not forceSUV4 and checkData(setup,r)

        doneCode=f'({globalDone}/{liverDone}/{liver1p5Done}/{suvMaxDone}/{suv4Done}/{suvMaxLesionsDone}/{suvMaxMetastasesDone})'
        codeDesc='(global/liver/liver1p5/suvMax/suv4/suvMaxLesions/suvMaxMetastases)'
        idVar='ParticipantId'
        pid=r[idVar]
        visitCode=r['visitCode']
        print(f'Done[{pid}/{visitCode}: {codeDesc}: {doneCode}')

        if globalDone and liverDone and suvMaxDone and \
                liver1p5Done and suv4Done and \
                suvMaxLesionsDone and suvMaxMetastasesDone:
            print('Skipping {} {}'.format(r['ParticipantId'],r['visitCode']))
            continue

        #PET
        for q in [petField]:
            localPath=statUtils.getImage(setup,r,q)
        if localPath=="NONE":
            continue
        pet=SimpleITK.ReadImage(localPath)

        #Seg
        segPaths=statUtils.getSegmentations(setup,r)
        if "NONE" in segPaths.values():
            os.remove(localPath)
            continue
        segKeys=list(segPaths.keys())
        for x in segPaths:
            print('Loaded {}/{}'.format(users[x],segPaths[x]))

        seg={x:SimpleITK.ReadImage(segPaths[x]) for x in segPaths}
        #check dimensions of seg and pet
        if not matchingDimensions(pet,seg):
            continue

        globalThreshold={x:0 for x in seg}
        if not globalDone:
            #get value for maximum in organs or liver mean and std
            setup['SUVdataset']='SUVanalysis'
            globalOutputs=thresholdAnalysis(setup,r,pet,seg,globalThreshold)
            if doUpload:
                uploadData(setup,r,globalOutputs)

        default={'SUVmax':0,'SUVmean':0,'SUVSD':0}
        setup['SUVdataset']='SUVanalysis'
        outputs=loadGlobals(setup,r,seg)

        #this are IDs from labkey
        liverId=1
        lesionId=4
        bmId=3
        metastasesId=6
    
        #liver threshold
        if not liverDone:
            liverThreshold={x:outputs[x].get(liverId,default)['SUVmean']
                +2*outputs[x].get(liverId,default)['SUVSD'] 
                for x in outputs}
            print(f'thr[liver]={liverThreshold}')
            setup['SUVdataset']='SUVanalysis_liver'
            liverOutputs=thresholdAnalysis(setup,r,pet,seg,liverThreshold)
            if doUpload:
                uploadData(setup,r,liverOutputs)

        #also for threshold=1.5
        if not liver1p5Done:
            liver1p5Threshold={x:1.5*outputs[x].get(liverId,default)['SUVmean']
                +2*outputs[x].get(liverId,default)['SUVSD'] 
                for x in outputs}
            print(f'thr[liver/1.5]={liver1p5Threshold}')
            setup['SUVdataset']='SUVanalysis_liver1p5'
            liver1p5Outputs=thresholdAnalysis(setup,r,pet,seg,liver1p5Threshold)
            if doUpload:
                uploadData(setup,r,liver1p5Outputs)

        #threshold=0.41 SUVmax
        if not suvMaxDone:
            idArr=[lesionId,bmId,metastasesId]
            fM=numpy.max
            suvMax={x:fM([outputs[x].get(y,default)['SUVmax'] for y in idArr])}
            suvMaxThreshold={x:0.41*suvMax[x] for x in suvMax}
            print(f'thr[SUVmax]={suvMaxThreshold}')

            setup['SUVdataset']='SUVanalysis_SUVmax'
            suvMaxOutputs=thresholdAnalysis(setup,r,pet,seg,suvMaxThreshold)
            if doUpload:
                uploadData(setup,r,suvMaxOutputs)

        #threshold=0.41 SUVMax (Lesion)
        if not suvMaxLesionsDone:
            suvMaxLesionsThreshold=\
                {x:0.41*outputs[x].get(lesionId,default)['SUVmax'] \
                    for x in outputs}
            print(f'thr[SUVmaxLesions]={suvMaxLesionsThreshold}')

            setup['SUVdataset']='SUVanalysis_SUVmaxLesions'
            suvMaxLesionsOutputs=thresholdAnalysis(setup,r,pet,seg,
                    suvMaxLesionsThreshold)
            if doUpload:
                uploadData(setup,r,suvMaxLesionsOutputs)


        #threshold=0.41 SUVMax (Metastases)
        if not suvMaxMetastasesDone:
            suvMaxMetastasesThreshold=\
                {x:0.41*outputs[x].get(metastasesId,default)['SUVmax'] \
                    for x in outputs}
            print(f'thr[SUVmaxMetastases]={suvMaxMetastasesThreshold}')

            setup['SUVdataset']='SUVanalysis_SUVmaxMetastases'
            suvMaxMetastasesOutputs=thresholdAnalysis(setup,r,pet,seg,\
                    suvMaxMetastasesThreshold)
            if doUpload:
                uploadData(setup,r,suvMaxMetastasesOutputs)




        if not suv4Done:
            suv4Threshold={x:4 for x in outputs}
            print(f'thr[SUV4]={suv4Threshold}')
            setup['SUVdataset']='SUVanalysis_SUV4'
            suv4Outputs=thresholdAnalysis(setup,r,pet,seg,suv4Threshold)
            if doUpload:
                uploadData(setup,r,suv4Outputs)

        
        #cleanup
        os.remove(localPath)

        for x in segPaths:
            os.remove(segPaths[x])


def loadGlobals(setup,r,seg):
    matchingVisits={
            'VISIT_0':['VISIT_0'],
            'VISIT_1':['VISIT_1','VISIT_11'],
            'VISIT_11':['VISIT_1','VISIT_11'],
            'VISIT_2':['VISIT_2','VISIT_12'],
            'VISIT_12':['VISIT_2','VISIT_12'],
            'VISIT_3':['VISIT_3','VISIT_13'],
            'VISIT_13':['VISIT_3','VISIT_13'],
            'VISIT_4':['VISIT_4','VISIT_14'],
            'VISIT_14':['VISIT_4','VISIT_14']}
    ids=[1,2,3,4,5,6]
    p=setup['project']
    d=setup['SUVdataset']
    idVar='ParticipantId'
    vCodes=';'.join(matchingVisits[r['visitCode']])
    vars=['MTV','TLG','SUVmean','SUVmax','SUVSD']
    outputs={}
    for x in seg:
        outputs[x]={}
        for i in ids:
            qFilter=[]
            qFilter.append({'variable':idVar,'value':r[idVar],'oper':'eq'})
            qFilter.append({'variable':'visitCode','value':vCodes,'oper':'in'})
            qFilter.append({'variable':'User','value':f'{x}','oper':'eq'})
            qFilter.append({'variable':'segment','value':f'{i}','oper':'eq'})
            ds=setup['db'].selectRows(p,'study',d,qFilter)
            outputs[x][i]={v:combineData(ds['rows'],v) for v in vars}
    return outputs

def combineData(rows,x):
    dt=numpy.array([float(r[x] or 0) for r in rows])
    #print('Combine data[{}] {}'.format(x,dt))
    r=0
    if x=='MTV' or x=='TLG':
        r=numpy.sum(dt)
    if x=='SUVmean':
        m=combineData(rows,'MTV')
        if m==0:
            r=0
        else:
            r=combineData(rows,'TLG')/m
    if x=='SUVmax':
        try:
            r=numpy.max(dt)
        except ValueError:
            #empty array
            r=0
    if x=='SUVSD':
        #this should be probably be done right, 
        #here I only want it to work also if one of the components is 0
        r=numpy.sqrt(numpy.sum(dt*dt))
    #print(f'Result: {r}')
    return r
    
def thresholdAnalysis(setup,r,pet,seg,thrs):
    #thresholds thrs are by segmentation author
    outputs={}
    for s in thrs:
        t=thrs[s]
        print(f'Thr[{s}]: {t}')
        #update radiomics setting
        setup['radiomics']['setting']['resegmentRange']=[thrs[s]]
        setup['radiomics']['setting']['resegmentShape']=True
        outputs[s]=getValuesForSegmentation(setup,r,pet,seg[s])
            
        _=[outputs[s][y].update({'threshold':thrs[s]}) for y in outputs[s]]
    print(outputs)
    return outputs


def getValues(setup,row,pet,seg):
    #seg is a list of segments

    rFile='radiomics.json'
    #short function names
    featureExtractor=radiomics.featureextractor.RadiomicsFeatureExtractor
    getStats=statUtils.getRadiomicsComponentStats

    with open(rFile,'w') as f:
        f.write(json.dumps(setup['radiomics']))

    setup['featureExtractor']=featureExtractor(rFile)


    #find labels associated with each (non-overlaping) segmentation
    ids=statUtils.getSegments(list(seg.values())[0])
    outputs={x:{} for x in seg}
    for x in seg:
        
        for id in ids:
            print('{} {}'.format(id,ids[id]))
            try:
                output=getStats(setup,pet,seg[x],ids[id])
            except ValueError:
                continue
            outputs[x][ids[id]]=output
        
    os.remove(rFile)
    return outputs

def getValuesForSegmentation(setup,row,pet,seg):

    #short function names
    featureExtractor=radiomics.featureextractor.RadiomicsFeatureExtractor
    getStats=statUtils.getRadiomicsComponentStats

    rFile='radiomics.json'

    with open(rFile,'w') as f:
        f.write(json.dumps(setup['radiomics']))

    setup['featureExtractor']=featureExtractor(rFile)


    #find labels associated with each (non-overlaping) segmentation
    ids=statUtils.getSegments(seg)
    outputs={}

    for id in ids:
        print('{} {}'.format(id,ids[id]))
        try:
            output=getStats(setup,pet,seg,ids[id])
        except ValueError:
            continue
        outputs[ids[id]]=output
        
    os.remove(rFile)
    return outputs


def uploadData(setup,r,outputs):
    baseVar=['ParticipantId','SequenceNum','patientCode','visitCode']
    p=setup['project']
    d=setup['SUVdataset']
    for x in outputs:
        for s in outputs[x]:
            output=outputs[x][s]
            output.update({x:r[x] for x in baseVar})
            output['User']=x
            output['segment']=s
            statUtils.updateDatasetRows(setup['db'],p,d,[output])
    
def checkData(setup,r):
    p=setup['project']
    d=setup['SUVdataset']
    qFilter=[]
    idVar='ParticipantId'
    qFilter.append({'variable':idVar,'value':r[idVar],'oper':'eq'})
    qFilter.append({'variable':'visitCode','value':r['visitCode'],'oper':'eq'})
    ds=setup['db'].selectRows(p,'study',d,qFilter)
    try:
        n=len(ds['rows'])
    except KeyError:
        #contraintuitevly return as if data were there
        #if ds has no rows, which typically happens
        #on failure to retrieve dataset, which is
        #related to missing datasets
        print(f'Dataset [{d}] not found. Returning OK')
        return True
    print('[{}:{}/{}] got {} rows.'.format(d,r[idVar],r['visitCode'],n))
    return n>0

def vEqual(a,b,eps=1e-8):
    delta=numpy.array(a)-numpy.array(b)
    aDelta=numpy.abs(delta)
    #print(f'Delta={delta} abs(Delta)={aDelta}')
    return all(aDelta<eps)

def matchingDimensions(pet,seg):
    #QA routine
    #check if dimensnions of pet and seg match
    #seg is a dict {'ID':sitkImage}
    
    
    #print('[PET] size={} origin={}'.format(pet.GetSize(),pet.GetOrigin()))
    sPET=pet.GetSize()
    ok=[vEqual(seg[x].GetSize(),sPET) for x in seg]
    if not all(ok):
        sz={x:seg[x].GetSize() for x in seg}
        print('Size mismatch PET[{}]: {}'.format(sPET,sz))
        return False
    return True

if __name__=='__main__':
    main(sys.argv[1])