pythonScripts/runStat.py

@@ -18,6 +18,12 @@ def main(parFile='../templates/statistics.json'):
     forceSUVMax=setup.get('forceSUVMax',False)
     forceLiver15=setup.get('forceLiver15',False)
     forceGlobal=setup.get('forceGlobal',False)
+    forceSUV4=setup.get('forceSUV4',False)
+    forceSUV4=setup.get('forceSUV4',False)
+    forceSUVMaxLesions=setup.get('forceSUVMaxLesions',False)
+    forceSUVMaxMetastases=setup.get('forceSUVMaxMetastases',False)
+
+    petField=setup.get('petField','petResampled')
 
     doUpload=setup.get('doUpload',True)
 
@@ -69,18 +75,33 @@ def main(parFile='../templates/statistics.json'):
         setup['SUVdataset']='SUVanalysis_SUVmax'
         suvMaxDone=not forceSUVMax and checkData(setup,r) 
 
+        setup['SUVdataset']='SUVanalysis_SUVmaxMetastases'
+        suvMaxMetastasesDone=not forceSUVMaxMetastases and checkData(setup,r)
+
+        setup['SUVdataset']='SUVanalysis_SUVmaxLesions'
+        suvMaxLesionsDone=not forceSUVMaxLesions and checkData(setup,r)
+
         setup['SUVdataset']='SUVanalysis_liver1p5'
         liver1p5Done=not forceLiver15 and checkData(setup,r)
 
-        doneCode=f'({globalDone}/{liverDone}/{liver1p5Done}/{suvMaxDone})'
-        print(f'Done: (global/liver/liver1p5/suvMax): {doneCode}')
+        setup['SUVdataset']='SUVanalysis_SUV4'
+        suv4Done=not forceSUV4 and checkData(setup,r)
 
-        if globalDone and liverDone and suvMaxDone and liver1p5Done:
+        doneCode=f'({globalDone}/{liverDone}/{liver1p5Done}/{suvMaxDone}/{suv4Done}/{suvMaxLesionsDone}/{suvMaxMetastasesDone})'
+        codeDesc='(global/liver/liver1p5/suvMax/suv4/suvMaxLesions/suvMaxMetastases)'
+        idVar='ParticipantId'
+        pid=r[idVar]
+        visitCode=r['visitCode']
+        print(f'Done[{pid}/{visitCode}: {codeDesc}: {doneCode}')
+
+        if globalDone and liverDone and suvMaxDone and \
+                liver1p5Done and suv4Done and \
+                suvMaxLesionsDone and suvMaxMetastasesDone:
             print('Skipping {} {}'.format(r['ParticipantId'],r['visitCode']))
             continue
 
         #PET
-        for q in ['petResampled']:
+        for q in [petField]:
             localPath=statUtils.getImage(setup,r,q)
         if localPath=="NONE":
             continue
@@ -96,7 +117,9 @@ def main(parFile='../templates/statistics.json'):
             print('Loaded {}/{}'.format(users[x],segPaths[x]))
 
         seg={x:SimpleITK.ReadImage(segPaths[x]) for x in segPaths}
-
+        #check dimensions of seg and pet
+        if not matchingDimensions(pet,seg):
+            continue
 
         globalThreshold={x:0 for x in seg}
         if not globalDone:
@@ -110,26 +133,18 @@ def main(parFile='../templates/statistics.json'):
         setup['SUVdataset']='SUVanalysis'
         outputs=loadGlobals(setup,r,seg)
 
-        #liver threshold
+        #this are IDs from labkey
         liverId=1
-        liverThreshold={x:outputs[x].get(liverId,default)['SUVmean']
-                +2*outputs[x].get(liverId,default)['SUVSD'] 
-                for x in outputs}
-        liver1p5Threshold={x:1.5*outputs[x].get(liverId,default)['SUVmean']
-                +2*outputs[x].get(liverId,default)['SUVSD'] 
-                for x in outputs}
         lesionId=4
         bmId=3
