limfomiPET/pythonScripts/statUtils.py

#load required libraries
import sys
import os
import SimpleITK
import json

#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')
sys.path.append(os.path.join(nixSuite,'wrapper'))
import nixWrapper
nixWrapper.loadLibrary('labkeyInterface')
import labkeyInterface
import labkeyDatabaseBrowser
import labkeyFileBrowser

def connectDB(server):
    #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):
        print('{} done'.format(localPath))
    else:
        if not fb.entryExists(urlPath):
            print('No file {}'.format(urlPath))
            return "NONE"
        fb.readFileToFile(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])