[vtkusers] Multiple issues with IntersectionPolyDataFilter

[Alex Goins]

I am trying to create multiple lines on a surface mesh by calculating the
intersection between the target mesh and a second surface/plane.  I am
having lots of issues with the IntersectionPolyDataFilter either not
finding the intersection or causing the program to crash.  If I move the
intersecting surface ever so slightly to one side, it seems to work fine.
Does anyone know what is causing it to crash or not return the
intersection?  In other cases, I am getting a lot of other warnings for "no
cell with correct orientation" and "edge not recovered, polygon fill
suspect" but I don't know if these are related to my previously mentioned
issues.  Here is some code to reproduce the issues.  Btw, I am running
version 7.1 from source (downloaded it about a month ago).  I have not been
able to create a minimum reproducible set of code to recreate the issue of
not finding intersections yet.  Any help is appreciated.

Thanks!
Alex

#include <vtkActor.h>
#include <vtkIntersectionPolyDataFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSmartPointer.h>

#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkMath.h>
#include <vtkSurfaceReconstructionFilter.h>
#include <vtkReverseSense.h>
#include <vtkContourFilter.h>
#include <vtkCellArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>


vtkSmartPointer<vtkPolyData>
createSurfaceFromSpline(vtkSmartPointer<vtkPolyData> line, double dist)
  {
    vtkSmartPointer<vtkPolyData> new_surface;
    vtkSmartPointer<vtkDataArray> normals =
line->GetPointData()->GetNormals();

    if(!normals)
    {
      cout << "no normals, cannot create surface\n";
      return new_surface;
    }

    new_surface = vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkCellArray> cells =
vtkSmartPointer<vtkCellArray>::New();
    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();

    // for each point, insert 2 points, one above and one below, to create
a new surface
    for(int i = 0; i < normals->GetNumberOfTuples(); ++i)
    {
      double* norm = normals->GetTuple(i);
      double* pt = line->GetPoints()->GetPoint(i);

      // insert the cell ids to create triangles
      if(i < normals->GetNumberOfTuples()-1)
      {
        vtkIdType start = i*2;

        cells->InsertNextCell(3);
        cells->InsertCellPoint(start);
        cells->InsertCellPoint(start+1);
        cells->InsertCellPoint(start+3);

        cells->InsertNextCell(3);
        cells->InsertCellPoint(start);
        cells->InsertCellPoint(start+3);
        cells->InsertCellPoint(start+2);
      }

      double new_pt[3];
      new_pt[0] = pt[0] + norm[0] * dist;
      new_pt[1] = pt[1] + norm[1] * dist;
      new_pt[2] = pt[2] + norm[2] * dist;

      points->InsertNextPoint( new_pt);

      new_pt[0] = pt[0] - norm[0] * dist;
      new_pt[1] = pt[1] - norm[1] * dist;
      new_pt[2] = pt[2] - norm[2] * dist;

      points->InsertNextPoint( new_pt);
    }

    new_surface->SetPolys(cells);
    new_surface->SetPoints(points);
    return new_surface;
  }

