CGAL users discussion list

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.
2) I obtain on average only 30 intersection points for about 30000
triangles. How to establish whether the number of intersections is
"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?

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