Constrained Delaunay2D doesn't cover the whole area

[MoonMarcher]

Hello VTK users,

I’ve been trying for quite some time to triangulate an area of points with a few polygons which must be hollowed up in the area. I’m using ConstrainedDelaunay2D as the basis of my implementation but the best I could reach so far is:

However what I’m looking for is something like this:

So far, I’ve tried different winding of points in polygons (CW/CCW), and changing “Alpha”, “Offset” and “Tolerance” (At this point I was shooting for anything!) but didn’t really help.

Any help would be highly appreciated.

[marcomusy]

if you start with polygons maybe you can shrink them with vtkShrinkFilter to make the delauney filter work better (?)
If you start from an image you can maybe threshold it into a mesh:

from vedo import *
pic = Picture('map.jpg')
msh = pic.tomesh()
arr1 = msh.pointdata['RGBA']
arr2 = np.sum((arr1-[50,50,50])**2, axis=1)
msh.threshold(arr2, above=40).triangulate()
msh = msh.extractLargestRegion().c('green3')
show(pic, msh, N=2, axes=1)

[MoonMarcher]

Thanks @marcomusy for your suggestion and your nice representation. However, the screenshot is taken from a Matlab visualization and I’m currently using C++ interface of VTK to implement a general visualization tool for our future simulations.

I’m using coordinates of points to create polygons and add field data in the next step to run delaunay on the area. Here’s what I’ve got so far:

    vtkNew<vtkPoints> points;
    points->SetNumberOfPoints(map.size());

    for(const auto &r : map)
    {
        points->SetPoint(r.first, 0.0, 0.0, 0.0);
    }

    vtkNew<vtkCellArray> aCellArray;
    int idOffset = points->GetNumberOfPoints();
    size_t nPoints = 0;

    for (size_t i = 0; i < blocks.size(); ++i)
    {
        vtkNew<vtkPolygon> aPolygon;
        idOffset += nPoints;
        nPoints = blocks[i]->outline->points.size();
        for (int i_p = 0; i_p < nPoints; ++i_p)
        {
            points->InsertNextPoint(blocks[i]->outline->points[nPoints-i_p-1]->x, blocks[i]->outline->points[nPoints-i_p-1]->y, 0.0);
            aPolygon->GetPointIds()->InsertNextId(idOffset+i_p);
        }
        aCellArray->InsertNextCell(aPolygon);
    }

    for(unsigned int t=0; t<N_t; t++)
    {
        for(const auto &r : map)
        {
            points->SetPoint(r.first, grid->x, grid->y, grid->F);
        }

        double bounds[6];
        points->GetBounds(bounds);

        vtkNew<vtkPolyData> aPolyData;
        aPolyData->SetPoints(points);

        vtkNew<vtkPolyData> boundary;
        boundary->SetPoints(aPolyData->GetPoints());
        boundary->SetPolys(aCellArray);

        vtkNew<vtkDelaunay2D> delaunay;
        delaunay->SetInputData(aPolyData);
        delaunay->SetSourceData(boundary);
        delaunay->Update();

        vtkNew<vtkDoubleArray> P;
        P->SetNumberOfComponents(1);
        P->SetName("P");

        for (vtkIdType i = 0; i < delaunay->GetOutput()->GetNumberOfPoints(); i++)
        {
            double val = delaunay->GetOutput()->GetPoint(i)[2];
            P->InsertNextValue(val);
        }

        delaunay->GetOutput()->GetPointData()->AddArray(P);

        vtkNew<vtkXMLPolyDataWriter> writer;
        writer->SetFileName((outputFileName + std::to_string(t) + ".vtp").c_str());
        writer->SetInputData(delaunay->GetOutput());
        writer->Write();
    }

First I fill the points container with my grid points- the ones which are going to have field data (initialized by zero) and add vertices of blocks (which I wish to cut them out from the final map) and then input field data at each time step in the loop. Finally I run delaunay and output the data.

It also throws this warning while writing data:

vtkDelaunay2D.cxx:1409  WARN| vtkDelaunay2D (0x39ac870): Edge not recovered, polygon fill suspect

Any suggestion?