c++vtkcsgvedo

How to compute boolean intersection between 3 geometries using vtk?

I'm trying to compute the intersection between 3 geometries using VTK.

I've started with the BooleanOperationPolyDataFilter example and made minor tweaks such as:

Adding a third sphere:

vtkNew<vtkSphereSource> sphereSource3;
    sphereSource3->SetCenter(0.25, 0, .5);
    sphereSource3->SetPhiResolution(21);
    sphereSource3->SetThetaResolution(21);
    sphereSource3->Update();
    
    input3 = sphereSource3->GetOutput();

Attempted to plug the result of the 1st boolean intersection (between 1st two spheres) and then as a test as Union with a third sphere:

booleanOperation->SetInputData(0, input1);
  booleanOperation->SetInputData(1, input2);
  
  vtkNew<vtkBooleanOperationPolyDataFilter> booleanOperation2;
  booleanOperation2->SetOperationToUnion();
  booleanOperation2->SetInputData(0, booleanOperation->GetOutput());
  booleanOperation2->SetInputData(1, input3);
 
  vtkNew<vtkPolyDataMapper> booleanOperationMapper;
  booleanOperationMapper->SetInputConnection(booleanOperation2->GetOutputPort());
  booleanOperationMapper->ScalarVisibilityOff();

Full code listing:

#include <vtkActor.h>
#include <vtkBooleanOperationPolyDataFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSmartPointer.h>
#include <vtkTriangleFilter.h>
 
#include <vtkBYUReader.h>
#include <vtkOBJReader.h>
#include <vtkPLYReader.h>
#include <vtkPolyDataReader.h>
#include <vtkSTLReader.h>
#include <vtkXMLPolyDataReader.h>
 
#include <vtkCylinderSource.h>
#include <vtkSphereSource.h>
#include <vtksys/SystemTools.hxx>
 
#include <chrono>
#include <ctime>
 
#include <vtkCamera.h>
namespace {
vtkSmartPointer<vtkPolyData> ReadPolyData(const char* fileName);
void PositionCamera(vtkRenderer* renderer, double* viewUp, double* position);
} // namespace
 
