|
@@ -0,0 +1,229 @@
|
|
|
+import nibabel
|
|
|
+import os
|
|
|
+import json
|
|
|
+import sys
|
|
|
+import numpy
|
|
|
+import matplotlib.pyplot
|
|
|
+#import chardet
|
|
|
+
|
|
|
+def buildPath(server,project,imageDir,patientCode,visitCode,imageType):
|
|
|
+ path='/'.join([server,'labkey/_webdav',project,'@files',imageDir,patientCode,visitCode])
|
|
|
+ tail='_notCropped_2mmVoxel'
|
|
|
+ if imageType=='Segm':
|
|
|
+ tail='_v5'
|
|
|
+ path+='/'+patientCode+'-'+visitCode+'_'+imageType+tail+'.nii.gz'
|
|
|
+ return path
|
|
|
+
|
|
|
+
|
|
|
+def getCOWAxis(seg,val,axis):
|
|
|
+#returns center of weight for segmentation image where val is selected
|
|
|
+ if axis==0:
|
|
|
+ #2,1 or 1,1
|
|
|
+ i1=2
|
|
|
+ i2=1
|
|
|
+ if axis==1:
|
|
|
+ #0,1 or 2,0
|
|
|
+ i1=2
|
|
|
+ i2=0
|
|
|
+ if axis==2:
|
|
|
+ #0,0 or 1,0
|
|
|
+ i1=1
|
|
|
+ i2=0
|
|
|
+
|
|
|
+ s=numpy.sum(numpy.sum(seg==val,i1),i2)
|
|
|
+ x=numpy.arange(len(s))
|
|
|
+ s0=0
|
|
|
+ for i in x:
|
|
|
+ if s[i]==0:
|
|
|
+ continue
|
|
|
+ s0=i
|
|
|
+ break
|
|
|
+ s1=len(s)
|
|
|
+ for i in numpy.arange(s0,len(s)):
|
|
|
+ if s[i]>0:
|
|
|
+ continue
|
|
|
+ s1=i
|
|
|
+ break
|
|
|
+ return [s0,numpy.average(x,weights=s),s1]
|
|
|
+
|
|
|
+def getGeometry(seg,val):
|
|
|
+#return center of weight of segmentation seg for segment val as a 3D vector
|
|
|
+ return [getCOWAxis(seg,val,x) for x in [0,1,2]]
|
|
|
+
|
|
|
+def getCOW(geom):
|
|
|
+ return [x[1] for x in geom]
|
|
|
+
|
|
|
+def getRange(geom):
|
|
|
+ return [[x[0],x[2]] for x in geom]
|
|
|
+
|
|
|
+
|
|
|
+def plot(imgs,t,val,tempBase):
|
|
|
+
|
|
|
+
|
|
|
+ segColors=[[0,0,0],[0.1,0.1,0.1],[0,0.2,0],[1,0,0],[0,0,1],[1,0,1]]
|
|
|
+
|
|
|
+#3-lungs 4-thyroid 5-bowel
|
|
|
+ delta=20
|
|
|
+ if val==4:
|
|
|
+ delta=40
|
|
|
+ window=350
|
|
|
+ level=40
|
|
|
+ geo=getGeometry(imgs['Segm'],val)
|
|
|
+ cowF=getCOW(geo)
|
|
|
+ rng=getRange(geo)
|
|
|
+ #print(rng)
|
|
|
+ cowI=[int(x) for x in cowF]
|
|
|
+ segment=imgs['Segm']==val
|
|
|
+
|
|
|
+ i0=rng[0][0]-delta
|
|
|
+ if i0<0:
|
|
|
+ i0=0
|
|
|
+ i1=rng[0][1]+delta
|
|
|
+ if i1>imgs['CT'].shape[0]:
|
|
|
+ i1=imgs['CT'].shape[0]
|
|
|
+ k0=rng[2][0]-delta
|
|
|
+ if k0<0:
|
|
|
+ k0=0
|
|
|
+ k1=rng[2][1]+delta
|
|
|
+ if k1>imgs['CT'].shape[2]:
|
|
|
+ k1=imgs['CT'].shape[2]
|
|
|
+
|
|
|
+ if t=='CT':
|
|
|
+ v0=level-0.5*window
|
|
|
+ v1=v0+window
|
|
|
+ matplotlib.pyplot.imshow(imgs['CT'][i0:i1,cowI[1],k0:k1].transpose(),cmap='gray',vmin=v0,vmax=v1)
|
|
|
+ if t=='PET':
|
|
|
+ matplotlib.pyplot.imshow(imgs['PET'][i0:i1,cowI[1],k0:k1].transpose(),cmap='inferno')
|
|
|
+ #blueish
|
|
|
+
|
|
|
+ if t=='CT':
|
|
|
+ rgb=segColors[val]
|
|
|
+ if t=='PET':
|
|
|
+ rgb=[1,1,1]
|
|
|
+
|
|
|
+ colors = [rgb+[c] for c in numpy.linspace(0,1,100)]
|
|
|
+
|
|
|
+ cmap = matplotlib.colors.LinearSegmentedColormap.from_list('mycmap', colors, N=5)
|
|
|
+
|
|
|
+ matplotlib.pyplot.imshow(segment[i0:i1,cowI[1],k0:k1].transpose(), cmap=cmap, alpha=0.2)
|
|
|
+ matplotlib.pyplot.gca().invert_yaxis()
|
|
|
+ outfile=os.path.join(tempBase,'slice{}_{}.png'.format(t,val))
|
|
|
+ matplotlib.pyplot.savefig(outfile)
|
|
|
+ return outfile
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+def main(parameterFile):
|
|
|
+#mask for segmentations
|
|
|
+
|
|
|
+ setupJSON=os.path.join(os.path.expanduser('~'),'.labkey','setup.json')
|
|
|
+ with open(setupJSON) 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(os.path.expanduser('~'),'.labkey','network.json')
|
|
|
+
|
|
|
+ net=labkeyInterface.labkeyInterface()
|
|
|
+ net.init(fconfig)
|
|
|
+ db=labkeyDatabaseBrowser.labkeyDB(net)
|
|
|
+ fb=labkeyFileBrowser.labkeyFileBrowser(net)
|
|
|
+
|
|
|
+ tempBase=os.path.join(os.path.expanduser('~'),'temp')
|
|
|
+
|
|
|
+ with open(parameterFile) as f:
|
|
|
+ pars=json.load(f)
|
|
|
+
|
|
|
+ project=pars['project']
|
|
|
+ dataset=pars['targetQuery']
|
|
|
+ schema=pars['targetSchema']
|
|
|
+
|
|
|
+ reportSchema=pars['reportSchema']
|
|
|
+ reportQuery=pars['reportQuery']
|
|
|
+ participantField=pars['participantField']
|
|
|
+
|
|
|
+ #all images from database
|
|
|
+ ds=db.selectRows(project,schema,dataset,[])
|
|
|
+
|
|
|
+ #input
|
|
|
+ imageResampledField={"CT":"ctResampled","PET":"petResampled","patientmask":"ROImask"}
|
|
|
+
|
|
|
+ rows=ds['rows']
|
|
|
+ rows=[ds['rows'][0]]
|
|
|
+
|
|
|
+ for r in rows:
|
|
|
+ print(r)
|
|
|
+ iTypes=['CT','PET','Segm']
|
|
|
+ needToCalculate=False
|
|
|
+ for t in ['CT','PET']:
|
|
|
+ idFilter={'variable':participantField,'value':r[participantField],'oper':'eq'}
|
|
|
+ visitFilter={'variable':'visitCode','value':r['visitCode'],'oper':'eq'}
|
|
|
+ verFilter={'variable':'version','value':pars['version'],'oper':'eq'}
|
|
|
+ typeFilter={'variable':'type','value':t,'oper':'eq'}
|
|
|
+
|
|
|
+ ds2=db.selectRows(project,reportSchema,reportQuery,[idFilter,visitFilter,verFilter,typeFilter])
|
|
|
+
|
|
|
+ if len(ds2['rows'])==0:
|
|
|
+ #skip if row is present
|
|
|
+ #there are in fact multiple rows for multiple organs...
|
|
|
+ needToCalculate=True
|
|
|
+ break
|
|
|
+
|
|
|
+ if not needToCalculate:
|
|
|
+ continue
|
|
|
+ imgs={}
|
|
|
+ for t in iTypes:
|
|
|
+
|
|
|
+ try:
|
|
|
+ imagePath=r['_labkeyurl_'+imageResampledField[t]]
|
|
|
+ except KeyError:
|
|
|
+ ds1=db.selectRows(project,pars['segmentationSchema'],pars['segmentationQuery'],\
|
|
|
+ [idFilter,visitFilter,verFilter])
|
|
|
+ imagePath=ds1['rows'][0]['_labkeyurl_segmentation']
|
|
|
+
|
|
|
+ localPath=os.path.join(tempBase,'image'+t+'.nii.gz')
|
|
|
+ if os.path.isfile(localPath):
|
|
|
+ os.remove(localPath)
|
|
|
+ fb.readFileToFile(imagePath,localPath)
|
|
|
+ img=nibabel.load(localPath)
|
|
|
+ imgs[t]=img.get_fdata()
|
|
|
+ print('Loading completed')
|
|
|
+ for t in ['CT','PET']:
|
|
|
+ for val in [3,4,5]:
|
|
|
+ outfile=plot(imgs,t,val,tempBase)
|
|
|
+ remoteDir=fb.buildPathURL(project,[pars['imageDir'],r['patientCode'],r['visitCode']])
|
|
|
+ imageFile=r['patientCode']+'-'+r['visitCode']+'_'+t+'_{}'.format(val)+'_'+pars['version']+'.png'
|
|
|
+ remoteFile='/'.join([remoteDir,imageFile])
|
|
|
+ fb.writeFileToFile(outfile,remoteFile)
|
|
|
+ print('Uploaded {}'.format(remoteFile))
|
|
|
+ os.remove(outfile)
|
|
|
+ organFilter={'variable':'organ','value':'{}'.format(val),'oper':'eq'}
|
|
|
+ typeFilter['value']=t
|
|
|
+ ds3=db.selectRows(project,reportSchema,reportQuery,\
|
|
|
+ [idFilter,visitFilter,verFilter,organFilter,typeFilter])
|
|
|
+ if len(ds3['rows'])>0:
|
|
|
+ mode='update'
|
|
|
+ frow=ds3['rows'][0]
|
|
|
+ else:
|
|
|
+ mode='insert'
|
|
|
+ frow={}
|
|
|
+ for f in [participantField,'patientCode','visitCode']:
|
|
|
+ frow[f]=r[f]
|
|
|
+ frow['organ']='{}'.format(val)
|
|
|
+ frow['type']=t
|
|
|
+ frow['version']=pars['version']
|
|
|
+ frow['file']=imageFile
|
|
|
+ db.modifyRows(mode,project,reportSchema,reportQuery,[frow])
|
|
|
+ print('Images uploaded')
|
|
|
+
|
|
|
+
|
|
|
+if __name__ == '__main__':
|
|
|
+ main(sys.argv[1])
|
|
|
+
|