CGAL users discussion list

[Sebastien Loriot <>]

Intersection of polygons is using the Arrangement package.
The projection traits is not a model of the traits concept required
by the arrangement.
So first you need a kernel with exact constructions for computing the
intersection (something like 
CGAL::Exact_predicates_exact_constructions_kernel).

What I would do:
- Pick a common plane
- use Plane_3::to_2d to produce 2d polygons
- compute the intersection
- use Plane_3::to_3d to turn the 2d polygon in a 3d one.
- compute the area of the 3d polygon (using Gauss formula with a point 
in the plane for exemple).

HTH,

Sebastien.



On 5/19/23 00:28, Maria Eduarda Veras Martins ( via 
cgal-discuss Mailing List) wrote:
> 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;
> #pragmaompparallelfor
> 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);
> }
> }
>
> #pragmaompcritical
> 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.
> image.png
>
> --
> You are currently subscribed to cgal-discuss.
> To unsubscribe or access the archives, go to
> https://sympa.inria.fr/sympa/info/cgal-discuss