Find the intersection of two models

Support
1

Hi there,

How can i find the points in which two models (a sort of cylinder and a bar) intersect one another?
i’ve tried messing up with numpy trying to search for intersections in data arrays but that doesn’t seem the right way.

Thanks!

image

2

Hi there,

A bit convoluted but here are some code examples that might get youon the right track useing the modeltomodel distance extension:

vtkOutput = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode")

parameters = {}
parameters["vtkFile1"] = preprocessedSurface
parameters["vtkFile2"] = modelNode1
parameters['distanceType'] = "absolute_closest_point"
parameters["vtkOutput"] = vtkOutput

cliNode = slicer.cli.runSync(slicer.modules.modeltomodeldistance, None, parameters)

VTKFieldData = vtkOutput.GetMesh().GetAttributesAsFieldData(0)
mesh = vtkOutput.GetMesh()

resultTableNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode", "defaultTable")
        
for j in range(0,VTKFieldData.GetNumberOfArrays()):
    #print(VTKFieldData.)
    resultTableNode.AddColumn(VTKFieldData.GetArray(j)) 
        
resultTableNode.AddColumn(mesh.GetPoints().GetData())
        
tableStorageNode = resultTableNode.CreateDefaultStorageNode()
tableStorageNode.SetFileName("TEST.csv")
tableStorageNode.WriteData(resultTableNode)

Haven’t tested the code in a while but this should generate a table of the output node, where you can look at the distances and see which are below a certain threshold and assume they are touching at those points.

Happy to answer FU quesions, let me know.

3

Thanks for the answer. I figured it out by myself just five minutes ago.

I’ll post the code, i used using ray-casting if anyone needs it.
Is it accurate enough or should i implement it your way? @nrex45.

#loading the model mesh
mesh = getNode(model_name).GetMesh()
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputData(mesh)

#defining the coords of the tube (could be any sort of line-like obj i think)
pSource = points[0]
pTarget = points[1]

#evalutaing the actual intersections
obbTree = vtk.vtkOBBTree()
obbTree.SetDataSet(mesh)
obbTree.BuildLocator()
pointsVTKintersection = vtk.vtkPoints()
code = obbTree.IntersectWithLine(pSource, pTarget, pointsVTKintersection, None)
pointsVTKIntersectionData = pointsVTKintersection.GetData()
noPointsVTKIntersection = pointsVTKIntersectionData.GetNumberOfTuples()

#dealing with the intersection points
pointsIntersection = [] #where the intersections will be stored for future usage
for index in range(noPointsVTKIntersection):
    point_of_intersection = pointsVTKIntersectionData.GetTuple3(index)
    pointsIntersection.append(point_of_intersection)
   (slicer.mrmlScene.AddNewNodeByClass( 'vtkMRMLMarkupsFiducialNode', ('Centro coronale ' +   str(index)))).AddControlPointWorld(point_of_intersection)

return pointsIntersection #note that this is present only because in my code is part of a function that returns the coords, the actual drawing of the markups happens even if you decide to don't store the point in any variable
1 Like
4

This method seems to find the intersection also.

# Expect console erors.
slicer.mrmlScene.RemoveNode(sphereModel)
slicer.mrmlScene.RemoveNode(cylinderModel)
slicer.mrmlScene.RemoveNode(intersectionModel_0)
slicer.mrmlScene.RemoveNode(intersectionModel_1)

cylinder = vtk.vtkCylinderSource()
cylinder.SetRadius(5.0)
cylinder.SetHeight(100.0)
cylinder.SetResolution(45)
cylinder.Update()
cylinderTriangulator = vtk.vtkTriangleFilter()
cylinderTriangulator.SetInputConnection(cylinder.GetOutputPort())
cylinderTriangulator.Update()
cylinderModel = slicer.modules.models.logic().AddModel(cylinderTriangulator.GetOutput())
cylinderModel.SetName("cylinderModel")
cylinderModel.GetDisplayNode().SetOpacity(0.2)
cylinderModel.GetDisplayNode().SetColor(1.0, 0.0, 0.0)


sphere = vtk.vtkSphereSource()
sphere.SetRadius(25.0)
sphere.SetPhiResolution(45)
sphere.SetThetaResolution(45)
sphere.Update()
sphereTriangulator = vtk.vtkTriangleFilter()
sphereTriangulator.SetInputConnection(sphere.GetOutputPort())
sphereTriangulator.Update()
sphereModel = slicer.modules.models.logic().AddModel(sphereTriangulator.GetOutput())
sphereModel.SetName("sphereModel")
sphereModel.GetDisplayNode().SetOpacity(0.2)
sphereModel.GetDisplayNode().SetColor(0.0, 1.0, 0.0)


intersector = vtk.vtkIntersectionPolyDataFilter()
intersector.SetInputConnection(0, cylinderTriangulator.GetOutputPort())
intersector.SetInputConnection(1, sphereTriangulator.GetOutputPort())
intersector.Update()

intersectionModel_0 = slicer.modules.models.logic().AddModel(intersector.GetOutputDataObject(0))
intersectionModel_0.SetName("intersectionModel_0")
intersectionModel_0.GetDisplayNode().SetColor(0.0, 0.0, 1.0)

"""
# Cylinder (first input) minus intersection
intersectionModel_1 = slicer.modules.models.logic().AddModel(intersector.GetOutputDataObject(1))
intersectionModel_1.SetName("intersectionModel_1")
intersectionModel_1.GetDisplayNode().SetColor(1.0, 1.0, 0.0)
"""
2 Likes