CGAL users discussion list

Dear Stephane,
some days ago you suggested that there are two kind of primitives  in
order to get back triangulation facets when combining a 3D mesh generation
procedure with an AABB tree.
I remind you that I had to compute the intersection points of a straight
line and a 3D object represented as a 3D complex (for example a
tetrahedron mesh in 3D mesh generation).  My query object is a line with
fixed direction and crossing a vertex of the mesh (I consider parallel
straight lines crossing each vertex for different intersection queries).
The second kind of primitives you mentioned such as Datum is Triangle_3
and Id is Tr::Facet_finite_iterator is especially useful for my purpose.
Actually I need to extract the information on both vertex coordinates of
the primitives and their indexes in the C3T3 data structure (I mean I need
to keep the map used in the mesh procedure).
As you can see in the attached code I get back the primitives from
Object_and_primitive_id as facet iterator.
The second constraint I have is to follow the ray from the "entrance" to
the "exit" point in the
object, that is set a new order, with respect to the tree, where the first
and the last intersections are on "opposite" sides of the surface of the
sphere (the attached code is for the simple test where straight lines are
parallel to the y axis, so entrance/exit points should have opposite sign
in y coordinate).
 So I fill a map container map<Point,Id> in order to sort the
intersections and avoid multiple counting of the primitives.
When looking at my output I observe that:
1) The sorting is not always successful.
2)There are a lot of facets with all three indexes equal to zero.
3) Different intersection points correspond to the same vertex index
(facets with two indexes equal to zero). I mean I have a different
accuracy inserting points in the map and associating indexes in the
triangulation.
Could it be related to facets not belonging the 3D complex?
Thanks in advance!



#include <list>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.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>

namespace CGAL {
        template <class Tr>
        class AABB_triangulation_triangle_primitive
        {
        public:
                // types
                typedef typename Tr::Finite_facets_iterator Id; // Id type
                typedef typename Tr::Geom_traits GeomTraits;
                typedef typename GeomTraits::Point_3 Point; // point type
                typedef typename GeomTraits::Triangle_3 Datum; // datum type

        private:
                // member data
                Id m_it; // iterator
                Datum m_datum; // 3D triangle

        public:
                AABB_triangulation_triangle_primitive() {}

                AABB_triangulation_triangle_primitive(Id it)
                : m_it(it)
                , m_datum(it->first->vertex((it->second+1)&3)->point(),
                                  
it->first->vertex((it->second+2)&3)->point(),
                                  
it->first->vertex((it->second+3)&3)->point()){}

        public:
                Id& id() { return m_it; }
                const Id& id() const { return m_it; }
                Datum& datum() { return m_datum; }
                const Datum& datum() const { return m_datum; }

                /// Returns a point on the primitive
                Point reference_point() const { return m_datum.vertex(0); }
        };
}  // end namespace CGAL

// 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 C3t3::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_triangulation_triangle_primitive<Tr> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef Tree::Object_and_primitive_id Object_and_primitive_id;
typedef Tree::Primitive_id Primitive_id;
typedef Tr::Finite_facets_iterator Id;
typedef K::Triangle_3 Datum;
// 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 sort_pred {
    bool operator()(const std::pair<Point,Id> &left, const
std::pair<Point,Id> &right) {
        return left.first < right.first;
    }
};

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

        // Criteria
        Facet_criteria facet_criteria(30, 0.9, 0.01); // angle, size, 
approximation
        Cell_criteria cell_criteria(3, 0.9); // 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);

        Tree
tree(c3t3.triangulation().finite_facets_begin(),c3t3.triangulation().finite_facets_end());
        Point a(0.0, 0.0, 0.0);
       Point b(0.0, 1.0, 0.0);
        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" << " "<<"crossed 
vertex"<<"
"<<c.x()<<" "<< c.y()<<" "<<c.z()<<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<Primitive_id> primitives;
                std::map<Point,Id> mymap;
                std::map<Point,Id>::iterator it;
                tree.all_intersected_primitives(line_vertex,
std::back_inserter(primitives));
                for(; it_intersections != intersections.end() ; 
++it_intersections)

                {
                        // get intersection object
                        Object_and_primitive_id op = *it_intersections;
                        CGAL::Object object = op.first;
                        Id id=op.second;


                        Point point;
                        Segment segment;
                        if(CGAL::assign(point,object))
                        {
                                it=mymap.begin();
                                mymap.insert (it, make_pair(point,id));
                                //coppia.unique();
                        }
                        else if(CGAL::assign(segment,object))
                                std::cout << "intersection object is a 
segment" <<std::endl;

                }
                int ipr=1;
                for ( it=mymap.begin() ; it != mymap.end(); it++ ){
                        std::cout << ipr++<<" "<<"NINT "<<"
"<<std::setprecision(25)<<(*it).first << " => " << std::endl;
                        C3t3::Cell_handle cell    = (*it).second->first;
                        int opposite_vertex_index = (*it).second->second;
                        for(int i = 0; i < 4; i++)
                                if(i != opposite_vertex_index) {
                                        std::cout << V[cell->vertex(i)]<< " 
"<< std::endl;
        }
                }


        }

        return 0;
}