#load required libraries import sys import os import SimpleITK import json def connectDB(server): #you should get nixSuite via git clone https://git0.fmf.uni-lj.si/studen/nixSuite.git #if you don't put it to $HOME/software/src/, you should update the path nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixSuite') if not os.path.isdir(nixSuite): nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixsuite') sys.path.append(os.path.join(nixSuite,'wrapper')) import nixWrapper nixWrapper.loadLibrary('labkeyInterface') import labkeyInterface import labkeyDatabaseBrowser import labkeyFileBrowser #check connectivity. This checks the configuration in $HOME/.labkey/network.json, #where paths to certificates are stored net=labkeyInterface.labkeyInterface() fconfig=os.path.join(os.path.expanduser('~'),'.labkey','{}.json'.format(server)) net.init(fconfig) #this reports the certificate used try: print('Using: {}'.format(net.connectionConfig['SSL']['user'])) except KeyError: pass #This gets a deafult CSRF code; It should report user name plus a long string of random hex numbers net.getCSRF() db=labkeyDatabaseBrowser.labkeyDB(net) fb=labkeyFileBrowser.labkeyFileBrowser(net) return db,fb def getUsers(db,project): ds=db.selectRows(project,'core','Users',[]) users={x['UserId']:x['DisplayName'] for x in ds['rows']} for u in users: print('{} {}'.format(u,users[u])) return users def getImage(setup, row, field, extraPath=None): fb=setup['fb'] pathList=[setup['imageDir'],row['patientCode'],row['visitCode']] if extraPath!=None: pathList.append(extraPath) pathList.append(row[field]) remotePath='/'.join(pathList) urlPath=fb.formatPathURL(setup['project'],remotePath) localPath=os.path.join(setup['localDir'],row[field]) if os.path.isfile(localPath): pass #print('{} done'.format(localPath)) else: if not fb.entryExists(urlPath): print('No file {}'.format(urlPath)) return "NONE" fb.readFileToFile(urlPath,localPath) print('Copying remote {} -> {} done'.format(urlPath,localPath)) return localPath def getSegmentations(setup,row): idField=setup['idField'] visitField=setup['visitField'] qFilter=[{'variable':x,'value':row[x],'oper':'eq'} for x in [idField,visitField]] dsSeg=setup['db'].selectRows(setup['project'],'study','Segmentations',qFilter) if len(dsSeg['rows'])<1: print('Failed to find segmentation for {}/{}'.format(row[idField],row[visitField])) return {0:"NONE"} return {r['User']:getImage(setup,r,'latestFile',extraPath='Segmentations') for r in dsSeg['rows']} def loadSetup(path): with open(path,'r') as f: setup=json.load(f) fC=[x for x in setup.keys() if x.find('Dir')>-1] for q in fC: setup[q]=setup[q].replace('_home_',os.path.expanduser('~')) return setup def getSegments(image): keys=image.GetMetaDataKeys() i=0 ids={} while True: id='Segment{}_ID'.format(i) value='Segment{}_LabelValue'.format(i) try: ids[image.GetMetaData(id)]=int(image.GetMetaData(value)) except RuntimeError: break i+=1 return ids def getStats(image): print(image.GetSize()) print(image.GetOrigin()) print(image.GetSpacing()) print(image.GetNumberOfComponentsPerPixel()) def getSegmentStats(pet,seg,label): o=getComponents(seg,label) cc=o['image'] shape_stats = SimpleITK.LabelShapeStatisticsImageFilter() #shape_stats.ComputeOrientedBoundingBoxOn() shape_stats.Execute(cc) intensity_stats = SimpleITK.LabelIntensityStatisticsImageFilter() intensity_stats.Execute(cc,pet) #output out=[(shape_stats.GetPhysicalSize(i), intensity_stats.GetMean(i), intensity_stats.GetStandardDeviation(i), intensity_stats.GetSkewness(i)) for i in shape_stats.GetLabels()] print(out) def getComponents(seg,label): cc=SimpleITK.Threshold(seg,lower=label,upper=label) ccFilter=SimpleITK.ConnectedComponentImageFilter() cc1=ccFilter.Execute(cc) return {'image':cc1, 'count':ccFilter.GetObjectCount()} def drop_array(v): return float(v) #if type(v) is numpy.ndarray: # return v[0] #return v def getSUVmax(vals): return drop_array(vals['original_firstorder_Maximum']) def getSUVmean(vals): return drop_array(vals['original_firstorder_Mean']) def getMTV(vals): try: return drop_array(vals['original_shape_MeshVolume']) except KeyError: return drop_array(vals['original_shape_VoxelVolume']) def getTLG(vals): V=vals['original_shape_VoxelVolume'] return V*getSUVmean(vals) def getCOM(vals): return vals['diagnostics_Mask-original_CenterOfMassIndex'] def getValue(vals,valueName): if valueName=='SUVmean': return getSUVmean(vals) if valueName=='SUVmax': return getSUVmax(vals) if valueName=='MTV': return getMTV(vals) if valueName=='TLG': return getTLG(vals) if valueName=='COM': return getCOM(vals) return 0 def getRadiomicsStats(setup,pet,seg,label): o=getComponents(seg,label) cc=o['image'] n=o['count'] output={} for i in range(n): output[i]=getRadiomicsComponentStats(setup,pet,seg,label) return output def getRadiomicsComponentStats(setup,pet,cc,label): vals=setup['featureExtractor'].execute(pet,cc,label=label) output={x:getValue(vals,x) for x in setup['values']} #for v in vals: # print('{}: {}'.format(v,vals[v])) for c in output: print('{}: {}'.format(c,output[c])) return output def findMatchingComponent(o,a,b,label): statFilter=SimpleITK.StatisticsImageFilter() overlap_measures_filter = SimpleITK.LabelOverlapMeasuresImageFilter() print('{}: [{}]:{} [{}]:{}'.format(label,a,o[a]['count'],b,o[b]['count'])) comps={v:{x+1:o[v]['image']==x+1 for x in range(o[v]['count'])} for v in [a,b]} stat={} for comp in comps: stat[comp]={} for x in comps[comp]: cc=comps[comp][x] statFilter.Execute(cc) stat[comp][x]={'sum':statFilter.GetSum()} for c in comps[a]: cc=comps[a][c] print('{}:{}'.format(c,stat[a][c]['sum'])) for d in comps[b]: cc1=comps[b][d] overlap_measures_filter.Execute(cc, cc1) print('\t {}:{} {}'.format(d,stat[b][d]['sum'],overlap_measures_filter.GetDiceCoefficient())) def evaluateByLesions(pet,seg,ids): for id in ids: print('{}: {}'.format(id,ids[id])) o={x:getComponents(seg[x],ids[id]) for x in seg} a=segKeys[0] for x in segKeys[1:]: findMatchingComponent(o,a,x,ids[id]) def updateDatasetRows(db,project,dataset,rows,filterVars=['ParticipantId','SequenceNum']): for r in rows: r['SequenceNum']+=0.01*r['segment'] qFilter=[{'variable':x,'value':''.format(r[x]),'oper':'eq'} for x in filterVars] ds=db.selectRows(project,'study',dataset,qFilter) if len(ds['rows'])>0: row=ds['rows'][0] row.update({x:r[x] for x in r if x not in filterVars}) db.modifyRows('update',project,'study',dataset,[row]) else: db.modifyRows('insert',project,'study',dataset,[r])