Segfault while finding intersection of vtkPolyData with python

Visualizing two polyData objects, I see that they overlap, however, I cannot calculate their intersection with VTK in python. What am I doing wrong??? I get the following error:

$ python test.py
VOLUME of probe1: 8.360612727741916 # of Polys: 344
VOLUME of probe2: 8.440116557580266 # of Polys: 480
Warning: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Filters/General/vtkIntersectionPolyDataFilter.cxx, line 1675
vtkIntersectionPolyDataFilter (0x7fe2c5463bb0): No cell with correct orientation found

ERROR: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Common/ExecutionModel/vtkExecutive.cxx, line 784
vtkCompositeDataPipeline (0x7fe2c54649e0): Algorithm vtkIntersectionPolyDataFilter(0x7fe2c5463bb0) returned failure for request: vtkInformation (0x7fe2c54654b0)
Debug: Off
Modified Time: 1314
Reference Count: 1
Registered Events: (none)
Request: REQUEST_DATA
FORWARD_DIRECTION: 0
ALGORITHM_AFTER_FORWARD: 1
FROM_OUTPUT_PORT: 0

ERROR: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Common/ExecutionModel/vtkExecutive.cxx, line 784
vtkCompositeDataPipeline (0x7fe2c545eee0): Algorithm vtkBooleanOperationPolyDataFilter(0x7fe2c545f530) returned failure for request: vtkInformation (0x7fe2c545fb50)
Debug: Off
Modified Time: 1069
Reference Count: 1
Registered Events: (none)
Request: REQUEST_DATA
FORWARD_DIRECTION: 0
ALGORITHM_AFTER_FORWARD: 1
FROM_OUTPUT_PORT: 0

0

My code is below, inputs are probe0.vtk and probe1.vtk
import vtk

def get_volume(vtk_poly):
    if vtk_poly.GetNumberOfPolys() == 0:
        return 0
    Mass = vtk.vtkMassProperties()
    Mass.SetInputData(vtk_poly)
    Mass.Update()
    return Mass.GetVolume()

def read_polydata(vtkfilename):
    reader_poly = vtk.vtkPolyDataReader()
    reader_poly.SetFileName(vtkfilename)
    reader_poly.Update()
    return reader_poly.GetOutput()

def get_mesh_intersection(mesh1, mesh2):
    '''
    Find the intersection volume between mesh vtk objects
    and return the volume of intersection
    '''
    booleanOperationFilter = vtk.vtkBooleanOperationPolyDataFilter()
    #booleanOperationFilter.GlobalWarningDisplayOff()
    booleanOperationFilter.SetInputData(0, mesh1)
    booleanOperationFilter.SetInputData(1, mesh2)
    booleanOperationFilter.SetOperationToIntersection()
    booleanOperationFilter.Update()
    intersection = booleanOperationFilter.GetOutput()
    return get_volume(intersection)



mesh_i = read_polydata('probe0.vtk')
mesh_j = read_polydata('probe1.vtk')
print('VOLUME of probe1:\t', get_volume(mesh_i), '# of Polys:\t', mesh_i.GetNumberOfPolys())
print('VOLUME of probe2:\t', get_volume(mesh_j), '# of Polys:\t', mesh_j.GetNumberOfPolys())

print(get_mesh_intersection(mesh_i, mesh_j))

This is very similar to How to extract surface from polydata

The vtk.vtkBooleanOperationPolyDataFilter is simply not producing a result which is causing the segfault down your pipeline. @cory.quammen would likely have more details on how to get the vtk.vtkBooleanOperationPolyDataFilter working

Also, note a lot of this code can be simplified with vtki. Perhaps if we figure out how to implement this, it could become a common routine in vtki

import vtki
import vtk

def get_mesh_intersection(mesh1, mesh2):
    '''
    Find the intersection volume between mesh vtk objects
    and return the volume of intersection
    '''
    alg = vtk.vtkBooleanOperationPolyDataFilter()
    alg.SetInputData(0, mesh2)
    alg.SetInputData(1, mesh1)
    alg.SetOperationToIntersection()
    alg.Update()
    intersection = vtki.wrap(alg.GetOutput())
    return intersection

