CGAL users discussion list

[Sebastien Loriot <>]

It's not a simple projection but a basis change, which means that 
distances are also changed. However to_3d restores that.

Best,

Sebastien.

On 5/19/23 20:14, Maria Eduarda Veras Martins ( via 
cgal-discuss Mailing List) wrote:
> plane_3::to_2d only removes the z coordinate. That way I could lose a 
> lot of information by removing this coordinate. In addition, I create 
> one polygon with segment_3, would that not be a problem in this 
> implementation?
>
> Thank you!
>
> Em sex., 19 de mai. de 2023 às 04:03, Sebastien Loriot 
> < <>> escreveu:
>
>     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
>     <https://sympa.inria.fr/sympa/info/cgal-discuss>
>      >
>
>     -- 
>     You are currently subscribed to cgal-discuss.
>     To unsubscribe or access the archives, go to
>     https://sympa.inria.fr/sympa/info/cgal-discuss
>     <https://sympa.inria.fr/sympa/info/cgal-discuss>
>
>
>
> --
> You are currently subscribed to cgal-discuss.
> To unsubscribe or access the archives, go to
> https://sympa.inria.fr/sympa/info/cgal-discuss