# Constrained Delaunay2D doesn't cover the whole area

**URL:** https://discourse.vtk.org/t/constrained-delaunay2d-doesnt-cover-the-whole-area/8132
**Category:** Support
**Created:** [March 24, 2022, 11:21am UTC](https://discourse.vtk.org/t/constrained-delaunay2d-doesnt-cover-the-whole-area/8132 "2022-03-24T11:21:33Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![MoonMarcher](https://discourse.vtk.org/user_avatar/discourse.vtk.org/moonmarcher/32/4879_2.png) [@MoonMarcher](https://discourse.vtk.org/u/MoonMarcher)
#### Post date: [March 24, 2022, 11:21am UTC](https://discourse.vtk.org/t/constrained-delaunay2d-doesnt-cover-the-whole-area/8132/1 "2022-03-24T11:21:33Z")

</div>

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](https://kitware.github.io/vtk-examples/site/Cxx/Filtering/ConstrainedDelaunay2D/) as the basis of my implementation but the best I could reach so far is:

 ![Screenshot from 2022-03-24 12-14-50](https://discourse.vtk.org/uploads/default/original/2X/9/9ff0afad7d69c8aabca3e9b7f7b84d046478a250.png)

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

 ![Screenshot from 2022-03-24 12-01-15](https://discourse.vtk.org/uploads/default/original/2X/2/226209feed208b0465f30d16841dafa01479b370.jpeg)

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.

---

<div class="post-metadata">

### Author: ![marcomusy](https://discourse.vtk.org/user_avatar/discourse.vtk.org/marcomusy/32/95_2.png) [@marcomusy](https://discourse.vtk.org/u/marcomusy)
#### Post date: [March 25, 2022, 12:13pm UTC](https://discourse.vtk.org/t/constrained-delaunay2d-doesnt-cover-the-whole-area/8132/2 "2022-03-25T12:13:25Z")

</div>

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:

```python
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)

```

 ![Screenshot from 2022-03-25 13-09-26](https://discourse.vtk.org/uploads/default/original/2X/7/779db474c5733010fc6e40decd42e09e54047b62.png)

---

<div class="post-metadata">

### Author: ![MoonMarcher](https://discourse.vtk.org/user_avatar/discourse.vtk.org/moonmarcher/32/4879_2.png) [@MoonMarcher](https://discourse.vtk.org/u/MoonMarcher)
#### Post date: [March 25, 2022, 1:03pm UTC](https://discourse.vtk.org/t/constrained-delaunay2d-doesnt-cover-the-whole-area/8132/3 "2022-03-25T13:03:19Z")

</div>

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:

```auto
    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:

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

```

Any suggestion?