CGAL users discussion list

[Maria Eduarda Veras Martins <>]

Hi, I have to find the polygons that were formed when I intersected planes with a tetrahedral mesh. I'm already getting a list with the segments resulting from the intersections but I don't know how to get the polygons from that.

here is my code so far:

#include <list>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Implicit_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/iterator.h>
#include <iostream>
#include <fstream>
#include <map>
#include <vector>
#include <string>
#include <CGAL/utility.h>
#include <omp.h>

// Domain
struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT (Function)(const Point&);
typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;

typedef C3t3::Triangulation Tr;
typedef C3t3::Cell_handle Cell_handle;
typedef Tr::Facet Facet;
typedef Tr::Finite_facets_iterator Finite_facets_iterator;
typedef Tr::Finite_vertices_iterator Finite_vertices_iterator;
typedef Tr::Finite_cells_iterator Finite_cells_iterator;
typedef Tr::Vertex_handle Vertex_handle;

typedef K::Plane_3 Plane;
typedef K::Line_3 Line;
typedef K::Direction_3 Direction;
typedef K::Triangle_3 Triangle;
typedef K::Segment_3 Segment;
typedef std::list<Triangle>::iterator Iterator;
typedef CGAL::AABB_triangle_primitive<K,Iterator> Primitive;
typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
typedef Tree::Object_and_primitive_id Object_and_primitive_id;
typedef Tree::Primitive_id Primitive_id;

// Function
FT sphere_function (const Point& p)
{
const FT x2=p.x()*p.x();
const FT y2=p.y()*p.y();
const FT z2=p.z()*p.z();

return x2+y2+z2-1;
}

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());
    }
};


int main() {
    // Domain
    Mesh_domain domain(sphere_function, K::Sphere_3(CGAL::ORIGIN, 1.));

    // Criteria
    Facet_criteria facet_criteria(30, 0.1, 0.01); // angle, size, approximation
    Cell_criteria cell_criteria(3, 0.1); // radius-edge ratio, size
    Mesh_criteria criteria(facet_criteria, cell_criteria);

    // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);

    // Output
    std::ofstream medit_file("outmatrix.mesh");
    c3t3.output_to_medit(medit_file);

    //std::cout<< "Number of verticesDT"<<" "<< c3t3.triangulation().number_of_vertices() << std::endl;
    //std::cout<< "Number of facetsDT"<<" "<< c3t3.triangulation().number_of_finite_facets() << std::endl;
    //std::cout<< "Number of facetsC3T3"<<" "<<c3t3.number_of_facets() << std::endl;
    //std::cout<< "Number of cellsC3T3"<<" "<<c3t3.number_of_cells() << std::endl;
    //std::cout<< "Number of cellsDT"<<" "<<c3t3.triangulation().number_of_finite_cells() << std::endl;


    // for each facet build a triangle
    std::list<Triangle> triangles;
    for(Finite_facets_iterator fi = c3t3.triangulation().finite_facets_begin(); fi != c3t3.triangulation().finite_facets_end(); fi++) {

        Cell_handle cell = fi->first;
        int opposite_vertex_index = fi->second;
        if(c3t3.is_in_complex(cell)) {
            Triangle t=c3t3.triangulation().triangle(cell,opposite_vertex_index);
            triangles.push_back(t);
        }
        else if (c3t3.is_in_complex(*fi)) {
            Triangle t=c3t3.triangulation().triangle(*fi);
            triangles.push_back(t);
        }

    }

    //build aabb tree
    Tree tree(triangles.begin(),triangles.end());

    //// set plane and find intersections
    //Point p(0, 0, 0);
    //Plane plane(p, K::Vector_3(1, 0, 0));
//
    //std::list <Object_and_primitive_id> intersections;
    //tree.all_intersections(plane, std::back_inserter(intersections));

    std::vector <Plane> planes;
    for (int i = 0; i < 10; i++) {
        double angle = i * 36.0 * (M_PI / 180);
        K::Vector_3 normal(cos(angle), sin(angle), 0);
        planes.emplace_back(Point (0, 0, 0), normal);
    }

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

    #pragma omp parallel for
    for (size_t i = 0; i < planes.size(); i++) {
        std::list <Object_and_primitive_id> intersections;
        tree.all_intersections(planes[i], std::back_inserter(intersections));
        #pragma omp critical
        planes_intersections[planes[i]] = intersections;
    }

    for (const auto& aux : planes_intersections) {
        std::cout << "Plane: " << aux.first << " Interseções: " << aux.second.size() << std::endl;
    }

    return 0;
}

Thank you!