Adding organ_percentile calculation, runs with segmentation.json

pythonScripts/organ_percentile.py

@@ -0,0 +1,153 @@
+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']
+    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,[])
+    
+    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)))
+    for r in rows:
+        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
+        remoteDir=fb.buildPathURL(project,[pars['imageDir'],\
+            r['patientCode'],r['visitCode']])
+        print('{}: {}'.format(petField,r[petField]))
+        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('Failed to get segmentation for {}/{}'.format(r[participantField],segVersion))
+
+        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
+        for level in [3,4,5]:
+            v=organ_percentile(niPET.get_fdata(),niSeg.get_fdata(),level,pv)
+            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'}
+                dsP=db.selectRows(project,schema,qQuery,[idFilter,sFilter,oFilter])
+                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
+                db.modifyRows(mode,project,schema,qQuery,[rowDSP])
+
+    
+
+    print('Done')
+
+
+if __name__ == '__main__':
+    main(sys.argv[1])
+
+
+