irAEMM/pythonScripts/generateNPZ.py

import os
import sys
import numpy
import random
import SimpleITK
import datetime

def initDB(site,returnFB=False):
    home=os.path.expanduser('~')
    print(home)
    nixSuite=os.path.join(home,'software','src','nixSuite')
    sys.path.append(os.path.join(nixSuite,'wrapper'))
    import nixWrapper
    nixWrapper.loadLibrary('labkeyInterface')
    import labkeyInterface
    import labkeyDatabaseBrowser
    import labkeyFileBrowser
    net=labkeyInterface.labkeyInterface()
    fconfig=os.path.join(home,'.labkey','{}.json'.format(site))
    net.init(fconfig)
    net.getCSRF()
    if not returnFB:
        return labkeyDatabaseBrowser.labkeyDB(net)
    return labkeyDatabaseBrowser.labkeyDB(net),\
        labkeyFileBrowser.labkeyFileBrowser(net)

def initOrthanc(site,returnFB=False):

    home=os.path.expanduser('~')
    nixSuite=os.path.join(home,'software','src','nixSuite')
    sys.path.append(os.path.join(nixSuite,'wrapper'))
    import nixWrapper
    nixWrapper.loadLibrary('orthancInterface')
    import orthancInterface
    import orthancDatabaseBrowser
    import orthancFileBrowser
    net=orthancInterface.orthancInterface()
    fconfig=os.path.join(home,'.labkey','{}.json'.format(site))
    net.init(fconfig)
    if not returnFB:
        return orthancDatabaseBrowser.orthancDB(net)
    return orthancDatabaseBrowser.orthancDB(net),\
        orthancFileBrowser.orthancFileBrowser(net)

def parsePos(x):
    return [float(y) for y in x.split('\\')]

def parseOrientation(x):
    z=[float(y) for y in x.split('\\')]
    O=[z[0:3],z[3:6]]
    O.append(list(numpy.cross(numpy.array(O[0]),numpy.array(O[1]))))
    return O

def getGeometry(code,odb,seriesId):
    #imagepositionpatient 0020-0032
    #imageorientationpatient 0020-0037
    #pixelSpacing 0028-0030
    #instance number 0020-0013
    #sliceThickness 0018-0050
    
    sd=odb.getData('series',seriesId)
    firstId=sd['Instances'][0] 
    
    ior=parseOrientation(odb.getDicomField(firstId,'0020-0037'))
    #spacing
    ips=parsePos(odb.getDicomField(firstId,'0028-0030'))
    ist=parsePos(odb.getDicomField(firstId,'0018-0050')) 
    #don't use slice thickness for pixel spacing
    #ips.append(ist[0])
    
    ipos=[parsePos(odb.getDicomField(y,'0020-0032')) for y in sd['Instances']] 
    #guess spacing from ipos
    m=[numpy.dot(numpy.array(ior[2]),numpy.array(x)) for x in ipos] 
    #sort array
    m1=numpy.sort(m)
    dm=m1[1:]-m1[:-1]
    zSpc=numpy.average(dm)
    zDev=numpy.std(dm)
    print(f'sliceThickness={zSpc} +- {zDev}')
    ips.append(zSpc)
    
    #direction, ie. scaled orientation
    #copy orientation to direction
    xdir=[y for y in ior]
    for i in range(3):
        xdir[i]=[ips[i]*w for w in xdir[i]]
    
    #find minimum from original array to idenfify lower slice
    im=numpy.argmin(numpy.array(m))
    print(im)
    org=ipos[im]
    print(org)
    print(xdir)
    return {f'{code}origin':org,f'{code}orientation':xdir}
    
        
        
def getFile(db,odb,project,outdir,n1=0,n2=1):
    schema='study'
    query='Imaging1'
    #coreDir=os.path.join(os.path.expanduser('~'),'temp')
    ds=db.selectRows(project,schema,query,[])
    rows=ds['rows']
    if n2>0: 
        rows=rows[n1:n2]
    imgs={'CT':'CT_orthancId','PET':'PETWB_orthancId'}
    files=[]
    for r in rows:
        code=r['pseudoCode']
        doUpdate=False
        if not code:
            code='{:16x}'.format(random.randrange(16**16))
            r['pseudoCode']=code
            doUpdate=True
        fname=os.path.join(outdir,f'{code}.npz')
        arr={}
        #debug  don't load array
        iSUV={x:getSUVfactor(odb,r[imgs[x]]) for x in imgs}
        arr.update({x:iSUV[x]*numpy.load(odb.getNumpy('series',r[imgs[x]])) for x in imgs})
        #debug don't calculate geometry
        igeo={x:getGeometry(x,odb,r[imgs[x]]) for x in imgs}
        #debug skip saving
        #continue
        _={arr.update(igeo[x]) for x in igeo}
        #print(arr)
        
        
        numpy.savez_compressed(fname,**arr)
        files.append(fname)
        if doUpdate:
            db.modifyRows('update',project,schema,query,[r])
    return files
        
def mergeDateAndTime(date,time):
    time=time.replace(' ','')
    dateTimeStr=' '.join([date,time])
    #print(f'[{date}] [{time}] [{dateTimeStr}]')
    dt=datetime.datetime.strptime(dateTimeStr,'%Y%m%d %H%M%S.%f')
    #print(dt)
    return dt
        
def getSUVfactor(odb,seriesId):
    sd=odb.getData('series',seriesId)
    orthancId=sd['Instances'][0] 
    
    MODALITY='0008-0060'
    PATIENT_WEIGHT='0010-1030'
    RADIOPHARMACEUTICAL_STARTTIME='0054-0016/0/0018-1072'
    RADIONUCLIDE_TOTALDOSE='0054-0016/0/0018-1074'
    RADIONUCLIDE_HALFLIFE='0054-0016/0/0018-1075'
    RADIOPHARMACEUTICAL_STARTDATETIME='0054-0016/0/0018-1078'
    SERIES_TIME='0008-0031'
    SERIES_DATE='0008-0021'
    modality=odb.getDicomField(orthancId,MODALITY)
    if modality!='PT':
        print('Returning 1')
        return 1
    weight = float(odb.getDicomField(orthancId,PATIENT_WEIGHT))
    #print(weight)
    halfLife = float(odb.getDicomField(orthancId,RADIONUCLIDE_HALFLIFE))#seconds
    #print(halfLife)
    injDose = float(odb.getDicomField(orthancId,RADIONUCLIDE_TOTALDOSE))#in Bq
    #print(injDose)
    #two options for time of activity measurement (full date time or time only, then match it with studyDate)
    #injDateTime = odb.getDicomField(orthancId,RADIOPHARMACEUTICAL_STARTDATETIME)
    #if troubles:
    injDateTime = mergeDateAndTime(odb.getDicomField(orthancId,SERIES_DATE),
            odb.getDicomField(orthancId,RADIOPHARMACEUTICAL_STARTTIME))
    scanDateTime=mergeDateAndTime(odb.getDicomField(orthancId,SERIES_DATE),
                                  odb.getDicomField(orthancId,SERIES_TIME))
    t=(scanDateTime-injDateTime).total_seconds()
    scanDose = injDose*numpy.exp(-(numpy.log(2)*t)/halfLife);
    SUV_factor = 1/((scanDose/weight)*0.001);
    print(SUV_factor)
    return SUV_factor