Merge branch 'master' of wiscigt.powertheword.com:oil/iraemm

pythonScripts/generateFigures.py

@@ -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])
+

pythonScripts/modifyDataset.py

@@ -0,0 +1,61 @@
+#a script to modify a patient from list. Current implementation deletes 
+#all data for patients identified by study regulatory number from study
+
+#basic python
+import os
+import subprocess
+import re
+import datetime
+import sys
+import json
+
+
+def main(parameterFile):
+    fhome=os.path.expanduser('~')
+
+    with open(os.path.join(fhome,".labkey","setup.json")) as f:
+        setup=json.load(f)
+
+    sys.path.insert(0,setup["paths"]["nixWrapper"])
+    import nixWrapper
+
+    nixWrapper.loadLibrary("labkeyInterface")
+
+    import labkeyInterface
+    import labkeyDatabaseBrowser
+    import labkeyFileBrowser
+
+
+    net=labkeyInterface.labkeyInterface()
+    net.init(os.path.join(fhome,'.labkey','network.json'))
+
+    db=labkeyDatabaseBrowser.labkeyDB(net)
+
+    #by default uses .labkey/Remote.json configuration
+    with open(parameterFile) as f:
+        pars=json.load(f)
+
+
+    project=pars['project']
+    schema=pars['schema']
+    query=pars['query']
+    
+
+    #study section ################
+
+
+
+#select patients enroled under regulatory number
+    
+    ds=db.selectRows(project,schema,query,[])
+
+    for r in ds['rows']:
+        r['View']='[VIEW]'
+        db.modifyRows('update',project,schema,query,[r])
+
+    print("Done")
+
+if __name__ == '__main__':
+    main(sys.argv[1])
+
+

pythonScripts/organ_percentile.py

@@ -0,0 +1,153 @@
+import numpy
+import sys
+import os
+import json
+import nibabel
+
+def organ_percentile(img, mask, level, p):
+#ORGAN_PERCENTILE - computes suv_x the pth percentile of the distribution of
+#   image intensity values in img from within some ROI mask
+#
+#   Inputs:
+#       img - image. ndarray
+#       mask - ROI mask. Must be same dimension as img. Must be binary ndarray
+#       p - percentile. Between 0-100
+#
+#   Outputs:
+#       suv_x - percentile of distribution. Defined by:
+#
+#                 suv_x
+#           p/100 = ∫ H(x) dx
+#                   0
+#
+#       where H(x) is the normalized distribution of image values within mask
+    #img = np.array(img)
+    #mask = np.array(mask)
+
+    h = img[mask == level]
+    suv_x = numpy.percentile(h, p)
+
+    return suv_x
+
+
+def main(parameterFile):
+    fhome=os.path.expanduser('~')
+    with open(os.path.join(fhome,".labkey","setup.json")) 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(fhome,'.labkey','network.json')
+
+
+    net=labkeyInterface.labkeyInterface()
+    net.init(fconfig)
+    db=labkeyDatabaseBrowser.labkeyDB(net)
+    fb=labkeyFileBrowser.labkeyFileBrowser(net)
+
+    with open(parameterFile) as f:
+        pars=json.load(f)
+
+    #segmentation layout
+    project=pars['project']
+    dataset=pars['targetQuery']
+    schema=pars['targetSchema']
+    segSchema=pars['segmentationSchema']
+    segQuery=pars['segmentationQuery']
+    qQuery=pars['percentileQuery']
+    segVersion=pars['version']
+    segVersionI=int(pars['versionNumber'])
+
+    tempBase=pars['tempBase']
+    if not os.path.isdir(tempBase):
+        os.makedirs(tempBase)
+    
+    participantField=pars['participantField']
+
+    #all images from database
+    ds=db.selectRows(project,schema,dataset,[])
+    
+    petField=pars['images']['PET']['queryField']
+    
+    rows=ds['rows']
+    #rows=[ds['rows'][0]]
+
+    pv=numpy.concatenate((numpy.linspace(10,50,5),
+            numpy.linspace(55,80,6),  
+            numpy.linspace(82,90,5),
+            numpy.linspace(91,100,10)))
+    for r in rows:
+        localPET=os.path.join(tempBase,'PET.nii.gz')
+        localSeg=os.path.join(tempBase,'Seg.nii.gz')
+        for f in [localPET,localSeg]:
+            if os.path.isfile(f):
+                os.remove(f)
+
+        #build image path
+        remoteDir=fb.buildPathURL(project,[pars['imageDir'],\
+            r['patientCode'],r['visitCode']])
+        print('{}: {}'.format(petField,r[petField]))
+        remotePET=remoteDir+'/'+r[petField]
+        print('{}:{}'.format(remotePET,fb.entryExists(remotePET)))
+
+        vFilter={'variable':'version','value':segVersion,'oper':'eq'}
+        idFilter={'variable':'patientCode','value':r['patientCode'], 'oper':'eq'}
+        visitFilter={'variable':'visitCode','value':r['visitCode'], 'oper':'eq'}
+
+        dsSeg=db.selectRows(project,segSchema,segQuery,[idFilter,visitFilter,vFilter])
+        if len(dsSeg['rows'])!=1:
+            print('Failed to get segmentation for {}/{}'.format(r[participantField],segVersion))
+
+        remoteSeg=remoteDir+'/'+dsSeg['rows'][0]['segmentation']
+
+        print('{}:{}'.format(remoteSeg,fb.entryExists(remoteSeg)))
+
+        fb.readFileToFile(remotePET,localPET)
+        fb.readFileToFile(remoteSeg,localSeg)
+        
+        niPET=nibabel.load(localPET)
+        niSeg=nibabel.load(localSeg)
+        #3 lungs
+        #4 thyroid
+        #5 bowel
+        for level in [3,4,5]:
+            v=organ_percentile(niPET.get_fdata(),niSeg.get_fdata(),level,pv)
+            for (x,y) in zip(pv,v):
+                #get existing entry
+                seqNum=r['SequenceNum']+0.0001*x+0.01*segVersionI
+                print('[{:.8f}] {}/{}: {}/{}'.format(seqNum,r['patientCode'],r['visitCode'],x,y))
+
+                sFilter={'variable':'SequenceNum','value':'{}'.format(seqNum),'oper':'eq'}
+                oFilter={'variable':'organ','value':'{}'.format(level),'oper':'eq'}
+                dsP=db.selectRows(project,schema,qQuery,[idFilter,sFilter,oFilter])
+                mode='update'
+                if len(dsP['rows'])==0:
+                    mode='insert'
+                    rowDSP={x:r[x] for x in [participantField,'patientCode','visitCode']}
+                    rowDSP['SequenceNum']=seqNum
+                    rowDSP['segmentationVersion']=segVersion
+                else:
+                    rowDSP=dsP['rows'][0]
+                rowDSP['percentile']=x
+                rowDSP['value']=y
+                rowDSP['organ']=level
+                db.modifyRows(mode,project,schema,qQuery,[rowDSP])
+
+    
+
+    print('Done')
+
+
+if __name__ == '__main__':
+    main(sys.argv[1])
+
+
+

