I'm trying to find the result of an intersection and difference between two triangles in CGAL.
My current approach is to convert the triangles to polygons, compute the intersection/difference and finally triangulate the result.
The following code gives me 'Variable used before being initialized (or CGAL bug)' in the intersection test. Any ideas what could be wrong?
#include <iostream>
#include <CGAL/basic.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/triangulate_polyhedron.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <vector>
#include <CGAL/Boolean_set_operations_2.h>
#include <list>
struct FaceInfo2
{
FaceInfo2(){}
int nesting_level;
bool in_domain(){
return nesting_level%2 == 1;
}
};
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
typedef CGAL::Polygon_2<Kernel> Polygon_2;
typedef CGAL::Polygon_with_holes_2<Kernel> Polygon_with_holes_2;
typedef std::list<Polygon_with_holes_2> Pwh_list_2;
typedef CGAL::Triangulation_vertex_base_2<Kernel> Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,Kernel> Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<Kernel,Fbb> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
typedef CGAL::Exact_predicates_tag Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel, TDS, Itag> CDT;
/*typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
typedef Kernel::Segment_2 Segment_2;
typedef CGAL::Polygon_2<Kernel> Polygon_2;
*/
typedef Kernel::Triangle_2 Triangle;
// from http://doc.cgal.org/latest/Triangulation_2/index.html#Subsection_37.8.2
void mark_domains(CDT& ct,
CDT::Face_handle start,
int index,
std::list<CDT::Edge>& border )
{
if(start->info().nesting_level != -1){
return;
}
std::list<CDT::Face_handle> queue;
queue.push_back(start);
while(! queue.empty()){
CDT::Face_handle fh = queue.front();
queue.pop_front();
if(fh->info().nesting_level == -1){
fh->info().nesting_level = index;
for(int i = 0; i < 3; i++){
CDT::Edge e(fh,i);
CDT::Face_handle n = fh->neighbor(i);
if(n->info().nesting_level == -1){
if(ct.is_constrained(e)) border.push_back(e);
else queue.push_back(n);
}
}
}
}
}
// from http://doc.cgal.org/latest/Triangulation_2/index.html#Subsection_37.8.2
//explore set of facets connected with non constrained edges,
//and attribute to each such set a nesting level.
//We start from facets incident to the infinite vertex, with a nesting
//level of 0. Then we recursively consider the non-explored facets incident
//to constrained edges bounding the former set and increase the nesting level by 1.
//Facets in the domain are those with an odd nesting level.
void
mark_domains(CDT& cdt)
{
for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){
it->info().nesting_level = -1;
}
std::list<CDT::Edge> border;
mark_domains(cdt, cdt.infinite_face(), 0, border);
while(! border.empty()){
CDT::Edge e = border.front();
border.pop_front();
CDT::Face_handle n = e.first->neighbor(e.second);
if(n->info().nesting_level == -1){
mark_domains(cdt, n, e.first->info().nesting_level+1, border);
}
}
}
std::vector<Triangle> triangulate(Pwh_list_2 & polyList ){
std::vector<Triangle> res;
for (Polygon_with_holes_2 polygon1 : polyList) {
CDT cdt;
auto outer = polygon1.outer_boundary();
cdt.insert_constraint(outer.vertices_begin(), outer.vertices_end(), true);
for (auto hole = polygon1.holes_begin(); hole != polygon1.holes_end(); hole++){
cdt.insert_constraint(hole->vertices_begin(), hole->vertices_end(), true);
}
for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin();
fit!=cdt.finite_faces_end();++fit)
{
std::cout <<"New triangle "<<std::endl;
if ( fit->info().in_domain() ) {
Triangle tri;
for (int i=0;i<3;i++){
auto point = fit->vertex(i)->point();
tri.vertex(i) =fit->vertex(i)->point();
std::cout <<point.x()<<" ";
}
res.push_back(tri);
std::cout <<std::endl;
}
}
}
return res;
}
void ProjectOtherTriangle(Triangle thisTri, Triangle tri,
std::vector<Triangle>& differenceTriangle, std::vector<Triangle>& intersectionTriangles){
Polygon_2 thisPoly;
for (int i=0;i<3;i++){
thisPoly.push_back(thisTri[i]);
}
Polygon_2 poly;
for (int i=0;i<3;i++){
poly.push_back(tri[i]);
}
// Compute the intersection of P and Q.
Pwh_list_2 intR;
CGAL::intersection (thisPoly, poly, std::back_inserter(intR));
// add into newTriangles
auto tris = triangulate(intR);
intersectionTriangles.insert(intersectionTriangles.end(), tris.begin(), tris.end());
// Find difference
CGAL::difference (thisPoly, poly, std::back_inserter(intR));
// add into currentTriangle
tris = triangulate(intR);
differenceTriangle.insert(differenceTriangle.end(), tris.begin(), tris.end());
}
void printTriangle(Triangle t){
using namespace std;
for (int i=0;i<3;i++){
cout << t[i] <<" ";
}
cout << endl;
}
int main ()
{
// Construct the two input polygons.
Triangle P;
P[0] = Point_2 (0, 0);
P[1] = Point_2 (2, 0);
P[2] = Point_2 (1, 1.5);
printTriangle(P);
Triangle Q;
Q[0] = Point_2 (0, 2);
Q[1] = Point_2 (1, 0.5);
Q[2] = Point_2 (2, 2);
printTriangle(Q);
std::vector<Triangle> currentTriangle;
std::vector<Triangle> newTriangle;
ProjectOtherTriangle(P,Q, currentTriangle, newTriangle);
using namespace std;
cout << "currentTriangle"<<endl;
for (auto t : currentTriangle){
printTriangle(t);
}
cout << "newTriangle"<<endl;
for (auto t : newTriangle){
printTriangle(t);
}
return 0;
}
You cannot modify a Triangle_2 object. The following is not correct:
Triangle P;
P[0] = Point_2 (0, 0);
P[1] = Point_2 (2, 0);
P[2] = Point_2 (1, 1.5);
printTriangle(P);
Triangle Q;
Q[0] = Point_2 (0, 2);
Q[1] = Point_2 (1, 0.5);
Q[2] = Point_2 (2, 2);
printTriangle(Q);
You should use the constructor from 3 points instead.