Compute the volume of the intersection between two vtkPolyData objects

[Robert_Sawko (Robert Manson-Sawko)]

Hello,

I would to compute a volume of the intersection between two vtkPolyData objects (or their bounding boxes). I have been mainly working with ParaView Python interface and created a programmable filter to display the intersection and print the volume

from vtk import vtkBooleanOperationPolyDataFilter
from vtk import vtkMassProperties

in1 = inputs[0].VTKObject
in2 = inputs[1].VTKObject
op = vtkBooleanOperationPolyDataFilter()
op.SetOperationToIntersection()
op.SetInputData(0, in1)
op.SetInputData(1, in2)
op.Update()
output.CopyStructure(op.GetOutput())
mass = vtkMassProperties()
mass.SetInputConnection(op.GetOutputPort())
mass.Update()
print(mass.GetVolume())

I started to test it and my initial euphoria quickly deflated when I realised that it doesn’t work for overlapping boxes:

(2781.218s) [paraview        ]    vtkPointLocator.cxx:845    ERR| vtkPointLocator (0x23010d70): No points to subdivide
(2781.219s) [paraview        ]vtkIntersectionPolyData:2410  WARN| No Intersection between objects 
(2781.219s) [paraview        ]vtkDistancePolyDataFilt:82     ERR| vtkDistancePolyDataFilter (0x2329da60): No points/cells to operate on
(2781.220s) [paraview        ]vtkDistancePolyDataFilt:82     ERR| vtkDistancePolyDataFilter (0x2329da60): No points/cells to operate on

I guess this is a limitation of the current implementation of the boolean filter? It seems to require an intersection to be nicely discretised. Is there a better way in VTK/ParaView to compute this quantity?

[lassoan (Andras Lasso)]

vtkBooleanOperationPolyDataFilter is one of the very, very few VTK filters that don’t really work. It may produce invalid output for simple, completely valid inputs. If you use it interactively and you find that if intersection computation fails or produces significant artifacts, then you can try to slightly change the pose of one input mesh. That sometimes fixes some issues. You can also try to change the mesh in other ways, subdivide, etc. to see if that helps.

However, there is a much better, more robust mesh Boolean operation implementation for vtk: vtkbool. At some point it was available in ParaView. It is also available in 3D Slicer (both in Python and C++). You can also build it from source. Hopefully in a couple of months it will be available as a VTK remote module that you can pip-install (@jcfr is working on the infrastructure for this).

[Robert_Sawko (Robert Manson-Sawko)]

Thanks for the pointer, @lassoan - I will definitely check it out. I tried a few workarounds with subdivision, but I couldn’t come up with anything that is robust enough for my target case (the above is my attempt at a minimal example).

Finally though, following the comment to a question from here I have implemented a simpler Monte Carlo estimation of this overlap which seems robust. I generate some points inside a reference cube, map them onto one of my cuboids using this expression and then use vtkExtractEnclosedPoints to check how many lie within within the other closed cuboid.

This is clearly a street-fighting mathematics solution, which I will aim to replace in the future - particularly if scalability becomes an issue.

[marcomusy]

the filter seems to work ok though, at least in this case, when you triangulate the meshes:

from vedo import Cube, show
c1 = Cube().triangulate().wireframe()
c2 = Cube().pos(0.5,0.4,0.3).rotateX(20).triangulate().wireframe()
cc = c1.boolean("intersect", c2)
print(cc.volume())
show(c1,c2,cc, axes=1)

[lassoan (Andras Lasso)]

Yes, in many cases the built-in VTK Boolean operation filter works. If you can afford to manually check the results and make manual adjustments as needed (e.g., slightly move meshes until there are no artifacts) then this performance may be sufficient.

[marcomusy]

true!! I just noticed by playing a bit with it that it doesn’t always work for all rotations, generating some artifacts:

It seems though that it’s because of the mesh resolution, if I apply the vtkLinearSubdivisionFilter then it’s ok:

from vedo import Cube, show
c1 = Cube()
c2 = Cube().pos(0.6,0.4,0.3).rotateX(30).rotateY(40)
c1.triangulate().subdivide(4, method=1).wireframe().alpha(0.1)
c2.triangulate().subdivide(4, method=1).wireframe().alpha(0.1)
cc = c1.boolean("intersect", c2).c("red5")
show(c1, c2, cc, f"Volume = {cc.volume()}", axes=1)

[Robert_Sawko (Robert Manson-Sawko)]

Thanks for this @marcomusy. So I think you’re saying that in my generic vtkPolyData case I simply need to triangulate my objects. In my target application I don’t really control their sizes so I may end up computing overlap between a very small object and a very large one so a fixed subdivision may need some additional logic to get enough points for the vtkBooleanOperationPolyDataFilter, but I think it’s doable.

Ultimately, I think these approaches (subdividing and Monte Carlo) can be made equivalent. I will try to abstract this steo in the target code and try both. Thanks.

[Freya_the_Goddess]

@Robert_Sawko

from vtk import

Hi Marco Musy,

I just download and try Vedo. I try to disable the internet connection and some examples can’t be run. Perhaps need to load data / object from directory, some people do not have stable internet connection. It is very great module.