pythonScripts/preprocess.py

@@ -18,7 +18,7 @@ def getStudyLabel(row,participantField='PatientId'):
     return getPatientLabel(row,participantField)+'-'+getVisitLabel(row)
 
 def runPreprocess_DM(matlab,generalCodes,niftiTools,studyDir):
-
+    print("Running matlab")
     #run after all directories have been assembled
     script="addpath('"+generalCodes+"');"
     script+="addpath('"+niftiTools+"');"
@@ -46,7 +46,7 @@ def getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
     if seriesId=="0":
         return False
 
-    print("{}: {}".format(im,seriesId))
+    print("getDicom: {}: {}".format(im,seriesId))
     fname=os.path.join(zipDir,\
             getStudyLabel(row,participantField)+'_'+im+".zip");
 
@@ -54,7 +54,7 @@ def getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
     if os.path.isfile(fname):
         print("Data already loaded. Skipping")
     else:
-        print("Loading data from orthanc")
+        print("getDicom: Loading data from orthanc")
         ofb.getZip('series',seriesId,fname)
 
     #unzip the zipped dicom series
@@ -111,6 +111,7 @@ def main(parameterFile):
     matlab=setup["paths"]["matlab"]
     generalCodes=setup["paths"]["generalCodes"]
     niftiTools=setup["paths"]["niftiTools"]
+    gzip=setup['paths']['gzip']
 
     net=labkeyInterface.labkeyInterface()
     net.init(fconfig)
@@ -148,7 +149,7 @@ def main(parameterFile):
 
     i=0
     for row in ds["rows"]:
