ENH: fixes #118

BoneReconstructionPlanner/BRPLib/helperFunctions.py

@@ -4,6 +4,7 @@
 
 from __main__ import vtk, slicer
 import numpy as np
+import logging
 
 def getIntersectionBetweenModelAnd1Plane(modelNode,planeNode,intersectionModel):
   plane = vtk.vtkPlane()
@@ -299,6 +300,8 @@ def createCylinder(name,R,H=50):
   #
   triangleFilter = vtk.vtkTriangleFilter()
   triangleFilter.SetInputConnection(tubeFilter.GetOutputPort())
+  # make strips valid
+  triangleFilter.PassLinesOff()
   #
   normalsFilter = vtk.vtkPolyDataNormals()
   normalsFilter.SetInputConnection(triangleFilter.GetOutputPort())
@@ -378,4 +381,233 @@ def getIntersectionPointsOfEachModelByMode(intersectionA,intersectionB,measureme
   if measurementMode == "proximal2proximal":
     return intersectionAtoBSimilarVectors[0][2], intersectionAtoBSimilarVectors[0][3]
   elif measurementMode == "distal2distal":
-    return intersectionAtoBSimilarVectors[-1][2], intersectionAtoBSimilarVectors[-1][3]
+    return intersectionAtoBSimilarVectors[-1][2], intersectionAtoBSimilarVectors[-1][3]
+
+class combineModelsRobustLogic:
+  def process(
+      inputModelA, 
+      inputModelB, 
+      outputModel, 
+      operation, 
+      numberOfRetries = 2, 
+      translateRandomly = 4, 
+      triangulateInputs = True
+    ):
+    """
+    Run the processing algorithm.
+    Can be used without GUI widget.
+    :param inputModelA: first input model node
+    :param inputModelB: second input model node
+    :param outputModel: result model node, if empty then a new output node will be created
+    :param operation: union, intersection, difference, difference2
+    :param numberOfRetries: number of retries if operation fails
+    :param translateRandomly: order of magnitude of the random translation
+    :param triangulateInputs: triangulate input models before boolean operation
+    """
+
+    if not inputModelA or not inputModelB or not outputModel:
+      raise ValueError("Input or output model nodes are invalid")
+
+    import time
+    startTime = time.time()
+    logging.info('Processing started')
+
+    import vtkSlicerCombineModelsModuleLogicPython as vtkbool
+
+    combine = vtkbool.vtkPolyDataBooleanFilter()
+
+    if operation == 'union':
+      combine.SetOperModeToUnion()
+    elif operation == 'intersection':
+      combine.SetOperModeToIntersection()
+    elif operation == 'difference':
+      combine.SetOperModeToDifference()
+    elif operation == 'difference2':
+      combine.SetOperModeToDifference2()
+    else:
+      raise ValueError("Invalid operation: "+operation)
+
+    transformToOutput = vtk.vtkGeneralTransform()
+    slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(inputModelA.GetParentTransformNode(), outputModel.GetParentTransformNode(), transformToOutput)
+    if transformToOutput is None:
+      transformToOutput = vtk.vtkTransform()
+    transformerA = vtk.vtkTransformPolyDataFilter()
+    transformerA.SetTransform(transformToOutput)
+    transformerA.SetInputData(inputModelA.GetPolyData())
+    transformerA.Update()
+    combine.SetInputData(0, transformerA.GetOutput())
+
+    transformToOutput = vtk.vtkGeneralTransform()
+    slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(inputModelB.GetParentTransformNode(), outputModel.GetParentTransformNode(), transformToOutput)
+    if transformToOutput is None:
+      transformToOutput = vtk.vtkTransform()
+    preTransformerB = vtk.vtkTransformPolyDataFilter()
+    preTransformerB.SetTransform(transformToOutput)
+    preTransformerB.SetInputData(inputModelB.GetPolyData())
+    preTransformerB.Update()
+    identityTransform = vtk.vtkTransform()
+    transformerB = vtk.vtkTransformPolyDataFilter()
+    transformerB.SetTransform(identityTransform)
+    transformerB.SetInputData(preTransformerB.GetOutput())
+    transformerB.Update()
+    combine.SetInputData(1, transformerB.GetOutput())
+
+    # first handle cases where inputs are not valid otherwise the boolean filter will crash
+    preBooleanOperationHandlingDone = False
+
+    modelAIsValid = transformerA.GetOutput().GetNumberOfPoints() != 0
+    modelBIsValid = transformerB.GetOutput().GetNumberOfPoints() != 0
+    modelAIsEmpty = not modelAIsValid
+    modelBIsEmpty = not modelBIsValid
+    bothModelsAreEmpty = modelAIsEmpty and modelBIsEmpty
+    onlyModelAIsValid = modelAIsValid and not modelBIsValid
+    onlyModelBIsValid = modelBIsValid and not modelAIsValid
+
+    if (
+      bothModelsAreEmpty or
+      ((operation == 'union') and onlyModelBIsValid) or
+      ((operation == 'union') and onlyModelAIsValid)
+    ):
+      appendFilter = vtk.vtkAppendPolyData()
+      appendFilter.AddInputData(transformerA.GetOutput())
+      appendFilter.AddInputData(transformerB.GetOutput())
+      combine = appendFilter
+      preBooleanOperationHandlingDone = True
+    elif (
+      (
+        (operation == 'intersection') and 
+        (onlyModelAIsValid or onlyModelBIsValid)
+      ) or
+      (
+        (operation == 'difference') and 
+        (onlyModelBIsValid)
+      ) or
+      (
+        (operation == 'difference2') and 
+        (onlyModelAIsValid)
+      )
+    ):
+      # return empty model
+      appendFilter = vtk.vtkAppendPolyData()
+      emptyPolyData = vtk.vtkPolyData()
+      appendFilter.AddInputData(emptyPolyData)
+      combine = appendFilter
+      preBooleanOperationHandlingDone = True
+    elif (
+      (operation == 'difference') and 
+      (onlyModelAIsValid)
+    ):
+      combine = transformerA
+      preBooleanOperationHandlingDone = True
+    elif (
+      (operation == 'difference2') and 
+      (onlyModelBIsValid)
+    ):
+      combine = transformerB
+      preBooleanOperationHandlingDone = True
+    # else:
+    #   combine is a vtkPolyDataBooleanFilter
+
+
+    # These parameters might be useful to expose:
+    # combine.MergeRegsOn()  # default off
+    # combine.DecPolysOff()  # default on
+    combine.Update()
+
+    collisionDetectionFilter = vtk.vtkCollisionDetectionFilter()
+    collisionDetectionFilter.SetInputData(0, transformerA.GetOutput())
+    collisionDetectionFilter.SetInputData(1, transformerB.GetOutput())
+    identityMatrix = vtk.vtkMatrix4x4()
+    collisionDetectionFilter.SetMatrix(0,identityMatrix)
+    collisionDetectionFilter.SetMatrix(1,identityMatrix)
+    collisionDetectionFilter.SetCollisionModeToFirstContact()
+
+    combineFilterSuccessful = combine.GetOutput().GetNumberOfPoints() != 0
+    if not combineFilterSuccessful and not preBooleanOperationHandlingDone:
+      for retry in range(numberOfRetries+1):
+        if (
+          operation == 'union'
+        ):
+          if retry == 0:
+            # check if the models are already intersecting
+            collisionDetectionFilter.Update()
+          if collisionDetectionFilter.GetNumberOfContacts() == 0:
+            # models do not touch so we append them
+            appendFilter = vtk.vtkAppendPolyData()
+            appendFilter.AddInputData(transformerA.GetOutput())
+            appendFilter.AddInputData(transformerB.GetOutput())
+            appendFilter.Update()
+            combine = appendFilter
+            break
+
+        if (
+          operation == 'intersection'
+        ):
+          if retry == 0:
+            # check if the models are already intersecting
+            collisionDetectionFilter.Update()
+          if collisionDetectionFilter.GetNumberOfContacts() == 0:
+            # models do not touch so we return an empty model
+            break
+
+        if (
+          operation == 'difference'
+        ):
+          if retry == 0:
+            # check if the models are already intersecting
+            collisionDetectionFilter.Update()
+          if collisionDetectionFilter.GetNumberOfContacts() == 0:
+            # models do not touch so we return modelA
+            combine = transformerA
+            break
+
+        if (
+          operation == 'difference2'
+        ):
+          if retry == 0:
+            # check if the models are already intersecting
+            collisionDetectionFilter.Update()
+          if collisionDetectionFilter.GetNumberOfContacts() == 0:
+            # models do not touch so we return modelB
+            combine = transformerB
+            break
+
+        if retry == 0 and triangulateInputs:
+          # in case inputs are not triangulated, triangulate them
+          triangulatedInputModelA = vtk.vtkTriangleFilter()
+          triangulatedInputModelA.SetInputData(inputModelA.GetPolyData())
+          # make strips valid
+          triangulatedInputModelA.PassLinesOff()
+          triangulatedInputModelA.Update()
+          transformerA.SetInputData(triangulatedInputModelA.GetOutput())
+          transformerA.Update()
+          triangulatedInputModelB = vtk.vtkTriangleFilter()
+          triangulatedInputModelB.SetInputData(inputModelB.GetPolyData())
+          # make strips valid
+          triangulatedInputModelB.PassLinesOff()
+          triangulatedInputModelB.Update()
+          preTransformerB.SetInputData(triangulatedInputModelB.GetOutput())
+          preTransformerB.Update()
+          transformerB.SetInputData(preTransformerB.GetOutput())
+
+        # retry with random translation if boolean operation fails
+        logging.info(f"Retrying boolean operation with random translation (retry {retry+1})")
+        transform = vtk.vtkTransform()
+        unitaryVector = [vtk.vtkMath.Random()-0.5 for _ in range(3)]
+        vtk.vtkMath.Normalize(unitaryVector)
+        import numpy as np
+        translationVector = np.array(unitaryVector) * (10**-translateRandomly)
+        transform.Translate(translationVector)
+        transformerB.SetTransform(transform)
+        transformerB.Update()
+        # recalculate the boolean operation
+        combine.SetInputData(1, transformerB.GetOutput())
+        combine.Update()
+
+    outputModel.SetAndObservePolyData(combine.GetOutput())
+    outputModel.CreateDefaultDisplayNodes()
+    # The filter creates a few scalars, don't show them by default, as they would be somewhat distracting
+    outputModel.GetDisplayNode().SetScalarVisibility(False)
+
+    stopTime = time.time()
+    logging.info('Processing completed in {0:.2f} seconds'.format(stopTime-startTime))

BoneReconstructionPlanner/BoneReconstructionPlanner.py

@@ -3443,7 +3443,7 @@ def makeBooleanOperationsToFibulaSurgicalGuideBase(self):
     biggerFibulaDentalImplantsCylindersModelsFolder = shNode.GetItemByName("Bigger Fibula Dental Implants Cylinders Models")
     biggerFibulaDentalImplantsCylindersModelsList = createListFromFolderID(biggerFibulaDentalImplantsCylindersModelsFolder)
 
-    combineModelsLogic = slicer.modules.combinemodels.widgetRepresentation().self().logic
+    combineModelsLogic = combineModelsRobustLogic
 
     surgicalGuideModel = slicer.modules.models.logic().AddModel(fibulaSurgicalGuideBaseModel.GetPolyData())
     surgicalGuideModel.SetName(slicer.mrmlScene.GetUniqueNameByString('FibulaSurgicalGuidePrototype'))
@@ -3823,7 +3823,7 @@ def makeBooleanOperationsToMandibleSurgicalGuideBase(self):
     biggerSawBoxesModelsFolder = shNode.GetItemByName("biggerSawBoxes Models")
     biggerSawBoxesModelsList = createListFromFolderID(biggerSawBoxesModelsFolder)
 
-    combineModelsLogic = slicer.modules.combinemodels.widgetRepresentation().self().logic
+    combineModelsLogic = combineModelsRobustLogic
 
     surgicalGuideModel = slicer.modules.models.logic().AddModel(mandibleSurgicalGuideBaseModel.GetPolyData())
     surgicalGuideModel.SetName(slicer.mrmlScene.GetUniqueNameByString('MandibleSurgicalGuidePrototype'))