int main(int, char *[])
{
  // Custom surface
  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();

  unsigned int gridSize = 10;
  for(unsigned int x = 0; x < gridSize; x++)
    {
    for(unsigned int y = 0; y < gridSize; y++)
      {
        points->InsertNextPoint(x  , y , 1.0 * cos(double(x)/2.0) - 1.0 *
sin(double(y)/2.0) + vtkMath::Random(0.0, 0.001));
      }
    }

  // Mesh custom surface
  vtkSmartPointer<vtkPolyData> mesh = vtkSmartPointer<vtkPolyData>::New();
  mesh->SetPoints(points);

  // surface reconstruction
  vtkSmartPointer<vtkContourFilter> cf;
  vtkSmartPointer<vtkSurfaceReconstructionFilter> surf =
  vtkSmartPointer<vtkSurfaceReconstructionFilter>::New();

  surf->SetInputData(mesh);
  surf->SetSampleSpacing(0.5);
  surf->SetNeighborhoodSize(5);
  surf->Update();

  cf = vtkSmartPointer<vtkContourFilter>::New();
  cf->SetInputConnection(surf->GetOutputPort());
  cf->SetValue(0, 0.0);
  cf->Update();

  vtkSmartPointer<vtkReverseSense> reverse =
vtkSmartPointer<vtkReverseSense>::New();
  reverse->SetInputConnection(cf->GetOutputPort());
  reverse->ReverseCellsOn();
  reverse->ReverseNormalsOn();
  reverse->Update();

  // create intersecting plane
  vtkSmartPointer<vtkPolyData> line = vtkSmartPointer<vtkPolyData>::New();
  vtkSmartPointer<vtkDoubleArray> norms =
vtkSmartPointer<vtkDoubleArray>::New();
  vtkSmartPointer<vtkPoints> line_points =
vtkSmartPointer<vtkPoints>::New();

  norms->SetNumberOfComponents(3);

  //////////////////////////////
  // NOTE: Change the first point (x,y) values to test various issues:
  // (0, -1, 0) -> code crash: "free(): invalid pointer" error
  // (1, -1, 0) -> No cell with correct orientation
  // (-2, -1, 0) -> code crash segmentation fault
  // remove the -1.0 in the below 'for' loop (see below NOTE); (-1, -1, 0)
-> lots of vtkDelaunay2D "Edge not recovered, polygon fill suspect" errors
before it finally crashes

  line_points->InsertNextPoint(0, -1, 0);
  /////////////////////////////

  line_points->InsertNextPoint(10, 10, 0);
  double n[3] = {0, 0, 1};

  norms->InsertNextTuple(n);
  norms->InsertNextTuple(n);
  line->SetPoints(line_points);
  line->GetPointData()->SetNormals(norms);

  vtkSmartPointer<vtkPolyData> mesh2 = vtkSmartPointer<vtkPolyData>::New();
  mesh2 = createSurfaceFromSpline(line, 5.0);


  vtkSmartPointer<vtkIntersectionPolyDataFilter> intersectionPolyDataFilter
=
    vtkSmartPointer<vtkIntersectionPolyDataFilter>::New();
    intersectionPolyDataFilter->SetInputData( 0, reverse->GetOutput() );

  for(int i = 0; i < 100; ++i)
  {
    cout << i << "\n";

    ////////////////////////////
    // NOTE: remove the -1.0 from the X component to replicate the
vtkDelaunay2D error (see NOTE above)
    double pt[3] = {float(i)/10.0 - 1.0, float(i)/10.0, 0};
    ////////////////////////////

    line->GetPoints()->SetPoint(1, pt);
    mesh2 = createSurfaceFromSpline(line, 5.0);

    intersectionPolyDataFilter->SetInputData( 1, mesh2 );
    intersectionPolyDataFilter->Update();
  }


  // Add displays
  vtkSmartPointer<vtkPolyDataMapper> mapper1 =
vtkSmartPointer<vtkPolyDataMapper>::New();
  vtkSmartPointer<vtkPolyDataMapper> mapper2 =
vtkSmartPointer<vtkPolyDataMapper>::New();

  // display 1
  mapper1->SetInputConnection( reverse->GetOutputPort() );
  mapper1->ScalarVisibilityOff();
  vtkSmartPointer<vtkActor> actor1 = vtkSmartPointer<vtkActor>::New();
  actor1->SetMapper( mapper1 );
  actor1->GetProperty()->SetOpacity(0.9);
  actor1->GetProperty()->SetColor(1,0,0);

  // display 2
  mapper2->SetInputData( mesh2);
  mapper2->ScalarVisibilityOff();
  vtkSmartPointer<vtkActor> actor2 = vtkSmartPointer<vtkActor>::New();
  actor2->SetMapper( mapper2 );
  actor2->GetProperty()->SetOpacity(0.9);
  actor2->GetProperty()->SetColor(0,1,0);

  vtkSmartPointer<vtkPolyDataMapper> intersectionMapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
  intersectionMapper->SetInputConnection(
intersectionPolyDataFilter->GetOutputPort() );
  intersectionMapper->ScalarVisibilityOff();

  vtkSmartPointer<vtkActor> intersectionActor =
vtkSmartPointer<vtkActor>::New();
  intersectionActor->SetMapper( intersectionMapper );

  vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();
  renderer->AddViewProp(actor1);
  renderer->AddViewProp(actor2);
  renderer->AddViewProp(intersectionActor);

  vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer( renderer );

  vtkSmartPointer<vtkRenderWindowInteractor> renWinInteractor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  renWinInteractor->SetRenderWindow( renderWindow );

  renderWindow->Render();
  renWinInteractor->Start();

  return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20170113/13a5d12c/attachment.html>