limfomiPET/pythonScripts/runStat.py

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

def main(parFile='../templates/statistics.json'):
    setup=statUtils.loadSetup(parFile)
    rFile='radiomics.json'

    #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'])
    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(setup['project'],'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']
    #params=os.path.join('..','templates','radiomics.yaml')
    setup['featureExtractor']=radiomics.featureextractor.RadiomicsFeatureExtractor(rFile)

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

    for r in rows:

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

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

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

        
        if liverDone and suvMaxDone and liver1p5Done:
            print('Skipping {} {}'.format(r['ParticipantId'],r['visitCode']))
            continue
        doneCode=f'({liverDone}/{liver1p5Done}/{suvMaxDone})'
        print(f'Done: (liver/liver1p5/suvMax): {doneCode}')
        #PET
        for q in ['petResampled']:
            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}


        try:
            thr=setup['radiomics']['setting']['resegmentRange'][0]
        except KeyError:
            thr=None

        setup['radiomics']['setting']['resegmentRange']=None        
        firstOrder=setup['radiomics']['featureClass']['firstorder']
        if 'Variance' not in firstOrder:
            firstOrder.append('Variance')
        #get value for maximum in organs or liver mean and std
        outputs=getValues(setup,r,pet,seg)
        setup['SUVdataset']='SUVanalysis'
        #uploadData(setup,r,outputs)

        print(outputs)
        
        default={'SUVmax':0,'SUVmean':0,'SUVSD':0}
        #liver threshold
        liverId=1
        liverThreshold={x:outputs[x].get(liverId,default)['SUVmean']
                +2*outputs[x].get(liverId,default)['SUVSD'] 
                for x in outputs}
        liver1p5Threshold={x:outputs[x].get(liverId,default)['SUVmean']
                +1.5*outputs[x].get(liverId,default)['SUVSD'] 
                for x in outputs}
        lesionId=4
        bmId=3
        suvMax={x:numpy.max([outputs[x].get(lesionId,default)['SUVmax'],
                    outputs[x].get(bmId,default)['SUVmax']])
                for x in outputs}
        suvMaxThreshold={x:0.41*suvMax[x] for x in suvMax}
    

        print('thr[liver]={} thr[liver/1.5]={} thr(suvmax)={}'.format(liverThreshold,liver1p5Threshold,suvMaxThreshold))
       

        if not liverDone:
            setup['SUVdataset']='SUVanalysis_liver'
            liverOutputs=thresholdAnalysis(setup,r,pet,seg,liverThreshold)
            uploadData(setup,r,liverOutputs)

        #also for threshold=1.5
        if not liver1p5Done:
            setup['SUVdataset']='SUVanalysis_liver1p5'
            liver1p5Outputs=thresholdAnalysis(setup,r,pet,seg,liver1p5Threshold)
            uploadData(setup,r,liver1p5Outputs)

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

        #skip threshold of 4
        doThreshold4=False
        if doThreshold4:
            #threshold of 4
            setup['radiomics']['setting']['resegmentRange']=[4]
            setup['radiomics']['setting']['resegmentShape']=True
            outputs4=getValues(setup,r,pet,seg)
            setup['SUVdataset']='SUVanalysis_SUV4'
            uploadData(setup,r,outputs4)

        
        #outputs=getValues(setup,users,r,pet)
        #uploadData(setup,r,outputs)

        #cleanup
        os.remove(localPath)

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


def thresholdAnalysis(setup,r,pet,seg,thrs):
    #thresholds thrs are by participant and region
    outputs={}
    for s in thrs:
        #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]]
    return outputs


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

    rFile='radiomics.json'

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

    setup['featureExtractor']=radiomics.featureextractor.RadiomicsFeatureExtractor(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=statUtils.getRadiomicsComponentStats(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):

    rFile='radiomics.json'

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

    setup['featureExtractor']=radiomics.featureextractor.RadiomicsFeatureExtractor(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=statUtils.getRadiomicsComponentStats(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']
    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'],setup['project'],setup['SUVdataset'],[output])
    
def checkData(setup,r):
    qFilter=[]
    qFilter.append({'variable':'ParticipantId','value':r['ParticipantId'],'oper':'eq'})
    qFilter.append({'variable':'visitCode','value':r['visitCode'],'oper':'eq'})
    ds=setup['db'].selectRows(setup['project'],'study',setup['SUVdataset'],qFilter)
    n=len(ds['rows'])
    print('[{}:{}/{}] got {} rows.'.format(setup['SUVdataset'],r['ParticipantId'],r['visitCode'],n))
    return n>0



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