CGAL users discussion list

[Maria Eduarda Veras Martins <>]

Hi!

I need to intersect several planes with a volumetric mesh. I'm making an aabb tree with the faces of the tetrahedron but I don't know how to efficiently intersect the planes.

Here is my code so far (only one plane):

#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>

// 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;

}

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

}

}

Tree tree(triangles.begin(),triangles.end());

// set plane and count 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));

//how to do it with many planes?

//I'm thinking of looping the plans, but I wanted something more efficient

return 0;

}