irAEMM/pythonScripts/organ_percentile.py

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'])

    participantField=pars['participantField']

    tempBase=pars['tempBase']
    if not os.path.isdir(tempBase):
        os.makedirs(tempBase)
    
    qFilter=[]
    try:
        vList=';'.join(pars['participants'])
        qFilter.append({'variable':participantField,'value':vList,'oper':'in'})
    except KeyError:
        pass
    try:
        vList=';'.join(pars['visits'])
        qFilter.append({'variable':'visitCode','value':vList,'oper':'in'})
    except KeyError:
        pass



    #all images from database
    ds=db.selectRows(project,schema,dataset,qFilter,view)
    
    petField=pars['images']['PET']['queryField']
    
    rows=ds['rows']
    print('Selected {} rows'.format(len(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])
        print('Deleting {} rows'.format(len(dsP['rows'])))
        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)
        print('Inserting {} rows'.format(len(dspRows)))
        db.modifyRows('insert',project,schema,qQuery,dspRows)

    

    print('Done')


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