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