-        suvMax={x:numpy.max([outputs[x].get(lesionId,default)['SUVmax'],
-                    outputs[x].get(bmId,default)['SUVmax']])
-                for x in outputs}
-        suvMaxThreshold={x:0.41*suvMax[x] for x in suvMax}
+        metastasesId=6
     
-        ftStr='thr[liver]={} thr[liver/1.5]={} thr(suvmax)={}'
-        print(ftStr.format(liverThreshold,liver1p5Threshold,suvMaxThreshold))
-       
-
+        #liver threshold
         if not liverDone:
+            liverThreshold={x:outputs[x].get(liverId,default)['SUVmean']
+                +2*outputs[x].get(liverId,default)['SUVSD'] 
+                for x in outputs}
+            print(f'thr[liver]={liverThreshold}')
             setup['SUVdataset']='SUVanalysis_liver'
             liverOutputs=thresholdAnalysis(setup,r,pet,seg,liverThreshold)
             if doUpload:
@@ -137,26 +152,65 @@ def main(parFile='../templates/statistics.json'):
 
         #also for threshold=1.5
         if not liver1p5Done:
+            liver1p5Threshold={x:1.5*outputs[x].get(liverId,default)['SUVmean']
+                +2*outputs[x].get(liverId,default)['SUVSD'] 
+                for x in outputs}
+            print(f'thr[liver/1.5]={liver1p5Threshold}')
             setup['SUVdataset']='SUVanalysis_liver1p5'
             liver1p5Outputs=thresholdAnalysis(setup,r,pet,seg,liver1p5Threshold)
             if doUpload:
                 uploadData(setup,r,liver1p5Outputs)
 
+        #threshold=0.41 SUVmax
         if not suvMaxDone:
+            idArr=[lesionId,bmId,metastasesId]
+            fM=numpy.max
+            suvMax={x:fM([outputs[x].get(y,default)['SUVmax'] for y in idArr])}
+            suvMaxThreshold={x:0.41*suvMax[x] for x in suvMax}
+            print(f'thr[SUVmax]={suvMaxThreshold}')
+
             setup['SUVdataset']='SUVanalysis_SUVmax'
             suvMaxOutputs=thresholdAnalysis(setup,r,pet,seg,suvMaxThreshold)
             if doUpload:
                 uploadData(setup,r,suvMaxOutputs)
 
-        #skip threshold of 4
-        doThreshold4=False
-        if doThreshold4:
-            #threshold of 4
-            setup['radiomics']['setting']['resegmentRange']=[4]
-            setup['radiomics']['setting']['resegmentShape']=True
-            outputs4=getValues(setup,r,pet,seg)
+        #threshold=0.41 SUVMax (Lesion)
+        if not suvMaxLesionsDone:
+            suvMaxLesionsThreshold=\
+                {x:0.41*outputs[x].get(lesionId,default)['SUVmax'] \
+                    for x in outputs}
+            print(f'thr[SUVmaxLesions]={suvMaxLesionsThreshold}')
+
+            setup['SUVdataset']='SUVanalysis_SUVmaxLesions'
+            suvMaxLesionsOutputs=thresholdAnalysis(setup,r,pet,seg,
+                    suvMaxLesionsThreshold)
+            if doUpload:
+                uploadData(setup,r,suvMaxLesionsOutputs)
+
+
+        #threshold=0.41 SUVMax (Metastases)
+        if not suvMaxMetastasesDone:
+            suvMaxMetastasesThreshold=\
+                {x:0.41*outputs[x].get(metastasesId,default)['SUVmax'] \
+                    for x in outputs}
+            print(f'thr[SUVmaxMetastases]={suvMaxMetastasesThreshold}')
+
+            setup['SUVdataset']='SUVanalysis_SUVmaxMetastases'
+            suvMaxMetastasesOutputs=thresholdAnalysis(setup,r,pet,seg,\
+                    suvMaxMetastasesThreshold)
+            if doUpload:
+                uploadData(setup,r,suvMaxMetastasesOutputs)
+
+
+
+
+        if not suv4Done:
+            suv4Threshold={x:4 for x in outputs}
+            print(f'thr[SUV4]={suv4Threshold}')
             setup['SUVdataset']='SUVanalysis_SUV4'