# Load meshes
mesh_i = vtki.read('probe0.vtk')
mesh_j = vtki.read('probe1.vtk')

print('VOLUME of probe1:\t', mesh_i.volume, '# of Polys:\t', mesh_i.n_faces)
print('VOLUME of probe2:\t', mesh_j.volume, '# of Polys:\t', mesh_j.n_faces)
intersection = get_mesh_intersection(mesh_i, mesh_j)
print('INTERSECTION VOLUME:\t', intersection.volume)

with output:

VOLUME of probe1:	 8.360612727741916 # of Polys:	 344
VOLUME of probe2:	 8.440116557580266 # of Polys:	 480
INTERSECTION VOLUME:	 0.0

Also, plotting is easy:

vtki.MultiBlock([mesh_i, mesh_j]).plot(show_edges=1, multi_colors=1)

or

p = vtki.Plotter()
p.add_mesh(mesh_i, color='g', show_edges=1)
p.add_mesh(mesh_j, color='b', show_edges=1)
p.show()

thank you so much! this is extremely helpful, and vtki makes working with vtk so much more pythonic. Thank you again!

I’ve used your exact example and I’m getting similar results but I also get errors. See below.

Could this be an issue with my installation? I’ve used pip install on a MacOS 10.14.2. to install both vtk and vtki

VOLUME of probe1: 8.360612727741916 # of Polys: 344
VOLUME of probe2: 8.440116557580266 # of Polys: 480
Warning: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Filters/General/vtkIntersectionPolyDataFilter.cxx, line 1675
vtkIntersectionPolyDataFilter (0x7ff71e259ff0): No cell with correct orientation found

ERROR: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Common/ExecutionModel/vtkExecutive.cxx, line 784
vtkCompositeDataPipeline (0x7ff71e25b4c0): Algorithm vtkIntersectionPolyDataFilter(0x7ff71e259ff0) returned failure for request: vtkInformation (0x7ff71e2caa90)
Debug: Off
Modified Time: 1787
Reference Count: 1
Registered Events: (none)
Request: REQUEST_DATA
FORWARD_DIRECTION: 0
ALGORITHM_AFTER_FORWARD: 1
FROM_OUTPUT_PORT: 0

ERROR: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Common/ExecutionModel/vtkExecutive.cxx, line 784
vtkCompositeDataPipeline (0x7ff71e2dfb00): Algorithm vtkBooleanOperationPolyDataFilter(0x7ff71e249650) returned failure for request: vtkInformation (0x7ff71bedff80)
Debug: Off
Modified Time: 1542
Reference Count: 1
Registered Events: (none)
Request: REQUEST_DATA
FORWARD_DIRECTION: 0
ALGORITHM_AFTER_FORWARD: 1
FROM_OUTPUT_PORT: 0

ERROR: In /Users/prabhu/src/git/VTKPythonPackage/standalone-build/VTK-source/Filters/Core/vtkMassProperties.cxx, line 95
vtkMassProperties (0x7ff71e259210): No data to measure…!

INTERSECTION VOLUME: 0.0

That’s not an issue with the installation or anything but is purely the vtk.vtkBooleanOperationPolyDataFilter (I am also getting that error). The more I read up on the vtk.vtkBooleanOperationPolyDataFilter, the less I understand how it works… but I think the cells in the two input meshes have to have specific orientations.

FYI: http://vtk.1045678.n5.nabble.com/How-to-use-vtkBooleanOperationPolyDataFilter-correctly-td5728839.html

why intersection.volume is 0.0 even they are joining?

Unfortunately, vtkBooleanOperationPolyDataFilter is very unreliable. It can give incorrect output for simple, completely valid input or can randomly crash. Classes in vtkbool work much better and its maintainer is usually willing to look into problems and fix them.