# Compute the volume of the intersection between two vtkPolyData objects

**URL:** https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338
**Category:** Support
**Tags:** python
**Created:** [April 20, 2022, 2:06pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338 "2022-04-20T14:06:42Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![Robert\_Sawko](https://discourse.vtk.org/user_avatar/discourse.vtk.org/robert_sawko/32/5048_2.png) [@Robert\_Sawko](https://discourse.vtk.org/u/Robert_Sawko)
#### Post date: [April 20, 2022, 2:06pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/1 "2022-04-20T14:06:42Z")

</div>

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](https://docs.paraview.org/en/latest/ReferenceManual/pythonProgrammableFilter.html) to display the intersection and print the volume

```python
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())

```

 ![two_spheres_intersection](https://discourse.vtk.org/uploads/default/original/2X/3/309ceebaa5c0b7842f9ffe5dee45716d9e5f689d.jpeg)

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

```auto
(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?

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [April 21, 2022, 1:08am UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/2 "2022-04-21T01:08:28Z")

</div>

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](https://github.com/zippy84/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).

---

<div class="post-metadata">

### Author: ![Robert\_Sawko](https://discourse.vtk.org/user_avatar/discourse.vtk.org/robert_sawko/32/5048_2.png) [@Robert\_Sawko](https://discourse.vtk.org/u/Robert_Sawko)
#### Post date: [April 27, 2022, 3:50pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/3 "2022-04-27T15:50:34Z")

</div>

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](https://math.stackexchange.com/questions/2778389/how-to-compute-volume-of-intersection-of-non-axis-aligned-cuboids-in-3d) 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](https://math.stackexchange.com/questions/2265255/mapping-a-3d-point-inside-a-hexahedron-to-a-unit-cube#2265507) 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.

---

<div class="post-metadata">

### Author: ![marcomusy](https://discourse.vtk.org/user_avatar/discourse.vtk.org/marcomusy/32/95_2.png) [@marcomusy](https://discourse.vtk.org/u/marcomusy)
#### Post date: [April 28, 2022, 11:04am UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/4 "2022-04-28T11:04:31Z")

</div>

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

```python
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)

```

 ![Screenshot from 2022-04-28 12-57-40](https://discourse.vtk.org/uploads/default/original/2X/0/0892f04b6fd6007c307654f1d1be094c50f84ebf.png)

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [April 28, 2022, 12:00pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/5 "2022-04-28T12:00:38Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![marcomusy](https://discourse.vtk.org/user_avatar/discourse.vtk.org/marcomusy/32/95_2.png) [@marcomusy](https://discourse.vtk.org/u/marcomusy)
#### Post date: [April 28, 2022, 12:08pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/6 "2022-04-28T12:08:06Z")

</div>

true!! I just noticed by playing a bit with it that it doesn’t always work for all rotations, generating some artifacts:  
 ![Screenshot from 2022-04-28 14-03-39](https://discourse.vtk.org/uploads/default/original/2X/3/3874ba74dee294919fb1306b50b77af2e6d44137.png)

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

```python
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)

```

 ![Screenshot from 2022-04-28 14-06-38](https://discourse.vtk.org/uploads/default/original/2X/c/c50facf9fd29347dceb40b383eb74bf899705f2d.png)

---

<div class="post-metadata">

### Author: ![Robert\_Sawko](https://discourse.vtk.org/user_avatar/discourse.vtk.org/robert_sawko/32/5048_2.png) [@Robert\_Sawko](https://discourse.vtk.org/u/Robert_Sawko)
#### Post date: [May 1, 2022, 4:01pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/7 "2022-05-01T16:01:08Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![Freya\_the\_Goddess](https://discourse.vtk.org/user_avatar/discourse.vtk.org/freya_the_goddess/32/4942_2.png) [@Freya\_the\_Goddess](https://discourse.vtk.org/u/Freya_the_Goddess)
#### Post date: [May 1, 2022, 4:19pm UTC](https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/8 "2022-05-01T16:19:36Z")

</div>

> [@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.
