Browse Source

Changing to Planar contour -> Ribbon model -> labelmap representatation conversion pathway as suggested in Slicer forums https://discourse.slicer.org/t/converting-segmentationnode-to-labelmap-takes-a-long-time-through-python/23696/6

Andrej Studen 2 months ago
parent
commit
e250127660
2 changed files with 133 additions and 13 deletions
  1. 2 1
      slicerScripts/config.json
  2. 131 12
      slicerScripts/parseDICOM.py

+ 2 - 1
slicerScripts/config.json

@@ -3,4 +3,5 @@
 	"schema":"study",
 	"schema":"study",
 	"query":"Imaging1",
 	"query":"Imaging1",
 	"baseDir":"/home/studen/temp",
 	"baseDir":"/home/studen/temp",
-"debug":0}
+   "debug":0
+}

+ 131 - 12
slicerScripts/parseDICOM.py

@@ -5,6 +5,8 @@ import numpy
 import zipfile
 import zipfile
 import shutil
 import shutil
 import os
 import os
+import signal
+from contextlib import contextmanager
 
 
 def main(configFile=None):
 def main(configFile=None):
 
 
@@ -17,6 +19,26 @@ def main(configFile=None):
    config.update(connectDB(config))
    config.update(connectDB(config))
    parseData(config,getMeanHeartDose)
    parseData(config,getMeanHeartDose)
 
 
+
+class TimeoutException(Exception): pass
+
+@contextmanager
+def time_limit(seconds):
+
+   def signal_handler(signum, frame):
+      raise TimeoutException("Timed out!")
+
+   signal.signal(signal.SIGALRM, signal_handler)
+   signal.alarm(seconds)
+
+   try:
+      yield
+   finally:
+      signal.alarm(0)
+
+
+
+
 def connectDB(setup):
 def connectDB(setup):
    nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixsuite')
    nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixsuite')
    sys.path.append(os.path.join(nixSuite,'wrapper'))
    sys.path.append(os.path.join(nixSuite,'wrapper'))
@@ -86,10 +108,17 @@ def getMeanHeartDose(config,r):
 #stores it as doseHeart to row r
 #stores it as doseHeart to row r
 #return True if r needs to be updated on the server
 #return True if r needs to be updated on the server
 #and False if r is unchanged
 #and False if r is unchanged
+
+   pid=r['ParticipantId']
    sid=r['orthancStudyId']
    sid=r['orthancStudyId']
    if not sid:
    if not sid:
-      print('No study for {}'.format(r['ParticipantId']))
+      print(f'No study for {pid}')
       return False
       return False
+
+   if r['comments']=='TIMEOUT':
+      print(f'Skipping {pid} - timeout')
+      return False
+   timeout=config.get('timeout',1200)
    doseHeart=r['doseHeart']
    doseHeart=r['doseHeart']
    if doseHeart:
    if doseHeart:
 #no need to update
 #no need to update
@@ -102,12 +131,19 @@ def getMeanHeartDose(config,r):
    msg=checkNodes(config,nodes)
    msg=checkNodes(config,nodes)
    if len(msg)>0:
    if len(msg)>0:
       r['comments']=msg
       r['comments']=msg
-   else:
-      r['doseHeart']=getMeanDose(nodes,'Heart')
+      clearDir(outDir)
+      removeNodes(nodes)
+      return True
+   print(f'Running with timeout={timeout} s')
+   r['doseHeart']=getMeanDoseRibbon(nodes,'Heart')
    clearDir(outDir)
    clearDir(outDir)
 #needs updating
 #needs updating
+#remove all created nodes from scene
+   removeNodes(nodes)
    return True
    return True
 
 
+
+
 def loadVolumes(dataDir):
 def loadVolumes(dataDir):
    nodeNames=[]
    nodeNames=[]
    with DICOMLib.DICOMUtils.TemporaryDICOMDatabase() as db:
    with DICOMLib.DICOMUtils.TemporaryDICOMDatabase() as db:
@@ -126,7 +162,8 @@ def loadVolumes(dataDir):
    nd=len(doseNodes)
    nd=len(doseNodes)
    print(f'vol:{nv} seg:{ns} dose: {nd}')
    print(f'vol:{nv} seg:{ns} dose: {nd}')
 
 
-   return {'vol':volumeNodes,'dose':doseNodes,'seg':segmentationNodes}
+   return {'all':nodes,'vol':volumeNodes,'dose':doseNodes,'seg':segmentationNodes}
+
 
 
 
 
 def checkNodes(config,nodes):
 def checkNodes(config,nodes):
@@ -141,26 +178,92 @@ def checkNodes(config,nodes):
       msg+=f'SEG[{nS}]'
       msg+=f'SEG[{nS}]'
    return msg
    return msg
 
 
+
+
+
+def removeNodes(nodes):
+   for node in nodes['all']:
+      slicer.mrmlScene.RemoveNode(node)
+
+
+
+
 def getMeanDose(nodes,target):
 def getMeanDose(nodes,target):
    segNode=nodes['seg'][0]
    segNode=nodes['seg'][0]
    seg=segNode.GetSegmentation()
    seg=segNode.GetSegmentation()
 
 
    segmentIds=seg.GetSegmentIDs()
    segmentIds=seg.GetSegmentIDs()
-#[seg.GetNthSegment(i) for i in range(seg.GetNumberOfSegments())]
 
 
    targetSegmentIds=[s for s in segmentIds if seg.GetSegment(s).GetName()==target]
    targetSegmentIds=[s for s in segmentIds if seg.GetSegment(s).GetName()==target]
-   #labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
-   #export=slicer.vtkSlicerSegmentationsModuleLogic.ExportSegmentsToLabelmapNode
-   #export(segNode, targetSegmentIds, labelmapVolumeNode, nodes['dose'][0])
-   #nodes.update({'labelMap':labelmapVolumeNode})
 
 
    doseNode=nodes['dose'][0]
    doseNode=nodes['dose'][0]
    doseArray = slicer.util.arrayFromVolume(doseNode)
    doseArray = slicer.util.arrayFromVolume(doseNode)
-   segmentArray = slicer.util.arrayFromSegmentBinaryLabelmap(segNode, targetSegmentIds[0], doseNode)
+   print('Dose array shape: {}'.format(doseArray.shape))
+   export=slicer.util.arrayFromSegmentBinaryLabelmap
+   segmentArray = export(segNode, targetSegmentIds[0], doseNode)
+   print('Segment array shape: {}'.format(segmentArray.shape))
+   doseVoxels = doseArray[segmentArray != 0]
+   meanDose=float(numpy.mean(doseVoxels))
+
+   print(f'Dose {meanDose}')
+#add a float() to avoid JSON complaining about float32 converion
+   return meanDose
+
+
+def getMeanDoseRibbon(nodes,target):
+   segNode=nodes['seg'][0]
+   seg=segNode.GetSegmentation()
+   seg.CreateRepresentation('Ribbon model')
+   seg.SetSourceRepresentationName('Ribbon model')
+
+   segmentIds=seg.GetSegmentIDs()
+
+   targetSegmentIds=[s for s in segmentIds if seg.GetSegment(s).GetName()==target]
+
+   doseNode=nodes['dose'][0]
+   doseArray = slicer.util.arrayFromVolume(doseNode)
+   print('Dose array shape: {}'.format(doseArray.shape))
+   export=slicer.util.arrayFromSegmentBinaryLabelmap
+   segmentArray = export(segNode, targetSegmentIds[0], doseNode)
+   print('Segment array shape: {}'.format(segmentArray.shape))
    doseVoxels = doseArray[segmentArray != 0]
    doseVoxels = doseArray[segmentArray != 0]
-   print(numpy.mean(doseVoxels))
+   meanDose=float(numpy.mean(doseVoxels))
+
+   print(f'Dose {meanDose}')
 #add a float() to avoid JSON complaining about float32 converion
 #add a float() to avoid JSON complaining about float32 converion
