Adding python routines for NPZ coding and decoding

pythonScripts/generateNPZ.py

@@ -0,0 +1,177 @@
+import os
+import sys
+import numpy
+import random
+import SimpleITK
+import datetime
+
+def initDB(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('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
+        codes.append(code)
+        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
+
+

pythonScripts/parseNPZ.py

@@ -0,0 +1,30 @@
+import SimpleITK
+import numpy
+
+def convertToITK(dt,label):
+    #convert to SimpleITKimage to display it
+    imgArray=dt[label]
+    org=dt['{}origin'.format(label)]
+    orient=dt['{}orientation'.format(label)]
+    spacing=[numpy.linalg.norm(x) for x in orient]
+    print(spacing)
+    orient=[orient[i]/spacing[i] for i in range(3)]
+    direction=numpy.reshape(orient,9)
+    print(direction)
+    img = SimpleITK.GetImageFromArray(imgArray)
+    img.SetOrigin(org)
+    img.SetDirection(direction)
+    img.SetSpacing(spacing)
+    return img
+
+def generateNiftyFromNPZ(infile):
+    #generate nifty files for PET and CT stored in NPZ
+    outdir=os.path.dirname(infile)
+    baseName=os.path.basename(infile).replace('.npz','')
+    dt=numpy.load(infile)
+    ct=parseNPZ.convertToITK(dt,'CT')
+    pet=parseNPZ.convertToITK(dt,'PET')
+    SimpleITK.WriteImage(ct,os.path.join(outdir,f'{baseName}_CT.nii.gz'))
+    SimpleITK.WriteImage(pet,os.path.join(outdir,f'{baseName}_PET.nii.gz'))
+    del pet
+    del ct