int main(int argc, char* argv[])
{
 
  const clock_t begin_time = clock();
 
  vtkSmartPointer<vtkPolyData> input1;
  vtkSmartPointer<vtkPolyData> input2;
  vtkSmartPointer<vtkPolyData> input3;
 
  std::string operation("intersection");
 
  if (argc == 4)
  {
    auto poly1 = ReadPolyData(argv[1]);
    vtkNew<vtkTriangleFilter> tri1;
    tri1->SetInputData(poly1);
    vtkNew<vtkCleanPolyData> clean1;
    clean1->SetInputConnection(tri1->GetOutputPort());
    clean1->Update();
    input1 = clean1->GetOutput();
 
    auto poly2 = ReadPolyData(argv[3]);
    vtkNew<vtkTriangleFilter> tri2;
    tri2->SetInputData(poly2);
    tri2->Update();
    vtkNew<vtkCleanPolyData> clean2;
    clean2->SetInputConnection(tri2->GetOutputPort());
    clean2->Update();
    input2 = clean2->GetOutput();
    operation = argv[2];
  }
  else
  {
    vtkNew<vtkSphereSource> sphereSource1;
    sphereSource1->SetCenter(0.25, 0, 0);
    sphereSource1->SetPhiResolution(21);
    sphereSource1->SetThetaResolution(21);
    sphereSource1->Update();
    input1 = sphereSource1->GetOutput();
 
    vtkNew<vtkSphereSource> sphereSource2;
    sphereSource2->Update();
    input2 = sphereSource2->GetOutput();
 
    vtkNew<vtkSphereSource> sphereSource3;
    sphereSource3->SetCenter(0.25, 0, .5);
    sphereSource3->SetPhiResolution(21);
    sphereSource3->SetThetaResolution(21);
    sphereSource3->Update();
    
    input3 = sphereSource3->GetOutput();
 
 
    if (argc == 2)
    {
      operation = argv[1];
    }
  }
 
  vtkNew<vtkNamedColors> colors;
 
  vtkNew<vtkPolyDataMapper> input1Mapper;
  input1Mapper->SetInputData(input1);
  input1Mapper->ScalarVisibilityOff();
  vtkNew<vtkActor> input1Actor;
  input1Actor->SetMapper(input1Mapper);
  input1Actor->GetProperty()->SetDiffuseColor(
      colors->GetColor3d("Tomato").GetData());
  input1Actor->GetProperty()->SetSpecular(0.6);
  input1Actor->GetProperty()->SetSpecularPower(20);
  input1Actor->SetPosition(input1->GetBounds()[1] - input1->GetBounds()[0], 0,
                           0);
  vtkNew<vtkPolyDataMapper> input2Mapper;
  input2Mapper->SetInputData(input2);
  input2Mapper->ScalarVisibilityOff();
  vtkNew<vtkActor> input2Actor;
  input2Actor->SetMapper(input2Mapper);
  input2Actor->GetProperty()->SetDiffuseColor(
      colors->GetColor3d("Mint").GetData());
  input2Actor->GetProperty()->SetSpecular(0.6);
  input2Actor->GetProperty()->SetSpecularPower(20);
  input2Actor->SetPosition(-(input1->GetBounds()[1] - input1->GetBounds()[0]),
                           0, 0);
  vtkNew<vtkBooleanOperationPolyDataFilter> booleanOperation;
  if (operation == "union")
  {
    booleanOperation->SetOperationToUnion();
  }
  else if (operation == "intersection")
  {
    booleanOperation->SetOperationToIntersection();
  }
  else if (operation == "difference")
  {
    booleanOperation->SetOperationToDifference();
  }
  else
  {
    std::cout << "Unknown operation: " << operation << std::endl;
    return EXIT_FAILURE;
  }
  booleanOperation->SetInputData(0, input1);
  booleanOperation->SetInputData(1, input2);
  
  vtkNew<vtkBooleanOperationPolyDataFilter> booleanOperation2;
  booleanOperation2->SetOperationToUnion();
  booleanOperation2->SetInputData(0, booleanOperation->GetOutput());
  booleanOperation2->SetInputData(1, input3);
 
  vtkNew<vtkPolyDataMapper> booleanOperationMapper;
  booleanOperationMapper->SetInputConnection(booleanOperation2->GetOutputPort());
  booleanOperationMapper->ScalarVisibilityOff();
 
  vtkNew<vtkActor> booleanOperationActor;
  booleanOperationActor->SetMapper(booleanOperationMapper);
  booleanOperationActor->GetProperty()->SetDiffuseColor(
      colors->GetColor3d("Banana").GetData());
  booleanOperationActor->GetProperty()->SetSpecular(.6);
  booleanOperationActor->GetProperty()->SetSpecularPower(20);
 
  //std::cout << float(clock() - begin_time) / CLOCKS_PER_SEC;
 
  vtkNew<vtkRenderer> renderer;
  renderer->AddViewProp(input1Actor);
  renderer->AddViewProp(input2Actor);
  renderer->AddViewProp(booleanOperationActor);
  renderer->SetBackground(colors->GetColor3d("Silver").GetData());
  vtkNew<vtkRenderWindow> renderWindow;
  renderWindow->AddRenderer(renderer);
  renderWindow->SetSize(640, 480);
  renderWindow->SetWindowName("BooleanOperationPolyDataFilter");
 
  double viewUp[3] = {0.0, 0.0, 1.0};
  double position[3] = {0.0, -1.0, 0.0};
  PositionCamera(renderer, viewUp, position);
  renderer->GetActiveCamera()->Dolly(1.4);
  renderer->ResetCameraClippingRange();
 
  vtkNew<vtkRenderWindowInteractor> renWinInteractor;
  renWinInteractor->SetRenderWindow(renderWindow);
 
  renderWindow->Render();
  renWinInteractor->Start();
 
  return EXIT_SUCCESS;
}
namespace {
vtkSmartPointer<vtkPolyData> ReadPolyData(const char* fileName)
{
  vtkSmartPointer<vtkPolyData> polyData;
  std::string extension =
      vtksys::SystemTools::GetFilenameExtension(std::string(fileName));
  if (extension == ".ply")
  {
    vtkNew<vtkPLYReader> reader;
    reader->SetFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else if (extension == ".vtp")
  {
    vtkNew<vtkXMLPolyDataReader> reader;
    reader->SetFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else if (extension == ".obj")
  {
    vtkNew<vtkOBJReader> reader;
    reader->SetFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else if (extension == ".stl")
  {
    vtkNew<vtkSTLReader> reader;
    reader->SetFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else if (extension == ".vtk")
  {
    vtkNew<vtkPolyDataReader> reader;
    reader->SetFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else if (extension == ".g")
  {
    vtkNew<vtkBYUReader> reader;
    reader->SetGeometryFileName(fileName);
    reader->Update();
    polyData = reader->GetOutput();
  }
  else
  {
    vtkNew<vtkSphereSource> source;
    source->Update();
    polyData = source->GetOutput();
  }
  return polyData;
}
 
void PositionCamera(vtkRenderer* renderer, double* viewUp, double* position)
{
  renderer->GetActiveCamera()->SetFocalPoint(0.0, 0.0, 0.0);
  renderer->GetActiveCamera()->SetViewUp(viewUp);
  renderer->GetActiveCamera()->SetPosition(position);
  renderer->ResetCamera();
  return;
}
} // namespace

The issue is I'm getting this error when performing the boolean operation between the sphere 1,2 results and sphere 3:

vtkDistancePolyDataFilt:82     ERR| vtkDistancePolyDataFilter (000001DE8F936820): No points/cells to operate on

I'm new to VTK so any hints/tips on what I could test/debug would be great !

Update:

Marco Musy's answer is correct (calling Update()). For reference, in my scenario this would be:

  booleanOperation->SetInputData(0, input1);
  booleanOperation->SetInputData(1, input2);
  // update 1st boolean operation
  booleanOperation->Update();

  vtkNew<vtkBooleanOperationPolyDataFilter> booleanOperation2;
  booleanOperation2->SetOperationToUnion();
  booleanOperation2->SetInputData(0, booleanOperation->GetOutput());
  booleanOperation2->SetInputData(1, input3);

The mapped might have done the Update() behind the scenes when I had one operation and I completely missed the fact that I need to manually update each boolean operation if I need to plug it's results into another.

His library vedo is really cool by the way: really easy and fun to use!

Solution

  • You might just be missing a call to Update() in the second step. The vtkBooleanOperationPolyDataFilter seems to be working fine, eg. (python):

    from vedo import Sphere
s1 = Sphere([0.25, 0, 0.0])
s2 = Sphere([0.00, 0, 0.0])
s3 = Sphere([0.25, 0, 0.5])
s1.boolean("+", s2).boolean("+", s3).show()

    enter image description here