CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

 wrote:
> Hi,
> I'm using the function unique for a list of intersections obtained by
> means of an AABB tree in order to avoid "multiple intersections". More
> specifically I'm searching for intersections with 3D triangles from a 3D
> mesh and my query object is a line with fixed direction and crossing a
> vertex of the mesh. In this case I obtain multiple intersections at the
> vertex site which are only partially removed by the unique function. As
> example, in the following code I use the mesh in
> examples/Mesh_3/mesh_implicit_sphere.cpp in order to create a list of all
> triangles belonging the 3D complex, (I want select both triangles on the
> surface and the volume), as set of primitives and then perform
> intersection queries "moving" the fixed direction on each vertex.
> I have some questions about the results:
> 1)Why unique function applied to the list doesn't remove all multiple
> intersections? It seems not linked with the choice of the geometrical
> kernel.
If you increase the precision of the output you will see that
the points are different.

> 2) I obtain on average only 30 intersection points for about 30000
> triangles. How to establish whether the number of intersections is
> "reasonable"?
What do you mean by "reasonable"?


> 3) I'm wrong expecting roughly uniform spacing for the intersection points
> when tuning radius bound of cells and boundary triangles to obtain a
> uniform mesh?
I think you are wrong. The fact that the mesh is uniform does not
imply that any line intersecting it has intersection points with 
tetrahedron uniformly spaced (the mesh is not made of axis aligned cube 
with the line being collinear to an axis).


S.
> Thanks in advance for any suggestion!
>
>
>
> #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::Ray_3 Ray;
> 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("out.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());
>         Point a(0.0, 0.0, 0.0);
>     Point b(0.7, 0.34, 0.61);
>
>         // counts #intersections
>     Ray ray_query(a,b);
>         Direction dir(ray_query);
>
>         std::map<Vertex_handle, int> V;
>         int inum = 1;
>         for( Finite_vertices_iterator vit =
> c3t3.triangulation().finite_vertices_begin();
>                 vit != c3t3.triangulation().finite_vertices_end();
>                 ++vit)
>         {
>                 V[vit] = inum++;
>                 Point c = vit->point();
>                 Line line_vertex(c,dir);
>
>     std::cout << tree.number_of_intersected_primitives(line_vertex)
>         << " intersections(s) with line query" << std::endl;
>
>         // computes all intersections with line query
>     std::list<Object_and_primitive_id> intersections;
>     tree.all_intersections(line_vertex,
> std::back_inserter(intersections));
>         std::list<Object_and_primitive_id>::iterator it_intersections=
> intersections.begin();
>     std::list<Point> points;
>         std::list<Point>::iterator it_points;
>
>         for(; it_intersections != intersections.end() ; ++it_intersections)
>
>     {
>         // get intersection object
>                 Object_and_primitive_id op = *it_intersections;
>                 CGAL::Object object = op.first;
>         Point point;
>                 Segment segment;
>         if(CGAL::assign(point,object))
>                 {
>                         
> points.push_back(Point(point.x(),point.y(),point.z()));
why not points.push_back(point);   ?
>                         //std::cout << "intersection object is a point" <<" "<< 
> point.x()<<"
> "<<point.y()<<" "<<point.z()<<std::endl;
>
>                 }
>                 else if(CGAL::assign(segment,object))
>                         std::cout << "intersection object is a segment" 
> <<std::endl;
>
>     }
>
>         //set order in intersection list
>         points.sort();
>                 int ip=1;
>                 for(it_points=points.begin(); it_points != points.end() ; 
> ++it_points) {
>                         std::cout <<  "myorderedlist contains:"<<" " << 
> ip++<<"
> "<<*it_points<<std:: endl;}
>                 std::cout <<"---" << std::endl;
>         //erase multiple intersections
>         points.unique();
>
>         int it=1;
>         for(it_points=points.begin(); it_points != points.end() ; 
> ++it_points) {
> std::cout <<  "myuniquelist contains:"<<" " << it++<<"
> "<<*it_points<<std:: endl;}
>
>         intersections.clear();
>         points.clear();
>
>         }
>
>         return 0;
> }