CGAL users discussion list

Hi!
In a recent mail I was wondering about the output obtained combining a 3D
mesh generation
procedure with an AABB tree.
I had to compute the intersection points of a straight line and a 3D
object represented as a 3D complex.  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). 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).
So a kind of primitives such as Datum is Triangle_3 and Id is
Tr::Facet_finite_iterator allows me to 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) there are a lot of facets with three or two indexes equal to "0", while
the vertex map in the triangulation starts with the index 1.
2)  the uniqueness of the map is not kept for this index, that is the
multiple index "0" is associated with different vertex coordinates.

I checked whether it was related to facets not belonging the 3D complex
(for both facets in the volume and on the surface by means of the function
c3t3.is_in_complex) but this doesn't  solve the problem.......

Thanks in advance for any suggestion or comments!



#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 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_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;
}

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

        // Criteria
        Facet_criteria facet_criteria(30, 0.8, 0.01); // angle, size, 
approximation
        Cell_criteria cell_criteria(3, 0.8); // 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.3, 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)
                << " intersected primitives" <<" "<< " with line query 
crossing the
vertex"<<" "<<std::setprecision(25)<<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));
                        }
                        else if(CGAL::assign(segment,object))
                                std::cout << "intersection object is a 
segment" <<std::endl;
                }
                //write intersection points coordinates/triangles indexes
                int ipr=1;
                for ( it=mymap.begin() ; it != mymap.end(); it++ ){
                        C3t3::Cell_handle cell    = (*it).second->first;
                        int opposite_vertex_index = (*it).second->second;
                if(c3t3.is_in_complex(cell)) {
                                std::cout << ipr++<<"th intersection"<<"
"<<std::setprecision(25)<<(*it).first << " => " << std::endl;
          for(int i = 0; i < 4; i++){
                          if(i != opposite_vertex_index) {
                          std::cout << "vertex index"<<" 
"<<V[cell->vertex(i)]<< " "<< std::endl;
                          std::cout <<" vertex coordinates:"<<" 
"<<cell->vertex(i)->point().x()
                          << "," << cell->vertex(i)->point().y()
                          << "," << cell->vertex(i)->point().z()<<std::endl;
                  }
                          }
                                }
                                else 
if(c3t3.is_in_complex(cell,opposite_vertex_index)) {
                                        std::cout << ipr++<<"th 
intersection"<<"
"<<std::setprecision(25)<<(*it).first << " => " << std::endl;

                                        for(int i = 0; i < 4; i++){
                                                if(i != 
opposite_vertex_index) {
                                                std::cout << "vertex 
index"<<" "<<V[cell->vertex(i)]<< " "<< std::endl;
                                                std::cout <<" vertex 
coordinates:"<<" "<<cell->vertex(i)->point().x()
                                                << "," << 
cell->vertex(i)->point().y()
                                                << "," << 
cell->vertex(i)->point().z()<<std::endl;
                                                std::cout<<"\n";
                                                }
                                                }
                                }

                        }

                }

                return 0;
        }