CGAL users discussion list

[Maria Eduarda Veras Martins <>]

Hi again, 

I need to find the intersection area of two polygons, the biggest problem is that they are in 3d. I tried to use projection traits but when I use some intersection function like do_intersect or CGAL::intersection, several errors occur.

Here is my code so far, im already able to create the polygons!

#include <iostream>

#include <list>

#include <map>

#include <vector>

#include <omp.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Boolean_set_operations_2.h>

#include <CGAL/Projection_traits_3.h>

#include <CGAL/Polygon_2.h>

#include <CGAL/Polygon_with_holes_2.h>

#include <CGAL/Polygon_2_algorithms.h>

#include <CGAL/AABB_tree.h>

#include <CGAL/AABB_traits.h>

#include <CGAL/Polyhedron_3.h>

#include <CGAL/AABB_face_graph_triangle_primitive.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Projection_traits_3 <K> GT;

typedef CGAL::Polygon_2 <GT> Polygon;

typedef CGAL::Polygon_with_holes_2 <GT> Polygon_with_holes;

typedef K::Point_3 Point;

typedef K::Plane_3 Plane;

typedef K::Vector_3 Vector;

typedef K::Segment_3 Segment;

typedef CGAL::Polyhedron_3<K> Polyhedron;

typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive_poly;

typedef CGAL::AABB_traits<K, Primitive_poly> Traits_poly;

typedef CGAL::AABB_tree<Traits_poly> Tree_poly;

typedef Tree_poly::Object_and_primitive_id Object_and_primitive_id_poly;

struct Plane_compare {

bool operator()(const Plane& p1, const Plane& p2) const {

// Use a combinação dos coeficientes dos planos para compará-los

return std::make_tuple(p1.a(), p1.b(), p1.c(), p1.d()) < std::make_tuple(p2.a(), p2.b(), p2.c(), p2.d());

}

};

struct Segment_compare {

bool operator() (const Segment& seg1, const Segment& seg2) const {

return seg1.source() < seg2.source() || (!(seg2.source() < seg1.source()) && seg1.target() < seg2.target());

}

};

Polygon seg_into_polygon (std::list <Segment>& segments, const Vector& plane_normal) {

std::map <Segment, std::pair <Point, Point>, Segment_compare> segments_map;

for (const auto& seg : segments) {

segments_map[seg] = {seg.source(), seg.target()};

}

Segment current_seg = segments.front();

GT gt {plane_normal};

Polygon polygon(gt);

polygon.push_back(current_seg.source());

polygon.push_back(current_seg.target());

segments_map.erase(current_seg);

while (!segments_map.empty()) {

if (segments_map.size() == 1) {

break;

}

bool find_flag = false;

Segment temp;

for (auto& it : segments_map) {

if (it.second.first == current_seg.target()) {

polygon.push_back(it.second.second);

temp = it.first;

current_seg = it.first;

find_flag = true;

break;

}

}

if (find_flag == true) {

segments_map.erase(temp);

}

else {

for (auto& it : segments_map) {

if (it.second.second == current_seg.target()) {

polygon.push_back(it.second.first);

temp = it.first;

current_seg = it.first;

find_flag = true;

break;

}

}

if (find_flag == true) {

segments_map.erase(temp);

}

}

}

//std::cout << "Poligono criado: " << polygon << std::endl;

return polygon;

}

int main()

{

Point p(1.0, 0.0, 0.0);

Point q(0.0, 1.0, 0.0);

Point r(0.0, 0.0, 1.0);

Point s(0.0, 0.0, 0.0);

Polyhedron polyhedron;

polyhedron.make_tetrahedron(p, q, r, s);

std::ofstream out ("tetrahedron.off");

out << polyhedron;

out.close();

// constructs AABB tree

Tree_poly tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);

//planes

std::vector <Plane> planes;

Point pE (0.58, 0, 0);

Point pF (0.58, 0, 0.42);

Point pG (0, 0.62, 0.38);

Plane p1 (pF, pG, pE);

planes.push_back(p1);

Point pH (0.26, 0, 0);

Point pI (0, 0, 0.5);

Point pJ (0, 0.32, 0);

Plane p2 (pH, pI, pJ);

planes.push_back(p2);

Plane p3 (pI, pG, pF);

planes.push_back(p3);

std::map <Plane, std::list <Segment>, Plane_compare> planes_intersections;

std::list <Polygon> polygon_teste;

#pragma omp parallel for

for (size_t i = 0; i < planes.size(); i++) {

std::list <Object_and_primitive_id_poly> intersections;

tree.all_intersections(planes[i], std::back_inserter(intersections));

std::list <Object_and_primitive_id_poly>::iterator it_intersections = intersections.begin();

std::list <Segment> segments;

for (; it_intersections != intersections.end(); ++it_intersections) {

Object_and_primitive_id_poly op = *it_intersections;

CGAL::Object obj = op.first;

Segment s;

if (CGAL::assign(s, obj)) {

segments.push_back(s);

}

}

#pragma omp critical

planes_intersections[planes[i]] = segments;

if (i == 2) {

polygon_teste.push_back(seg_into_polygon(segments, planes[i].orthogonal_vector()));

}

else

Polygon temp = seg_into_polygon(segments, planes[i].orthogonal_vector());

}

//prints for debugging

for (const auto& aux : planes_intersections) {

std::cout << "Plano:" << aux.first << std::endl;

std::cout << "Segmentos: " << std::endl;

std::list <Segment> list_seg = aux.second;

for (const auto& item : list_seg) {

std::cout << item << std::endl;

}

std::cout << std::endl;

}

//intersection with coplanar polygon

//testing with plane (-0.05x - 0.07y - 0.36z = -0.18)

//polygon c = (0, 0, 0.5) (0, 0.62, 0.38) (0.58, 0, 0.42)

//frature cell = (1.07, -0.66, 0.48) (0.27, 0.14, 0.44) (0.75, 0.39, 0.32)

std::list <Segment> temp_cel;

Segment k (Point(1.07, -0.66, 0.48), Point(0.27, 0.14, 0.44));

temp_cel.push_back(k);

Segment l (Point(0.75, 0.39, 0.32), Point(1.07, -0.66, 0.48));

temp_cel.push_back(l);

Segment m (Point(0.27, 0.14, 0.44), Point(0.75, 0.39, 0.32));

temp_cel.push_back(m);

//green polygon

Polygon frature_cel = seg_into_polygon(temp_cel, p3.orthogonal_vector());

//red polygon

Polygon polygon_c = polygon_teste.front();

std::cout << "frature_cel: " << frature_cel << std::endl;

std::cout << "Area = " << frature_cel.area() << std::endl;

std::cout << "polygon c: " << polygon_c << std::endl;

std::cout << "Area = " << polygon_c.area() << std::endl;

return 0;

}

to make it clearer, the purple area is the one I need to find.