-   return float(numpy.mean(doseVoxels))
+   return meanDose
+
+
+
+def getMeanDoseSlicerRT(nodes,target):
+#use SlicerRT BatchProcessing example to extract LabelMaps
+#alternatove to getMeanDose since it was taking forever to build label map
+#no success, still takes forever
+   segNode=nodes['seg'][0]
+   seg=segNode.GetSegmentation()
+
+   segmentIds=seg.GetSegmentIDs()
+
+   targetSegmentIds=[s for s in segmentIds if seg.GetSegment(s).GetName()==target]
+
+
+   doseNode=nodes['dose'][0]
+   labelMapNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
+   export=slicer.vtkSlicerSegmentationsModuleLogic.ExportSegmentsToLabelmapNode
+   export(segNode, targetSegmentIds, labelMapNode, doseNode)
+
+   doseArray = slicer.util.arrayFromVolume(doseNode)
+   print('SlicerRT: Dose array shape: {}'.format(doseArray.shape))
+   segmentArray = slicer.util.arrayFromVolume(labelMapNode)
+   print('SlicerRT: Segment array shape: {}'.format(segmentArray.shape))
+   doseVoxels = doseArray[segmentArray != 0]
+#add a float() to avoid JSON complaining about float32 converion
+   meanDose= float(numpy.mean(doseVoxels))
+   print(f'SlicerRT: Dose {meanDose} Gy')
+   slicer.mrmlScene.removeNode(labelMapNode)
+   return meanDose
+
+
 
 
 def getDicomInstances(config,sid):
 def getDicomInstances(config,sid):
    odb=config['odb']
    odb=config['odb']
@@ -168,8 +271,11 @@ def getDicomInstances(config,sid):
    sd=odb.getStudyData(sid)
    sd=odb.getStudyData(sid)
    series=sd['Series']
    series=sd['Series']
    instances=[]
    instances=[]
+   validModalities=['RTSTRUCT','RTDOSE']
    for s in series:
    for s in series:
       sed=odb.getSeriesData(s)
       sed=odb.getSeriesData(s)
+      if sed['MainDicomTags']['Modality'] not in validModalities:
+         continue
       instances.extend(sed['Instances'])
       instances.extend(sed['Instances'])
    #download instances one by one
    #download instances one by one
    baseDir=config.get('baseDir',os.path.join(os.path.expanduser('~'),'temp'))
    baseDir=config.get('baseDir',os.path.join(os.path.expanduser('~'),'temp'))
@@ -182,6 +288,8 @@ def getDicomInstances(config,sid):
    return outDir
    return outDir
 
 
 
 
+
+
 def getDicomZip(config,sid):
 def getDicomZip(config,sid):
    ofb=config['ofb']
    ofb=config['ofb']
    baseDir=config.get('baseDir',os.path.join(os.path.expanduser('~'),'temp'))
    baseDir=config.get('baseDir',os.path.join(os.path.expanduser('~'),'temp'))
@@ -195,6 +303,9 @@ def getDicomZip(config,sid):
    os.remove(path)
    os.remove(path)
    return outDir
    return outDir
 
 
+
+
+
 def extractZip(config,fname):
 def extractZip(config,fname):
 #flattens the zip files in the baseDir/bname directory 
 #flattens the zip files in the baseDir/bname directory 
 #where bname is the basename of the file without the .zip suffix
 #where bname is the basename of the file without the .zip suffix
@@ -215,10 +326,18 @@ def extractZip(config,fname):
 
 
    return outDir
    return outDir
 
 
+
+
+
 def clearDir(outDir):
 def clearDir(outDir):
    if os.path.isdir(outDir):
    if os.path.isdir(outDir):
       shutil.rmtree(outDir)
       shutil.rmtree(outDir)
 
 
+
+
+
+
+
 if __name__=='__main__':
 if __name__=='__main__':
    try:
    try:
       main(sys.argv[1])
       main(sys.argv[1])