-
+        print("Starting row id:{} seq:{}".format(row[participantField],row['SequenceNum']))
         #interesting files are processedDir/studyName_CT_notCropped_2mmVoxel.nii
         #asn processedDir/studyName_PET_notCropped_2mmVoxel.nii
         volumeFileNames={im:\
@@ -181,6 +182,8 @@ def main(parameterFile):
     
         #setup the directory structure for preprocess_DM
         studyDir=os.path.join(tempBase,getStudyLabel(row,participantField))
+        print("Making local directories in {}".format(studyDir))
+
         if not os.path.isdir(studyDir):
             os.mkdir(studyDir)
 
@@ -224,7 +227,7 @@ def main(parameterFile):
 
             for f in volumeFiles.values():
                 print("Running gzip {}".format(f))
-                outText=subprocess.check_output(["/bin/gzip",f])
+                outText=subprocess.check_output([gzip,f])
                 print(outText.decode('utf-8'))
 
         #upload local files to remote

pythonScripts/runSegmentationnnUNet.py

@@ -6,7 +6,6 @@ import nibabel
 import shutil
 import sys
 import pathlib
-import SimpleITK
 import numpy
 
 #nothing gets done if you do import
@@ -198,13 +197,14 @@ def doSegmentation(parameterFile):
 
  
     i=0
+    #for debugging
     rows=[ds['rows'][0]]
+    #production mode
     rows=ds['rows']
     for row in rows:
        
 
-        #check if file is already there
-        #dummy tf to get the suffix
+        #build file name 
         sfx=pars['images']['segmentation']['suffix']
         outpath=fb.buildPathURL(pars['project'],[pars['imageDir'],row['patientCode'],row['visitCode']])
         outName=addVersion(\
@@ -248,7 +248,7 @@ def doSegmentation(parameterFile):
         outRow['Version']= pars['version']
         outRow['Segmentation']= outName
         print(db.modifyRows(mode,pars['project'],pars['segmentationSchema'],pars['segmentationQuery'],[outRow]))
-        #pull results back to LabKey
+        #push results back to LabKey
     print("Done")
 
 

templates/segmentationIRAEMM_ONKO.json

@@ -0,0 +1,87 @@
+{
+ "setVariables":["__tempBase__","__segBase__","__roiFile__","__petFile__","__ctFile__","__segFile__","__modelName__"],
+ "setVariablesComment":"this variables will get updated with local values like home and can be used to set variables further on",
+ "__tempBase__":"__home__/temp/iraemm",
+ "__segBase__":"/home/studen/software/src/iraemm/segmentation",
+ "__roiFile__":"testMask.nii.gz",
+ "__ctFile__":"testCT.nii.gz",
+ "__petFile__":"testPET.nii.gz",
+ "__segFile__":"segmentation.nii.gz",
+ "__modelName__":"DM_defaults.DM_train_VISCERAL16_Fold1.final.2019-10-01.07.46.19.932616.model.ckpt",
+ "tempBase":"__tempBase__",
+ "model":"__model__",
+ "project":"IPNUMMprospektiva/Study",
+ "targetSchema":"study",
+ "targetQuery":"Imaging1",
+ "viewName":"segmentationReview",
+ "participantField":"ParticipantId",
+ "segmentationSchema":"study",
+ "segmentationQuery":"Segmentations",
+ "reportQuery":"reportImages",
+ "reportSchema":"lists",
+ "percentileQuery":"SUVQuantiles",
+ "imageDir":"preprocessedImages",
+ "version":"v5",
+ "versionNumber":"5",
+ "images":{
+	"CT":{
+		"queryField":"ctResampled",
+		"tempFile":"__ctFile__"},
+	"PET":{
+		"queryField":"petResampled",
+		"tempFile":"__petFile__"},
+	"patientmask":{
+		"queryField":"ROImask",
+		"tempFile":"__roiFile__"},
+	"segmentation":{
+		"suffix":".nii.gz"
+	}
+ },
+ "replacePattern":{
+	 "__workDir__":"__tempBase__",
+	 "__roi__":"__tempBase__/__roiFile__",
+	 "__pet__":"__tempBase__/__petFile__",
+	 "__ct__":"__tempBase__/__ctFile__",
+	 "__seg__":"__tempBase__/__segFile__",
+	 "__model__":"__modelName__"
+ },
+ "nnUNet":{
+	 "ModelId":"501",
+	 "configuration":"3d_fullres",
+	 "env":{
+		"nnUNet_raw_data_base":"__tempBase__",
+		"nnUNet_preprocessed":"__tempBase__",
+		"RESULTS_FOLDER":"/home/studen/software/src/iraemmsegmentationmodels"
+	 }
+ },
+ "deepmedic": {
+	 "config":{
+		 "model":{
+		 	"template":"__segBase__/model/modelConfig.cfg.template",
+		 	"out":"__tempBase__/modelConfig.cfg"
+	 	},
+	 	"test":{
+			"template":"__segBase__/test/testConfig.cfg.template",
+		 	"out":"__tempBase__/testConfig.cfg"
+	 	},
+		"predictions":{
+			"template":"__segBase__/test/testNamesOfPredictions.cfg.template",
+			"out":"__tempBase__/testNamesOfPredictions.cfg"
+		},
+		"CT":{
+			"template":"__segBase__/test/testChannels_CT.cfg.template",
+			"out":"__tempBase__/testChannels_CT.cfg"
+		},
+		"ROI":{
+			"template":"__segBase__/test/testRoiMasks.cfg.template",
+			"out":"__tempBase__/testRoiMasks.cfg"
+		}
+
+
+
+	 }
+ }
+
+
+
+}