import numpy import sys import os import json import nibabel def organ_percentile(img, mask, level, p): #ORGAN_PERCENTILE - computes suv_x the pth percentile of the distribution of # image intensity values in img from within some ROI mask # # Inputs: # img - image. ndarray # mask - ROI mask. Must be same dimension as img. Must be binary ndarray # p - percentile. Between 0-100 # # Outputs: # suv_x - percentile of distribution. Defined by: # # suv_x # p/100 = ∫ H(x) dx # 0 # # where H(x) is the normalized distribution of image values within mask #img = np.array(img) #mask = np.array(mask) h = img[mask == level] suv_x = numpy.percentile(h, p) return suv_x def main(parameterFile): fhome=os.path.expanduser('~') with open(os.path.join(fhome,".labkey","setup.json")) as f: setup=json.load(f) sys.path.insert(0,setup["paths"]["nixWrapper"]) import nixWrapper nixWrapper.loadLibrary("labkeyInterface") import labkeyInterface import labkeyDatabaseBrowser import labkeyFileBrowser fconfig=os.path.join(fhome,'.labkey','network.json') net=labkeyInterface.labkeyInterface() net.init(fconfig) db=labkeyDatabaseBrowser.labkeyDB(net) fb=labkeyFileBrowser.labkeyFileBrowser(net) with open(parameterFile) as f: pars=json.load(f) #segmentation layout project=pars['project'] dataset=pars['targetQuery'] schema=pars['targetSchema'] view=pars['viewName'] segSchema=pars['segmentationSchema'] segQuery=pars['segmentationQuery'] qQuery=pars['percentileQuery'] segVersion=pars['version'] segVersionI=int(pars['versionNumber']) tempBase=pars['tempBase'] if not os.path.isdir(tempBase): os.makedirs(tempBase) participantField=pars['participantField'] #all images from database ds=db.selectRows(project,schema,dataset,[],view) petField=pars['images']['PET']['queryField'] rows=ds['rows'] #rows=[ds['rows'][0]] pv=numpy.concatenate((numpy.linspace(10,50,5), numpy.linspace(55,80,6), numpy.linspace(82,90,5), numpy.linspace(91,100,10))) #debug, select #rows1=[r for r in rows \ # if r['patientCode']=='NIX-LJU-D2002-IRAE-A010' \ # and r['visitCode']=='VISIT_6'] #rows=rows1 #rows=[r for r in rows if r[participantField]=='956/18'] for r in rows: #is segmentation available localPET=os.path.join(tempBase,'PET.nii.gz') localSeg=os.path.join(tempBase,'Seg.nii.gz') for f in [localPET,localSeg]: if os.path.isfile(f): os.remove(f) #build image path try: remoteDir=fb.buildPathURL(project,[pars['imageDir'],\ r['patientCode'],r['visitCode']]) except TypeError: print('[{}/{}]: Skipping - Failed to build remote dir'.\ format(r[participantField],r['SequenceNum'])) continue print('{}: {}'.format(petField,r[petField])) if r[petField]==None: print('[{}/{}]: Skipping - Missing PET'.\ format(r[participantField],r['SequenceNum'])) continue remotePET=remoteDir+'/'+r[petField] print('{}:{}'.format(remotePET,fb.entryExists(remotePET))) vFilter={'variable':'version','value':segVersion,'oper':'eq'} idFilter={'variable':'patientCode','value':r['patientCode'], 'oper':'eq'} visitFilter={'variable':'visitCode','value':r['visitCode'], 'oper':'eq'} dsSeg=db.selectRows(project,segSchema,segQuery,[idFilter,visitFilter,vFilter]) if len(dsSeg['rows'])!=1: print('[{}/{}] Skipping - missing segmentation, version [{}]'.\ format(r[participantField],r['SequenceNum'],segVersion)) continue remoteSeg=remoteDir+'/'+dsSeg['rows'][0]['segmentation'] print('{}:{}'.format(remoteSeg,fb.entryExists(remoteSeg))) fb.readFileToFile(remotePET,localPET) fb.readFileToFile(remoteSeg,localSeg) niPET=nibabel.load(localPET) niSeg=nibabel.load(localSeg) #3 lungs #4 thyroid #5 bowel dsP=db.selectRows(project,schema,qQuery,[idFilter,visitFilter,vFilter]) db.modifyRows('delete',project,schema,qQuery,dsP['rows']) dspRows=[] for level in [1,2,3,4,5,6,7,8]: try: v=organ_percentile(niPET.get_fdata(),niSeg.get_fdata(),level,pv) except IndexError: print('Error for {}/{}: {}'.format(r['patientCode'],r['visitCode'],level)) continue for (x,y) in zip(pv,v): #get existing entry seqNum=r['SequenceNum']+0.0001*x+0.01*segVersionI #print('[{:.8f}] {}/{}: {}/{}'.format(seqNum,r['patientCode'],r['visitCode'],x,y)) # sFilter={'variable':'SequenceNum','value':'{}'.format(seqNum),'oper':'eq'} # oFilter={'variable':'organ','value':'{}'.format(level),'oper':'eq'} # mode='update' # if len(dsP['rows'])==0: # mode='insert' rowDSP={x:r[x] for x in [participantField,'patientCode','visitCode']} rowDSP['SequenceNum']=seqNum rowDSP['segmentationVersion']=segVersion # else: # rowDSP=dsP['rows'][0] rowDSP['percentile']=x rowDSP['value']=y rowDSP['organ']=level dspRows.append(rowDSP) db.modifyRows('insert',project,schema,qQuery,dspRows) print('Done') if __name__ == '__main__': main(sys.argv[1])