How to correctly use VTK ConstrainedDelaunay2D?

Question

I've started from the VTK ConstrainedDelaunay2D example and added my own points:

#include <vtkSmartPointer.h>
#include <vtkDelaunay2D.h>

#include <vtkCellArray.h>
#include <vtkProperty.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPolygon.h>
#include <vtkMath.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkNamedColors.h>
#include <vtkVersionMacros.h> // For version macros

int main(int, char *[])
{
  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
  
  int ptsHeight = 400;
  std::vector<std::vector<int>> pts{ {166, 127},{103, 220},{166, 190},{174, 291},{189, 226},{227, 282},{213, 187},{242, 105},{196, 131},{182, 83} };
  for (size_t i = 0; i < pts.size(); i++)
  {
      // !important: flip y
      int x = pts[i][0];
      int y = ptsHeight - pts[i][1];
      points->InsertNextPoint(x, y, 0);
  }

  vtkSmartPointer<vtkPolyData> aPolyData = vtkSmartPointer<vtkPolyData>::New();
  aPolyData->SetPoints(points);

  // Create a cell array to store the polygon in
  vtkSmartPointer<vtkCellArray> aCellArray = vtkSmartPointer<vtkCellArray>::New();

  // Define a polygonal hole with a clockwise polygon
  vtkSmartPointer<vtkPolygon> aPolygon = vtkSmartPointer<vtkPolygon>::New();

  for (unsigned int i = 0; i < pts.size(); i++)
  {
      aPolygon->GetPointIds()->InsertNextId(i); 
  }

  aCellArray->InsertNextCell(aPolygon);

  // Create a polydata to store the boundary. The points must be the
  // same as the points we will triangulate.
  vtkSmartPointer<vtkPolyData> boundary =
   vtkSmartPointer<vtkPolyData>::New();
  boundary->SetPoints(aPolyData->GetPoints());
  boundary->SetPolys(aCellArray);

  // Triangulate the grid points
  vtkSmartPointer<vtkDelaunay2D> delaunay =
   vtkSmartPointer<vtkDelaunay2D>::New();
  delaunay->SetInputData(aPolyData);
  delaunay->SetSourceData(boundary);

  // Visualize
  vtkSmartPointer<vtkPolyDataMapper> meshMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  meshMapper->SetInputConnection(delaunay->GetOutputPort());
  
  vtkSmartPointer<vtkNamedColors> colors =
    vtkSmartPointer<vtkNamedColors>::New();

  vtkSmartPointer<vtkActor> meshActor =
    vtkSmartPointer<vtkActor>::New();
  meshActor->SetMapper(meshMapper);
  meshActor->GetProperty()->EdgeVisibilityOn();
  meshActor->GetProperty()->SetEdgeColor(colors->GetColor3d("Peacock").GetData());
  meshActor->GetProperty()->SetInterpolationToFlat();
  meshActor->GetProperty()->SetBackfaceCulling(true);

  // Create a renderer, render window, and interactor
  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);

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

  // Add the actor to the scene
  renderer->AddActor(meshActor);
  //renderer->AddActor(boundaryActor);
  renderer->SetBackground(colors->GetColor3d("Mint").GetData());

  // Render and interact
  renderWindow->SetSize(640, 480);
  renderWindow->Render();
  renderWindowInteractor->Start();

  return EXIT_SUCCESS;
} 

I'm experiencing two issues:

1. I get different results if I flip the Y coordinates: why is that ?
2. Why are there faces pointing in the wrong direction (flipped normal / wrong winding )?

Here's what I mean by the 1st issue:

If I don't flip the Y coordinates I get this:

I get the same effect if I don't flip the Y axis but insert the boundary polygon in reverse order:

for (unsigned int i = 0; i < pts.size(); i++)
  {
      aPolygon->GetPointIds()->InsertNextId(pts.size() - 1 - i); 
  }

I don't think I fully understand how the boundary/constraint works. I thoght that the same points should produce the same triangulation wether the vertices are flipped vertically or not. (I suspect the order of indices changes then ?)

Regarding the second issue (unpredictable flipped faces) I'm not sure what the best way forward is. I had a look at the vtkDelaunay2D class and couldn't find anything related. (I've tried setting projection plane mode to VTK_DELAUNAY_XY_PLANE, but it didn't seem to affect the output)

I've also tried to use vtkPolyDataNormals but got no output:

vtkSmartPointer<vtkPolyDataNormals> normalGenerator = vtkSmartPointer<vtkPolyDataNormals>::New();
normalGenerator->SetInputData(delaunay->GetOutput());
normalGenerator->ComputePointNormalsOff();
normalGenerator->ComputeCellNormalsOn();
normalGenerator->FlipNormalsOn();
normalGenerator->Update();

(normalGenerator's output has 0 cells and points)

Is there a way to compute constrained delaunay triangulation for a list of 2d points and ensure all the faces point the same way ? (If so, how ? Would it be possible to do this with the vtkDelaunay2D class alone or is it necessary to use other filters?)

Any hints/tips are more than welcome :)

I'm using VTK 8.2 by the way.

Answer 1

the flipping in y effectively reverses the faces orientation (what is clockwise becomes anti-clockwise, like in a mirror). I'm not sure I can reproduce your example above. A quick test in python seems to give the expected behavior, maybe you can start from this and map it to your c++ version:

import vedo

pts = [
    [166, 127],
    [103, 220],
    [166, 190],
    [174, 291],
    [189, 226],
    [227, 282],
    [213, 187],
    [242, 105],
    [196, 131],
    [182, 83],
]

ids = [[2,4,6], [0,2,8]]  # faces to erase by pt-index (clockwise)
dly = vedo.delaunay2D(pts, mode='xy', boundaries=ids)
dly.c('grey5').lc('red4').lw(2)

labels = vedo.Points(pts).labels('id').z(1)

vedo.show(labels, dly, axes=1)