CGAL users discussion list

[Bob Bill <>]

Let's say I have two polygons with non-empty intersection. I computed their intersection using CGAL::intersection free function. Then I wanted to mesh the obtained intersection. My approach has been the following: 

- compute the intersection

- loop over vertices and store coordinates

- create a CDT with those vertices, where the last one coincides with the first one. 

The problem is that the output says the resulting number of faces is 6, instead of 5. By printing the coords for each face, I noticed that one face is made by the vertices 

(0 1), (0.6 0), (0.257248 0.571254) which are the red dots in the attachment. Why is this so? And also, how can I avoid that this is considered as a face?

The code is the following:

#include <iostream>

#include <vector>

#include <CGAL/Boolean_set_operations_2.h>

#include <CGAL/Constrained_Delaunay_triangulation_2.h>

#include <CGAL/Delaunay_mesh_face_base_2.h>

#include <CGAL/Delaunay_mesh_size_criteria_2.h>

#include <CGAL/Delaunay_mesher_2.h>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/IO/write_VTU.h>

#include <CGAL/Polygon_2.h>

#include <CGAL/Polygon_with_holes_2.h>

#include <CGAL/Triangulation_2.h>

#include <CGAL/draw_polygon_2.h>

#include <CGAL/draw_polygon_with_holes_2.h>

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;

typedef CGAL::Polygon_2<Kernel> Polygon;

typedef CGAL::Polygon_with_holes_2<Kernel> Polygon_with_holes_2;

typedef Polygon::Point_2 Point;

typedef Polygon::Segment_2 Segment;

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Triangulation_vertex_base_2<K> Vb;

typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;

typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;

typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;

typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;

typedef CDT::Vertex_handle Vertex_handle;

int main() {

// define first polygon

const std::vector<Point> pts{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}};

const Polygon p(pts.cbegin(), pts.cend());

// define second polygon

const std::vector<Point> pts2{{0, 1}, {0.6, 0}, {2, 1}, {0, 2}};

const Polygon p2(pts2.cbegin(), pts2.cend());

std::vector<Polygon_with_holes_2> poly_list;

auto oi = CGAL::intersection(p, p2, std::back_inserter(poly_list));

Polygon intersect = poly_list[0].outer_boundary();

// Create the triangulation

CDT cdt;

Vertex_handle first_vertex =

cdt.insert(CDT::Point(CGAL::to_double(intersect.vertices_begin()->x()),

CGAL::to_double(intersect.vertices_begin()->y())));

std::vector<Vertex_handle> vec_vertices;

for (auto it = intersect.vertices_begin(); it != intersect.vertices_end();

++it) {

// Now I have coords for each vertex, I can start building the triangulation

Vertex_handle va = cdt.insert(

CDT::Point(CGAL::to_double(it->x()), CGAL::to_double(it->y())));

vec_vertices.push_back(va);

}

vec_vertices.push_back(first_vertex); // so to close the polygon

for (unsigned int i = 0; i < vec_vertices.size() - 1; ++i) {

cdt.insert_constraint(vec_vertices[i], vec_vertices[i + 1]);

}

CGAL::refine_Delaunay_mesh_2(cdt, Criteria(0.125, 0.8));

std::cout << "Number of vertices: " << cdt.number_of_vertices() << '\n';

std::cout << "Number of faces: " << cdt.number_of_faces() << '\n';

std::ofstream output("intersected_mesh.vtu");

CGAL::write_vtu(output, cdt);

return 0;

}

Attachment: mesh_intersection.png
Description: PNG image