[vtkusers] Multiple issues with IntersectionPolyDataFilter

[Cory Quammen]

Hi Alex,

Yeah, that filter, and the boolean operations filter, fails on some
configurations of input geometry, and it can be difficult to track
down what is going on. That's why they are not more robust. I
apologize for that, but it is a tricky operation to get right.

At the moment, there are not resources to put into making the filter
more robust, but we welcome community contributions to that end.

Thanks,
Cory

On Fri, Jan 13, 2017 at 5:46 PM, Alex Goins <alexgoins301 at gmail.com> wrote:
> 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;
> }
>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Search the list archives at: http://markmail.org/search/?q=vtkusers
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/vtkusers
>



-- 
Cory Quammen
Staff R&D Engineer
Kitware, Inc.