CGAL users discussion list

Hi,
I'm sorry for posting again this message but I'm not able to understand
the result which I obtain and where the procedure doesn't work.
I computed the intersection points of a ray and a 3D object represented as
a 3D complex (for example a tetrahedron mesh in 3D mesh generation).  I
combined a 3D mesh generation procedure with an AABB tree. My query object
is a line with fixed direction and crossing a vertex of the mesh (I
consider parallel rays crossing each vertex for different intersection
queries) and I use a kind of AABB_primitive, such as Datum is Triangle_3
and Id is Tr::Facet_finite_iterator. Thus I can get back the triangulation
facet from the intersected Object_and_primitive_id.
For each line which crosses a vertex, I need to keep/control just the
intersection point "before" this vertex. The first step is to define the
point "before" with respect to the vertex of interest. My idea is to
select the two cells incident to the vertex (crossed by a fixed line) and
having as facet opposite to this vertex the facet intersected by the line.
In order to find this cells and the intersection point "before", I stored
a vector pair of  intersection points/intersected primitives  which are
facets of the triangulation. Then I run through the facets and I select
just the facet which have their "opposite vertex_index" equal to the index
of the vertex of interest (crossed by a fixed line).
 I would expect to find just two incident cells and "compare" the two
intersection points lying in the facets opposite to the vertex.
Conversely I detect a lot of facets where the condition to be facets
opposite to the vertex is fulfilled.

So I'm wondering whether:

1) Is my criterion sufficient to choose the two opposite facet? Or I have
to select all the cells incident to the vertex before setting the
condition on the "opposite vertex_index"? (In this case I have an error in
the first argument of the function
c3t3.triangulation().incident_cells(vh,std::back_inserter(cells)), where I
set std::vector<Cell_handle> cells. I see that vh=Vertex_handle but I
don't know how set vh to the vertex pointed by the finite_vertex iterator
I'm using).

2) How  can I compare the two intersection points lying on these facets in
order to establish the point on the line "before" and "after" the vertex
under study? Is a sort function a good choice or I have to evaluate the
sign of the scalar product between the direction of the oriented line and
a vector connecting the two points?

3) Moreover I'm setting the condition c3t3.is_in_complex(cell) to be
assured of getting facets belonging the 3D complex. Does it correspond to
get facets both inside the volume and on the surface or surface facets are
excluded?

I'm sending just the part where my criterion is used. I hope it can be
understood:

....DEFINE CLASS AABB PRIMITIVES (Tr::Finite_facets_iterator Id)
....DEFINE TYPES

int main()
{....MESH A SPHERE AND SET THE TREE
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 = 0;
for( Finite_vertices_iterator vit =
c3t3.triangulation().finite_vertices_begin();
        vit != c3t3.triangulation().finite_vertices_end();
        ++vit)
{
        V[vit] = inum++;
}

for( Finite_vertices_iterator vit =
c3t3.triangulation().finite_vertices_begin();
        vit != c3t3.triangulation().finite_vertices_end();
        ++vit)
{
       Point c = vit->point();
        Line line_vertex(c,dir);
        ......COMPUTE INTERSECTIONS AND STORE A COUPLE
std::vector<std::pair<Point,Id> >myvector

        for ( it=myvector.begin() ; it != myvector.end(); it++ ){
                C3t3::Cell_handle cell    = (*it).second->first;
                int opposite_vertex_index = (*it).second->second;
           if(opposite_vertex_index= V[vit]) {
                           if(c3t3.is_in_complex(cell)) {
std::cout << " intersection/coordinates"<<"
"<<std::setprecision(15)<<(*it).first << " => " << std::endl;
                  for(int i = 0; i < 4; i++){
                                          const Vertex_handle& vh =
(*((*it).second)).first->vertex(i);
std::cout<<" vertex index:"<<" " <<V[vh] << " "<<std::endl;
std::cout <<" vertex coordinates:"<<" "<<vh->point().x()<<"
"<<vh->point().y()<<" "<<vh->point().z()<<std::endl;
                                             }
   }
  }
 }
}
        return 0;
}

Thank for any suggestions!