-            uploadData(setup,r,outputs4)
+            suv4Outputs=thresholdAnalysis(setup,r,pet,seg,suv4Threshold)
+            if doUpload:
+                uploadData(setup,r,suv4Outputs)
 
         
         #cleanup
@@ -168,6 +222,7 @@ def main(parFile='../templates/statistics.json'):
 
 def loadGlobals(setup,r,seg):
     matchingVisits={
+            'VISIT_0':['VISIT_0'],
             'VISIT_1':['VISIT_1','VISIT_11'],
             'VISIT_11':['VISIT_1','VISIT_11'],
             'VISIT_2':['VISIT_2','VISIT_12'],
@@ -316,11 +371,38 @@ def checkData(setup,r):
     qFilter.append({'variable':idVar,'value':r[idVar],'oper':'eq'})
     qFilter.append({'variable':'visitCode','value':r['visitCode'],'oper':'eq'})
     ds=setup['db'].selectRows(p,'study',d,qFilter)
-    n=len(ds['rows'])
+    try:
+        n=len(ds['rows'])
+    except KeyError:
+        #contraintuitevly return as if data were there
+        #if ds has no rows, which typically happens
+        #on failure to retrieve dataset, which is
+        #related to missing datasets
+        print(f'Dataset [{d}] not found. Returning OK')
+        return True
     print('[{}:{}/{}] got {} rows.'.format(d,r[idVar],r['visitCode'],n))
     return n>0
 
+def vEqual(a,b,eps=1e-8):
+    delta=numpy.array(a)-numpy.array(b)
+    aDelta=numpy.abs(delta)
+    #print(f'Delta={delta} abs(Delta)={aDelta}')
+    return all(aDelta<eps)
 
+def matchingDimensions(pet,seg):
+    #QA routine
+    #check if dimensnions of pet and seg match
+    #seg is a dict {'ID':sitkImage}
+    
+    
+    #print('[PET] size={} origin={}'.format(pet.GetSize(),pet.GetOrigin()))
+    sPET=pet.GetSize()
+    ok=[vEqual(seg[x].GetSize(),sPET) for x in seg]
+    if not all(ok):
+        sz={x:seg[x].GetSize() for x in seg}
+        print('Size mismatch PET[{}]: {}'.format(sPET,sz))
+        return False
+    return True
 
 if __name__=='__main__':
     main(sys.argv[1])

pythonScripts/statUtils.py

@@ -65,10 +65,11 @@ def getImage(setup, row, field, extraPath=None):
     return localPath
 
 def getSegmentations(setup,row):
+    segName=setup.get('segmentationDataset','Segmentations')
     idField=setup['idField']
     visitField=setup['visitField']
     qFilter=[{'variable':x,'value':row[x],'oper':'eq'} for x in [idField,visitField]]
-    dsSeg=setup['db'].selectRows(setup['project'],'study','Segmentations',qFilter)
+    dsSeg=setup['db'].selectRows(setup['project'],'study',segName,qFilter)
     if len(dsSeg['rows'])<1:
         print('Failed to find segmentation for {}/{}'.format(row[idField],row[visitField]))
         return {0:"NONE"}

templates/statistics.json

@@ -2,9 +2,11 @@
  "project":"limfomiPET/Study",
  "imagingDataset":"Imaging1",
  "SUVdataset":"SUVanalysis",
+ "segmentationDataset":"Segmentations",
  "localDir":"_home_/temp/limfomiPET",
  "idField":"ParticipantId",
  "visitField":"visitCode",
+ "petField":"petResampled",
  "commentParticipants":["1487/17"],
  "commentVisits":["VISIT_2"],
  "n":20,