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